A conservative scheme for solving coupled surface-bulk convection-diffusion equations with an application to interfacial flows with soluble surfactant

(1)

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

A conservative scheme for solving coupled surface-bulk

convection–diffusion equations with an application to

interfacial ﬂows with soluble surfactant

Kuan-Yu Chen

a

, Ming-Chih Lai

b

,

∗

a_{Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 30010, Taiwan}

b_{Center of Mathematical Modeling and Scientiﬁc Computing & Department of Applied Mathematics, National Chiao Tung University, 1001,}

Ta Hsueh Road, Hsinchu 30010, Taiwan

a r t i c l e

i n f o

a b s t r a c t

Article history:

Received 29 November 2011

Received in revised form 17 August 2013 Accepted 3 October 2013

Available online 10 October 2013

Keywords:

A conservative scheme Interfacial ﬂow Soluble surfactant

Adsorption–desorption process Immersed boundary method

Many physical problems arising in biological or material sciences involve solving partial differential equations in deformable interfaces or complex domains. For instance, the surfactant (an amphiphilic molecular) which usually favors the presence in the fluid interface may couple with the surfactant soluble in one of bulk domains through adsorption and desorption processes. Thus, it is important to accurately solve coupled surface-bulk convection–diffusion equations especially when the interface is moving. In this paper, we first rewrite the original bulk concentration equation in an irregular domain (soluble region) into a regular computational domain via the usage of the indicator function so that the concentration flux across the interface due to adsorption and desorption processes can be termed as a singular source in the modified equation. Based on the immersed boundary formulation, we then develop a new conservative scheme for solving this coupled surface-bulk concentration equations which the total surfactant mass is conserved in discrete sense. A series of numerical tests has been conducted to validate the present scheme. As an application, we extend our previous work [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] to the soluble case. The effects of solubility of surfactant on drop deformations in a quiescent and shear flow are investigated in detail.

1. Introduction

Many problems in biological, physical and material sciences involve solving partial differential equations in complex domains or deformable interfaces. In particular, the underlying material quantities in the bulk domain may couple with the one in the interface through adsorption and desorption processes. Meanwhile, the concentration of surface quantities might change the physical behavior of the interface through the modiﬁcations of interfacial forces. For instance, surfactant molecules typically consisting of a hydrophilic head and a hydrophobic tail may adsorb and desorb between bulk ﬂuids and the interface so that the interfacial tension can be reduced. Meanwhile, this non-uniform distribution of surfactant molecules produce extra force (Marangoni force) along the tangential direction to affect the dynamics. In practice, the surfactant may be soluble only to some portion of bulk domain enclosed by the interface where the interface and the soluble region are evolving simultaneously. In order to simulate this problem, we have to introduce two surfactant concentrations in the

*

Corresponding author.

(2)

system; namely, the surface concentration along the interface, and the volume concentration in the bulk region. Thus, one need to solve a coupled system of surface-bulk convection–diffusion equations[18,32,25].

Another example comes from cell biology applications where proteins inside the cell can diffuse and bind to the mem-brane whereas memmem-brane-bound proteins can dissociate and diffuse to the inner cytoplasm[19]. To simulate this problem, one need to solve a coupled system of surface-volume reaction–diffusion equations. Many other examples in physical and biological systems that have the similar adsorption or desorption mechanisms in the dynamics can be found in Ref.[25].

It is known that solving a coupled system of surface-bulk equations in complex domains or deformable interfaces nu-merically is quite challenging especially when the interface (or the interior boundary of domains) is moving. Even in the case of only surface material (without bulk coupling), developing numerical methods for convection–diffusion equations on an evolving interface is still of major interest in scientific computing community. These methods include surface element method [5,6,9,10], level set method [2,7,29,24,17], and phase field method [23,25,11], just to name a few. Front tracking method is typically more accurate for one-dimensional interface (a curve in 2D space), but it can be more complicated to implement for two-dimensional surface especially the implementation involves surface mesh distortion or even topological changes. In[14], we have successfully developed a mass conservative scheme for convection–diffusion equation on moving interface and applied to simulate the interfacial flows with insoluble surfactant[14–16]. A recent work of Khatri and Torn-berg[13] used segment projection method to represent the interface and solve the surfactant equation. More up-to-dated numerical methods for solving Navier–Stokes flows with insoluble surfactant can be found in[13]as well.

In this paper, we shall extend our previous work to soluble surfactant case. However, as a very first step, we need to develop a numerical scheme for solving coupled surface-bulk convection–diffusion equations. There are at least three major numerical issues from our point of view. (1) How to handle the adsorption and desorption between the interface and the bulk accurately? (2) How to maintain the total surfactant mass conserved during the evolution? (3) Since the surfactant might be soluble to only one of buck fluid, how to avoid the surfactant being present in other bulk regions via either convection or diffusion mechanism? Here, we formulate the coupled surface-bulk convection–diffusion equations in the immersed boundary framework so that the adsorption and desorption processes can be termed as a singular source in the bulk equation. Moreover, by using the indicator function, we can embed the bulk equation into the whole computational domain so that regular Eulerian finite difference scheme can be applied without handling the complicated moving irregular domain. We develop a new conservative scheme for solving the coupled bulk-surface concentration equations which the total surfactant mass can be conserved exactly in discrete sense. By introducing the indicator function and solving the bulk equation in the regular computational domain, one can avoid evaluating the surfactant flux across the interface due to adsorption and desorption processes.

The present formulation is similar to other front tracking approaches such as in[32,18]but differs from their numerical computations. For instance, in order to let the surfactant be depleted from only one bulk phase, some one-sided discretized delta functions were used in[32]which results the numerical integration of the discrete function does not yield the exact value of unity. The authors have tried different forms of one-sided delta function and the mass error is within 1%. Here, we use the traditional discrete delta function for the spreading and interpolating operators in the immersed boundary method so that the surfactant mass leaking error is much smaller compared with[32]. There are other numerical methods in literature for interfacial ﬂows with soluble surfactant dynamics such as in[1,4,26,31].

The rest of the paper is organized as follows. In Section2, we present a coupled surface-bulk concentration model and discuss their conservation property. By employing the indicator function, we then embed the bulk equation into a regu-lar Cartesian computational domain. Based on our immersed boundary formulation, we develop a conservative scheme for solving the coupled surface-bulk equations in Section 3. As an application, we apply the present scheme to solve Navier– Stokes ﬂow with soluble surfactant in Section 4. In Section5, a detailed numerical tests have been conducted to validate our present scheme and study the effect of soluble surfactant on drop deformations in a quiescent and shear ﬂow.

