7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Numerical evidence for a novel nonaxisymmetric bubble shape regime
in square channel Taylor flow
Martin Wojrner
Karlsruhe Institute of Technology (KIT), Institute of Nuclear and Energy Technologies,
P.O. Box 3640, 76021 Karlsruhe, Germany
martin.n\ come nl6 krIil idu
Keywords: Square channel, Taylor flow, transition to annular flow, numerical simulation
Abstract
In this paper we present numerical investigations of cocurrent downward Taylor flow in a square vertical minichannel by a
volumeoffluid method. The focus is on the shape of the front and rear meniscus at different values of the capillary number,
which is here in the range 0.04 I Ca I 0.66. At low values of Ca the bubble tip and rear adopt a hemispherical shape. With
increase of Ca the curvature of the front meniscus increases while that of the rear meniscus decreases. This behavior is
qualitatively similar to that in circular capillaries. For the highest value of Ca we find a novel regime in a square channel,
where the bubble tip and body are axisymmetric, whereas the rear meniscus is not axisymmetric but shows symmetry with
respect to the midplanes and diagonals of the channel. We show that this break of symmetry is related to the pressure
distribution in the liquid phase, which is nonuniform in crosssections close to the bubble rear with the highest values
occurring in the channels comers. Results for the evolution of the bubble shape during the transition from Taylor flow to
annular flow are also presented.
In the present study we investigate the shape of the front
and rear meniscus for cocurrent downward Taylor flow in
a square channel for varying values of the capillary
number in more detail. At the highest value of Ca, we find
a novel steady bubble shape regime, where the bubble tip
and body are axisymmetric while the rear of the bubble is
not. In an axial crosssection close to the rear meniscus, the
bubble shape adopts a square with rounded corners. To the
best of our knowledge, this regime has not been found in
experiments or numerical simulations so far. Here, we
show that this break of symmetry from an axisymmetric
state to one with midplane and diagonal symmetry is
related to the pressure distribution in the liquid film
surrounding the bubble and has its origin in a substantial
crosssectional variation of the pressure in the film and
corner region of the square channel. Finally, we investigate
in this paper the transient evolution of the bubble shape
that occurs in a square minichannel when there is a
transition from Taylor flow to annular flow.
Computational setup
In this section we give a short description of the numerical
method and the computational setup. The timedependent
threedimensional computations are performed with an
inhouse computer code, called TURBITVOF. This code
solves the NavierStokes equation with surface tension
term in nondimensional single field formulation for two
incompressible Newtonian fluids with constant viscosity
and coefficient of surface tension on a regular staggered
Cartesian grid by a finite volume method. All spatial
derivatives are approximated by central differences. Time
Introduction
Taylor flow is a special kind of slug flow in small channels,
where the liquid slugs which separate the elongated
bulletshaped bubbles (Taylor bubbles) are free from gas
entrainment. Taylor flow occurs in microfluidic devices
for applications in life sciences, material synthesis and
chemical process engineering (Giinther & Jensen, 2006;
Kreutzer et al., 2005). The hydrodynamics of Taylor flow
is mainly governed by two nondimensional groups, the
capillary number Ca E puL Ug / o and the Reynolds
number Re a pt~ DI, LL / pUL. Here, LL is the magnitude of
the bubble velocity, puL and pt~ are the liquid viscosity and
density, respectively, a is the coefficient of surface tension
and DI, is the channel hydraulic diameter.
In vertical round channels, the bubble shape is always
axisymmetric, whereas for cocurrent upward flow in
square vertical channels the bubble shape is axisymmetric
for Ca > 0.04 and nonaxisvmmetric for Ca < 0.04
(Thulasidas et al., 1995). Recently, we used an inhouse
computer code based on the volumeoffluid method and
performed a computational study of the cocurrent
downward Taylor flow of nitrogen bubbles in a viscous
liquid (squalane) in a square minichannel for capillary
numbers in the range 0.04 I Ca I 0.66 (Keskin et al.,
2010). We compared the computed steady bubble shape
with experimental flow visualizations obtained by a high
speed CCD camera and found good agreement. The
experimental and numerical results show both an
increase/decrease of the curvature at the front/rear
meniscus as the capillary number increases.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Results and Discussion
Bubble shape in numerical simulations
In Figure 1 we show the steady bubble shape at different
values of Ca as computed by Keskin et al. (2010). The
displayed bubble shape does not correspond to a certain
isosurface of the liquid volume fraction but is obtained as
follows. For each mesh cell that contains both phases the
center of area of the plane which represents the interface in
this mesh cell is computed. The centers of area of such
planes in neighboring mesh cells that contain an interface
are then connected to form triangles or quadrangles so that
a closed surface is obtained.
From Figure 1 we observe that with increase of the
capillary number the front of the bubble becomes more
pointed while the bubble rear becomes more flat (see also
Keskin et al., 2010). This result is in agreement with the
numerical study of Taha & Cui (2006) who found that in a
square channel a single Taylor bubble acquires spherical
ends at low values of Ca, while the back of the bubbles
changes from convex to concave at high Ca. This behavior
is similar to circular channels (Martinez & Udell, 1989;
Giavedoni, & Saita, 1999).
integration is performed by an explicit third order
RungeKutta method. A divergence free velocity field at
the end of each time step is enforced by a projection
method, in which the resulting Poisson equation is solved
by a conjugate gradient technique. The dynamic evolution
of the interface is computed by an unsplit volumeoffluid
method with piecewise planar interface reconstruction. For
further details about the governing equations and the
numerical method we refer to Sabisch et al. (2001) and
Ojztaskin et al. (2009).
For the computational setup, we follow the procedure of
our previous papers and consider one unit cell, which
consists of one gas bubble and one liquid slug. We use in
axial (vertical) direction periodic boundary conditions to
mimic the influence of the trailing and leading bubble in
Taylor flow. Noslip boundary conditions are applied at the
four lateral walls of the square channel. In accordance with
the experiments, the inner dimensions of the square
channel are 1 mm x 1 mm. The axial length of the unit cell
L,, is either 4 mm or 6 mm. In the former case the uniform
grid consists of 80 x 320 x 80 cubic mesh cells, while in
the latter case the number of grid cells is 80 x 480 x 80. So
in total up to about 3106 meSh cells are used and typically
a few 10,000 up to 100,000 time steps are computed.
In accordance to the experimental results reported in Bauer
(2007) and Keskin et al. (2010) we use as continuous
liquid phase squalane (C3 H62) While the disperse gas phase
is nitrogen. These experiments were performed at a
pressure of 20 bar. The corresponding fluid properties of
nitrogen used in the simulations are PG 23.6 kg/m3 and
pG~=0.01804 mPas. Since the physical properties of
squalane at a pressure of 20 bar are not available to our
knowledge, we use the known (constant) properties at
standard conditions which are pe~ 802 kg/m3, Lu=0.029
Pas, and 0 = 0.0286 N/m. Thus, the viscosity of squalane
is about 30 times higher than that of water.
The initial phase distribution of the simulations was
defined by placing an elongated axisymmetric bubble of
given volume on the channel axis. By this procedure the
gas holdup a is fixed. The initial velocity field for both
phases was given by fluids at rest, or, to save CPU time, by
a constant axial velocity, or by a parabolic axial velocity
profile within the channel crosssection (which was axially
uniform). Starting from these initial conditions, the flow is
driven by a prescribed constant source term in the axial
momentum equation which corresponds to the axial
pressure drop along the unit cell. In the course of the
simulation, the evolution from the initial velocity field and
prescribed bubble shape toward a fully developed Taylor
flow is computed. The latter is assured by recording the
mean axial gas and liquid velocities in the computational
domain and continuing the simulation till both velocities
approach constant terminal values. This final steady bubble
shape corresponds to a certain (constant) bubble velocity
and capillary number Ca. The Reynolds number Re is
related to the capillary number by Re = LaCa, where La 
oplDh L2 iS the Laplace number, which is constant here
(La = 27.27). Furthermore, the Weber number is given by
We = CaRe = pLDhUB2Io
Figure 1: Side view (top) and perspective view (bottom)
of computed steady bubble shape for Lee = 6 mm and a =
0.4 for five different values of the capillary number. From
left to right the values of (Ca; Re) are (0.045; 1.22),
(0.117; 3.19), (0.17; 4.64), (0.25; 6.81), (0.491; 13.4). (For
further details about the simulations see Keskin et al.,
2010). The solid lines indicate the size of the
computational domain.
7
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Thus, by virtue of Eq. (1) at any point on the bubble
surface, the local curvature of the interface, K, is related to
the pressure p and normal viscous stress r' on both sides
of the interface. For a Newtonian fluid, the normal viscous
stress at the liquid and gas side of the interface is given by
Taylor (1961) showed for a circular pipe that at small
values of the capillary number the streamlines in the liquid
slug show a recirculation pattern in a frame of reference
moving with the bubble, while at large Ca socalled bypass
flow occurs. Thulasidas et al. (1997) found experimentally
that for cocurrent upward flow in a square channel this
transition to bypass flow occurs for Ca 0.47. The highest
capillary number in the simulations of Keskin et al. (2010)
is obtained for case 4 20 B (in this case denomination the
first and second numbers indicate the axial length and the
gasholdup while the letter refers to different values of the
pressure drop). For this case it is Lee = 4 mm, a = 0.2, Ca =
0.655 and Re = 17.86. It is interesting to note that thought
this value is higher than 0.47, there is now bypass flow
because it is ry UB tot = 1.96 < 2.096 (see discussion in
the next subsection). This suggests that the transition from
recirculation flow to complete bypass flow occurs in
downward flow at higher values of Ca than in upward flow.
This conjecture is supported by Fig. 5 in Kececi et al.
(2009), where ry is displayed as function of Ca and
numerical results for upward and downward flow are
compared. This figure shows that for a given value of Ca,
ryis higher for upward flow than for downward flow.
In Figure 2 we show the bubble shape for case 4_20_B.
This figure shows that the rear of the bubble is flat. At the
axial position of the bubble, where the liquid film is
thinnest and the bubble diameter is largest, the bubble
shape is axisymmetric. At the rear of the bubble, however,
this symmetry is lost. Figure 2 b) shows that the red line
which denotes the position of the interface in a
crosssection close to the bubble rear is not a circle but
adopts the shape of a square with rounded corners. Thus, at
the rear of the bubble the interface is not axisvmmetric.
Instead, there is symmetry with respect to the two
midplanes and diagonals of the square channel. To the
authors knowledge this behavior has not been observed so
far neither in experiments nor in numerical simulations.
While the measurements of the threedimensional shape of
the rear meniscus of a moving Taylor bubble is a difficult
subject by itself, it is further complicated by the fact that
the axial extension of the nonaxisymmetric shape reported
here is only about 50 pLm.
Theoretical analysis of local interface curvature
In this subsection we investigate the reasons for the shape
of the rear meniscus as displayed in Figure 2 b). To this
end we consider the dynamic boundary condition at the
interface, i.e. the local force balance between pressure,
viscous stresses and surface tension. Here, we are only
interested in the projection of this dynamic boundary
condition in direction of the unit vector n,, which is normal
to the interface and points into the continuous liquid phase.
Then it is
pL,1 G,1 + LI G~,1 = ~ 1
with Kbeing twice the mean local interface curvature
S= 2H =1/ Rmin+1/ Rmax (2)
1 T
LI VvL +(VVL 1
T MG VVG +(VVG
G,1 1
Figure 2: Bubble shape for case 4_20_B in a) perspective
view and b) view from behind (only one half of the bubble
is shown). The two lines indicate the bubble shape in an
axial crosssection in the middle of the bubble (blue line)
and very close to the bubble rear (red line).
An evaluation of the local pressure field shows that the
pressure in the bubble is constant (see also Figure 3).
Therefore, we assume per = PB = COnst. In the present
simulations, the ratio between the gas and liquid viscosity
takes a value of about 0.0006. Therefore, we assume that
the normal viscous stress at the gas side of the interface
can be neglected. Then Eq. (1) yields the following
relation for the nondimensional local interface curvature,
K rDh =L
B, L ~
When the bubble moves either in positive or negative
vdirection (unit vector e ), then the dyadic product at the
bubble front and rear becomes n,n, = eye, and we have with
VL = (UL vLWL)T for the normal viscous stress on the liquid
side of the interface at the bubble tip and rear the relation
Li,tip rea = r VL, + (VVL) TI~ ti rear
c1L (6)
= 24u ~~_
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The velocity gradient in ydirection in the liquid phase at
the bubble tip and rear may be approximated as (see Fig. 4)
Here, 1101l is the magnitude of the axial velocity in the
liquid slug on the channel centerline and Lnwbl iS the axial
distance across which this velocity change occurs. For a
fully developed liquid slug the axial centerline velocity is
given by
K(xl ) =", (pB L (X
Thus, the knowledge of the pressure field on the liquid side
of the interface is the key to understand the bubble shape.
Since variations of the pressure in the channel
crosssection are related to the velocity field, we
investigate in the next subsection both, the local pressure
field and the local velocity field.
Analysis of local pressure and velocity field
Figure 3 shows for case 4 20 B of Keskin et al. (2010) the
axial profile of the pressure for different positions within
the channel crosssection. While in the governing
equations adding a constant to the pressure field is without
influence as only the pressure gradient matters, in Figure 3
for clarity of presentation a constant has been added to the
pressure values so that the minimum pressure at the bottom
of the channel (v = 0) is zero. From Figure 3 it is evident
that in the liquid slug the pressure within an axial
CTOsssection is uniform as it is typical for laminar single
phase flow in a straight channel. Furthermore, in the liquid
Slug the axial pressure drop is constant. Figure 3 also
shows that the pressure in the bubble is about constant (see
profile for channel centerline in the range y = 2 4 mm). It
further indicates that a notable variation of the pressure
within an axial crosssection occurs only at the bubble rear.
There, the pressure difference between the bubble and the
liquid is clearly larger in the lateral film region than in the
corner flow region. At the rear edge of the interface the
normal viscous stresses may be neglected: therefore for
this region Eq. (12) may be used to estimate the local
interface curvature. As Figure 3 shows the difference
between the pressure in the bubble and on the liquid side
of the interface is finite in lateral direction, whereas it is
almost zero in diagonal direction. According to Eq. (12),
therefore, the local interface curvature close to the center
of the channel side walls is finite (and positive) while it is
about zero in diagonal direction. This explains the
crosssectional interface shape as given by the red line in
Figure 2b).
Figure 3 shows that the axial pressure drop in the liquid
phase across the bubble length occurs to the main part over
the first third of the bubble. In the region, where the
thickness of the lateral liquid film is about constant, the
pressure is nearly constant too, and no significant axial
pressure drop exists.
It is obvious that in the bubble region the pressure field is
closely related to the velocity field. Hazel & Heil (2002)
analyzed the steady propagation of a semiinfinite bubble
into a channel of rectangular crosssection by solving the
freesurface Stokes equations and investigated in detail the
pressure and velocity field in the liquid phase mn axial
crosssections at various distances from the bubble tip.
They found that in square channels the fluid particles tend
to move towards the corners, which offer less resistance to
the flow than the thinner regions along the sides of the
channel. The ensuing transverse flows induce a transverse
pressure gradient that lowers the fluid pressure in the
E L,cl = Cotot
2.0962Jtot
Here, Jtor & a B + (1 8) UL iS the total superficial
velocity. 11 denotes the magnitude of the mean axial liquid
velocity within the unit cell. Using Eq. (8) and defining r
a UB/J we can write
Dh PB L,tlp) 2Dh Cn C
Ktip = Lml,
Lnb,tip l
Lnyohbl,rea Cr
Dh PB L,rear)
Ker =
Thus, the effect of the normal viscous stress on the
curvature of the bubble front and rear depends on the ratio
C. / ry For recirculation flow it is ry < C. and the normal
viscous stress acts to increase the curvature of the bubble
tip and to decrease the curvature of the bubble rear. For r
> C. we have bypass flow and Eqs. (9) and (10) suggest
that the curvature of the bubble tip is decreasing with
increase of Ca, while that of the bubble rear is increasing.
Here, we have always ry< C. and the second term on the
r.h.s. of Eq. (9) is positive.
Combining Eqs. (9) and (10) yields
Dh PL,rear PL,tip )
Ktip Krear
+2 r Lbl,tlp 1b,rear (I
Thus, the difference in curvature of the front and rear
meniscus depends on the pressure drop in the liquid phase
over the length of the bubble. For a given pressure
difference pB Lp,tp, this difference in curvature increases
with increase of Ca for ry< C. and decreases with increase
of Ca for ry> C, where ryitself is a function of Ca.
With exception of the bubble front and rear stagnation
point, the normal viscous stresses may be neglected. Then,
the following relation for the local curvature at an arbitrary
point x, lying on the interface should approximately hold
ur , , , .
0 1 2 3
y [mm]
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.001m(812Pa 541Pa)
Ktl
0.0286N/m
22.096 (3
+0.655 1(3
0.2 1.96
= 9.475 +0.455 = 9.93
and
0.001m(812Pa 770Pa)
K =
rear 0.0286N/m
0.655 1.9 (14)
0.16 1.96
= 1.468 0.568 = 0.9
While these values should be considered as estimates of
the magnitudes of the relevant forces, they nevertheless
show that in the present case 4 20 B the curvature of the
bubble tip is one order of magnitude larger than the one at
the bubble rear. Furthermore it is evident that at the bubble
tip the normal viscous stresses are one order of magnitude
smaller than the pressure forces, while at the bubble rear
the magnitude of the normal viscous stresses is about 38%
of that of the pressure forces. Since at the bubble rear the
normal viscous stresses on the liquid side of the interface
act opposite to the pressure forces for ry < C., they are
definitely important and should not be neglect at high
values of the capillary number.
corners, where the gasliquid interface moves radial
outwards. In the region behind the bubble tip, surface
tension acts to restore the gasliquid interface to an
axisymmetric shape, and the nonuniform pressure
distribution continues to drive fluid into the corner until the
bubble diameter in lateral and diagonal direction are equal.
As concerns the case 4 20 B, Figure 3 shows that in axial
regions which are closer to the bubble tip than to the
bubble rear, the pressure within a certain axial
crosssection is indeed slightly lower in the corners than at
the channel sides. However, in axial crosssections which
are closer to the bubble rear than to the bubble tip, we find
the opposite behavior, i.e. higher pressure values in the
channel corners than at the channel sides.
*Inn
centerline (itk)
 (40;40)
lateral (i;k) diagonal (i;k)
(3;40)  (7;7)
(6;40) (14;14)
 600
400
200 
Figure 3: Axial profile of pressure at different positions mn
the channelcross section for case 4 20 B in Keskin et al.
(2010). The values (i,k) denote the mesh cell index within
the channel crosssection (1 I i I 80; 1 I k & 80).
In Figure 4 we show axial profiles of the liquid axial
velocity at different positions within the channel
crosssection. We find that the axial liquid velocity is, aS
expected, in general higher in the corner region than in the
lateral liquid film region. However, at a position close to
the bubble rear and close to the interface the liquid axial
velocity is about the same in lateral and diagonal direction
(compare in Figure 4 the profile for i = k =14 with the one
for i = 7, k =40 at y 4mm). A comparison of the velocity
components within the channelcross section (not shown
here) indicates that the magnitude of the liquid velocity in
the crosssection is higher in the channel comers than at
the channel sides. Thus, in a crosssection close to the
bubble rear the liquid moves faster inward in the channel
corners than at the channel sides. This inward directed
liquid flow is associated with curved streamlines, which
result in the presence of nonnegligible inertial forces in a
nonuniform crosssectional pressure distribution which
has, in the present case, maxima in the channel corners.
However, further investigations are needed to elucidate the
intimate relation between the local pressure and velocity
fields and the interface curvature at the bubble rear.
From Figure 4 we can evaluate the length Lilvbl and find
Lilvbl,tip / Dh = 0.2 and Lbl,rear / Dh = 0.16, respectively. With
these data, the numerical values of the different terms in
Eqs. (9) and (10) become
anr
~ *C~ Y

0.2

 0.4 
0.6
centerline (i;k)
(4;40)
lateral (i;k)
(7;40)
diagonal (i~k)
(14;14)
0 1 2 3 4
Figure 4: Axial profile of liquid axial velocity at different
positions in the channelcross section for case 4 20 B in
Keskin et al. (2010). The values (i,k) denote the mesh cell
index within the channel crosssection (1 I i I 80; 1 I k I
80).
Transition from slug to annular flow
In this subsection we present and discuss results of a new
transient simulation for a square minichannel, where the
pressure gradient is so high that cocurrent downward
Taylor flow is not stable and a transition from slug flow to
annular flow occurs. In this simulation it is Le = 4 mm and
a = 0.4: we denote this case by 4 40 C. The initial bubble
shape is displayed in Figure 6. The initial velocity is the
 J L
ot [slu
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
liquid side of the interface in the lateral liquid film and the
liquid corner flow diminishes. Due to this, in the time
period t = 5 8 ms the rear meniscus is almost
axisymmetric again. For t > 8 ms the liquid slug length is
zero and the tip of the trailing bubbles penetrates in the
concave region of the rear meniscus of the trailing bubble.
At this stage the present simulation breaks down. The
reason is that during the process of the merger of the
trailing and leading interfaces a large number of small gas
and liquid fragments are produced. The curvature of these
tiny structures cannot adequately be resolved by the grid.
same in both phases and is constant, namely T>L =vG = 0.3
m/s. The prescribed mean pressure difference between the
top and bottom crosssection of the channel is for this case
550.3 Pa. This value is about 27% lower than for case
4 20 B, where it is 754.4 Pa. However, the hydrostatic
pressure head that must be overcome to drive the flow
downward is much lower for case 4_40_C. This is because
the gas holdup is 40% instead of 20% for case 4 20 B.
Figure 5 shows the temporal evolution of Jrot, LL and El
during the course of the simulation. As can be seen, the
downward velocity of the bubble strongly increases in time
whereas the magnitude of the mean liquid velocity in the
domain is slightly decreasing. At about 5 ms the slope of
LE versus t decreases somewhat before it is strongly
increasing for t > 8 ms. Thus, the bubble velocity
approaches no constant value: instead the bubble is always
accelerating.
I!
II
Figure 6: Visualization of instantaneous bubble shape for
case 4 40 C during the transition from Taylor flow to
annular flow for different instants in time. Top row (from
left to right): 0 ms, 0.27 ms, 0.53 ms, 0.67 ms, 1.93 ms,
3.72 ms: bottom row (from left to right): 5.05 ms, 6.65 ms,
7.45 ms, 7.99 ms, 8.42 ms, 8.62 ms.
Figure 5: Transition from Taylor flow to annular flow:
temporal evolutions of total superficial velocity, bubble
velocity, mean liquid velocity and liquid slug length for
case 4 40 C.
Also shown in Figure 5 is the time history of the liquid
slug length L,1. For each instant in time, this length is
evaluated as the axial distance along which in the entire
channel crosssection only liquid is present. The initial
liquid slug length is about 1 mm. At t = 5 ms this value is
reduced to about 0.3 mm while for t > 8 ms its value is
zero.
Figure 6 shows the bubble shapes for case 4_40_C at 12
different instants in time. As time proceeds we observe that
the rear meniscus of the bubble quickly deforms from the
initially hemispherical shape to a flat interface, which is
reached at t 2 ms. Characteristic for the time period from
t = 2 5 ms is a concave rear meniscus on the channel axis
which goes along with a non axisymmetric bubble shape
in some radial distance from the channel centerline. This
bubble shape will be discussed below. As time proceeds
further, the diameter of the bubble body continues to
decrease while the length of the bubble increases and
accordingly the length of the liquid slug decreases. As a
consequence of the decreasing bubble diameter, the flow
along the circumferential of the bubble becomes more
uniform because the difference in the axial velocity on the
Figure 7: Bubble shape for case 4_40_C at t = 3.72 ms in
a) perspective view and b) view from front (only one half
of the bubble is shown). The three lines indicate the
bubble shape in an axial crosssection in the middle of the
bubble (black line) and close to the bubble rear (blue and
red line).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
In Figure 7 we investigate the shape of the bubble rear for t
= 3.72 ms in more detail. Similar to Figure 2 we display
the shape of the interface in the channel crosssection at
certain axial positions. It is evident that the body of the
bubble is axisynunetric (black line in Fig. 7 b). Close to
the bubble rear we find again the characteristic shape of a
square with rounded corners (blue line) that was already
found for case 4 20 B, see Fig. 2 b). However, as
indicated by the small blue circle, in the present case the
interface at the rear meniscus is not flat but slightly
concave. In an axial position slightly more upstream we
find four closed isolines of the bubble interface (red lines).
These represent cusps that fonn where the bubble rear is
close to the middle of the lateral four walls of the vertical
channel. A detailed inspection of Fig. 7b) shows that the
red and blue lines are not perfectly synunetric with respect
to the vertical channel midplanes or the channel diagonals
(note that our simulation is full 3D and no kind of
synunetry is assumed). Instead, the respective patterns
slightly change in time. We speculate that these interface
defonnations and fluctuations are due to the continuous
acceleration of the bubble. It is an open question if such a
shape of the rear meniscus may persist under steady bubble
motion. To investigate this issue, a new simulation is
currently under way where we take the final time step of
case 4 20 B as initial condition, but increase the driving
pressure difference by about 10%.
Conclusions
In this paper we used results of threedimensional
numerical simulations of viscous cocurrent downward
Taylor flow in a square vertical minichannel to investigate
the shape of the front and rear meniscus at different values
of the capillary number. It is found that the curvature of the
front meniscus increases with increase of Ca while that of
the rear meniscus decreases. This behavior is in qualitative
agreement with that in circular channels. For the highest
capillary number, Ca = 0.65, we find a novel
nonaxisymmetric bubble shape regime which has to the
author's knowledge not been reported so far. This new
regime is characterized by a bubble shape which is
axisynunetric at the tip and body of the bubble but adopts
the shape of a square with rounded corners at the bubble
rear. We showed that this symmetry breaking is related to a
nonunifonn pressure distribution in crosssections close to
the rear meniscus, with local pressure maxima in the
channel corners.
We also investigated the evolution of the bubble shape
during the transition from Taylor flow to annular flow in a
square channel. During this transient simulation, the
bubble is always accelerating and the diameter of the
axisymmetric bubble body decreases in time while the
length of the bubble increases. In certain stages of this
simulation the rear meniscus shows the novel shape
described above and even fonns cusps close to the middle
of the four channel side walls. However, it is unclear if the
latter shape exists only when the bubble accelerates or also
under steady bubble motion.
Nomenclature
Ca
Dh
ey
g
H
Jtot
K
k
La
Lnyb1
Lee
nx
Rmm
Re
v
x
We
y
x
Capillary number ()
hydraulic diameter (m)
unit nonnal vector in axial direction ()
gravitational constant (m1 i'!
mean interface curvature (1/m)
mesh cell index in xdirection ()
total superficial velocity (m/s)
mesh cell index in ydirection ()
nondimensional interface curvature ()
mesh cell index in zdirection ()
Laplace number ()
thickness of nonnal viscous boundary layer (m)
length of the liquid slug (m)
Length of the unit cell (m)
unit nonnal vector to the interface ()
pressure ('j nr1)
minimum principal radius of curvature (m)
maximum principal radius of curvature (m)
Reynolds number ()
magnitude of bubble velocity (m/s)
magnitude of mean liquid velocity (m/s)
velocity field (m/s)
wallnonnal coordinate (m)
Weber number ()
axial coordinate (positive in upw. direction) (m)
Cartesian coordinates, x = (x,y,z)T m
wallnonnal coordinate (m)
Greek letters
a gas holdup in computational domain ()
K interface curvature (1/m)
pu dynamic viscosity (Pa s)
0 coefficient of surface tension (N/m)
ry = Gk/Jtot()
nonnal viscous stress (Pa)
Subscripts
B bubble
G gas phase
i interface
L liquid phase
nvb1 nonnal viscous boundary layer
rear rear meniscus on channel axis
tip front meniscus on channel axis
uc unit cell
References
Bauer, T. Experimental and theoretical investigations of
monolithic reactors for threephase catalytic reactions.
PhD thesis, Technical University Dresden, 2007.
Giavedoni, M.D., Saita, F.A. The rear meniscus of a long
bubble steadily displacing a Newtonian liquid in a
capillary tube. Phys. Fluids 11 (1999) 786.
Giinther, A., Jensen, K.F. Multiphase microfluidics: from
flow characteristics to chemical and material synthesis.
Lab Chip 6 (2006) 1487.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Ojztaskin, M.C., Wdmer, M., Soyhan, H.S. Numerical
investigation of the stability of bubble train flow in a
square minichannel. Phys. Fluids 21 (2009) 042108.
Sabisch, W., Wdmer, M., Grijtzbach, G, Cacuci, D.G. 3D
volumeoffluid simulation of a wobbling bubble in a
gasliquid system of low Morton number. 4th Int. Conf.
Multiphase Flow, ICMF2001, New Orleans, LA, USA,
May 27 June 1, 2001, CDROM.
Taha, T., Cui, Z.F.. CFD modelling of slug flow inside
square capillaries. Chem. Eng. Sci. 61 (2006) 665.
Taylor, G.I. Deposition of a viscous fluid on the wall of a
tube. J. Fluid Mech. 10 (1961) 161.
Thulasidas, T.C., Abraham, M.A., Cerro, R.L. Bubbletrain
flow in capillaries of circular and square cross section.
Chem. Eng. Sci. 50 (1995) 183.
Thulasidas, T.C., Abraham, M.A., Cerro, R.L. Flow
patterns in liquid slugs during bubbletrain flow inside
capillaries. Chem. Eng. Sci. 52 (1997) 2947.
Hazel, A.L., Heil, M. The steady propagation of a
semiinfinite bubble into a tube of elliptical or rectangular
crosssection. J. Fluid Mech. 470 (2002) 91.
Kececi, S., Wdmer, M., Onea, A., Soyhan, H.S.
Recirculation time and liquid slug mass transfer in
cocurrent upward and downward Taylor flow. Catalysis
Today 147S (2009) S125.
Keskin, Oj., Wdmer, M., Soyhan, H.S., Bauer, T.,
Deutschmann, O., Lange, R. Viscous cocurrent downward
Taylor flow in a square minichannel. AIChE J., in press,
doi: 10.1002/aic.12113.
Kreutzer, M.T., Kapteijn, F., Moulijn, J.A., Heiszwolf, J.J.
Multiphase monolith reactors: chemical reaction
engineering of segmented flow in microchannels. Chem.
Eng. Sci. 60 (2005) 5895.
Martinez, M.J., Udell, K.S. Boundary integral analysis of
the creeping flow of long bubbles in capillaries. J. Appl.
Mech. 56 (1989) 211.
