An immersed boundary method for interfacial flows with insoluble surfactant

(1)

An immersed boundary method for interfacial ﬂows

with insoluble surfactant

Ming-Chih Lai

a,*

, Yu-Hau Tseng

a

, Huaxiong Huang

b

a

Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan

b

Department of Mathematics and Statistics, York University, Toronto, Ontario, Canada M3J 1P3 Received 15 September 2007; received in revised form 1 February 2008; accepted 16 April 2008

Available online 26 April 2008

Abstract

In this paper, an immersed boundary method is proposed for the simulation of two-dimensional fluid interfaces with insoluble surfactant. The governing equations are written in a usual immersed boundary formulation where a mixture of Eulerian flow and Lagrangian interfacial variables are used and the linkage between these two set of variables is pro-vided by the Dirac delta function. The immersed boundary force comes from the surface tension which is affected by the distribution of surfactant along the interface. By tracking the interface in a Lagrangian manner, a simplified surfactant transport equation is derived. The numerical method involves solving the Navier–Stokes equations on a staggered grid by a semi-implicit pressure increment projection method where the immersed interfacial forces are calculated at the begin-ning of each time step. Once the velocity value and interfacial configurations are obtained, surfactant concentration is updated using the transport equation. In this paper, a new symmetric discretization for the surfactant concentration equa-tion is proposed that ensures the surfactant mass conservaequa-tion numerically. The effect of surfactant on drop deformaequa-tion in a shear flow is investigated in detail.

Keywords: Immersed boundary method; Interfacial ﬂow; Navier–Stokes equations; Surfactant

1. Introduction

In this paper, we propose an immersed boundary method for the simulation of two-dimensional fluid inter-faces with insoluble surfactant. Surfactant are surface active agents that adhere to the fluid interface and affect the interface surface tension. Surfactant play an important role in many applications in the industries of food, cosmetics, oil, etc. For instance, the daily extraction of ore rely on the subtle effects introduced by the presence of surfactant[5]. In a liquid–liquid system, surfactant allow small droplets to be formed and used as an emul-sion. Surfactant also play an important role in water purification and other applications where micro-sized bubbles are generated by lowering the surface tension of the liquid–gas interface. In microsystems with the

* _{Corresponding author. Tel.: +886 3 5731361; fax: +886 3 5724679.}

E-mail addresses:mclai@math.nctu.edu.tw(M.-C. Lai),yhtseng@math.nctu.edu.tw(Y.-H. Tseng),hhuang@yorku.ca(H. Huang). Journal of Computational Physics 227 (2008) 7279–7293

(2)

presence of interfaces, it is extremely important to consider the effect of surfactant since in such cases the cap-illary effect dominates the inertia of the fluids[20].

The immersed boundary (IB) method proposed by Peskin[14], has been applied successfully to blood–valve interaction and other biological problems. The IB formulation employs a mixture of Eulerian and Lagrangian variables, where the immersed boundary is represented by a set of discrete Lagrangian markers embedding in the Eulerian fluid domain. Those markers can be treated as force generators to the fluid while being carried by the fluid motion. The interaction between the Lagrangian force generators (markers) and the fluid motion, described by variables defined on the fixed Eulerian grid, is linked by a properly chosen discretized delta func-tion. Most IB applications in the literature belong to the fluid–structure problems, and they can be found in a recent review of Peskin[15]. However, there is comparatively less work on the application of the IB method to viscous, incompressible multi-phase flow problems. Perhaps the most successful one is the front-tracking method proposed by Tryggvason et al. [21,22] which uses an approach similar to the immersed boundary method.

In the case of interfacial flows with surfactant, Ceniceros [4]used a hybrid level set and front tracking approach to study the effects of surfactant on the formation of capillary waves. Lee and Pozrikids[12]used Peskin’s immersed boundary idea to study the effects of surfactant on the deformation of drops and bubbles in Navier–Stokes flows. The surfactant convection–diffusion equation in these papers is based on the formulation proposed by Wong et al.[23], and the conservation of total mass of surfactant on the interface has not been rigorously investigated numerically.

James and Lowengrub[9] have proposed a surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant. Instead of solving the surfactant concentration equation based on Stone’s derivation [19] directly, the authors relate the surfactant concentration to the ratio of the surfactant mass and surface area so that they are tracked independently. The method has been applied to study the axis-symmetric drop deformation in extensional flows. Recently, Xu et al. [25] develop a level-set method for interfacial Stokes flows with surfactant. Their method couples surfactant transport, solved in an Eulerian domain [26] with Stokes flow field, solved by the immersed interface method [11] with jump conditions across the interface. However, the method does not conserve the mass automatically and numerical scaling is used to enforce the conservation of surfactant on the interface numerically. Recently, Muradoglu and Tryggvason [13] have proposed a front-tracking method for computation of interfacial flows with soluble surfactant. They consider the axis-symmetric motion and deformation of a viscous drop moving in a circu-lar tube.