2. A coupled surface-bulk concentration model

As in [25], we consider the same coupled bulk-surface material (or surfactant) concentration model in which the ad-sorption and dead-sorption can be occurred on the moving deformable interface. Consider a domain

Ω

in R2 _{and there is an}

interface

Σ, which is a simple closed curve immersed in

Ω. The interior of the interface is

Ω

0, and the exterior is

Ω

1

so that

Ω

= Ω

0

∪ Ω

1, see the illustration of these domains in Fig. 1. The interface is represented by a Lagrangian form

X

(

α

,

t

),

0

α

Lb, where

α

is the Lagrangian material coordinate attached to the interface which is not necessarily to be

the arc-length parameter. The unit tangent vector of the interface can be written as

τ

=

∂X

∂α

/

|

∂X

∂α

|

; thus, the unit outward normal vector n pointing into

Ω

1 can be deﬁned accordingly. In addition, the interface

Σ

is moving with a given velocity

ﬁeld u

= (

u

,

v

)

in

Ω; that is,

∂

X

(

α

,

t

)

∂

t

=

U

(

α

,

t

)

=

Ω

u

(

x

,

t

)δ

x

−

X

(

α

,

t

)

dx

,

(1)

where

δ(

x

)

= δ(

x

)δ(

y

)

is the two-dimensional Dirac delta function. We use the above usual delta function formulation in the immersed boundary method[20]to represent the interpolation of the velocity ﬁeld into the interface. Here we assume

(3)

Fig. 1. Illustration of domains.

the velocity ﬁeld is incompressible (

∇ ·

u

=

0) in

Ω

and no ﬂow boundary condition (u

·

n1

=

0) is imposed on

∂Ω

= ∂Ω

1.

Notice that, in later section, the velocity ﬁeld can be obtained by solving the Navier–Stokes equations.

It is assumed that the surfactant exists on the interface as a monolayer and is adsorbed from or desorbed into the bulk ﬂuid in

Ω

1; that is, the surfactant is soluble in the exterior bulk

Ω

1 but not in the interior one

Ω

0. Therefore, we have to

introduce two surfactant concentrations in the system; namely, the surface concentration

Γ (

α

,

t

)

along the interface

Σ, and

the bulk concentration C(x

,

y

,

t

)

in the region

Ω

1. By taking the adsorption and desorption of bulk surfactant into account,

the dimensionless surface concentration equation can be modiﬁed as

∂Γ

∂

t

+ (∇

s

·

u

)Γ

=

1 Pes

∇

2 s

Γ

+ (

Sa

/λ)

Cs

(

1

− Γ ) −

Sd

Γ,

(2)

where

∇

s

= (

I

−

n

⊗

n

)

∇

and

∇

s2

= ∇

s

· ∇

s are the surface gradient and surface Laplacian operators, respectively. The

dimensionless number Pes is the surface Peclet number, Sa and Sd are the absorption and desorption Stanton number,

respectively, and

λ

is the dimensionless adsorption depth. Those parameters are deﬁned as Pes

=

U∞R

/

Ds

,

Sa

=

ka

/

U∞

,

Sd

=

kdR

/

U∞

,

λ

= Γ∞

/(

C∞R

)

where R,U_∞

, Γ

∞

,

C_∞ are the reference values for the length, ﬂow velocity, the surface and bulk concentration, and ka

,

kd

are the absorption and desorption coeﬃcients. Cs is the bulk surfactant concentration adjacent to the interface which can

be deﬁned later. The above non-dimensionalization process can be found in[32,18,14]. Notice that, as in[14], the interface is tracked in Lagrangian manner and the surface concentration is deﬁned at the material point, so the time derivative in Eq.(2)has the meaning of the material derivative naturally.

The dimensionless bulk concentration in the exterior region

Ω

1 [18,25,26,32]can be written as

∂

C

∂

t

+

u

· ∇

C

=

1 Pe

∇

2_C

_,

₍₃₎ 1

λ

Pe

∂

C

∂

n

Σ

= (

Sa

/λ)

Cs

(

1

− Γ ) −

Sd

Γ

∂

C

∂

n1

∂Ω1

=

0

,

(4)

where Pe is the Peclet number, n is the unit normal vector on

Σ

pointing into

Ω

1 and n1 is the unit outward normal to

the boundary

∂Ω

1

= ∂Ω

.

Eqs.(2)–(4)describe the present coupled surface-bulk concentration equations. Since the ﬂuid is incompressible and no ﬂow velocity boundary condition is imposed on

∂Ω

1, one can conclude that the total surfactant mass (the surfactant mass

on the interface

Σ

and the mass in the bulk region

Ω

1) must be conserved. The conservation property can be proved easily

as follows.

By ﬁrst taking the integration of

Γ

over the interface, then applying the time derivative and using Eq.(2), the rate of change of surfactant mass in the interface

Σ

can be written as

d dt

Σ

Γ

dl

=

Σ

∂Γ

∂

t

+ (∇

s

·

u

)Γ

dl

=

Σ

1 Pes

∇

2 s

Γ

+ (

Sa

/λ)

Cs

(

1

− Γ ) −

Sd

Γ

dl

=

Σ

(

Sa

/λ)

Cs

(

1

− Γ ) −

Sd

Γ

dl

,

(5)

(4)

where dl

= |

∂_∂X_α

|

d

α

is the arc-length element. The last equality is obtained by using the fact that

Σ

is a closed interface. Similarly, the rate of change of surfactant mass in the bulk region

Ω

1 can be obtained by taking the integration of C over

the region

Ω

1, applying the time derivative, and then using Eq.(3)as d dt

Ω1 C dx

=

Ω1 DC Dt dx

=

Ω1 1 Pe

∇

2_{C dx}

₌

∂Ω1 1 Pe

∂

C

∂

n1 dl

−

Σ 1 Pe

∂

C

∂

ndl

= −

Σ

SaCs

(

1

− Γ ) −

Sd

λΓ

dl

,

(6)

where the last equality is obtained by using the boundary conditions of Eq. (4). One can immediately lead to the total surfactant mass conservation by summing Eqs.(5)and(6)so we have

d dt

Ω1 C dx

+ λ

Σ

Γ

dl

=

0

.

(7)

2.1. An embedding bulk concentration equation in a regular Cartesian domain

As mentioned before, solving the bulk concentration equation involves solving a convection–diffusion equation in an evolving irregular domain

Ω

1. In order to describe the solution in a regular Cartesian domain

Ω

= Ω

0

∪ Ω

1, we introduce

the indicator function H deﬁned as H

(

x

,

t

)

=

Ω1

