Journal of Computational Physics 230 (2011) 4466–4487
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Efficient numerical methods for multiple surfactant-coated bubbles in a two-dimensional stokes flow Mary Catherine A. Kropinski a,⇑,1, Enkeleida Lushi b,2 a b
Department of Mathematics, Simon Fraser University, Burnaby, British Columbia, Canada V5A 1S6 Courant Institute of Mathematical Sciences, New York University, New York, NY 10012, USA
a r t i c l e
i n f o
Article history: Received 5 January 2010 Received in revised form 7 February 2011 Accepted 14 February 2011 Available online 2 March 2011 Keywords: Fluid interface Stokes flow Insoluble surfactant Integral equations Bubbles Fast-multipole method
a b s t r a c t We present efficient and highly accurate numerical methods to compute the deformation of surfactant-coated, two-dimensional bubbles in a slow viscous flow. Surfactant acts to locally alter the surface tension and thereby change the nature of the interface motion. In this paper, we restrict our attention to the case of a dilute insoluble surfactant. The convection–diffusion equation for the surfactant concentration on the interface is coupled with the Stokes equations in the fluid domain through a boundary condition based on the Laplace-Young condition. The Stokes equations are first recast as an integral equation and then solved using a fast-multipole accelerated iterative procedure. The computational cost per time-step is only O(N log N) operations, with N being the number of discretization points on the interface. The bubble interfaces are described by a spectral mesh and is advected according to the fluid velocity in such a manner so as to preserve equal arc length spacing of marker points. This equal arc length framework has the dual advantage of dynamically maintaining the spatial mesh and allowing efficient, implicit treatment of the stiffest terms in the dynamics. Several phenomenologically different examples are presented. Crown Copyright Ó 2011 Published by Elsevier Inc. All rights reserved.
1. Introduction Surface active agents (surfactants), whether present as impurities or deliberately added, critically alter the dynamics of multiphase flow systems. Surfactant molecules typically consist of hydrophilic heads and hydrophobic tails, and thus they tend to adhere to, and accumulate on, interfaces between fluids. Their impact is to locally lower the surface tension and thereby alter the dynamics, often critically. Important phenomenological effects that may be influenced by surfactant include bubble or droplet creation by tip streaming, cusp formation, and drop coalescence or breakup. Not surprisingly, surfactant can play an important role in many industrial and biomedical applications such as drug delivery, hydro-desulfurization of crude oil, polymer blending, paints and plastic production, emulsifications, or pulmonary functioning. Computationally modeling these flows is very challenging. The interfacial surface tension depends on the surfactant concentration through a possibly-nonlinear equation of state. Nonuniform capillary (normal) and Marangoni (tangential) stresses are induced in the fluid, which significantly affects the deformation of the interface. The evolution equation for the surfactant concentration is stiff; surfactant diffuses by molecular mechanisms and is transported along the interface by fluid convection. There may be adsorption/desorption of interfacial surfactant to/from the bulk of the fluid. This evolution ⇑ Corresponding author. 1 2
E-mail addresses:
[email protected] (M.C.A. Kropinski),
[email protected] (E. Lushi). Supported in part by the Natural Sciences and Engineering Research Council of Canada. Supported in part by the Applied Mathematical Sciences Program of the US Department of Energy under Contract DEFG0200ER25053.
0021-9991/$ - see front matter Crown Copyright Ó 2011 Published by Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.02.019
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4467
equation is coupled with the equations governing the fluid motion, and this coupling is highly nonlinear. The shape of the interface can become very complex: regions of high curvature might develop and/or topological changes may occur. Numerically investigating the often-subtle phenomenological changes caused by the presence of surfactant (for example, whether a cusp forms in finite time, when and where surfactant caps form, the conditions required for tip streaming) clearly require a high degree of spatial and temporal accuracy. However, depending on the methods used, requiring higher accuracy may be achievable only at the expense of the size of problem that can be investigated. Even with significant advances in computer architecture, this trade-off continues to be negotiated. Recent examples in the literature support this point. Xu et al. [37] present a level set/immersed interface method for two-dimensional interfaces with insoluble surfactants. They present examples of up to four drops in a shear flow, with both a linear and nonlinear equation of state. Their methods include novel techniques to correct the characteristic loss of mass often seen in level set methods. However, even with these corrections, the change in drop area or surfactant mass can be up to 1%, with computational costs preventing the refinement needed for further accuracy. Integral equation methods are more amenable to numerical methods that can achieve high spatial accuracy, however, standard implementations [3,26] can be very costly to compute. Not surprisingly, interfaces are often left under-resolved. For example [26] reports numerical oscillations in the interface shape due to inadequate spatial resolution and in [3], quantities were conserved to only 0.5%. The linearity of the Stokes equations can be exploited by a variety of analytic techniques. Seigel and co-authors [31] use slender-body theory to examine the conditions required for tip streaming in a single axisymmetric bubble, presenting analytical results valid in certain regimes, only. In this paper, we hope to provide motivation that when suitably chosen numerical methods are coupled with modern, fast, algorithms, it is possible to retain a high degree accuracy without sacrificing efficiency. Here, we present highly-accurate and efficient numerical methods for computing the motion of inviscid interfaces (bubbles) in a two-dimensional Stokes flow in the presence of an insoluble surfactant. We employ integral equation methods based on the complex variable theory for the biharmonic equation and further developed by Kropinski [16,17] for studying interfacial motion in a Stokes flow. The discretization of the integral equations is spectrally-accurate, and the iterative solution is accelerated by using the fast multipole method (FMM) [5,27] to compute the matrix–vector products. With N points in the discretization of the boundary, the integral-equation solve requires only O(N) operations, versus standard implementations of iterative schemes (for example [20]) which require O(N2) operations. For the motion of the interface, we employ the ideas of Hou et al. [13] and recast the evolution equations into an equal arc length frame. In this manner, we are able to dynamically maintain marker points at equal arc length intervals. More importantly, we can easily remove the stiffness caused by the diffusion term in the surfactant transport equation, since an implicit treatment of this term becomes explicit in Fourier space. The dynamics, then, are amenable to efficient treatment by high-order, semi-implicit methods such as the IMEX Runge–Kutta methods discussed in [2]. There are a variety of other analytical and computational approaches, each having their unique advantages and disadvantages. Integral equation methods are a natural choice for Stokes flow, and we are by no means the first to adopt this approach (a non-exhaustive list includes [3,9,20,25,39]). In most of these studies, a cubic spline or a lower-order representation is used to describe the interface. Siegel [29,30] exploits the complex-variable theory for the biharmonic to derive analytical and semi-analytical methods (largely limited to studying a single bubble); these provide useful examples to verify other numerical methods. Methods that naturally handle topological changes (for example, drop breakup or coalescence) include level set and volume of fluid methods. For example, Xu et al. [37] use at a level set method to look at insoluble surfactants on droplets. Yang et al. [38] develop a hybrid level-set/front-tracking method to research interface problems for unstructured triangular grids. These methods are excellent for handling complicated, multi-component environments, but they tend to be low-order accurate. There has been recent work using the full Navier–Stokes equations [8,18,22]: however, these nascent studies are low-order accurate and simulations typically involve only a small number of bubbles or drops. There are a number of limitations in the physical model we investigate here: It is in two-dimensions. We do not consider the case of droplets containing fluid of a different viscosity. We consider the case, only, of insoluble surfactant with a linear constitutive relation between surfactant concentration and surface tension (although a nonlinear relationship can easily be substituted for the linear one). We cannot deal with topological changes. While the physical model is limited, we believe our methods are useful for two reasons. First, our techniques are highly accurate and they will serve as a useful benchmark for other methods that are developed to solve more complicated models. Second, we feel this paper serves as a preliminary step in demonstrating the philosophy that by coupling modern fast algorithms with well-conditioned integral equations, it is possible to develop highly accurate and efficient tools for investigation problems of scientific interest. While the problems presented in this paper may be of limited scientific scope, currently, the methods can and will be extended to include, for example, three-dimensions (the fast multipole method has already been adapted to handle the potentials associated with the Stokes equations in three-dimensions [35]) and/or drops with different viscosity ratios ([16] includes this effect, albeit in the clean flow case, only). We begin in the next section by outlining the governing equations for inviscid interfaces in a Stokes flow and the convection–diffusion equation for an insoluble surfactant. We show how the two are coupled through the Laplace-Young condition and an equation of state that relates surface tension to surfactant concentration. In Section 3, we outline the
4468
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
complex variable theory for the biharmonic equation, and present the resulting integral equations for the problems at hand. The equal arc length framework is discussed in Section 4, and we discuss the interface and surfactant dynamics in this framework. Our numerical methods are presented in Section 5. In Section 6 we use a semi-analytical solution of Siegel’s [29] to verify our numerical methods. Further examples demonstrating a variety of different phenomena are given in Section 7. In Section 8, we state our conclusions. 2. The governing equations 2.1. The fluid dynamics We consider the motion of a collection of two-dimensional closed fluid interfaces immersed in a slow viscous flow. The fluid domain is multiply connected, the interfaces enclose regions of an inviscid fluid, e.g. air, and are thus excluded from the flow domain, and is unbounded, or wall-bounded in extent (see Fig. 1). We assume that the Reynolds number is small, thus the equations governing the fluid motion are the Stokes equations. Presented in their non-dimensional form, these are:
r2 u ¼ rp;
r u ¼ 0;
x 2 X;
ð1Þ
where u = (u, v, 0) is the fluid velocity and p is the pressure. Here, the variables have been scaled with respect to suitable characteristic values: length is the radius a of a typical bubble, velocity as r0/l where r0 is the characteristic surface tension (the surface tension without surfactant) and l is the kinematic viscosity (assumed to be constant), time as a l/r0, and pressure as r0/a. The boundary condition on the interface is given by the Laplace-Young condition,
ðp pk Þn þ 2En ¼ rjn þ rs r þ Bng x;
x 2 @ Xk ;
ð2Þ
where pk is the internal pressure inside oXk (effectively, this internal pressure ensures the area of each bubble remains constant, c.f. [17] for details). Here, E is the rate-of-strain tensor
Eij ¼
1 @ui @uj ; þ 2 @xj @xi
where the indices i and j take on the values 1 or 2 corresponding to the x or y directions, respectively, r is the (variable) surface tension and j is the local curvature of oX (see Fig. 1 for orientation). The rsr term in (2) gives the tangential stress (Marangoni force), which results from the dependence of the surface tension on a nonuniform surfactant concentration. Gravity effects are included through the term containing the bond number B, which is defined by
B¼
Dðq qb Þga2
r0
;
where qb is the density of the inviscid fluid/gas in the bubble. Here, ng is the unit normal pointing in the direction of gravity. We also require that the following far-field conditions be satisfied:
u ! u1 ;
p ! p1 ðtÞ;
as jxj ! 1:
We assume each oXk is parametrized by a 2 [0, 2p]. For x(a, t) 2 oXk, the motion of the interface is determined by integrating the following in time:
dxða; tÞ ¼ Un þ Ts; dt
x 2 @ X:
ð3Þ
Fig. 1. A collection of bubbles in an unbounded viscous fluid domain X with boundary @ X. The component boundaries are denoted by @ X1, @ X2, . . . , @ XM. The unit normal n points out of X and h is the tangent angle to C. The local curvature is given by j = hs, where s is arc length increasing in a clockwise direction.
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4469
Here, the tangential velocity T is
Tða; tÞ ¼ Tð0; tÞ
Z
a
ha0 Uda0 þ
0
a 2p
Z 2p
ha0 Uda0 ;
ð4Þ
0
where h is the tangent angle to oXk (again, see Fig. 1 for orientation). This choice for T maintains the equal arc length frame: sa(t) = jxaj = Lk(t)/2p, where Lk(t) is the total length of oXk at time t. As discussed in [13,16,17], the equal arclength framework eliminates the need to post process the discretization points along the interface, and it prevents the stiffness exhibited when these points are allowed to cluster together. As we will see in the following section, it will also enable efficient temporal integration of the surfactant convection–diffusion equation. 2.2. The surfactant convection–diffusion equation The non-uniformity of the surface tension r arises from its dependence on the surface surfactant concentration C. This is given by an equation of state of the form r = r(C). The relation between surface tension and surfactant concentration can be described by a logarithmic dependence called the Langmuir equation of state [9,30,37]. When the actual surfactant concentration is much lower than the closest packing case, i.e. the dilute limit, the equation of state can be approximated by a linear relationship (c.f. discussions in [24,26,37,39]). Written in non-dimensional form, this relationship is
r ¼ ð1 bCÞ;
ð5Þ
where b is a non-dimensional constant expressing the sensitivity of the surface tension to the surfactant concentration; it is often called the elasticity parameter. We will use the linear equation here but note that the logarithmic one can be easily substituted in with no additional computational complication. As shown in [32,36], the evolution equation for an insoluble surfactant on an interface is
@C 1 dx rs C rs ðCðu sÞsÞ þ Cðrs nÞU; ¼ r2s C þ Pe dt @t
x 2 @ Xk :
ð6Þ
The surfactant diffusion is controlled through the surface Peclet number,
Pe ¼
r0 a ; lDs
where Ds is the surfactant diffusivity constant. We simplify (6) by observing that in two-dimensions, rs = s@ s. Let S be the tangential component of fluid velocity on oX, S = u s, and recall from (3) that s dx/dt = T. Employing the Frenet formulae @s/@s = jn, and @n/@s = js we obtain
@C 1 @2C @ C @ðCSÞ jU C; þT ¼ Pe @s2 @s @t @s
on @ Xk :
Finally, rewriting this equation in terms of a gives
@C 1 @2C T @C 1 @ðCSÞ þ jU C; ¼ 2 @t Pesa ðtÞ @ a2 sa ðtÞ @ a sa ðtÞ @ a
on @ Xk :
ð7Þ
3. The Sherman–Lauricella integral equations Following the work in [16,17], we reformulate the Stokes equations as integral equations based on complex variable theory for the biharmonic equation. The integral equations are essentially identical to those for the clean-flow problem, with minor modification to take into account variable surface tension. Below is a summary of the necessary details, and we refer the reader to [21,23] for the development of this complex variable theory, and to [16,17] for the application to the problems at hand. The two-dimensional Stokes Eq. (1) can be simplified by introducing a stream function W(x, y) which satisfies the relations u = Wy, v = Wx. The Stokes equations are then replaced by the biharmonic equation for W(x, y), D2 W(x, y) = 0. Any plane biharmonic function W(x, y) can be expressed by Goursat’s formula as
Wðx; yÞ ¼ Reðz/ðzÞ þ w0 ðzÞÞ; where / and w are analytic functions of the complex variable z = x + i y, and Re (f) denotes the real part of the complex-valued function f. All vectors, n, s, and ng are replaced by their counterparts in complex variables n, s, and ng, respectively. The functions /(z) and w(z) are known as Goursat functions, and all of the physical variables can be expressed in terms of these analytic functions. For completeness, we list them here (c.f. [19]):
iðu þ iv Þ ¼ /ðzÞ þ z/0 ðzÞ þ wðzÞ; f þ ip ¼ 4/0 ; E11 þ iE12 ¼ E22 þ iE21 ¼ iðz/00 þ w0 Þ;
4470
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
where f = DW is the vorticity. In addition, for a point s 2 oXk, the outward normal in the complex plane is given by n = iss. An expression for the stress acting on the interface is
ðp pk Þn þ 2En 2
@ p / s/0 w þ k s : @s 2
Substituting the above and jn = sss, s = ss into the interface condition (2) yields
2
@ p @; ss Þ ¼ ðrss Þ iBReðng s ss Þ; / z/0 w þ k s ¼ rsss þ rs ss iB Reðng s @s @s 2
s 2 @ Xk :
Integrating the above with respect to s gives our final stress condition on the boundary:
þ pk s ¼ 1 rss þ i B / z/0 w 2 2 2
Z
s
ss Þds: Reðng s
ð8Þ
0
We note for future reference that the net force Fk acting on oXk is related to the jump of the right-hand side of (8) about oXk. Specifically, this value is
F k ¼ Bng Ak ;
ð9Þ
where Ak is the area bounded by oXk (c.f. [17] for more details). The suitable far-field expressions for the Goursat functions are
i 1 / /1 ðzÞ p1 ðtÞz þ GðtÞ; 4 2
1 w w1 ðzÞ GðtÞ; 2
as z ! 1;
ð10Þ
where the functions p1(t) and G(t) are determined as part of the solution (G(t) must be included to ensure that / and w can be uniquely described). The functions /1 and w1 are obtained from the far-field velocity:
/1 ðzÞ þ z/01 ðzÞ þ w1 ðzÞ ¼ iðu1 þ iv 1 Þ: The integral equations are derived by finding suitable representations for /(z) and w(z). This is discussed extensively in [17], and we summarize these representations here:
/ðzÞ ¼
1 2p i
Z @X
xðn; tÞ nz
i dn p1 ðtÞz þ /1 ðzÞ þ GðtÞ 4
M Bng X Am logðz zm Þ; 8p m¼1 Z Z nxðn; tÞ 1 xðn; tÞdn þ xðn; tÞdn 1 dn þ w1 ðzÞ wðzÞ ¼ 2p i @ X nz 2pi @ X ðn zÞ2 M B X zm : GðtÞ þ i ng Am logðz zm Þ þ ng Am 8p m¼1 z zm
i
In the above, logarithmic singularities have been placed at the centres zk of each oXk. These singularities are Stokeslet singularities, and their strength is related to the total force acting on each bubble which, in turn, is given by (9). If B – 0, Stokes paradox results in a logarithmic growth of the velocity field. In order to eliminate the consequences of Stokes paradox, we follow the same approach that was taken in [10,15,17,24] and assuming that the bubbles are embedded in a half space S. By convention, we assume that S is the upper half plane so that a solid wall is placed at y = 0, and on this wall, the fluid velocity is zero. By changing the direction of gravity, ng, it is trivial to change the relative position of the wall. The boundary condition on the wall is satisfied by using the method of images and formulae for the reflected sources are given in [15]. Effectively, the reflected Stokeslets act to cancel out the logarithmic growth of the velocity field. Substituting these representations into (8), letting z tend to a point s on the contour oX and using the classical formulae for the limiting values of Cauchy-type integrals, we obtain the Sherman–Lauricella integral equation,
Z Z 1 ns 1 ns i ðn; tÞd ~k xðn; tÞd ln þ x ðp pk Þs a n s 2p i @ X 2 1 2pi @ X ns Z a rðCÞ @ s B sa0 Þda0 ¼ Reðng s /1 ðsÞ þ s/01 ðsÞ þ w1 ðsÞ þ i 2 @s 2 0 M B X s zm þ 2ng Am argðs zm Þ þ ing ; s 2 @ Xk : 8p m¼1 s zm
xðs; tÞ þ
ð11Þ
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4471
~k ðtÞ ak ðtÞ 2GðtÞ; p1 ðtÞ, and pk(t) are determined as part of the solution through Here, a
~k ¼ a
Z
xðn; tÞds;
@ Xk
p1 ðtÞ ¼
2
p
Re
Z @ X1
2
pk ðtÞ ¼ p1 ðtÞ þ
p
Re
xðn; tÞ dn; ðn z1 Þ2 Z xðn; tÞ ðn zk Þ2
@ Xk
dn:
4. Dynamics on the interface We now consider the motion of the interface. In order to reduce notational clutter, we assume below that s 2 oXk, and that all variables and parameters are assumed to reside on the component curve under consideration. Much of what follows is discussed in depth in [16,17] with minor modifications required to include variable surface tension. Rewriting (3) in terms of complex variables yields
ds sa sa ¼ Ui þ T ; dt sa sa
s 2 @ Xk ;
ð12Þ
g and where U ¼ Refðu þ iv Þn
u þ iv ¼
1 2p
I
xðn; tÞ
@X
Z dn dn 1 ns þ þ u1 þ iv 1 þ ðu þ iv Þs : þ xðn; tÞd ns ns 2p @ X ns
ð13Þ
H In the above, x(s, t) is found from the solution to (11), and denotes a principal value integral. If gravity effects are included, s it is understood that (u + i v) takes on the value from the appropriate singular contributions. Following the work of [13,16], a small-scale decomposition of (3) yields
ds 1 H½rðCÞsa ; 2 sa ðtÞ dt where H is the Hilbert transform. For explicit time-stepping methods, this will typically lead to a low-order stability constraint, Dt = O(saDa); however, it is clear that problems will arise when surface tension becomes small, which may happen in regions of large surfactant concentration. The diffusion term in (7) leads to a higher-order stability constraint, which is diagonalizable by the Fourier transform,
b ðkÞ @C k b b ðkÞ þ RðkÞ; C ¼ @t Pes2a ðtÞ 2
ð14Þ
^ where R contains the non-diffusive terms in (7), and RðkÞ is its Fourier transform. A similar approach was taken in [6] (albeit with a finite difference discretization) to study the effects of surfactant on capillary waves; to our knowledge it has not been applied to surfactant-laden interface motion in a Stokes flow. Since the surfacant is insoluble, the total amount of surfacant on the interface must remain fixed,
d dt
Z 2p
Cða; tÞsa ðtÞda ¼ 0:
ð15Þ
0
This is a constraint that is used to check the accuracy of our numerical solutions. 5. Numerical methods Using the methods outlined in [13,16], we obtain a reparametrization of each component curve, s 2 @ Xk, so that we achieve an initial equal arc length distribution of Nk marker points. We again label a 2 [0, 2p) as the parameter; s is uniformly discretized in a, where hk = 2p/Nk is the mesh spacing. All derivatives with respect to a are calculated spectrally via FFTs, and Nk is increased/decreased so as to maintain spectral resolution. At every time step, then, we are required to solve for N comP plex-valued unknowns, where N ¼ M 1 N k . Below, we discuss the time integration schemes and in Section 5.2 we outline the spatial discretization. 5.1. Time integration schemes The system of evolution equations can be written as
" # 0 E @ s 1 2 þ ; ¼ @t C Pe s12a @@ aC2 R
4472
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
where E is the right-hand side of (12). We use IMEX (Implicit – Explicit) Runge–Kutta schemes that treat the first term on the right-hand side in the above implicitly, and the remainder explicitly. This is a natural choice, since the methods are self-starting and we anticipate that the diffusion term will be stiff. For most examples, we use the second-order Midpoint method found in [2] and which we denote by ‘‘RK2’’. For comparison, we also use a third-order scheme, which we call ‘‘RK3’’. More details on RK2 and RK3 can be found in [2]. Let Dt be the time step size and (.)n be the approximate solution at time t = nDt. To implement RK2, then, we first obtain intermediate values
Dt n E ; 2 " # 2 Dt k bn; b ¼ C b n þ Dt R C 1þ k 2 Peðsa Þ2 k 2 k
s ¼ sn þ
ð16Þ
where
sa ¼
1 2p
Z 2p 0
jsa jda:
Final values are obtained through
snþ1 ¼ sn þ DtE ; Cnþ1 ¼ Cnþ1 þ
@ 2 C þ DtR : Peðsa Þ2 @ a2
Dt
ð17Þ
Aliasing errors are suppressed using a high-order Krasny filter (c.f. [16,17] for a more careful discussion on this point). 5.2. Discretization of the integral equation In order to solve (11), we use the Nyström discretization based on the trapezoidal rule since it achieves super-algebraic convergence for smooth data on smooth boundaries. Associated with each such point aj = (j 1)hk, j = 1, . . . , Nk on oXk, are the discrete variables skj , xkj and Ckj . The derivative values sa at aj will be denoted by mj. The discretization of the Sherman–Lauricella Eq. (11) is
xkj þ
Nm M X X
m K 1 ðskj ; sm i Þ xi þ
m¼1 i¼1
Nm M X X
m K 2 ðskj ; sm i Þxi þ hk
m¼1 i¼1
Nk X i¼1
xki þ i
Nk hk skj X
p
i¼1
" Re
xki mki
ðski zk Þ2
# ¼ g kj ;
ð18Þ
where
g kj
Z r Ckj k B jhk k k 0 k k ¼ m / ð s Þ þ s / ð s Þ þ w ð s Þ þ i Reðng sk ska0 Þda0 1 1 j j 1 j j 2 0 2jmkj j j ( ) M skj zm B X k þ 2ng Am argðsj zm Þ þ ing k : 8p m¼1 sj zm
The discrete kernels K1 and K2 are given by
0 1 mmi A hm @ mm i ; K1 s s ¼ k 2pi sm smi skj i sj 0 1 m k m m m s s i j i m h C m B i K 2 ðskj ; sm @ 2 A: i Þ ¼ 2pi sm sk smi skj i j
When
k j;
m i
skj ¼ smi , K1 and K2 are replaced by the appropriate limits
hk k k j jm j; 2p j j 2 h mkj k k k k jj k ; K 2 sj ; sj ¼ 2p jmj j
K1
skj ; skj ¼
where jkj denotes the curvature at the point skj . As was done in [16,17], the linear Eq. (18) are solved iteratively using a FMM-accelerated generalized minimum residual method GMRES [28]. Since the integral equations are Fredholm equations of the second kind, the number of iterations required to solve the system to a fixed precision is bounded; thus the total computational cost per time step is O(N). Since the FFT is used to evaluate derivatives on the interface, the total cost is O(N log N).
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4473
Once the solution xj has been computed, we calculate the velocity at the grid points through discretizing (13): k
ukj þ iv j ¼
hk
p
ðxa Þkj þ
Nm M X X
m K 01 ðskj ; sm i Þ xi þ i
m¼1 i¼1
Nm M X X
s
k k s m K 2 ðskj ; sm i Þxi þ ðu þ iv Þj þ ðu1 þ iv 1 Þj
ð19Þ
m¼1 i¼1
where K2 is defined previously and K 01 is
0 1 mmi A hm @ m m i s s ¼ þ k 2p sm smi skj i sj ( ) ðma Þkj hk Re K 01 skj ; skj ¼ : 2p mkj
K 01
k j;
m i
This quadrature rule is spectrally accurate and is based on subtracting off the singularity in the principal value integral. Again, FMM is used to evaluate (19). CLOSE CONTACT: As we will show in the following section, the above techniques are both efficient and accurate when the component curves @ Xk are well separated from each other. However, when distinct interfaces come into close contact, the integral operators in (11) become close to singular. As discussed in [11,12], the kernels can become under-resolved numerically even when the density and, in our case, the geometry of the interface are well resolved. Therefore, in addition to requiring an increase in Nk as the shape of the interface becomes more complex, we need to detect when component curves come into close contact and adapt accordingly. In [11,12], an adaptive Gaussian quadrature was implemented to locally refine the mesh in regions of close contact between distinct component curves. We have found that in the context of moving interfaces, it is more efficient to adopt the procedure discussed in [15], which is to use the trapezoid rule but increase the number of discretization points sufficiently so that the kernels of the integral operators remain well resolved. We now summarize this procedure. Let dkm measure the distance of closest approach between bubble oXk and oXm. As discussed in [15], in order to adequately resolve the interactions between these two bubbles, the following heuristic relationship between the mesh spacing and dkm should be satisfied:
maxðdsk ; dsm Þ <
1 dkm ; 4
ð20Þ
where ds is the mesh spacing measured in arc length, dsk = hk sa. At every time step, Nk is increased to ensure that the above is satisfied. Care, of course, must be taken in order to ensure a particular configuration of bubbles satisfies (20). Finding dkm between all particles by directly comparing the distances between all points would require O(N2) operations. Instead, we first embed the configuration of bubbles in a hierarchy of grid refinement. The first level of refinement is a single box which contains all of the bubbles. The next level in the hierarchy is achieved by subdividing each box in the current level into four regions; this subdivision continues until the width of a box at the finest level of refinement is just larger than 4maxkdsk. All of the discretization points are then assigned to the boxes at this finest level. Determining whether or not (20) is satisfied reduces to checking this constraint in each non-empty fine box containing discretization points from different component curves. The computational cost for this is O(N). 6. Semi-analytical solutions and method verification In [33], complex-variable theory for the biharmonic equation is exploited in order to derive analytical solutions for a class of polynomial-shaped time-evolving bubbles with constant surface tension. In [29,30], Siegel modifies these solutions to include surfactant effects. In this section, we summarize these semi-analytic solutions and use them to check the accuracy of our methods. We focus on the simple case of an evolving elliptical bubble in a pure strain flow of the form u ? (Qx, Qy) as jzj ? 1. 0 0 (Here, Q is the dimensionless strain rate, Q = Q l/r0, where Q is the dimensioned strain rate. The initial bubble shape is described by
zð#; 0Þ ¼ að0Þei# þ bð0Þei# ;
# 2 ½0; 2p;
where a(0) and b(0) are real and positive. The bubble area remains constant, A(t) = p[a2(t) b2(t)] A(0). As argued in [29], the bubble will evolve according to
zð#; tÞ ¼ aðtÞei# þ bðtÞei# ;
# 2 ½0; 2p:
The evolution of the parameters a(t) and b(t) are found through
dðabÞ ¼ 2abI0 ða; bÞ þ Qa2 ; dt
ð21Þ
4474
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
where
I0 ða; bÞ ¼
1 4p
Z 2p
rðCð#0 ; tÞÞ
0
jz#0 ð#0 ; tÞj
d#0 :
Eq. (21) is coupled to the surfactant concentration through the integral I0. In this framework, the convection–diffusion equation for the surfactant becomes:
@C 1 @ C# dz C# 1 @ 1 z## þ Re ImðPÞ; ReðPÞ þ ¼ Im Pejz# j @# jz# j dt z# jz# j @# jz# j @t z#
ð22Þ
where
Pð#; tÞ ¼
u z# C; jz# j
and
2 3 0 I2p dz rðCð#0 ; tÞÞ # # 05 2i# 4 1 cot d# u ¼ z# e dt 4p 2 jz# ð#0 ; tÞj 0 dz e2i# rðCÞ H : ¼ þ z# dt jz# j 2i The evolution of the surfactant coated bubble requires numerically integrating the evolution Eqs. (21) and (22) given an initial bubble profile and distribution of surfactant, and for given values of b, Q and Pe. The surfactant concentration is represented spectrally with enough points to ensure full resolution, and we use a high-order Runge–Kutta method [4] for the time integration. Examples of initially circular bubbles in straining flow with two different strain rates are shown in Figs. 2 and 3. The lower strain rate, Q = 0.15, gives rise to a steady state, while a steady state does not exist for the higher strain rate (c.f.
Fig. 2. Bubble profile and surfactant concentration for Siegel’s semi-analytic solutions to a bubble in a straining flow, Q = 0.15, Pe = 103, and b = 0.1.
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4475
Fig. 3. Bubble profile and surfactant concentration for Siegel’s semi-analytic solutions to a bubble in a straining flow, Q = 0.30, Pe = 103, and b = 0.1.
[29,30] for a more thorough discussion on the phenomena resulting from different parameter values). We use Siegel’s semianalytic solutions to check the temporal errors in our numerical methods and the results are shown in Fig. 4. In these simulations, N is large enough to ensure full spatial resolution. Two time integrators are used: RK2 and RK3 (c.f. Section 5.1), and the errors in the total surfactant concentration and in the bubble geometry are shown in Fig. 4. With a time step Dt = 0.01, RK2 keeps all errors to O(105) or below. Furthermore, errors in mass and surfactant conservation remain on or below O(104). For RK3 and with an initial time step of 0.005, the errors are O(107), with mass and surfactant being conserved to, at worst, O(106)%. We also use this example to show the linear scaling of the computational costs with N in Table 1. 7. Numerical results In this section we illustrate the strengths of our numerical methods with a variety of examples demonstrating different phenomenological effects. In all examples, Nk is kept sufficiently large to ensure full spatial resolution. The algorithms have been implemented in Fortran and the tolerance for GMRES is set to at most 1010. As a measure of error, we track any loss in area or surfactant concentration. Unless otherwise stated, these errors are below 103. All timings cited are for a single processor on a Mac Pro with a quad-core Intel Xeon processor. 7.1. Bubbles in a linear straining flow It has been widely studied, both experimentally and mathematically (e.g. in [34,29]), that for a range of parameters b, Pe, and Q, a bubble placed in a linear strain flow will reach an elliptical steady-state profile. In this steady-state regime, it is possible for the surfactant to be entirely swept to the bubble tips and eventually form stagnant surfactant caps. Examples of steady surfactant-covered bubbles at two different strain rates are shown in Fig. 5. Fig. 6 shows an example of two bubbles in a linear straining flow. The two bubbles are initially placed so that the are separated by a distance of 0.4. As the bubbles evolve, a dimpling region forms between the two bubbles. As shown in Fig. 7, the size of this dimpled region is influence by surfactant. It is interesting to note that as b increases, not only does the height of
4476
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
Fig. 4. Errors for the RK2 (dashed lines) and RK3 (solid lines). A and B are the ellipses’ major and minor axes, respectively, and jmax is the maximum curvature in absolute value. The figure on the left shows cumulated errors for Q = 0.15 (N = 256, Dt = 0.01, CPU time is 197 s for RK2, 301 s for RK3). The figure on the right is for Q = 0.3 (N = 2048, Dt = 0.005, CPU time is 1265 s for RK2, 1895 s for RK3). In both cases, Pe = 103, b = 0.1.
Table 1 The CPU time in seconds required to compute the bubble profile and surfactant concentration for Siegel’s semi-analytic solutions for a bubble in a straining flow, Q = 0.15, with Pe = 103 and b = 0.1, and for increasing values of N. In each, Dt = 0.01, and the final time was t = 2.0. N
CPU time
128 256 512 1024 2048
7.16 12.76 22.06 39.45 74.26
the dimpling region increase, but so does the minimum separation between the two interfaces. A summary of this effect is shown in Table 2. 7.2. Bubbles in a cubic straining flow In an experimental setup with a bubble placed in the middle of a four roller mill, G.I. Taylor witnessed that at a critical rotation rate of the rollers, the bubble achieved pointed ends [34]. This phenomenon was further investigated by De Bruijn
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4477
Fig. 5. Stable surfactant caps for Q = 0.15 and Q = 0.2, Pe = 103, b = 0.1. The surfactant caps are shown clearly at the bubble tips.
Fig. 6. Two bubbles in a linear straining flow, Q = 0.2, at t = 0, t = 20, and t = 60. The plot on the left is the clean flow example (CPU time to t = 20 was 52.2 min, to t = 60 was 58 h, N = 8192 at t = 60), and on the right, Pe = 103, b = 0.15 (CPU time to t = 20 was 12.7 min, to t = 60 was 186.0 min, N = 2048 at t = 60).
[7], who showed that jets or small bubbles may be emitted from these pointed tips in a process called tip streaming. De Bruijn’s experiments suggest that tip streaming occurs when gradients in surface tension develop.
4478
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
Fig. 7. Two bubbles in a linear straining flow, Q = 0.2, at t = 20 for b = 0, 0.05, 0.10, 0.15 (least through most deformed, respectively).
Table 2 Separation distances for the example of two bubbles in a shear flow; Dymin measures the minimum distance between the two interfaces near the tips and Dyd measures the height of the dimpled region.
Dymin
b
0.00 0.05 0.10 0.15
D yd
t = 20
t = 60
t = 20
t = 60
0.029 0.110 0.127 0.137
0.014 0.032 0.040 0.053
0.368 0.398 0.472 0.514
0.463 0.614 0.643 0.659
Antanovskii [1] showed that the four-roller-mill plane flow employed by Taylor is described, to leading order, by the following far-field velocity and pressure fields:
u1 ¼ Q ðx; yÞ þ Q ðx3 þ 3xy2 ; 3x2 y y3 Þ;
P1 ¼ p1 þ 6Q ðx2 y2 Þ:
These correspond to the following far-field expressions for the Goursat functions (See Table 3)
i /1 ¼ Q z3 ; 2
w1 ¼ iQz:
Hear, Q is the linear nondimensional strain rate, and Q is a low-order, additional, rotational strain rate. Due to the presence of a cubic term in /1, this type of flow is often called cubic-strain flow. The stability of the resulting spindle-like bubble shapes are studied in [25,30]. For a range of parameters the surface tension force is inadequate to balance the outward motion of the bubble, hence the profile becomes unstable. Here we present numerical examples to demonstrate this. The response to different values of b for a single strain rate Q = 0.3 is shown in Fig. 8. This figure also shows the corresponding surface tension on the interface. (To make a better comparison, we follow
Table 3 The CPU time in minutes and the final value of N (N = 1024 initially) required to compute the cubic straining flow examples up to t = 2.5. Here, Q = 0.3, Pe = 103. In each, the initial time step is Dt = 0.005. This time step is cut in have each time N is doubled (N is increased so as to maintain spectral resolution). b
maxN
CPU time
0.00 0.05 0.10 0.15
2048 4096 8192 32768
2.6 5.9 18.9 406.3
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4479
Fig. 8. On the left, the evolution up to t = 2.5 of initially circular bubbles with different values of b (Pe = 103), in a cubic strain flow with Q = 0.3, = 0.1. On the right, we show the rescaled surface tension, which is lowered considerably at the bubble tip as b increases.
[25] and rescale the non-dimensional surface tension by r = (1 bC)/ (1 b). This ensures that for all values of b, r(1) = 1.) The surfactant has accumulated at the bubble tips, lowering the surface tension there.
4480
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
Further details of the time evolution of different quantities are shown in Fig. 9. In this figure, we plot the evolution of the deformation of the bubble, defined by:
D¼
Rmax Rmin ; Rmax þ Rmin
where Rmin, Rmax are the minimum and maximum values of the bubble radius, respectively, as well as maximum normal velocity of the fluid, the curvature at the tip, and the ratio of maximum to minimum surfactant concentration. Finally, a comparison of the tip at different values of b is shown in Fig. 10. The above examples demonstrate that our numerical methods can track the time-evolution of such shapes far into cusp formation (we have accurately tracked increases in tip curvature up to O(103)). A detailed study of the stability of the cusped bubbles for different parameters Q, Pe, b, and , or the role of a linear versus nonlinear equation of state is not in the scope of this paper, but may be considered in future work.
7.3. Bubbles in a shear flow In this next example, we consider 5 bubbles with an initial configuration shown in Fig. 11 in a shear flow, u ? (0, Cay) as x ? 1, where Ca is the dimensionless shear rate, or Capillary number. (The capillary number expresses the ratio of viscous stresses due to the shear flow and surface tension.) We use this example to demonstrate the capabilities of our method to simulate a larger number of closely interacting interfaces. Fig. 12 shows the interacting bubbles at two separate times; at t = 10, the 5 bubbles are tightly clustered together and at t = 20, they have become well separated. It is interesting to note that more points are required to spatially resolve the clean flow situation, the likely reason being that the bubbles come into closer contact than when surfactant is present. A very similar example was presented in [37] using a level set method. These authors report that their computations break down soon after the time when tight clustering occurs, which they believe to be a result of inadequate spatial resolution in the near-contact regions.
Fig. 9. The above plots show the evolution of the deformation, fluid normal velocity at the tip, the maximum curvature and ratio of maximum to minimum surfactant concentration for different values of b. The clean case b = 0 appears to be approaching steady state.
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4481
Fig. 10. A close up of the tip region for the bubbles in a cubic strain flow at t = 2.6. The curvature values at the tips are 23.34, 56.27, 183.74, and 2221.07.
Fig. 11. Initial bubble configuration for the shear flow example.
7.4. Buoyancy effects In this example, we consider the buoyancy-driven motion of a bubble under the influence of surfactant. In particular, we are interested in the case of a bubble falling towards a plane wall. This example was considered in [24], but without the effects of surfactant. The rising case (with surfactants) was studied in [14]. As indicated in [24], these falling bubbles can become susceptible to gravitational Rayleigh–Taylor instabilities. Depending on the value of B, as the bubble falls close to the wall, it may start to spread and may develop a dimple at the midpoint on the leading edge. A layer of fluid can become trapped between the interface and the wall, and the dimple may continue to grow. It would appear that for larger values of B, the situation is more susceptible to instabilities. The role surfactant might play remains in question.
4482
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
Fig. 12. Bubbles in a shear flow with Ca = 0.6. On the left, a comparison of the clean flow (solid black line) versus the presence of surfactant (Pe = 103, b = 0.10) at time t = 10. The number of points on each of the clean bubbles is 2048 by the end of the simulation and 1024 on the dirty bubbles. The plot on the rights shows the bubbles at t = 20 (total CPU time was 17 h).
Fig. 13. An initially circular bubble with centre (0, 5) is released and falls towards a plane wall with B ¼ 1. The clean interface is shown in black. For the bubbles with surfactant, Pe = 103. The bubble lagging furthest behind has b = 0.4, the other has b = 0.2.
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4483
Fig. 14. An initially circular bubble with centre (0, 5) is released and falls towards a plane wall with B ¼ 10. The clean interface is shown in black. For the bubble with surfactant, Pe = 103, b = 0.2.
Fig. 15. Velocity and vorticity plots at t = 1.5, B ¼ 10.
The results of simulations at two different values of B are shown in Figs. 13 and 14. For B ¼ 1, the evolution of the clean bubble has reached a near steady state by t = 100. As the bubble becomes close to the wall, it develops a slight dimple. As the dirty bubbles fall, surfactant is swept away from the leading edge towards the trailing edge. Surfactant acts to retard the drop slightly; the clean drop reaches the wall first, while the bubble with the highest value of b is last. Interestingly enough, the
4484
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
Fig. 16. Velocity and vorticity plots at t = 5.92, B ¼ 10.
Fig. 17. A close up of the dimple region: velocity and vorticity plots at t = 5.92, B ¼ 10.
Fig. 18. Initial configuration of 144 bubbles.
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4485
Fig. 19. Final configuration of 144 bubbles at t = 1.5. Initially, Nk = 64 on all of the bubbles (i.e. the total number of points is 9216). In the clean flow case (the solid curves), the number of points on the bubbles in the final configuration ranged from 128 to 1024 for a total of 108032 points, and the computation required almost 19 h. For b = 0.1, Pe = 103 (the dashed curves), the number of points on the bubbles in the final configuration ranged from 128 to 4096 for a total of 141056 points, and the computation required almost 32 h. For b = 0.2, Pe = 103 (the dotted curves), the number of points on the bubbles in the final configuration ranged from 128 to 4096 for a total of 190848 points, and the computation required almost 84 h. In all cases, the area of each bubble is conserved up to O(104).
Fig. 20. The surfactant concentration on one of the 144 rising bubbles at t = 1.5. The clean interface is shown in black. The most deformed interface occurs when b = 0.2.
b = 0.4 bubble forms a dimple before the b = 0.2 bubble. We did not evolve the dirty bubbles further to determine if they had reached steady state or not. For B ¼ 10, we were unable to evolve the clean bubble much past t = 5.92 before the calculation became unstable. By all appearances, this instability is a physical one, not a numerical one, but this requires more careful investigation. We were able to evolve the dirty bubble with b = 0.2 slightly further in time, but again, the calculation becomes unstable. As can be seen in Fig. 14, at the earlier stages in the evolution, there is a small, but visible, difference on the trailing edges between the clean and dirty bubbles. (Velocity and vorticity plots of the flow around the dirty bubble at t = 1.5 is shown in Fig. 15.) This difference diminishes at later times. Velocity and vorticity plots of the flow around the dirty bubble at t = 5.92 is shown in Fig. 16, with a close up of the dimpling region in Fig. 17. The velocity plots at both times shows a recirculation zone just trailing the bubble tip. This was also observed in [24].
4486
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
We finish off our examples with a larger-scale simulation of 144 rising bubbles. The initial distribution of bubbles is shown in Fig. 18. During the course of the simulation, the number of points on individual bubbles increases due to deformation and/or close approach to other bubbles in the assemblage. We computed the solution up to t = 1.5 (see Fig. 19 for a comparison of the clean flow with b = 0.1 and b = 0.2.) Fig. 20 shows a more detailed comparison of one of the bubbles in the assemblage.
8. Conclusion We have presented efficient and highly-accurate numerical methods to compute the motion of multiple surfactantcoated, two-dimensional interfaces in a slow viscous flow. These methods have the ability to maintain a high degree of accuracy in space and time while remaining highly efficient. This efficiency is measured by the low operation count per time step (O(N) or O(N log N), where N is the number of spatial nodes in the discretization), and the low-order stability constraint. A spectral mesh is used to represent the interface, and by working in the equal arc length frame, the equations become amenable to high-order semi-implicit time-stepping methods. We have demonstrated the efficacy of our techniques on a number of phenomenologically different examples, including some relatively large-scale examples that were computed using only modest computational resources. In a number of cases (such as flows including gravitational effects or bubbles in a cubic straining flow), we have indicated that a more in depth and careful parametric investigation may be warranted. We believe we have developed methods that will serve as a useful benchmark for other, less accurate methods which are aimed at more complicated physical models. Immediate future work will extend these methods to three-dimensions and to consider drops containing fluid of different viscosity ratios. References [1] L.K. Antanovskii, Formation of a pointed drop in Taylor’s four-roller mill, J. Fluid Mech. 327 (1996) 325–342. [2] U.M. Ascher, S.J. Ruuth, R.J. Spiteri, Implicit-explicit Runge–Kutta methods for time-dependent partial differential equations, Appl. Numer. Math. 25 (1997) 151–167. [3] I. Bazhlekov, P.D. Anderson, H.E.H. Meijer, Numerical investigation of the effect of insoluble surfactants on drop deformation and breakup in simple shear flow, J. Colloid Interf. Sci. 298 (2006) 369–394. [4] R.W. Brankin, I. Gladwell, L.F. Shampine, RKSUITE release 1.0, 1991.
. [5] J. Carrier, L. Greengard, V. Rokhlin, A fast adaptive multipole algorithm for particle simulations, SIAM J. Sci. Stat. Comput. 9 (1988) 669–686. [6] H.D. Ceniceros, The effects of surfactants on the formation and evolution of capillary waves, Phys. Fluids 15 (2003) 245–256. [7] R.A. De Bruijn, Tipstreaming of drops in simple shear flows, Chem. Eng. Sci. 48 (1993) 277–284. [8] M.A. Drumright-Clarke, Y. Renardy, The effect of insoluble surfactant at dilute concentration on drop breakup under shear with inertia, Phys. Fluids 16 (2004) 14–21. [9] C.D. Eggleton, T.M. Tsai, K.J. Stebe, Tip streaming from a drop in the presence of surfactants, Phys. Rev. Lett. 87 (4) (2001) 048302. [10] L. Greengard, M.C. Kropinski, A. Mayo, Integral equation methods for Stokes flow and isotropic elasticity in the plane, J. Comput. Phys. 125 (1996) 403– 414. [11] Johan Helsing, Thin bridges in isotropic electrostatics, J. Comput. Phys. 127 (1996) 142–151. [12] Johan Helsing, Leslie Greengard, On the numerical evaluation of elastostatic fields in locally isotropic two-dimensional composites, J. Mech. Phys. Solids 46 (8) (1998) 1441–1462. [13] T.Y. Hou, J.S. Lowengrub, M.J. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys. 114 (1994) 312–338. [14] R.A. Johnson, A. Borhan, Stability of the shape of a surfactant-laden drop translating at low Reynolds number, Phys. Fluids 12 (2000) 773–784. [15] M.C.A. Kropinski, Integral equation methods for particle simulations in creeping flows, Comput. Math. Appl. 38 (1999) 67–87. [16] M.C.A. Kropinski, An efficient numerical method for studying interfacial motion in two-dimensional creeping flows, J. Comput. Phys. 171 (2001) 479– 508. [17] M.C.A. Kropinski, Numerical methods for multiple inviscid interfaces in creeping flows, J. Comput. Phys. 180 (2002) 1–24. [18] M.C. Lai, Y.H. Tseng, H. Huang, An immersed boundary method for interfacial flows with insoluble surfactant, J. Comput. Phys. 227 (2008) 7279–7293. [19] W.E. Langlois, Slow Viscous Flow, Macmillan Co., New York NY, 1964. [20] X. Li, C. Pozrikidis, Effect of surfactants on drop deformation and on the rheology of dilute emulsions in Stokes flow, J. Fluid Mech. 341 (1997) 165–194. [21] S.G. Mikhlin, Integral Equations, Pergamon, London, 1957. [22] M. Muradoglu, G. Tryggvason, A front-tracking method for computation of interfacial flows with soluble surfactants, J. Comput. Phys. 227 (2008) 2238– 2262. [23] N.I. Muskhelishveli, Some Basic Problems of the Mathematical Theory of Elasticity, Noordhoff, Groningen, 1953. [24] C. Pozrikidis, The deformation of a liquid drop moving normal to a plane wall, J. Fluid Mech. 215 (1990) 331–363. [25] C. Pozrikidis, Numerical studies of cusp formation at fluid interfaces in Stokes flow, J. Fluid Mech. 357 (1998) 29–57. [26] C. Pozrikidis, Numerical investigation of the effect of surfactants on the stability and rheology of emulsions and foam, J. Eng. Math. 41 (2001) 2637– 2658. [27] V. Rokhlin, Rapid solution of integral equations of classical potential theory, J. Comput. Phys. 60 (1985) 187–207. [28] Y. Saad, M.H. Schultz, GMRES: a generalized minimum residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. [29] M. Siegel, Influence of surfactant on rounded and pointed bubbles in two-dimensional Stokes flow, SIAM J. Appl. Math. 59 (1999) 1998–2027. [30] M. Siegel, Cusp formation for time-evolving bubbles in two-dimensional Stokes flow, J. Fluid Mech. 412 (2000) 227–257. [31] M. Siegel, M. Booty, Steady deformation and tip-streaming of a slender bubble with surfactant in an extensional flow, J. Fluid Mech. 544 (2005) 243– 275. [32] H. Stone, A simple derivation of the time-dependent convection–diffusion equation for surfactant transport along a deforming interface, Phys. Fluids A 2 (1990) 111–112. [33] S. Tanveer, G.L. Vasconcelos, Time-evolving bubbles in two-dimensional Stokes flow, J. Fluid Mech. 301 (1995) 325–344. [34] G.I. Taylor, The formation of emulsions in definable fields of flow, Proc. R. Soc. Lond. A 146 (1934) 501–523. [35] A. K Tornberg, L. Greengard, A fast multipole method for the three-dimensional Stokes equations, J. Comput. Phys. 227 (2008) 1613–1619. [36] H. Wong, D. Rumschitzki, C. Maldarelli, On the surfactant mass balance at a deforming fluid interface, Phys. Fluids A 8 (1996) 3203–3204. [37] J.J. Xu, J. Lowengrub, H.K. Zhao, A level-set method for interfacial flows with surfactant, J. Comput. Phys. 212 (2006) 590–616.
M.C.A. Kropinski, E. Lushi / Journal of Computational Physics 230 (2011) 4466–4487
4487
[38] X. Yang, A.J. James, J. Lowengrub, X. Zheng, V. Cristini, An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids, J. Comput. Phys. 217 (2006) 364–394. [39] S. Yon, C. Pozrikidis, A finite-volume/boundary element method for flow past interfaces in the presence of surfactants, with applications to shear flow past a viscous drop, Comput. Fluids 27 (1998) 879–902.