In this paper, we propose an immersed boundary method to simulate the interfacial problems with insol-uble surfactant. By tracking the interface in a Lagrangian manner, the surfactant concentration equation becomes much simpler than the one in[23]. Our numerical method involves solving the Navier–Stokes equa-tions on a staggered grid by a semi-implicit pressure increment projection method where the immersed inter-facial forces are calculated at the beginning of each time step. A new symmetric discretization for the surfactant concentration equation is proposed so that the total mass of surfactant is conserved numerically. The eﬀect of surfactant on drop deformation in a shear ﬂow is then investigated in detail.

The rest of the paper is organized as follows. In Section 2, we present the governing equations which includes the immersed boundary formulation and the surfactant concentration equation in Lagrangian coor-dinates on the interface. The numerical method is described in Section3which includes an algorithm of solv-ing the Navier–Stokes equations and a conservative scheme for the surfactant equation. The eﬀect of surfactant on drop deformation in a shear ﬂow is investigated numerically in Section 4. Some concluding remarks and brief discussion on future directions are given in Section5.

2. The governing equations

Consider an incompressible two-phase flow problem consisting of fluids 1 and 2 in a fixed two-dimensional square domain X¼ ½a; b ½c; d ¼ X1[ X2where an interface R separates X1 from X2. Here, we assume the

interface is a simple closed curve immersed in the fluid domain, and is contaminated by the surfactant so that the distribution of the surfactant changes the surface tension accordingly. In each fluid region, the Navier– Stokes equations are satisfied as

(3)

qi oui ot þ ðui rÞui ¼ r Tiþ qig; in Xi; ð1Þ r ui¼ 0; in Xi; ð2Þ u¼ ub; inoX; ð3Þ

where for i¼ 1; 2 in each ﬂuid domain, Ti¼ piIþ liðruiþ ruTiÞ is the stress tensor, piis the pressure, uiis

the ﬂuid velocity, q_iis the density, l_i is the viscosity, and g is the gravitational constant. It is well-known that, across the interface R, the velocity is continuous

½u_R¼ uj_R;2 uj_R;1¼ 0 ð4Þ

and the normal stress jump is balanced by the interfacial force F (deﬁned only on R) as

½Tn_Rþ F ¼ 0; ð5Þ

where n is the unit normal vector on R directed towards fluid 2. Since it is not easy to solve the Navier–Stokes equations(1) and (2)in X with jump conditions(4) and (5)on R, especially when the interface is moving. In order to formulate the problem using the immersed boundary approach, we simply treat the interface as an immersed boundary that exerts force F to the fluids and moves with local fluid velocity. In this paper, we con-sider the case of equal viscosity l1¼ l2¼ l, equal density q1¼ q2¼ q, and neglect gravity. However, the

cur-rent formulation can be extended straightforwardly to general two phase flow with diffecur-rent density and viscosity. The present interfacial force term in delta function formulation and the surfactant concentration equation are the same as the single phase problem. The major difference comes from the Navier–Stokes for-mulation and their numerics.

2.1. Immersed boundary formulation

Throughout this paper, the interface R is represented by a parametric form ðX ðs; tÞ; Y ðs; tÞÞ; 0 6 s 6 Lb,

where s is the parameter of the initial conﬁguration of the interface, which is not necessarily the arc-length. Using the non-dimensionalization process[9,25], we can write down our governing equations in the usual immersed boundary formulation as follows.

ou otþ ðu rÞu þ rp ¼ 1 Rer 2_u_þ 1 ReCaf; ð6Þ r u ¼ 0; ð7Þ fðx; tÞ ¼ Z R Fðs; tÞdðx Xðs; tÞÞ ds; ð8Þ oXðs; tÞ ot ¼ uðXðs; tÞ; tÞ ¼ Z X uðx; tÞdðx Xðs; tÞÞ dx; ð9Þ Fðs; tÞ ¼ o osðrðs; tÞsðs; tÞÞ; ð10Þ sðs; tÞ ¼ oX os oX os : ð11Þ

The dimensionless numbers are the Reynolds number (Re) describing the ratio between the inertial force and the viscous force, and the capillary number (Ca) describing the strength of the surface tension. Eqs.(8) and (9)

represent the interaction between the immersed interface and the fluids. In particular, Eq. (8) describes the force (f) acting on the fluid due to the interfacial force (F), which is defined only on the interface and must be balanced by the normal stress as shown in Eq.(5). Here, r is the surface tension, and s is the unit tangent vector on the interface. Eq.(9)states that the interface moves with the fluid velocity which is consistent with

(4). The present formulation employs a mixture of Eulerian (x) and Lagrangian (X) variables which are linked by the two-dimensional Dirac delta function dðxÞ ¼ dðxÞdðyÞ.

(4)

The interfacial force F arises from the surface tension and its form is derived from Laplace–Young condi-tion[10]. One can further take derivatives explicitly so that

Fðs; tÞ ¼ o osðrsÞ ¼ or ossþ r os os¼ or ossþ rjn oX os ; ð12Þ

where j is the curvature of the interface and n is the unit outward normal. The first term on the right-hand side of Eq.(12)is the Marangoni force (the tangential force) and the second one is the capillary force (the normal force). (Note that, we have different sign convention in the capillary term since the sign of curvature is different from that in the literature[8,9,12,25]. For circular interface, the present curvature is negative.) Note also that if the surface tension is a constant, then the force only exerts in the normal direction. However, when the inter-face is contaminated by the surfactant, the distribution of the surfactant changes the surinter-face tension accord-ingly. Generally speaking, the higher the surfactant concentration, the less the surface tension. The relation between surface tension and surfactant concentration can be described by the Langmuir equation of state

[17]. As in[4], the following linear approximation of Langmuir equation is used

rðCÞ ¼ rcð1 bCÞ; ð13Þ

where C is the surfactant concentration, rc is the surface tension of a clean interface, and b satisfying

0 6 b < 1 is a dimensionless number that measures the sensitivity of surface tension to changes in surfactant concentration.

In order to close the system, we still need one more equation for surfactant concentration evolution. As mentioned before, surfactant are insoluble to the buck fluids so they are simply convected and diffused along the interface. Since there is no exchange between the interface and the bulk fluids, the total mass of the sur-factant must be conserved. The equation of sursur-factant concentration is derived in next subsection.

2.2. Surfactant concentration equation

The basic equation for surfactant transport equation along a deforming interface has been derived by Scri-ven[18], Aris[2], and Waxman[24]. All three papers derived the surfactant equation relying heavily on dif-ferential geometry. Stone [19], however, presented a simple derivation of the time-dependent convective– diffusion equation for surfactant transport along a deforming interface. In this subsection, we present a slightly different derivation from Stone for the surfactant transport equation which will be used as one of our governing equations for numerical computation. Our derivation is in the same spirit of the immersed boundary approach. A more detailed derivation for surfactant concentration equation along a two-dimen-sional parametric deforming surface in three-dimentwo-dimen-sional fluid domains can be found in our recent work[7]. Let LðtÞ be an interfacial segment where the surfactant concentration (the mass of the surfactant per unit length) is defined. Since the surfactant remain on the material element and do not transport or diffuse to the surrounding bulk fluids, the mass on the segment is conserved

d dt

Z

LðtÞ

Cðl; tÞ dl ¼ 0; ð14Þ

where dl is the arc-length element. To apply the time derivative more easily, we rewrite the above equation in terms of the initial parameter s as

d dt Z Lð0Þ Cðs; tÞoX os ds ¼ 0: ð15Þ

By taking the time derivative inside the integral, we obtain Z Lð0Þ oC ot oX os þ C o ot oX os ds¼ 0: ð16Þ

Note that, in our present formulation, both the interface and surfactant concentration are tracked in a Lagrangian manner. Thus, the time derivative of the ﬁrst term in Eq.(16)is exactly the material derivative

(5)

of Stone’s derivation[19]. The time derivative of the second term is due to interface stretching. Now we need to compute the rate of the stretching factor, and using Eq.(9), we have

o ot oX os ¼ oX os o os oX ot þoY os o os oY ot joX osj ¼ oX os ou osþ oY os ov os joX osj ¼ oX os ru oX os þoY os rv oX os joX osj ¼ ou os s _oX os ¼ ðrs uÞ oX os : ð17Þ

Here, the notation rs u means the surface divergence which is used commonly in the literature. Since the

material segment is arbitrary, we thus have oC

ot þ ðrs uÞC ¼ 0: ð18Þ

If we allow surfactant diﬀusion along the interface, we obtain the surfactant transport-diﬀusion equation as oC ot þ ðrs uÞC ¼ 1 Pes o os oC os oX os oX os ; ð19Þ

where Pes is the surface Peclet number [9]. We note that surface diﬀusion is also written in terms of initial

parameter s.

Let us summarize this section by pointing out the differences and similarities between our present surfactant equation(19)and the ones derived in the literature[19,23]. As we discussed before, the present time derivative is exactly the material derivative with the material parameter s fixed, while the time derivative used in[19]is keeping the material coordinates X fixed. Wong et al.[23]argued that the time derivative term in Stone’s sur-factant equation causes ambiguity in numerical discretization since the material coordinates is time-dependent as well. Wong et al.[23]provide an alternative derivation for the surfactant equation, where the concentration time derivative is applied by keeping the material parameter s fixed. This is exactly what we have done here. It is interesting (but not surprising) to conclude that the surfactant concentration equation in[23]can be simpli-fied to our present form(19)by substituting Eq.(9) into their equation.

3. Numerical method

In this paper, the fluid flow variables are defined on a staggered marker-and-cell (MAC) mesh introduced by Harlow and Welsh [6]; that is, the pressure is defined on the grid points labelled as x¼ ðxi; yjÞ ¼

ðði 1=2Þh; ðj 1=2ÞhÞ for i; j ¼ 1; 2 . . . ; N , the velocity components u and v are deﬁned at ðxi1=2; yjÞ ¼

ðði 1Þh; ðj 1=2ÞhÞ and ðxi; yj1=2Þ ¼ ðði 1=2Þh; ðj 1ÞhÞ, respectively, where the spacing h ¼ Dx ¼ Dy.

For the immersed interface, we use a collection of discrete points sk¼ kDs; k ¼ 0; 1; . . . M such that the

Lagrangian markers are denoted by Xk ¼ XðskÞ ¼ ðXk; YkÞ. The surfactant concentration Ck, surface tension

rk are deﬁned at the ‘‘half-integer” points given by skþ1=2¼ ðk þ 1=2ÞDs. Without loss of generality, for any

function deﬁned on the interface /ðsÞ, we approximate the partial derivativeo/ os by

Ds/ðsÞ ¼

/ðs þ Ds=2Þ /ðs Ds=2Þ

Ds : ð20Þ

By using this ﬁnite diﬀerence convention, the interface stretching factor can be approximated byj DsXkj, and

thus the unit tangent vector sk are also deﬁned at the ‘‘half-integer” points.

Let Dt be the time step size, and n be the superscript time step index. At the beginning of each time step, e.g., step n, the variables Xn

k ¼ Xðsk; nDtÞ, Cnk ¼ Cðskþ1=2; nDtÞ, u

n_{¼ uðx; nDtÞ, and p}n1=2_{¼ pðx; ðn 1=2ÞDtÞ are all}

given. The details of the numerical time integration are as follows. 1. Compute the surface tension and unit tangent on the interface as

rn k ¼ rcð1 bCnkÞ; ð21Þ sn_k¼ DsX n k jDsXnkj; ð22Þ

(6)

both of which hold for skþ1=2¼ ðk þ 1=2ÞDs. Then we deﬁne the interface force as

Fn_k ¼ Dsðrn_ksn_kÞ; ð23Þ

at point Xk.

2. Distribute the force from the markers to the ﬂuid by fnðxÞ ¼X

k

Fn_kdhðx XnkÞDs; ð24Þ

where the smooth version of Dirac delta function in[15]is used.

3. Solve the Navier–Stokes equations. This can be done by the following second-order accurate projection method[3], where the nonlinear term is approximated by the Adams–Bashforth scheme and the viscous term is approximated by the Crank–Nicholson scheme.

ðu rhÞunþ1=2¼ 3 2ðu n_r hÞun 1 2ðu n1_r hÞun1; ð25Þ u_un Dt þ ðu rhÞu nþ1=2_{¼ rp}n1=2_þ 1 2Rer 2 hðu _{þ u}n_{Þ þ} f n ReCa; ð26Þ u¼ ub; onoX; ð27Þ r2 h/ nþ1 _¼rh u Dt ; o/ on ¼ 0; on oX; ð28Þ unþ1¼ u_Dtr h/nþ1; ð29Þ pnþ1=2¼ pn1=2þ /nþ1rh u 2Re : ð30Þ

Hererhis the standard centered diﬀerence operator on the staggered grid. One can see that the above

Na-vier–Stokes solver involves solving two Helmholtz equations for velocity u _{and one Poisson equation for}

pressure. These elliptic equations are solved using the fast Poisson solver provided by the public software package Fishpack[1].

4. Interpolate the new velocity on the ﬂuid lattice points onto the marker points and move the marker points to new positions. Unþ1_k ¼X x unþ1dhðx XnkÞh 2_; ð31Þ Xnþ1_k ¼ Xn kþ DtU nþ1 k : ð32Þ

5. Update surfactant concentration distribution Cnþ1_k . Since the surfactant is insoluble, the total mass on the interface must be conserved. Thus, it is important to develop a numerical scheme for the surfactant concentration equation to preserve the total mass. This can be done as follows.

Firstly, let us rewrite the surfactant concentration equation(19) by multiplying the stretching factor on the both sides of the equation as

oC ot oX os þ ðrs uÞ oX os C ¼Pe1s o os oC os oX os : ð33Þ

Then substitute Eq.(17)of rate of stretching factor into the above equation, we have oC ot oX os þoto oX os C ¼Pe1s o os oC os oX os : ð34Þ

(7)

Now we discretize the above equation by the Crank–Nicholson scheme in a symmetric way as Cnþ1_k Cn k Dt DsXnþ1k þ DsXnk 2 þ jDsXnþ1k j jDsXnkj Dt Cnþ1_k þ Cn k 2 ¼ 1 2Pes 1 Ds ðCnþ1 kþ1 Cnþ1k Þ=Ds ðjDsXnþ1kþ1j þ jDsXnþ1k jÞ=2 ðC nþ1 k C nþ1 k1Þ=Ds ðjDsXnþ1k j þ jDsXnþ1k1jÞ=2 ! þ 1 2Pes 1 Ds ðCn kþ1 C n kÞ=Ds ðjDsXn_kþ1j þ jDsXn_kjÞ=2 ðC n k C n k1Þ=Ds ðjDsXn_kj þ jDsXn_k1jÞ=2 : ð35Þ

Since the new interface marker location Xnþ1_k is obtained in the previous step, the above discretization re-sults in a symmetric tri-diagonal linear system which can be solved easily. More importantly, the total mass of surfactant is conserved numerically; that is,

X k Cnþ1_k jDsXnþ1k jDs ¼ X k Cn_kjDsXnkjDs: ð36Þ

(Note that, the summation is exactly the mid-point rule discretization for the integral in Eq. (15).) The above equality can be easily derived by taking the summation of both sides of Eq.(35)and using the peri-odicity of those quantities.

4. Numerical results

The effect of surfactant on the deformation of a drop is of considerable interest in polymer and emulsion industries. It is also a good theoretical model for illustrating subtle physics in viscous interfacial flow. In this section, the immersed boundary method is applied to study the effect of surfactant on drop deformation in Navier–Stokes flows.

Following the set up in [25], we consider a computational domain X¼ ½5; 5 ½2; 2 where a circular drop of radius one is initially located at the center of the domain. We apply a steady shear ﬂow to the drop; that is, we set the boundary condition ub¼ ð0:5y; 0Þ, for 2 6 y 6 2. For comparison purposes, both clean

(without surfactant) and contaminated (with surfactant) drops are used in these computations. Using the equation of state given by Eq.(13), b¼ 0 implies no contamination, in which case we do not need to solve the surfactant equation(19). Throughout this paper, we set rc ¼ 1 so the clean interface has a uniform surface

tension r¼ rc. For the contaminated case, the initial surfactant concentration is uniformly distributed along

the interface such that Cðs; 0Þ ¼ 1. Unless otherwise, we set the Reynolds number Re ¼ 10, the capillary num-ber Ca¼ 0:5, the surface Peclet number Pes¼ 10, and the parameter b ¼ 0:25.

4.1. Convergence test of ﬂuid velocity and surfactant concentration

Before we proceed, we ﬁrst carry out the convergence study of the present method. Here, we perform dif-ferent computations with varying Cartesian mesh h¼ Dx ¼ Dy ¼ 0:04; 0:02; 0:01; 0:005. The Lagrangian mesh is chosen as Ds h=2 and the time step size is Dt ¼ h=8. The solutions are computed up to time T ¼ 1.

Since the analytical solution is not available in these simulations, we choose the results obtained from the ﬁnest mesh as our reference solution and compute the L2error between the reference solution and the solution

obtained from the coarser grid. Table 1 shows the mesh reﬁnement analysis of the velocity u, v, and the

Table 1

The mesh reﬁnement analysis of the velocity u, v, and the surfactant concentration C

h ku urefk2 Rate kv vrefk2 Rate kC Crefk2 Rate

0.04 4.9739E03 – 4.1656E03 – 1.4551E02 –

0.02 2.1476E03 1.21 1.8169E03 1.20 6.3542E03 1.20

(8)

surfactant concentration C. One can see that the error decreases substantially when the mesh is refined, and the rate of convergence is about 1.5. Notice that, the fluid variables are defined at the staggered grid and the surfactant concentration is defined at ‘‘half-integer” grid, so when we refine the mesh, the numerical solutions will not coincide with the same grid locations. In these runs, we simply use a linear interpolation to compute the solutions at the desired locations. We attribute this is part of the reason why the rate of convergence behaves less than second-order.

4.2. Clean vs. contaminated interface

To examine the effect of the surfactant on interfacial dynamics, we compare a drop with and without sur-factant in a steady shear flow. When the sursur-factant are present in the interface, the surface tension can be reduced significantly, cf. equation of state(13). Throughout the rest of this paper, we use a uniform Cartesian mesh h¼ Dx ¼ Dy ¼ 0:02, and a Lagrangian grid with size Ds h=2. The time step size is set to be Dt ¼ h=8.

Fig. 1shows the time evolution plots of drop deformation in a steady shear flow field. Here, we consider three different values of b in Eq.(13); namely, b¼ 0 (dotted, clean interface), b ¼ 0:25 (dash-dotted), and b¼ 0:5 (solid). As expected, the magnitude of drop deformation increases when the value of b increases, as in the case of Stokes flow[25].Fig. 2shows the vorticity plot for the drop with surfactant near the left and the right tips. One can see that two vortices with positive and negative signs are generated near the drop tips. During the drop deformation, the Lagrangian markers will gradually sweep into the tips and cause clus-tered distribution near the tips. If the markers become too crowdedly or too coarsely distributed, it will affect the numerical accuracy. Thus, in order to maintain the numerical stability and accuracy, we need to perform grid redistribution if necessary. The detail is given as follows.

In each time step, we compute the distance between two adjacent markers. If the distance is within an inter-val½0:25h; h, then we basically keep the original resolution. However, if the distance is smaller than 0:25h, then we remove some of the markers. Similarly, when the distance is larger than h, we add more points between these two markers. In general, we just keep the distance between two adjacent markers in a reasonable range. One important thing during the grid redistribution process is to keep the mass conservation of the sur-factant. This can be done in a local way. For instance, in the segment of adding more grid points, we simply

–3 –2 –1 0 1 2 3 –2 –1 0 1 2 –3 –2 –1 0 1 2 3 –2 –1 0 1 2 –3 –2 –1 0 1 2 3 –2 –1 0 1 2 –3 –2 –1 0 1 2 3 –2 –1 0 1 2

Fig. 1. The time evolution of a drop in a shear ﬂow with clean (b¼ 0, ‘’) and contaminated interface (b ¼ 0:25, ‘- - - -’, b ¼ 0:5, ‘–––––’).

(9)

distribute the surfactant mass into those points uniformly. On the other hand, in the segment of removing grid points, we add up those surfactant mass to be a new surfactant concentration in the new combining segment. Thus, the overall surfactant mass is conserved exactly without any scaling.

Plots of the corresponding surfactant concentration (left column) and surface tension (right column) vs. arc-length are given inFig. 3. For the surfactant concentration plot, we omit the case of clean interface since

vorticity –3.2 –3 –2.8 –2.6 –2.4 –2.2 –2 –1.2 –1 –0.8 –0.6 –0.4 –0.2 0 vorticity 2 2.5 3 0 0.2 0.4 0.6 0.8 1 1.2 1.4

Fig. 2. The vorticity plot for the drop with surfactant near the left and right tips (b¼ 0:5; T ¼ 12).

0 5 10 0.2 0.4 0.6 0.8 1 1.2 0 5 10 0.4 0.6 0.8 1 0 5 10 0.2 0.4 0.6 0.81 1.2 0 5 10 0.5 1 0 5 10 0.2 0.4 0.6 0.81 1.2 0 5 10 0.5 1 0 5 10 0.2 0.4 0.6 0.81 1.2 0 5 10 0.5 1

Fig. 3. Distributions of the surfactant concentration (left) and the corresponding surface tension (right). Notations and parameters are same as inFig. 1.

(10)

the concentration is zero everywhere on the interface. It can be seen from this figure, the drop is elongated by the shear flow so that the total length of the interface is increased from the rest state. Since there is no surfac-tant transferred between the interface and the fluid, the surfacsurfac-tant concentration is diluted on a portion of the interface, partly due to the elongation of the interface, but mainly because it is swept to the drop tips. As a result, the smallest surface tension occurs at the drop tips. One can also see that the value of b affects the sur-factant concentration by shifting the distributions slightly along the drop length. Once again, we confirm the same qualitative behavior as in[25].

InFig. 4, the corresponding capillary (deﬁned as rjjoX

osj=ðReCaÞ, left column) and the Marangoni forces

(deﬁned asor

os=ðReCaÞ, right column) are plotted vs. the arc-length for diﬀerent cases of b. Since the capillary

force depends on the curvature and surface tension, we see that the largest capillary force occurs at the drop tips due to the high curvature there. For clean interface, the Marangoni force is obviously zero.

InFig. 5, we present four different plots: namely, (a) total mass of the surfactant; (b) the error of total mass, mðtÞ mð0Þ; (c) total area of the drop; (d) total length of the drop interface. Clearly, the present method pre-serves the total surfactant mass and the errors reach machine precision. However, there is a slight area losing or fluid leakage in the drop as shown inFig. 5c. It seems that the drop without surfactant has a more serious leakage than the ones with surfactant. It is well-known that the fluid leakage often appears in the simulation of immersed boundary method. In[16], Peskin and Printz proposed an improved volume (area in 2D) conserva-tion scheme for the immersed boundary method by constructing a discrete divergence operator based on the interpolation scheme. Here, however, the area loss is not that significant, thus no modification is applied. Once again, we can see fromFig. 5d that the drop with surfactant has larger deformation than the one without sur-factant due to the increase of total length of the interface.

4.3. Linear vs. nonlinear equation of state

In this test, we use the same set up as in the previous one except that a simpliﬁed form of nonlinear Lang-muir equation of state rðCÞ ¼ rcð1 þ lnð1 bCÞÞ is used and compared with the results of the linear equation

of the state. InFig. 6, the evolution of the drop under steady shear ﬂow is shown at diﬀerent times using the linear (dotted) and nonlinear (solid) equations of state with b¼ 0:5. Once again, our results are consistent

0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1 0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1 0 5 10 – 0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1 0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1

(11)

with those in[25], i.e., drop deformation increases when the nonlinear equation of state is used. The corre-sponding surfactant concentrations and surface tensions are shown inFig. 7. One can easily see that the non-linear equation of state generates smaller surface tension at drop tips which leads to a larger deformation. As shown in Fig. 8, the capillary forces are roughly similar but the Marangoni force for the nonlinear case is slightly larger at the drop tips. The four diﬀerent plots for both linear and nonlinear cases including the total mass of the surfactant, the error of total mass, the total area of the drop, and the total length of the drop are shown inFig. 9. 0 2 4 6 8 10 12 6.2832 6.2832 6.2832 6.2832 6.2832 6.2832 time

total mass of surfactant

0 2 4 6 8 10 12 –3 –2 –1 0 1 2x 10 –14 time

error of surfactant mass

0 2 4 6 8 10 12 3.135 3.14 3.145 3.15 time

total area of the drop

0 2 4 6 8 10 12 6 7 8 9 10 11 12 13 time

total length of the drop

Fig. 5. (a) Total mass of the surfactant. (b) Time plot of mðtÞ mð0Þ. (c) Total area of the drop. (d) Total length of the drop interface. Notations and parameters are same as inFig. 1.

–3 –2 –1 0 1 2 3 –2 –1 0 1 2 –3 –2 –1 0 1 2 3 –2 –1 0 1 2 –3 –2 –1 0 1 2 3 –2 –1 0 1 2 –3 –2 –1 0 1 2 3 –2 –1 0 1 2

(12)

4.4. Eﬀect of capillary number on drop deformation

As the last test, we perform the study on how different capillary numbers affect the drop deformation. Here, we fix the Reynolds number Re¼ 10 and the surface Peclet number Pes¼ 10. We vary the capillary number as

0 5 10 0.2 0.4 0.6 0.8 1 0 5 10 0.4 0.6 0.8 0 5 10 0.2 0.4 0.6 0.8 1 0 5 10 0.4 0.6 0.8 0 5 10 0.2 0.4 0.6 0.8 1 0 5 10 0.4 0.6 0.8 0 5 10 0.2 0.4 0.6 0.8 1 0 5 10 0.4 0.6 0.8

Fig. 7. Distributions of the surfactant concentration (left) and the corresponding surface tension (right). Notations and parameters are same as inFig. 6. 0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1 0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1 0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1 0 5 10 –0.8 –0.6 –0.4 –0.2 0 0 5 10 –0.1 0 0.1

(13)

Ca¼ 0:05; 0:25; 0:5; 1:0 and perform our runs up to time T ¼ 4. As conﬁrmed in previous literature such as

[12], a larger capillary number means a smaller surface tension (with the viscosity ﬁxed) so the drop under shear ﬂow can deform more substantially. This is exactly what we see in our simulation as illustrated in

Fig. 10. We also make runs by varying the different surface Peclet number while keeping the Reynolds and capillary numbers fixed. However, the effect of surface Peclet number is not as significant as the effect of the capillary number on drop deformations, so we omit the results here.

5. Conclusion

In this paper, we have developed an immersed boundary method for two-dimensional ﬂuid interfacial prob-lems with insoluble surfactant. The governing equations are formulated in a usual immersed boundary

frame-0 2 4 6 8 10 12 6.2832 6.2832 6.2832 6.2832 6.2832 6.2832 time

total mass of surfactant

0 2 4 6 8 10 12 –4 –3 –2 –1 0 1 2x 10 –14 time

error of surfactant mass

0 2 4 6 8 10 12 3.135 3.14 3.145 3.15 time

total area of the drop

0 2 4 6 8 10 12 6 8 10 12 14 time

total length of the drop

Fig. 9. (a) Total mass of the surfactant. (b) Time plot of mðtÞ mð0Þ. (c) Total area of the drop. (d) Total length of the drop interface. Notations and parameters are same as inFig. 6.

–3 –2 –1 0 1 2 3 –1.5 –1 –0.5 0 0.5 1 1.5

Fig. 10. The eﬀect of capillary number Ca on the drop deformation (Ca¼ 0:05: ‘’, Ca ¼ 0:25: ‘– – – – –’, Ca ¼ 0:5: ‘- - - -’, Ca ¼ 1:0: ‘––––’).

(14)

work where a mixture of Eulerian fluid and Lagrangian interfacial variables are used, with the linkage between those two different variables is provided by Dirac delta function. The immersed boundary force comes from the surface tension which is affected by the distribution of surfactant along the interface. By tracking the inter-face in a Lagrangian manner, a simplified surfactant concentration equation can be obtained. The numerical method involves solving the Navier–Stokes equations on a staggered grid by a semi-implicit pressure incre-ment projection method where the immersed interfacial forces are calculated at the beginning of each time step. Once the velocity values and interfacial configurations are obtained, a new symmetric discretization for the surfactant concentration equation is used such that the total mass of surfactant is conserved numerically.

As a next step, we will generalize the present algorithm to simulate two phase ﬂows with distinct densities and viscosities. In particular, we plan to study the eﬀect of soluble surfactant on drop detachment from a solid surface, i.e., a problem with moving contact points/lines. Finally, we plan to generalize the current work to 3D simulations.

Acknowledgments

M.-C. Lai is supported in part by National Science Council of Taiwan under Research Grant NSC-95-2115-M-009-010-MY2 and MoE-ATU project. H. Huang is supported by grants from the Natural Science and Engineering Research Council (NSERC) of Canada and the Mathematics of Information Technology and Complex Systems (MITACS) of Canada. We thank Dr. Y.-N. Young at NJIT for useful discussions. References

[1] J. Adams, P. Swarztrauber, R. Sweet, Fishpack – a package of Fortran subprograms for the solution of separable elliptic partial diﬀerential equations, 1980. <http://www.netlib.org/ﬁshpack>.

[2] R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Prentice-Hall, Englewood Cliﬀs, NJ, 1962.

[3] D.L. Brown, R. Cortez, M.L. Minion, Accurate projection methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 168 (2001) 464–499.

[4] H.D. Ceniceros, The eﬀects of surfactants on the formation and evolution of capillary waves, Phys. Fluids 15 (1) (2003) 245–256. [5] P.G. De Gennes, F. Brochard, D. Quere, Gouttes, Bulles, Perles Ondes, Edition Berlin, 2002.

[6] F.H. Harlow, J.E. Welsh, Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with a free surface, Phys. Fluids 8 (1965) 2181–2189.

[7] H. Huang, M.-C. Lai, H.-C. Tseng, A parametric derivation of the surfactant transport equation along a deforming ﬂuid interface, submitted for publication.

[8] M. Hameed, M. Siegel, Y.-N. Young, J. Li, M.R. Booty, D.T. Papageorgiou, Inﬂuence of insoluble surfactant on the deformation and breakup of a bubble or thread in a viscous ﬂuid, J. Fluid Mech. 594 (2008) 307–340.

[9] A.J. James, J.S. Lowengrub, A surfactant-conserving volume-of-ﬂuid method for interfacial ﬂows with insoluble surfactant, J. Comput. Phys. 201 (2004) 685–722.

[10] L.D. Landau, E.M. Lifshitz, Fluid Mechanics, Pergamon Press, New York, 1958.

[11] R. LeVeque, Z. Li, The immersed interface method for elliptic equations with discontinuous coeﬃcients and singular sources, SIAM J. Numer. Anal. 31 (1994) 1019–1044.

[12] J. Lee, C. Pozrikidis, Eﬀect of surfactants on the deformation of drops and bubbles in Navier–Stokes ﬂow, Comput. Fluids 35 (2006) 43–60.

[13] M. Muradoglu, G. Tryggvason, A front-tracking method for computations of interfacial ﬂows with soluble surfactants, J. Comput. Phys. 227 (2008) 2238–2262.

[14] C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252. [15] C.S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 1–39.

[16] C.S. Peskin, B.F. Printz, Improved volume conservation in the computation of ﬂows with immersed elastic boundaries, J. Comput. Phys. 105 (1993) 33–36.

[17] Y. Pawar, K.J. Stebe, Marangoni eﬀects on drop deformation in an extensional ﬂow: the role of surfactant physical chemistry, I. Insoluble surfactants, Phys. Fluids 8 (1996) 1738.

[18] L.E. Scriven, Dynamics of a ﬂuid interface, Chem. Eng. Sci. 12 (1960) 98.

[19] H.A. Stone, A simple derivation of the time-dependent convective–diﬀusion equation for surfactant transport along a deforming interface, Phys. Fluids A 2 (1) (1990) 111–112.

[20] P. Tabeling, Introduction to Microﬂuidics, Oxford University Press, 2005.

[21] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase ﬂow, J. Comput. Phys. 169 (2001) 708–759.

(15)

[22] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous incompressible multi-ﬂuid ﬂows, J. Comput. Phys. 100 (1992) 25– 37.

[23] H. Wong, D. Rumschitzki, C. Maldarelli, On the surfactant mass balance at a deforming ﬂuid interface, Phys. Fluids 8 (11) (1996) 3203–3204.

[24] A.M. Waxman, Dynamics of a couple-stress ﬂuid membrane, Stud. Appl. Math. 70 (1984) 63.

[25] J.-J. Xu, Z. Li, J.S. Lowengrub, H.-K. Zhao, A level-set method for interfacial ﬂows with surfactant, J. Comput. Phys. 212 (2006) 590–616.

[26] J.-J. Xu, H.-K. Zhao, An Eulerian formulation for solving partial diﬀerential equations along a moving interface, J. Sci. Comput. 19 (2003) 573–593.