δ(

x

− ˜

x

)

dx

˜

=

1 if x

∈ Ω1

,

0 if x

∈ Ω0

.

(8)

The above indicator function is nothing but the Heaviside function across the interface. By taking the gradient and then divergence operators on both sides, we have

∇

H

(

x

,

t

)

=

Σ n

δ(

x

−

X

)

dl

,

∇

2H

(

x

,

t

)

= ∇ ·

Σ n

δ(

x

−

X

)

dl

.

(9)

Thus, the indicator function can be obtained by solving Poisson equation with a singular source term[28].

By introducing the indicator function, one can rewrite the bulk concentration equation (3)in

Ω

1 with the absorption

and desorption on the interface described in Eq.(4)into one equation in the whole domain

Ω

as

∂

t

(

H C

)

+ ∇ · (

u H C

)

=

1 Pe

∇ · (

H

∇

C

)

−

Σ

SaCs

(

1

− Γ ) −

Sd

λΓ

δ(

x

−

X

)

dl

.

(10)

Here, we rewrite the convection term in a divergence form since the velocity is incompressible. Notice that, in the do-main

Ω

1 where the indicator function H

=

1 so the above equation recovers to the original bulk surfactant equation(3).

Moreover, since H

=

0 in the domain

Ω

0, we have H C

=

0 no matter what are the values of C which restricts that the

surfactant is insoluble in

Ω

0. One should mention that the above bulk concentration equation(10)has the similar form as

the diffuse-interface approach proposed by Teigen et al.[25]except the expression of last term. By taking the integration of Eq.(10)over the domain

Ω, we have

d dt

Ω H C dx

=

Ω

∂

t

(

H C

)

dx

= −

Ω

∇ · (

u H C

)

dx

+

Ω 1 Pe

∇ · (

H

∇

C

)

dx

−

Ω

Σ

SaCs

(

1

− Γ ) −

Sd

λΓ

δ(

x

−

X

)

dl dx

= −

Σ

SaCs

(

1

− Γ ) −

Sd

λΓ

dl

,

where the last equality is obtained by using the no ﬂow velocity boundary condition u

·

n1

=

0, the zero surfactant ﬂux

∂C

∂n1

|

∂Ω1

=

0, and the integral identity of the Dirac delta function

Ω

δ(

x

−

X

)

dx

=

1

.

(11)

One can immediately see that the above rate of change of bulk surfactant mass in the domain

Ω

is exactly the same as the one in the domain

Ω

1as shown in Eq.(6). Based on our new formulation, the bulk surfactant adjacent to the interface can

(5)

Cs

(

α

,

t

)

=

Ω

H C

δ

x

−

X

(

α

,

t

)

dx

.

(12)

3. A conservative scheme for solving the coupled surface-bulk concentration equations

In this section, we describe a numerical scheme to solve the coupled surface-bulk concentration equations(2)and(10). For simplicity, let us assume the computational domain to be a rectangular domain

Ω

= [

a

,

b

] × [

c

,

d

]

. Within this domain, a uniform lattice grid with mesh width h is employed where the given velocity components u and v are deﬁned at usual staggered grid[12]at

(

xi−1/2

,

yj

)

= (

a

+ (

i

−

1)h

,

c

+ (

j

−

1/2)h

)

and

(

xi

,

yj−1/2

)

= (

a

+ (

i

−

1/2)h

,

c

+ (

j

−

1)h

), respectively.

However, the discrete indicator function Hi,j, and the bulk surfactant concentration Ci,j are both deﬁned at the cell center

labeled as x

= (

xi

,

yj

)

= (

a

+ (

i

−

1/2)h

,

c

+ (

j

−

1/2)h

). For the immersed interface, we use a collection of discrete points

α

k

=

k

α

,

k

=

0,1, . . .M such that the Lagrangian markers are denoted by Xk

=

X

(

α

k

)

= (

Xk

,

Yk

). The discrete stretching

factor

|

Dα Xk

|

and the unit tangent vector

τ

k are also deﬁned at the “half-integer” points

α

k+1/2

= (

k

+

1/2)

α

, where the

unit tangent can be approximated by

τ

k

=

DαXk

/

|

DαXk

| =

Xk+1

−

Xk

α

Xk+1

−

Xk

α

.

(13)

Once we have deﬁned the unit tangent vector on the interface, the unit normal vector nkcan be calculated straightforwardly.

The surface concentration

Γ

kis also deﬁned at the “half-integer” points.

Let

t 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

_{, Γ}

n _{and C}n _{are all given.}

Step 1: Compute the new interface position. Since the velocity field u is given through the computation, the first step is to interpolate the new velocity at the fluid lattice points into the marker points and to move the marker points to new positions. Un_k+1

=

i j un_{i j}+1

δ

h

xi j

−

Xnk

h2

,

Xn_k+1

=

Xn_k

+

tUn_k+1

,

(14) where

δ

h is a two-dimensional discrete delta function used in the immersed boundary method such as the one in[20,30].

Step 2: Compute the indicator function. Based on the new interface position Xn_k+1, one can compute the corresponding indicator function Hn+1 by numerically solving Eq.(9)as

∇

2 hH n+1 i j

= ∇

h

·

k nn_k+1

δ

h

xi j

−

Xnk++11/2

DαX n+1 k

α

,

(15)

where X_kn₊+1₁_/₂

= (

X_kn+1₊₁

+

Xn_k+1

)/2. The above discrete Poisson equation can be solved eﬃciently by a fast direct solver

provided by the popular software package FISHPACK[3]. The detailed accuracy issue about solving Eq.(15)can be found in our recent work[8].

Step 3: Compute the surface concentration. As in[14], we ﬁrst multiply Eq.(2)by the interface stretching factor and rewrite the equation in terms of the material coordinate explicitly to obtain

∂Γ

∂

t

∂

X

∂

α

+ (∇

s

·

u

)

∂

X

∂

α

Γ =

_Pe1_s

∂

α

∂Γ

∂

α

∂

X

∂

α

+

qs

∂

X

∂

α

where qs

= (

Sa

/λ)

Cs

(1

− Γ ) −

Sd

Γ

represents the ﬂux between surface and bulk concentration. Using the identity of

∂ ∂t

|

∂X

∂α

| = (∇

s

·

u

)

|

∂X

∂α

|

, the above equation becomes

∂Γ

∂

t

∂

_∂

_α

X

+

_∂

∂

_t

∂

_∂

X

_α

Γ =

_Pe1_s

_∂

∂

_α

∂Γ

_∂

_α

∂

_∂

X

_α

+

qs

∂

_∂

X

_α

.

(16)

We then discretize the above equation by the following Crank–Nicholson scheme.

Γ

_kn+1

− Γ

_kn

t

|

DαXn_k+1

| + |

DαXn_k

|

2

+

|

DαXn_k+1

| − |

DαXn_k

|

t

Γ

_kn+1

+ Γ

_kn 2

=

1 2Pes 1

α

_(Γ

n+1 k+1

− Γ

n+1 k

)/

α

(

|

DαXn_k₊₁+1

| + |

DαXn_k+1

|)/

2

−

(Γ

n+1 k

− Γ

n+1 k−1

)/

α

(

|

DαXn_k+1

| + |

DαXn_k₋₁+1

|)/

2

+

1 2Pes 1

α

_(Γ

n k+1

− Γ

kn

)/

α

(

|

DαXn_k₊₁

| + |

DαXn_k

|)/

2

−

(Γ

_kn

− Γ

_kn₋₁

)/

α

(

|

DαXn_k

| + |

DαXn_k₋₁

|)/

2

+

((

Sa

/λ)

C∗k

(

1

− Γ

n+1 k

)

−

Sd

Γ

n+1 k

)

|

DαX n+1 k

| + ((

Sa

/λ)

Ckn

(

1

− Γ

kn

)

−

Sd

Γ

kn

)

|

DαXnk

|

2

.

(17)

(6)

The adjacent bulk concentration C_k∗and C_kn in last term can be obtained through Eq.(12)as C_k∗

=

i j Hn_{i j}+1Cn_{i j}

δ

h

xi j

−

Xn_k+₊1₁_/₂

h2

,

C_kn

=

i j Hn_{i j}Cn_{i j}

δ

h

xi j

−

Xn_k₊₁_/₂

h2

,

where X_kn₊+1₁_/₂

= (

X_kn₊+1₁

+

Xn_k+1

)/2 since the surface concentration is deﬁned at the material coordinate

α

k+1/2. Since the new

interface marker location Xn_k+1and the indicator function Hn+1are both obtained in previous steps, the above discretization results in a symmetric tri-diagonal linear system which can be solved easily. One should notice that the above surface concentration equation (16) and its numerical scheme (17) can be modiﬁed accordingly if the Lagrangian markers are required to be equally distributed. The detail of an equi-distributed technique for Lagrangian markers by adding an artiﬁcial tangential velocity can be found in our recent work[15].

Step 4: Compute the bulk concentration. The last step is to update the bulk concentration C . We discretize the bulk con-centration equation(10)by the following Crank–Nicholson scheme

(

H C

)

n_{i j}+1

− (

H C

)

n_{i j}

t

+ ∇

h

· ₍

_{uH C}

₎

n+1 i j

+ (

uH C

)

n i j 2

=

1 Pe

∇

h

· (

H

∇

hC

)

_{i j}n+1

+ ∇

h

· (

H

∇

hC

)

n_{i j} 2

−

Qi j (18)

where Qi j is the discrete version of the singular integral in Eq.(10)as

Qi j

=

1 2

k

SaCk∗

1

− Γ

_kn+1

−

Sd

λΓ

_kn+1

DαXn_k+1

δ

h

xi j

−

Xn_k₊₁+1_/₂

α

+

1 2

k

SaCnk

1

− Γ

_kn

−

Sd

λΓ

kn

DαXnk

δ

h

xi j

−

Xnk+1/2

α

.

The difference operator

∇

h

= (

Dx_h

,

D_hy

)

is the regular centered difference approximation on the staggered grid to the gradient

operator. For instance, Dx_h

(

u H C

)

i j

=

ui+1/2,jHi+1/2,jCi+1/2,j

−

ui−1/2,jHi−1/2,jCi−1/2,j

h

,

Dx_h

H D_hxC

_{i j}

=

Hi+1/2,j

(

Ci+1,j

−

Ci j

)/

h

−

Hi−1/2,j

(

Ci,j

−

Ci−1,j

)/

h

,

where Hi+1/2,j

,

Ci+1/2,j are the approximate values deﬁned at the cell edges that can be evaluated as the average of two

neighboring values. One can deﬁne the difference approximation in the y direction D_hyin the similar fashion.

One should notice that, since H

=

0 in

Ω

0, to avoid the division by zero in above scheme, we regularize the indicator

function H by using

√

H2

₊

2 _{instead of H itself where}

_{is chosen about 10}−6_{. This regularization is commonly adopted}

in literature such as in[25].

By taking the summation of Eq. (17)over the interface and using the periodicity of those quantities due to the fact of the closed interface, it leads to the following discrete rate of change of surface concentration

λ

t

k

Γ

_kn+1

DαXn_k+1

α

−

k

Γ

_kn

DαXnk

α

=

k

(

SaCk∗

(

1

− Γ

n+1 k

)

−

Sd

λΓ

kn+1

)

|

DαX n+1 k

| + (

SaCkn

(

1

− Γ

kn

)

−

Sd

λΓ

kn

)

|

DαXnk

|

2

α

.

Similarly, by taking the summation of Eq. (18)over the whole computational domain and applying the discrete boundary conditions of velocity and bulk concentration (u

·

n1

=

0,

∇

hC

·

n1

=

0 on

∂Ω), we then obtain

1

t

i j

(

H C

)

n_{i j}+1h2

−

i j

(

H C

)

n_{i j}h2

= −

k

(

SaC_k∗

(

1

− Γ

_kn+1

)

−

Sd

λΓ

_kn+1

)

|

DαXkn+1

| + (

SaCkn

(

1

− Γ

kn

)

−

Sd

λΓ

kn

)

|

DαXnk

|

2

α

where we use the discrete analogue of the delta function integral identity

i j

δ

h

(

xi j

−

Xk+1/2

)

h2

=

1

.

(19)

By combining the above two summations, we have proved the following discrete conservation for the total surfactant mass as

(7)

i j

(

H C

)

n_{i j}+1h2

+ λ

k

Γ

_kn+1

DαXnk+1

α

=

i j

(

H C

)

n_{i j}h2

+ λ

k

Γ

_kn

DαXnk

α

.

(20)

4. Navier–Stokes ﬂow with soluble surfactant

Consider an incompressible flow problem consisting of two-phase fluids in a fixed two-dimensional square domain

Ω

=

Ω

0

∪ Ω

1 where an interface

Σ

separates

Ω

0 from

Ω

1 as illustrated in Fig. 1. As in previous section, it is assumed that

the surfactant exists on the interface as a monolayer and is adsorbed from or desorbed to the bulk ﬂuid in

Ω

1; that is,

the surfactant is soluble in the exterior bulk but not in the interior one. The interface is contaminated by the surfactant so that the distribution of the surfactant changes the surface tension accordingly. In order to formulate the problem using the immersed boundary approach, we simply treat the interface as an immersed boundary that exerts force to the fluids and moves with local fluid velocity. For simplicity, we assume equal viscosity and density for both fluids, and neglect the gravity. Certainly, the present Navier–Stokes solver can be replaced by the one with different density and viscosity ratios.

As in[14], the non-dimensional Navier–Stokes ﬂow in the usual immersed boundary formulation can be written as

∂

u

∂

t

+ (

u

· ∇)

u

+ ∇

p

=

1 Re

∇

2_u

₊

f ReCa

,

(21)

∇ ·

u

=

0

,

(22) f

(

x

,

t

)

=

Σ F

(

α

,

t

)δ

x

−

X

(

α

,

t

)

d

α

,

(23)

∂

X

(

α

,

t

)

∂

t

=

U

(

α

,

t

)

=

Ω u

(

x

,

t

)δ

x

−

X

(

α

,

t

)

dx

,

(24) F

(

α

,

t

)

=

∂

α

σ

(

α

,

t

)

τ

(

α

,

t

)

,

(25)

where u is the ﬂuid velocity and p is the pressure. The dimensionless numbers are the Reynolds number (Re

=

ρ

U_∞R

/

μ

) describing the ratio between the inertial force and the viscous force, and the Capillary number (Ca

=

μ

U_∞

/

σ

_∞) describing the strength of the surface tension. The presence of surfactant will reduce the surface tension of the interface by the Langmuir equation of state[22]

σ

=

1

+

El ln

(

1

− Γ ),

(26)

where

σ

is the surface tension, and El is the elasticity number measuring the sensitivity of the surface tension to the surfactant concentration. Since the surfactant is soluble in

Ω

1, we need to solve the coupled surface-bulk concentration

equations(2)–(4)to close the system.

In the following, we describe how to march one time step for the solutions. At the beginning of each time step, the interface position, the ﬂuid velocity, the surface and bulk concentrations must be given. The numerical algorithm is as follows.

1. Compute the surface tension and unit tangent on the interface.

2. Distribute the interfacial force from the Lagrangian markers into the ﬂuid. 3. Solve the Navier–Stokes equations by the projection method.

4. Interpolate the new velocity on the ﬂuid lattice points onto the marker points and move the marker points to new positions as shown in Eq.(14).

5. Compute the indicator function as Eq.(15).

6. Compute the surface surfactant concentration using the scheme of Eq.(17). 7. Compute the bulk surfactant concentration using the scheme of Eq.(18).

Note that, the detailed numerical implementation of ﬁrst four steps is quite standard in immersed boundary method and can be found in any related literature. Here, we just use the same solver as in our previous work[14]. The last four steps are exactly the same four steps shown in previous section.

5. Numerical results

In order to validate our numerical scheme described in previous sections, we perform a series of tests to check if our code is implemented correctly. Since we have detailed discussions on the numerical conservation of surface concentration in our previous work[14], here we focus on the numerical solution of bulk concentration and its coupling with the surface. Throughout this section, we choose the dimensionless adsorption depth

λ

=

1 which refers approximately to a drop of micron in glycerol/water solutions of the polyethoxylate surfactant C12E6 [21]. The computational domain is chosen as

(8)

Table 1

The L2and L1errors and their convergent rates for the bulk diffusion with a ﬁxed interface at T=0.5.

N×N (H C)2N− (H C)N 2 rate (H C)2N− (H C)N 1 rate 64×64 2.0344E−02 – 1.0424E−02 – 128×128 1.6218E−02 0.3270 4.7816E−03 1.1243 256×256 1.0940E−02 0.5679 2.4963E−03 0.9376 512×512 3.8809E−03 1.4951 2.9607E−04 3.0757 Table 2

The convergent study of numerical leakage relative error ML

M0 at T=0.5 for three different regularization parameter.

N×N =10−6 _rate ₌₀_._1h _rate ₌₀_._1h2 _rate

64×64 1.7175E−05 – 1.4376E−04 – 2.1261E−05 –

128×128 1.2707E−05 0.4346 8.1861E−05 0.8124 1.3810E−05 0.6225 256×256 7.0408E−06 0.8518 4.3141E−05 0.9241 7.3280E−06 0.9142 512×512 3.3182E−06 1.0853 2.1641E−05 0.9953 3.3904E−06 1.1119

5.1. Bulk diffusion with a ﬁxed interface

As the ﬁrst test, we consider the diffusion of surfactant in the exterior bulk phase

Ω

1. The interface (the inner boundary

of

Ω

1) is ﬁxed and chosen as a circle centered at

(

x0

,

y0

)

= (

0,0)with radius r0

=

0.3. To consider only the bulk diffusion,

we set the velocity to be zero and turn off the bulk and surface coupling; that is, we solve Eqs. (3)–(4)with u

=

0 and Sa

=

Sd

=

0. We also set the Peclet number Pe

=

100. The initial bulk concentration is set as

C

(

x

,

y

,

0

)

=

₀

_.

₅

₊

₀

_.

_{4 cos}

₍

_π

_x

₎

_cos

₍

_π

_y

₎

_{if r}

_>

₂

_.

_5r 0

,

¯

C

+ (

0

.

4 cos

(

π

x

)

cos

(

π

y

)

+

0

.

5

− ¯

C

)

w

(

r

)

if r0

r

2

.

5r0

,

0 otherwise

.

(27) where the function

w

(

r

)

=

1 2

1

−

cos

(

r

−

r0

)

π

1

.

5r0

(28) with r

=

(

x

−

x0

)

2

+ (

y

−

y0

)

2, andC is the concentration value at the interface. Here, we choose

¯

C as the average value

¯

of the function 0.5

+

0.4 cos(

π

x

)

cos(

π

y

)

along the interface r

=

r0. Note that, by the choice of w(r

), one can immediately

check that this initial condition satisﬁes the boundary condition Eq.(4)since w(r0

)

=

w

(

r0

)

=

0.

We ﬁrst carry out the convergence study of the present scheme. Here, we perform different computations with vary-ing grid numbers N

×

N

,

N

=

64,128,256,512,1024 in

Ω

= [−

1,1

] × [−

1,1

]

so the spatial mesh width is h

=

2/N. The Lagrangian marker width is chosen as

α

≈

h

/2 and the time step size is

t

=

h

/8. The solutions are computed at time

T

=

0.5. Since the analytical solution is not available in these simulations, we compute the error between two successive grids denoted by

(

H C

)

2N

− (

H C

)

N

, so the rate of convergence can be computed as rate

=

log2 (

H C)N−(H C)N/2

(H C)2N−(H C)N . Table 1

shows the L2 and L1errors and their convergent rates for the bulk concentration. One can see that the rate of convergence

is averagely ﬁrst-order which is a typical accuracy behavior using the immersed boundary method[20].

Fig. 2shows the bulk concentration along the horizontal line y

=

0 at different times using the grid number N

=

256. One can see that inside the interface, the bulk concentration is almost zero which reﬂects the insolubility to the region

Ω

0.

Meanwhile, the surfactant diffuses gradually as time evolves in these plots. To see if the total surfactant mass is conserved in the domain, we plot the relative error deﬁned by Mt−M0

M0 , where M0 is the initial total surfactant mass while Mt is the total mass computed at time t. The upper panel of Fig. 3shows the time evolutionary plot of the total mass relative error. One can see the error is in the magnitude of 10−13 which indicates that the present scheme has excellent conservation property.

Although in the plots of Fig. 2, the bulk concentration is indistinguishable zero inside the interface; however, there is still an O

(

)

amount of mass leaking into the region

Ω

0. In the lower panel ofFig. 3, we show the time evolutionary plot

of the relative error for the leaking mass inside the interface. Here, the relative error is deﬁned as ML

M0 in which the leaking mass ML is computed by ML

=

i j

(

H C

)

i jh2where Hi j

0.01. Note that, the present relative error is within 10−5 which is

only about 0.001%.

We attribute this mass leakage to the regularization of the indicator function. To further illustrate the dependence of regularization parameter

on the effect of numerical leakage error, we test three different choices of

=

10−6

,

0.1h and 0.1h2 _{as shown in}_{Table 2. The errors are all computed up to time T}

₌

_{0.5. One can see that all three cases show similar}

ﬁrst-order convergence as the mesh is reﬁned. Meanwhile, the errors of

=

10−6 and

=

0.1h2 _{are both comparable, and}

(9)

Fig. 2. The bulk concentration along the horizontal line y=0 at different times.

Fig. 3. Upper panel: The time evolutionary plot of total mass relative error. Lower panel: The time evolutionary plot of leaking mass relative error inside the interface.

(10)

Fig. 4. The bulk concentration at different times. Left column: The bulk concentration contour plots and the interface positions. Right column: The bulk concentration plots along the horizontal line y=0. The dashed line in the ﬁrst right plot denotes the initial bulk concentration along y=0.

5.2. Bulk convection–diffusion with a moving interface

In second test, we use the same initial setup as the first test but add the convection effect into the equation. The flow and the interface are now moving with a prescribed incompressible velocity field u

= (

u

,

v

)

as

u

= −

1 2

1

+

cos

(

π

x

)

sin

(

π

y

),

v

=

1 2

1

+

cos

(

π

y

)

sin

(

π

x

).

(29)

The center of the circular interface is shifted to the point of

(

x0

,

y0

)

= (

0.1,0)initially so that the asymmetric ﬂow can be

developed.

Table 3shows the mesh reﬁnement analysis of the bulk concentration in

Ω. Again, one can see that the rate of

conver-gence is ﬁrst-order even with a moving interface.

Fig. 4 shows the bulk concentration at different times. The left column shows the contour plots of bulk concentration and the instantaneous interface positions. The right column of the ﬁgure shows the concentration plots along the horizontal

(11)

Table 3

The L2and L1errors and their convergent rates for the bulk convection–diffusion with a moving interface at T=0.5.

N×N (H C)2N− (H C)N 2 rate (H C)2N− (H C)N 1 rate

64×64 2.3101E−02 – 1.1570E−02 –

128×128 1.5788E−02 0.5491 5.3424E−03 1.1149

256×256 1.1002E−02 0.5210 2.7364E−03 0.9652

512×512 3.9145E−03 1.4909 3.2451E−04 3.0759

line y

=

0. One can see that the ﬂow is asymmetric and both convection and diffusion indeed have apparent effects on the distribution of bulk concentration.Fig. 5shows the relative errors for total mass in

Ω

and leaking mass into

Ω

0. Again, the

total surfactant mass (upper panel) is conserved almost as well as the ﬁrst test and the leaking mass error is still controlled within 0.002%. Note that, the oscillatory behavior of the latter error is due to the convection of the ﬂow.

5.3. Surface-bulk coupling with a moving interface

In this test, we consider the surface-bulk coupling of the system with a moving interface; that is, we need to solve the coupling equations (2)–(4) with the prescribed velocity u deﬁned in Eq. (29). Here, we choose the Peclet numbers Pe

=

Pes

=

100 and the adsorption and desorption Stanton numbers as Sa

=

Sd

=

1.0. As previously, the initial interface is a

circle with the center located at

(

x0

,

y0

)

= (

0.1,0)and the radius r0

=

0.3. The initial bulk concentration is deﬁned as C

(

x

,

y

,

0

)

=

0

.

5

(

1

−

x2

₎

2 _{if r}

_>

₂

_.

_5r 0

,

0

.

5

(

1

−

x2

)

2w

(

r

)

if r0

r

2

.

5r0

,

0 otherwise (30) where the function w(r

)

is deﬁned as same as Eq.(28). We also set the initial surface concentration to be identically zero (i.e.

Γ (

α

,

0)

=

0) so that signiﬁcant surface absorption can be observed near the interface. As before, by the choice of w

(

r

), the initial bulk and surface concentrations satisfy the boundary condition Eq.

(4).Table 4shows the mesh reﬁnement analysis of the bulk and surface concentrations in L2norm. One can see that the rate of convergence is roughly ﬁrst-order

for both bulk and surface concentrations.

Fig. 6shows the bulk concentration at different times based on the results of grid number N

=

256. The upper panel shows the contour plots of bulk concentration and the instantaneous interface positions. The lower left column of the ﬁgure shows the bulk concentration plots along the horizontal line y

=

0, while the lower right column shows the surface concentration along the interface. Since the total mass is conserved, the bulk concentration decreases while the surface concentration increases signiﬁcantly due to the surface adsorption process. The relative errors for total mass in

Ω

and leaking mass into

Ω

0 are shown inFig. 7. Again, the total surfactant mass is conserved as well as before and the leaking

(12)

Fig. 6. Upper panel: The bulk concentration contour plots and the interface positions. Lower left: The bulk concentration plots along the horizontal line y=0. The dashed line in the ﬁrst plot denotes the initial bulk concentration. Lower right: The surface concentration along the interface is plotted in counter-clockwise direction starting from the point marked by “◦”.

Table 4

The L2errors and their convergent rates for the bulk and surface concentrations at T=0.5.

N×N (H C)2N− (H C)N 2 rate Γ2N− ΓN 2 rate

64×64 1.9479E−03 – 1.9494E−04 –

128×128 8.5982E−04 1.1798 5.2554E−05 1.8911

256×256 4.3928E−04 0.9689 1.5999E−05 1.7158

(13)

Table 5

The L2errors and their convergent rates for the bulk and surface surfactant concentrations, and the ﬂuid velocity ﬁeld at T=0.5.

N×N (H C)2N− (H C)N 2 rate Γ2N− ΓN 2 rate 64×64 1.9723E−02 – 1.4433E−03 – 128×128 1.2555E−02 0.6516 4.7029E−04 1.6178 256×256 9.2176E−03 0.4457 1.9054E−04 1.3034 N×N u2N−uN 2 rate v2N−vN 2 rate 64×64 4.3642E−03 – 4.0541E−03 – 128×128 1.2665E−03 1.7848 1.2908E−03 1.6511 256×256 4.0467E−04 1.6460 4.3488E−04 1.5695

5.4. A freely oscillating drop

We now consider a freely oscillating drop immersed in a quiescent ﬂow domain

Ω

= [−

1,1

] × [−

1,1

]

as our first numerical test for Navier–Stokes flow with soluble surfactant. Since the flow is quiescent initially, the initial velocity of Navier–Stokes is zero everywhere, and the no-slip boundary conditions are imposed on the computational boundary. The initial drop interface is an ellipse defined as X(

α

,

0)

= (

0.6 cos

α

,

0.15 sin

α

), 0

α

2

π

. Unlike the previous tests, we now set the initial bulk concentration in

Ω

1 (outside of the ellipse) to be a constant C

¯

s, i.e. C

(

x

,

y

,

0)

= ¯

Cs in

Ω

1. The bulk

surfactant adjacent to the interface Cs

(

α

)

is naturally equal toC

¯

swhich leads to zero Neumann condition ∂_∂C_n

=

0 along the

interface. We then choose the initial surface concentration

Γ (

α

,

0)to be another constant

Γ

¯

such that

(

Sa

/λ) ¯

Cs

(1

− ¯

Γ )

=

Sd

Γ

¯

. Based on the above choices of the initial bulk and surface concentrations, one can immediately see those initial

conditions satisfy the boundary conditions in Eq. (4). Here, we choose the dimensionless absorption length as

λ

=

1, and the absorption and desorption numbers as Sa

=

3 and Sd

=

1, respectively. We also choose the initial bulk constantC

¯

s

=

0.5

which leads the initial surface concentration constant to be

Γ

¯

=

0.6. Other dimensionless numbers are chosen as Pe

=

10, Pes

=

100, Re

=

10, Ca

=

2, and El

=

0.5.

As before, we ﬁrst carry out a convergence study of the present numerical scheme described in Section 4. Here, we perform four different computations with varying grid numbers N

=

64,128,256,512. The Lagrangian marker width is chosen as

α

≈

h

/2 and the time step size is

t

=

h

/8. Since the analytical solution is not available in these simulations,

we compute the L2 error between two successive grids and the rate of convergence is computed as in previous examples.

Table 5shows the mesh reﬁnement analysis of the bulk and surface concentration, and the velocity ﬁeld at T

=

0.5. One can see that the rate of convergence still behaves like ﬁrst-order. The time evolutionary plots of relative errors for the total mass and the leaking mass inside the interface are similar to those in previous tests so we omit here.

Next, we make the comparison of a freely oscillating drop in a quiescent ﬂow with a soluble and insoluble surfactant. Note that, for the insoluble surfactant case, there is no surface-bulk coupling thus no need to solve the bulk concentration Eq.(3).Fig. 8shows the comparison between insoluble (denoted by dash-dotted line) and soluble (denoted by solid line)

(14)

Fig. 8. The comparison of insoluble (denoted by “-·”) and soluble (denoted by “-”) interfacial ﬂows for a freely oscillating drop. The dotted lines in the ﬁrst plots of lower panel denote the initial concentrations. The surface concentrations along the interfaces are plotted in counter-clockwise direction starting from the point marked by “◦”.

cases for a freely oscillating drop based on the results of grid number N

=

256. The upper panel shows the interface positions at different times, while the lower left and right show the bulk concentration along y

=

0 (soluble case only) and surface concentrations, respectively. It is important to mention that due to the shrinkage of the interface (the total perimeter of the interface becomes smaller), the surface concentration for the insoluble case is generally higher than the initial constant concentration

Γ

¯

. Here, the surface concentrations are plotted in terms of initial conﬁguration coordinate

α

rather than the arc-length coordinate.

Due to the surface tension force, the drop tends to oscillate until it reaches to a stationary circular shape. For the soluble case, the surface surfactant will be desorbed into the neighboring bulk region (see the lower left panel) so the bulk concentration is higher near the interface while the surface concentration is lower than the one of the insoluble case (see

(15)

Fig. 9. The interface positions of clean (denoted by “-·”) and soluble (denoted by “-”) interfacial ﬂows for a drop under shear ﬂow.

the lower right panel). As a result, the surface tension will not be reduced as much as in the insoluble case so the drop with soluble surfactant tends to oscillate faster than the insoluble case. Our numerical simulation conﬁrms the above drop behaviors qualitatively.

5.5. A drop under shear ﬂow

As a ﬁnal test, we put a circular drop centered at the origin with radius r0

=

0.3 initially in the domain

Ω

= [−

1,1

] ×

[−

1,1

]

under shear ﬂow. The initial and boundary conditions for the velocity used in this simulation are based on the setup of Tornberg et al. in[27]. That is, the velocity boundary conditions are u

= ±

1, v

=

0 at y

= ±

1, respectively and u,v are both periodic at x

= ±

1. As also suggested in[27], we set the initial velocity as the quiescent flow to avoid unrealistic flow interior to the drop during the transition caused by using the linear velocity initially. That is, the flow is set to motion simply by the boundary shear. The boundary conditions for the bulk concentration equation are chosen as ∂_∂C_y

=

0 at y

= ±

1

(16)

Fig. 10. Left column: The evolutionary bulk concentration along y=0 for soluble case. Right column: The evolutionary surface concentration for soluble (denoted by “-”) case. The surface concentrations along the interface are plotted in counter-clockwise direction starting from the point marked by “◦” in Fig. 9.

and are periodic at x

= ±

1 to be consistent with the ﬂow conditions. Therefore, the total surfactant mass is still conserved in this case. The initial bulk concentration is deﬁned as

C

(

x

,

y

,

0

)

=

sin2

(

π

x

)

if r

>

2

.

5r0 sin2

(

π

x

)

w

(

r

)

if r0

r

2

.

5r0 0 otherwise (31)

where the function w

(

r

)

is deﬁned as same as Eq. (28). We also set the initial surface concentration to be identically zero (i.e.

Γ (

α

,

0)

=

0) so that signiﬁcant surface absorption can be expected near the interface. As before, by the choice of w(r

), the initial bulk and surface concentrations satisfy the boundary condition Eq.

(4)at the interface r

=

r0 because

(17)

Table 6

The L2errors and their convergent rates for the bulk and surface surfactant concentrations, and the ﬂuid velocity ﬁeld at T=0.5.

N×N (H C)2N− (H C)N 2 rate Γ2N− ΓN 2 rate 64×64 6.2645E−03 – 1.1654E−03 – 128×128 3.8501E−03 0.7022 3.9683E−04 1.5542 256×256 2.3775E−03 0.6954 1.2670E−04 1.6471 N×N u2N−uN 2 rate v2N−vN 2 rate 64×64 5.8396E−04 – 3.7335E−04 – 128×128 2.7426E−04 1.0903 1.8524E−04 1.0111 256×256 1.4308E−04 0.9386 9.1926E−05 1.0109

of w(r0

)

=

w

(

r0

)

=

0. Other dimensionless numbers are

λ

=

1, Sa

=

3, Sd

=

1, Pe

=

10, Pes

=

100, Re

=

10, Ca

=

4/3, and

El

=

0.5.

As in previous test, we ﬁrst perform four different computations by varying grid number N

=

64,128,256,512 with the associated mesh width h

=

2/N. The Lagrangian marker width is chosen as

α

≈

h

/2 and the time step size is

t

=

h

/8.

Table 6shows the mesh reﬁnement analysis of the bulk and surface concentrations, and the velocity ﬁeld at T

=

0.5. One can see that again the rate of convergence behaves like ﬁrst-order in general.

Fig. 9shows the evolutionary interface positions of clean (denoted by dash-dotted line) and soluble surfactant (denoted by solid line) cases for a drop under shear ﬂow based on the results of grid number N

=

256. The clean drop bears no surfactant along the interface throughout the evolution so no bulk and surface surfactant equations are needed to be solved and the surface tension remains to be a constant

σ

=

1. (Note that, we use the clean drop as a comparison simply because of zero initial surface concentration is chosen in present setting.) Due to shear stresses, both drops will be elongated and gradually aligned with the ﬂow directions. For the soluble case, the interface will start to absorb the bulk surfactant so the bulk concentration decreases while the surface concentration increases in the beginning, see Fig. 10in detail. Later, both absorption and desorption processes become more balanced so the bulk and surface concentrations become quite steady. As expected, the largest surface concentration appears to occur at the drop tips after the drop aligned with the ﬂow. The drop with soluble surfactant has smaller surface tension than the clean drop so the deformation tends to be larger. One can see fromFig. 9that the clean drop approaches to a steady state shape after T

=

9.0 while the soluble surfactant drop continues to deform slightly afterwards. Our numerical results are physically reasonable and qualitatively consistent with those obtained in other literature such as in[26].

Acknowledgements

The authors thank the referee for pointing out the consistent usage of initial and boundary conditions for bulk and surface concentrations in simulations. We also thank the editor for pointing out Ref.[27]of using the initial quiescent ﬂow in the simulations of a drop under shear ﬂow. This work was supported in part by National Science Council of Taiwan under research grants NSC-98-2115-M-009-014-MY3, NSC-100-2911-I-009-504, and the support of NCTS in Taiwan.

References

[1]S. Adami, X.Y. Hu, N.A. Adams, A conservative SPH method for surfactant dynamics, J. Comput. Phys. 229 (2010) 1909–1926.

[2]D. Adalsteinsson, J. Sethian, Transport and diffusion of material quantities on propagating interfaces via level set methods, J. Comput. Phys. 185 (2003) 271–288.

[3] J. Adams, P. Swarztrauber, R. Sweet, Fishpack – A package of Fortran subprograms for the solution of separable elliptic partial differential equations, available athttp://www.netlib.org/ﬁshpack, 1980.

[4]M.R. Booty, M. Siegel, A hybrid numerical method for interfacial fluid flow with soluble surfactant, J. Comput. Phys. 229 (2010) 3864–3883. [5]M. Burger, Numerical simulation of anisotropic surface diffusion with curvature-dependent energy, J. Comput. Phys. 203 (2005) 602–625. [6]E. Bänsch, P. Morin, R.H. Nochetto, A finite element method for surface diffusion: the parametric case, J. Comput. Phys. 203 (2005) 321–343. [7]M. Bertalmio, L.-T. Cheng, S.J. Osher, G. Sapiro, Variational problems and partial differential equations on implicit surfaces, J. Comput. Phys. 174 (2001)

759–780.

[8]K.-Y. Chen, K.-A. Feng, Y. Kim, M.-C. Lai, A note on pressure accuracy in immersed boundary method for Stokes ﬂow, J. Comput. Phys. 230 (2011) 4377–4383.

[9]G. Dziuk, C.M. Elliott, Finite element on evolving surfaces, IMA J. Numer. Anal. 27 (2007) 262–292. [10]G. Dziuk, C.M. Elliott, Surface ﬁnite elements for parabolic equations, J. Comput. Math. 25 (2007) 385–407.

[11]C.M. Elliott, B. Stinner, V. Styles, R. Welford, Numerical computation of advection and diffusion on evolving diffuse interfaces, IMA J. Numer. Anal. 31 (2011) 786–812.

[12]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.

[13]S. Khatri, A.-K. Tornberg, A numerical method for two phase ﬂows with insoluble surfactant, Comput. Fluids 49 (2011) 150–165.

[14]M.-C. Lai, Y.-H. Tseng, H. Huang, An immersed boundary method for interfacial ﬂows with insoluble surfactant, J. Comput. Phys. 227 (2008) 7279–7293. [15]M.-C. Lai, Y.-H. Tseng, H. Huang, Numerical simulation of moving contact lines with surfactant by immersed boundary method, Commun. Comput.

Phys. 8 (2010) 735–757.

[16]M.-C. Lai, C.-Y. Huang, Y.-M. Huang, Simulating the axisymmetric interfacial ﬂows with insoluble surfactant by immered boundary method, Int. J. Numer. Anal. Model. 8 (2011) 105–117.