## 國

## 立

## 交

## 通

## 大

## 學

### 應用數學系

### 碩

### 士

### 論

### 文

### 模擬帶有表面活性劑之軸對稱界面流問題

### Simulating Interfacial Flows with Surfactant – Axisymmetric Case

### 研 究 生：黃義閔

### 指導教授：賴明治 教授

### 模擬帶有表面活性劑之軸對稱界面流問題

### Simulating Interfacial Flows with Surfactant – Axisymmetric case

### 研 究 生：黃義閔 Student：Yi-Min Huang

### 指導教授：賴明治 Advisor：Ming - Chih Lai

### 國 立 交 通 大 學

### 應 用 數 學 系

### 碩 士 論 文

A Thesis

Submitted to Department of Applied Mathematics College of Science

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Applied Mathematics

June 2009

Hsinchu, Taiwan, Republic of China

**模擬帶有表面活性劑之軸對稱界面流問題 **

### 學生：

黃義閔### 指導教授

：賴明治### 國立交通大學應用數學學系﹙研究所﹚碩士班

### 摘

### 要

### 本論文把 Immersed boundary method 發展到軸對稱圓柱座標系中。為方便分

### 佈界面所產生的力至流體以及內插流體流速至界面的移動速度，須藉由 Dirac

### delta function 將 Eulerian 流體和 Lagrangian 界面變數合併。我們使用

### semi-implicit second-order projection method 來解 Navier-Stokes 方程式。把液滴

### 放入不同的流場中，觀察有加表面活性劑與沒加表面活性劑的差異性。表面活性

### 劑會分佈在界面中，並且影響表面張力的大小。在解完上述的方程式後，我們加

### 入人工速度去調整界面上網格的位置，使其分佈均勻。然後比較有加人工速度與

### 沒加人工速度之間的差異。最後我們把液滴放到固體平面板上，測試 moving

### contact line 的問題。

**Simulating Interfacial Flows with Surfactant – **

**Axisymmetric Case **

### student：

Yi-Min Huang### Advisors：Dr.

Ming-Chih Lai### Department﹙Institute﹚of

Applied Mathematics### National Chiao Tung University

### ABSTRACT

### In this paper, we develop an immersed boundary method in the axisymmetric

### cylindrical coordinates. Coupling the Eulerian flow and Lagrangian interfacial

### variables by the Dirac delta function is convenient to spread the interface force and

### interpolate the interfacial velocities. We solve the Navier-Stokes equations by

### semiimplicit second-order projection scheme. Test a variety of cases that the drops

### deform in the di_erent velocity field and observe the di_erence of the clean drop and

### the drop contaminated with the surfactant. The surface tension is a_ected by the

### distribution of surfactant along the interface. After those steps, we impose the

### artificial velocity to control the distribution of marker positions and compare this with

### the original manner without imposing artificial velocity. Finally, we put a drop on a

### solid to observe the contact line moving.

### 誌

### 謝

### 本篇論文的完成，首先要感謝我的指導教授 ─ 賴明治老師。在老師的細心

### 指點之下，讓我一點一滴熟悉數值偏微分方程以及科學計算的方法。我也從自己

### 的研究中找到對數值偏微分方程以及科學計算這塊領域的興趣與成就感。除了指

### 導教授外，在我們科學計算研究室中的其他學長們也給了我很多的幫助。感謝陳

### 冠羽學長在我剛踏入這塊領域時，指導我撰寫程式語言時所遇到的困難以及課業

### 上的問題。感謝曾昱豪學長與黃仲尹學長很有耐心地教我流體力學的相關理論以

### 及使用 Immersed boundary method 時會碰到的疑難雜症。

### 在論文口試期間，承蒙楊肅煜教授與郭鴻基教授費心審閱並提供許多寶貴意

### 見，使得本論文變得更加完備，學生永遠銘記在心。

### 平常除了做數值計算研究外，我也要感謝同一屆的同學們與研究室的學長

### 們，他們陪伴我閒話家常，一起去運動，交換各種休閒興趣的心得，讓我的研究

### 生生活多采多姿，處處充滿著樂趣。

### 最後，我要好好地感謝我的家人，是他們鼓勵並支持著我做研究的興趣，讓

### 我可以無所顧忌地完成我的學業。希望能與他們，以及在我身邊所有關心我的

### 人，一起分享這篇論文完成的喜悅與榮耀。

### 目

### 錄

### 中文提要

### ………

### i

### 英文提要

### ………

### ii

### 誌謝

### ………

### iii

### 目錄

### ………

### Iv

### 一、

### Introduction ………

### 1

### 二、

### The governing equations ………

### 2

### 2.1

### Immersed boundary formulation ………

### 3

### 2.2

### Surfactant concentration equation ………

### 5

### 2.3

### Controlling the interfacial markers uniformly …

### 6

### 三、

### Numerical method ………

### 7

### 四、

### Numerical results ………

### 11

### 4.1

### Oscillating drop ………

### 11

### 4.2

### Convergence test of fluid velocity and

### surfactant concentration ………

### 12

### 4.3

### Clean vs. contaminated interface ………

### 13

### 4.3.1 The steady state of the clean interface in the

### extensional flow ………

### 13

### 4.3.2 The effect for the controlled markers on the

### interface ………

### 15

### 4.4

### Forming the dumbbell - shaped drop in another

### fluid velocity field ………

### 18

### 4.5

### Moving contact line on the slip boundary ………

### 20

### 4.5.1 Convergence test of fluid velocity and

### surfactant concentration ………

### 20

### 4.5.2 Clean vs. contaminated interface with

### θ

s### =60˚ …

### 21

### 4.5.3 Clean vs. contaminated interface with

_{θ}

_{s}

_{=120˚… }

### 23

### 4.5.4 Imposing initial velocity on our moving contact

### line case ………

### 26

### 五、

### Conclusion ………

### 29

### Appendix

### ………

### 29

### A

### Equations in cylindrical coordinates ………

### 29

### B

### Artificial velocity in Eq.(2.19) ………

### 30

### C

### Artificial velocity in Eq.(2.29)………

### 31

### 1 Introduction

The main purpose of this paper is to develop an immersed boundary method for Navier-Stokes equations in cylindrical geometries, and to use this solver to simulate the three-dimensional fluid interfaces with insoluble surfactant. We shall restrict ourselves in this paper to the axisymmetric case. Surfactant are surface active agents such as soap, deter-gent and food grease, which form a kind of skin on the surface of liquid and affect the interface surface tension. Effects of the surfactant are also important in bio-mechanical flows. For example, some insects drift on the water by injecting a chemical at the rear. This chemical reduces the surface tension behind their bodies so that the insects are pulled forward.

We refer to [5] and [6] in order to transform the immersed boundary formulation [7] in three-dimensional coordinates into the cylindrical coordinates. For our axisymmetric

*case, we must set the boundary conditions on r = 0 such as ur* *= 0 and ∂uz/∂r = 0 [5].*

In [7], the surfactant concentration equation in Lagrangian coordinates on the interface is presented, but it is in Cartesian coordinates. For testing a variety of axisymmetric cases, we should transform the equation in Cartesian coordinates into cylindrical coordinates according to [8].

Since our initial velocity field has to satisfy Navier-Stoke equations and the axisym-metric boundary conditions, the uniaxial extensional flow [9] is a good option. Surfactant on the interface drifts along the direction of the fluid flowing. Similarly, the markers in Lagrangian coordinates on the interface are affected by the flow. The following section will introduce a artificial method about marker redistribution.

For testing the case about moving contact line, put a drop on a solid on the slip bound-ary and control its equilibrium contact angle. Then, we impose the surfactant on the interface to observe that the drop slips with the time evolution. This simulation result is similar to that we use the detergent to remove grease. According to [12], we have Robin boundary condition on the slip boundary. We have to modify our solver for Navier-Stokes equations since this special boundary condition.

The remaining sections of this paper is organized as follows. In Section 2, we present the governing equations which includes the immersed boundary formulas and the sur-factant concentration equation in Lagrangian coordinates on the interface. Note that we impose artificial velocities on some equations to modify the distance between markers. The numerical method is described in Section 3 which includes an algorithm of solving the Navier-Stokes equations and a mass conservative scheme for the surfactant equation. We show our results about the effect of surfactant on drop deformation in some velocity fields in Section 4. In Section 5, we summarize some conclusions and our future work.

### 2 The governing equations

In a fixed three-dimensional cubic domain, we consider an incompressible two-phase flow problem and assume the interface is a simple closed surface immersed in the fluid

domain. The interface of the drop Σ separates two parts, Ω1 inside the interface and Ω2

*outside the interface from our domain. i.e. Ω = [a, b] × [c, d] × [e, f ] = Ω*1

S

Ω2. As

the interface is contaminated by surfactant, the distribution of the surfactant changes the
surface tension accordingly. In each fluid region, the Navier-Stokes equations are satisfied
as
ρ*i*
∂u*i*
*∂t* + (u*i*· ∇) u*i*
!
= ∇ · T*i*+ ρ*i*g, *in Ωi*, (2.1)
∇ · u*i* = 0, *in Ωi*, (2.2)
u = u*b*, *in ∂Ω,* (2.3)

*where for i = 1, 2 in each fluid domain, Ti* *= −pi*I + µ*i*(∇u*i* + ∇u*Ti*) is the stress tensor,

*pi* is the pressure, u*i* is the fluid velocity, ρ*i* is the density, µ*i* is the viscosity, and g is the

gravitational constant.

The velocity is continuous across the interface Σ. Thus,

[u]Σ = u|Σ,2− u|Σ,1 = 0 (2.4)

and the normal stress jump is balanced by the interfacial force F (defined only on Σ) as

[Tn]Σ+ F = 0 (2.5)

where n is the unit normal vector on Σ directed towards fluid 2. In this paper, we use the immersed boundary method to approach the velocity solutions of the Navier-Stokes equations. The immersed boundary ( the interface ) exerts force F to the fluids and moves with local fluid velocity. For computing easily, we first consider the case of equal viscosity

µ1 = µ2 = µ, equal density ρ1 = ρ2 = ρ, and neglect gravity. Hence, the interfacial force

term in delta function formulation and the surfactant concentration equation are the same problem in the single phase.

To be convenient to compute our numerical experiment in the single phase, we first non-dimensionlize the following three dimensional equations,

ρ ∂u
*∂t* + (u · ∇) u
!
*= −∇p + ∇ · (2µE) + f,* (2.6)
∇ · u = 0, (2.7)
∂Γ
*∂t* + (∇*s· u) Γ = Ds*∇
2
*s*Γ, (2.8)

*where u is the velocity, p is the pressure, ρ is the density, µ is the viscosity, Γ is the*

*surfactant concentration, Dsis diffusivity, and E = (∇u+∇uT*)/2 is the rate of deformation

*Assume the dimensionless numbers, x = r*0*x*∗*, u = U*∞u∗*, t = (r*0*/U*∞*)t*∗*, p = ρU*∞2 *p*∗,

*f = f*0f∗, and Γ = Γ∞Γ∗. Then, the non-dimensional equations are presented as follows.

∂u∗
*∂t*∗ + (u
∗_{· ∇}∗_{) u}∗ _{= −∇}∗* _{p}*∗

_{+}1

*Re*∇ ∗

*∗*

_{· (2E}_{) +}

*r*0

*f*0

*ρU*2 ∞ f∗

_{,}

_{(2.9)}∇∗

_{· u}∗

_{= 0,}

_{(2.10)}∂Γ∗

*∂t*∗ + ∇ ∗

*s*· u∗ Γ∗

_{=}1

*Pes*∇∗2

*s*Γ∗, (2.11)

*where Re = ρU*∞*r*0*/µ and Pes* *= r*0*U*∞*/Ds*. In the axisymmetric cylindrical coordinates,

*we spread the interfacial force into fluid by Dirac delta function δ(x) = δ(r)δ(z) as follows.*
*f (x, t) =*

Z

Σ

∂στ

*∂s* *δ (x − X (s, t)) ds,*

where x is the fluid position, X is the interfacial position, σ is the surface tension, τ is the
*tangential vector on the interface, s is the interfacial parameter, and Σ is the domain of*
parameter.

Assume the dimensionless number σ = σ0σ∗, then we have

Z
Σ
∂στ
*∂s* *δ (x − X (s, t)) ds =*
σ0
*r*2
0
Z
Σ
∂σ∗_{τ}
*∂s* δ (x
∗_{− X}∗_{) ds.}

*We suppose Ca = µU*∞/σ0and use ∇ · u = 0, then

∂u∗
*∂t*∗ + (u
∗_{· ∇}∗_{) u}∗ _{= −∇}∗* _{p}*∗

_{+}1

*Re*∆ ∗

_{(u}∗

_{) +}1

*ReCa*Z Σ ∂σ∗

_{τ}

*∂s*δ (x ∗

_{− X}∗

_{) ds,}_{(2.12)}∇∗

_{· u}∗

_{= 0,}

_{(2.13)}∂Γ∗

*∂t*∗ + ∇ ∗

*s*· u∗ Γ∗

_{=}1

*Pes*∇∗2

*s*Γ∗. (2.14)

*In fact, the characteristic velocity U*∞ *= Gr*0*where G is the principle strain rate and 2r*0is

*the equivalent diameter of the bubble. Hence, we have Re = ρ(Gr*0*)r*0*/µ, Ca = µGr*0/σ0,

*and Pes* *= r*0*U*∞*/Ds*. The following subsections will introduce the operator forms of

governing equations in the axisymmetric cylindrical coordinates.

### 2.1 Immersed boundary formulation

Throughout this paper, we consider the Navier-Stokes equations in cylindrical
coordi-nates instead of Cartesian coordicoordi-nates, and the domain is axisymmetric. It is convenient
*to compute our data in a cylindrical coordinate system (r, θ, z). Due to the azimuthal *
*sym-metry, the flow depends spatially on only two cylindrical coordinates(r, z). The interface*

*Σ is represented by a parametric form (r(s, t), z(s, t)), 0 ≤ s ≤ Lb, where s is the parameter*

of the initial configuration of the interface, which is not necessarily the arc-length. By using the non-dimensionlization process, we can write down our governing equa-tions in the usual immersed boundary formulation as follows.

*∂ur*
*∂t* *+ ur*
*∂ur*
*∂r* *+ uz*
*∂ur*
*∂z* +
*∂p*
*∂r* =
1
*Re*
∂2_{u}*r*
*∂r*2 +
1
*r*
*∂ur*
*∂r* +
∂2_{u}*r*
*∂z*2 −
*ur*
*r*2
!
+ 1
*ReCafr*, (2.15)
*∂uz*
*∂t* *+ ur*
*∂uz*
*∂r* *+ uz*
*∂uz*
*∂z* +
*∂p*
*∂z* =
1
*Re*
∂2_{u}*z*
*∂r*2 +
1
*r*
*∂uz*
*∂r* +
∂2_{u}*z*
*∂z*2
!
+ 1
*ReCafz*, (2.16)
1
*r*
*∂rur*
*∂r* +
*∂uz*
*∂z* = 0, (2.17)
*f (x, t) =*
Z
Σ
*F (s, t) δ (x − X (s, t)) ds,* (2.18)
*∂X (s, t)*
*∂t* *= U (s, t) = u (X (s, t) , t) =*
Z
Ω
*u (x, t) δ (x − X (s, t)) dx,* (2.19)
*F (s, t) =* ∂
*∂s(σ (s, t) τ (s, t)) ,* (2.20)
*τ (s, t) =*
∂X
*∂s*
∂X
*∂s*
, (2.21)

*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. The Eq.(2.18) represents the interfacial force F spreads to

*the fluid such that we get the force f = ( fr, fz*) on the fluid, and the force F on the interface

*is balanced by the normal stress. Here, u = (ur, uz) is the velocity ( ur*is the radial velocity

*and uz*is the axial velocity ), σ is the surface tension, and τ is the unit tangent vector on

the interface. The Eq.(2.19) represents the interface moves with the fluid velocity. The
Eulerian (x) and Lagrangian (X) variable are used in the three-dimensional Dirac delta
*function δ(x) = δ(r)δ(z).*

The interfacial force F arises from the surface tension and its form is derived from
Laplace-Young condition. If σ is not a constant, then we further take the derivatives in
Eq.(2.20) as follows.
*F (s, t) =* ∂
*∂s*(στ) =
∂σ
*∂s*τ + σ
∂τ
*∂s* =
∂σ
*∂s*τ + σκn
∂X* _{∂s}*
, (2.22)

where κ is the curvature of the interface and n is the unit outward normal. In Eq.(2.22),
*(∂σ/∂s) τ is the Marangoni force (the tangential force) and σκn |∂X/∂s| is the capillary*
force (the normal force). Note that if the surface tension is a constant, then the force
only exerts in the normal direction. Let the interface be contaminated by the surfactant

which reduces the surface tension, we can find the phenomenon that the surface tension is smaller as the surfactant concentration is higher. The relation between surface tension and surfactant concentration can be described by the non-linear Langmuir equation. The following is the approximation of non-linear Langmuir equation.

σ (Γ) = σ*c*(1 + ln (1 − βΓ)) , (2.23)

where Γ is the surfactant concentration, σ*c*is the surface tension of a clean interface, and

β is a dimensionless number that measures the sensitivity of surface tension to changes in surfactant concentration.

Since we want to close the system, the surfactant is insoluble to the buck fluids such that it is convected and diffused along the interface. Hence, there is no exchange between the interface and the bulk fluids. By those reasons, the total mass of the surfactant must be conserved. The equation of surfactant concentration is derived in next subsection.

### 2.2 Surfactant concentration equation

In order to make the total mass of the surfactant conserved, we present a slightly different derivation from Stone [10] for the surfactant transport equation in the three-dimensional cylindrical coordinates.

In axisymmetric cylindrical coordinates, we can regard the surface as the curve in r-z
*coordinates. Thus, let L(t) be an interfacial segment where the surfactant concentration*
(the mass of the surfactant per unit surface area) 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)*

*Γ (l, t) rdl = 0,* (2.24)

*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)*
*Γ (s, t)*_{}∂X
*∂s*
* rds = 0.* (2.25)

By taking the time derivative inside the integral, we obtain
Z
*L(0)*
∂Γ
*∂t*
∂X_{∂s}* r + Γ _{∂t}*∂
∂X

_{∂s}*r*! !

*ds = 0.*(2.26)

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

exactly the material derivative of Stones derivation. Now we need to compute the rate of the stretching factor, and using Eq.(2.19), we have

∂
*∂t*
∂X_{∂s}* r*
!
= *∂r*
*∂t*
∂X_{∂s}* + r _{∂t}*∂
∂X

*=*

_{∂s}*ur*

*r*

*r*∂X

_{∂s}*+ r*

*∂r*

*∂s∂s*∂

*∂r*

*∂t*+

*∂z*

*∂s∂s*∂

*∂z*

*∂t*∂X

*∂s*=

*ur*

*r*

*r*∂X

_{∂s}*+ r*

*∂r*

*∂s*

*∂ur*

*∂s*+

*∂z*

*∂s*

*∂uz*

*∂s*∂X

*∂s*(2.27) =

*ur*

*r*

*r*∂X

_{∂s}*+ r*

*∂r*

*∂s*

*∇ur*· ∂X

*+*

_{∂s}*∂z*

*∂s*

*∇uz*· ∂X

*∂X*

_{∂s}*∂s*=

*ur*

*r*

*r*∂X

*+ ∂u*

_{∂s}_{∂τ}· τ !

*r*

_{}∂X

*∂s*= (∇

*s· u) r*∂X

*.*

_{∂s}Here, the notation ∇*s·u = ur/r+(∂u/∂τ)·τ means the surface divergence in axisymmetric*

cylindrical coordinates. Since the material segment is arbitrary, we thus have
*r*_{}∂X
*∂s*
∂Γ_{∂t}*+ r*
∂X* _{∂s}*
( ∇

*s*· u ) Γ = 0. (2.28)

If we allow surfactant diffusion along the interface, we obtain the surfactant
transport-diffusion equation as
*r*_{}∂X
*∂s*
∂Γ_{∂t}*+ r*
∂X* _{∂s}*
( 5

*s*· u ) Γ = 1

*Pes*∂

*∂s*"

*r*∂Γ

*∂s*! , ∂X

*# (2.29)*

_{∂s}*where Pesis the surface Peclet number and the boundary condition of Γ is ∂Γ/∂s = 0 on*

*∂L (t) . We note that surface diffusion is also written in terms of initial parameter s.*

### 2.3 Controlling the interfacial markers uniformly

When there is the velocity field in the domain, the interface may be deformed ,and our computational markers on the interface concentrate on some part of interface with the direction of the fluid flowing. This phenomenon makes the larger error for the drop volume than the one which let the interfacial markers be distributed uniformly.

Since drop deformation is only related to the force in the normal direction, we impose an artificial velocity on Eq.(2.19) to redistribute the interfacial markers uniformly. Hence,

substitute Eq.(2.19) into
*∂X (s, t)*
*∂t* =
Z
Ω
*u (x, t) δ (x − X (s, t)) dx + UA(s, t) τ*
*= u (X (s, t) , t) + UA(s, t) τ*
*= U (s, t) + UA(s, t) τ,* (2.30)

*where the scalar UA* _{is the artificial velocity, and it can modify the marker positions such}

*that they avoid being moved by the tangential force. This condition ∂ |∂X/∂s| /∂s = 0*
means that the tangential displacement of markers remains unchanged. According to
those conditions, we obtain

*UA _{(s, t) − U}A_{(0, t) =}*

*s*2π Z

_{2π}0 ∂U

*∂s*0 · τ 0

*0*

_{ds}_{−}Z

*0 ∂U*

_{s}*∂s*0 · τ 0

*0*

_{ds}_{,}

_{(2.31)}

*where UA _{(0, t) = 0, since we should fix one marker to modify their relative positions.}*

After adding the artificial velocity term, the new marker positions have changed, so Eq.(2.29) has to add some term to preserve that the equation holds. The following is the modified surfactant transport-diffusion equation,

*r*_{}∂X
*∂s*
∂Γ_{∂t}*+ r*
∂X* _{∂s}*
( 5

*s*· u ) Γ − ∂

*∂s*

*rUA*Γ= 1

*Pes*∂

*∂s*"

*r*∂Γ

*∂s*! , ∂X

*# . (2.32)*

_{∂s}### 3 Numerical method

In this paper, the fluid flow variables are defined on a center of the cell which is divided
by uniform grid mesh; that is, the pressure is defined on the grid points labelled as x =
*(ri−1/2, zj−1/2) = ((i − 1/2)h, ( j − 1/2)h) for i, j = 1, 2, . . . , N, the velocity components*

*ur* *and uz* *are defined at (ri−1/2, zj−1/2) = ((i − 1/2)h, ( j − 1/2)h), where the spacing h =*

*∆r = ∆z . For the immersed interface, we use a collection of discrete points sk* *= k∆s, k =*

*1, 2, . . . , M such that the Lagrangian markers are denoted by X = X(sk) = (rk, zk*). The

surfactant concentration Γ*k* and surface tension σ*k* are defined at the half-integer points

*given by sk−1/2* *= (k − 1/2)∆s. Without loss of generality, for any function defined on the*

*interface φ(s), we approximate the partial derivative ∂φ/∂s by*

*Dsφ (s) =*

*φ (s + ∆s/2) − φ (s − ∆s/2)*

*∆s* . (3.1)

*Since |Ds*X*k*| can approximate the interface stretching factor by using this finite difference

convention, the unit tangent vector τ*k*are defined 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*

*k* *= X (sk, n∆t) , Γn _{k}*

*= Γ sk−1/2, n∆t*

, u*n* _{=}

*u (x, n∆t), and pn* _{= p (x, n∆t) are all given. The details of the numerical time integration}

1. Compute the surface tension and unit tangent vector on the interface as
σ*n _{k}* = σ

*c*1 + ln 1 − βΓ

*n*, (3.2) τ

_{k}*n*

*k*=

*Ds*X

*nk*

*Ds*X

*nk*, (3.3)

*both of which hold for sk−1/2* *= (k − 1/2)∆s. Then we define the interfacial force as*

F*n*

*k* *= Ds* σ*nk*τ*nk*

, (3.4)

at point X*k*.

2. Distribute the force from the markers to the fluid by

f*n*(x) =X

*k*

F*n _{k}*δ

*h*x − X

*nk*

*∆s,* (3.5)

where the smooth version of Dirac delta function,

δ*h(x) =*
( _{1}
*4h*
1 + cos*πx*
*2h*
*, if −2h ≤ x ≤ 2h*
0, otherwise, (3.6)
is used.

3. Solve the Navier-Stokes equations. This can be done by the following semi-implicit second-order projection scheme, where the nonlinear term is approximated by a second-order extrapolation to avoid solving a nonlinear system at each time step.

3 ˜u*n+1*_{− 4u}*n*_{+ u}*n−1*
*2∆t* +
2u*n*_{· ˜∇}
*h*
u*n*_{−}_{u}*n−1*_{· ˜∇}
*h*
u*n−1*_{+ ˜∇p}n_{=} 1
*Re*˜∆˜u
*n+1*_{+} 1
*ReCa*f
*n*_{,} _{(3.7)}
˜u = u*b, on ∂Ω,* (3.8)
˜∇*h*· u*n+1*= 0, (3.9)
˜∇2
*h*φ*n+1*=
3 ˜∇*h*· ˜u
*2∆t* ,
∂φ
∂n *= 0, on ∂Ω* (3.10)
˜∇*h*φ*n+1*=
3˜u*n+1*_{− u}*n+1*
*2∆t* , (3.11)

*pn+1 _{= p}n*

_{+ φ}

*n+1*

_{−}1

*Re*˜∇*h*· ˜u

*n+1*_{.} _{(3.12)}

*The boundary conditions of the velocities are ur* *= 0 and ∂uz/∂r = 0 on r = 0. Here*

the gradient operator ˜∇*h* *= (∂/∂r, ∂/∂z) , the divergent operator ˜∇h*· = (1* _{r}_{∂r}*∂

*r,*∂)· and

_{∂z}the Laplace operator
˜
∆*h*˜u*n+1*=
∂2_{u}n+1*r*
*∂r*2 +
1
*r*
*∂un+1*
*r*
*∂r* +
∂2_{u}n+1*r*
*∂z*2 −
*un+1*
*r*
*r*2 ,
∂2_{u}n+1*z*
*∂r*2 +
1
*r*
*∂un+1*
*z*
*∂r* +
∂2_{u}n+1*z*
*∂z*2
!
˜∇2
*h*φ*n+1* =
∂2_{φ}*n+1*
*∂r*2 +
1
*r*
∂φ*n+1*
*∂r* +
∂2_{φ}*n+1*
*∂z*2

are the standard centered difference operator on the grid which is defined at the center of the uniform grid mesh. One can see that the above Navier-Stokes solver

involves solving two Helmholtz equations for velocity ˜u*n+1*_{and one Poisson }

equa-tion for pressure. These elliptic equaequa-tions are solved by using the fast Poisson solver which is provided by the public software package Fishpack.

4. Interpolate the new velocity on the fluid lattice points onto the marker points and move the marker points to new positions.

U*n+1*
*k* =
X
x
u*n+1*_{δ}
*h*(x − X*nk)h*2, (3.13)
X*n+1 _{k}* = X

*n*

_{k}*+ ∆tUn+1*. (3.14)

_{k}In order to redistribute the interfacial markers more uniformly, we should add the
artificial velocity,
*UAn+1*
*k* =
*k∆s*
2π
*M*
X
*k*0_{=0}
*Ds*U*n+1k* · τ*n+1k* *∆s −*
*k*
X
*k*0_{=0}
*Ds*U*n+1k* · τ*n+1k* *∆s.* (3.15)

From the above equations, we have the modified marker positions,

X*n+1*

*k* = X*nk* *+ ∆tUn+1k* *+ ∆tUAn+1k* τ*n+1k* . (3.16)

5. Update surfactant concentration distribution Γ*n+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. Substitute Eq.(2.27) of rate of stretching factor into the equation(2.29), we have

*r*_{}∂X
*∂s*
∂Γ* _{∂t}* +

*∂ ∂X*

_{∂t}

_{∂s}*r*! Γ = 1

*Pes*∂

*∂s*"

*r*∂Γ

*∂s*! , ∂X

*# (3.17)*

_{∂s}Now we discretize the above equation by the Crank-Nicholson scheme in a
sym-metric way as
Γ*n+1*
*k* − Γ*nk*
*∆t*
*rn+1*
*k*
*Ds*X*n+1k* * + rnk*
*Ds*X*nk*
2 +
*rn+1*
*k*
*Ds*X*n+1k* * − rnk*
*Ds*X*nk*
*∆t*
Γ*n+1*
*k* + Γ*nk*
2
= 1
*Pes*
1
*∆s*
*r*
*n+1*
*k+1/2*
*|Ds*X*n+1k+1| + |Ds*X*n+1k* |
Γ*n+1*
*k+1*− Γ*n+1k*
*∆s* −
*rn+1*
*k−1/2*
*|Ds*X*n+1k* *| + |Ds*X*n+1k−1*|
Γ*n+1*
*k* − Γ*n+1k−1*
*∆s*
+ *r*
*n*
*k+1/2*
*Ds*X*nk+1* +*Ds*X*nk*
Γ*n*
*k+1*− Γ*nk*
*∆s* −
*rn*
*k−1/2*
*Ds*X*nk* +*Ds*X*nk−1*
Γ*n*
*k*− Γ*nk−1*
*∆s*
(3.18)
*Where k is “half-integer” index, i.e. k = 0.5, 1.5, . . . , M − 1/2. According to the*

new interfacial marker position X*n+1*

*k* obtained in the previous step and the boundary

condition (Γ*n+1*

*k* − Γ*n+1k−1)/∆s = 0, where k = 0.5 or M + 1/2, on the endpoints*

of the interface, Eq.(3.18) results 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*
Γ*n+1*
*k* *|Ds*X*n+1k* *|rn+1k−1/2∆s =*
X
*k*
Γ*n*
*k|Ds*X*nk|rnk−1/2∆s* (3.19)

The summation in Eq.(3.19) is the mid-point rule discretization for the integral in Eq.(2.25).The thing that the total mass is conserved can be derived by taking the summation of both sides of Eq.(3.18) and using the boundary condition.

Note that, when the artificial velocity is imposed on our equations like Eq.(3.16), Eq.(3.18) also have to impose the artificial velocity as follows.

Γ*n+1*
*k* − Γ*nk*
*∆t*
*rn+1*
*k*
*Ds*X*n+1k* * + rnk*
*Ds*X*nk*
2 +
*rn+1*
*k*
*Ds*X*n+1k* * − rnk*
*Ds*X*nk*
*∆t*
Γ*n+1*
*k* + Γ*nk*
2
−
*r*
*n+1*
*k+1/2UAn+1k+1/2*
Γ*n+1*
*k+1* + Γ*n+1k*
*/2 − rn+1*
*k−1/2UAn+1k−1/2*
Γ*n+1*
*k* + Γ*n+1k−1*
/2
*2∆s*
+*r*
*n*
*k+1/2UAnk+1/2*
Γ*n*
*k+1*+ Γ*nk*
*/2 − rn*
*k−1/2UAnk−1/2*
Γ*n*
*k* + Γ*nk−1*
/2
*2∆s*
= 1
*Pes*
1
*∆s*
*r*
*n+1*
*k+1/2*
*|Ds*X*n+1k+1| + |Ds*X*n+1k* |
Γ*n+1*
*k+1*− Γ*n+1k*
*∆s* −
*rn+1*
*k−1/2*
*|Ds*X*n+1k* *| + |Ds*X*n+1k−1*|
Γ*n+1*
*k* − Γ*n+1k−1*
*∆s*
+ *r*
*n*
*k+1/2*
*Ds*X*n _{k+1}* +

*Ds*X

*n*Γ

_{k}*n*

*k+1*− Γ

*nk*

*∆s*−

*rn*

*k−1/2*

*Ds*X

*n*+

_{k}*Ds*X

*n*Γ

_{k−1}*n*

*k*− Γ

*nk−1*

*∆s* .(3.20)

### 4 Numerical results

### 4.1 Oscillating drop

First, it is important to check how much the drop volume leaks or bulks up since
it affects the computational accuracy. We test a simple case that a spheroid drop in the
static liquid will be shrunk and turn into a sphere drop. Consider a computational domain
Ω = [0, 1] × [−1, 1] in the axisymmetrical cylindrical coordinates, and put a spheroid drop
*with the semi-major axis whose length is a = 0.75 and with the semi-minor axis whose*
*length is b = 1/3 into our domain. Then, given the computational mesh h = ∆r = ∆z,*
*∆t = h/8, Reynolds number Re = 25, and Capillary number Ca = 0.04 for our fluid.*

0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 0 , (V
a − Ve)/Ve = 0 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 0.25 , (V
a − Ve)/Ve = 0.12587 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 0.5 , (V
a − Ve)/Ve = 0.26542 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 0.75 , (V
a − Ve)/Ve = 0.36258 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 1 , (V_{a} − V_{e})/V_{e} = 0.35525 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 1.25 , (V_{a} − V_{e})/V_{e} = 0.40839 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 1.5 , (V_{a} − V_{e})/V_{e} = 0.50401 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 1.75 , (V_{a} − V_{e})/V_{e} = 0.6027 %

0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 2 , (V_{a} − V_{e})/V_{e} = 0.68024 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 2.25 , (V_{a} − V_{e})/V_{e} = 0.75274 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 2.5 , (V_{a} − V_{e})/V_{e} = 0.83533 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
r−axis
z−axis
Time = 2.75 , (V_{a} − V_{e})/V_{e} = 0.92554 %

Figure 4.1: In the fluid velocity field, we observe the time evolution of a spheroid drop in
*a static liquid with Re = 25, Ca = 0.04 and ∆t = h/8.*

In figure 4.1, the oscillation of the clean drop depends on the fluid velocity field and
*the incompressibility of the drop. The spheroid drop with a = 0.75 and b = 1/3 oscillates*
*from the time T = 0 to T = 2, and its steady state is the sphere drop with the radius*
*r 0.4368. After T = 2, the amplitude of vibration is smaller than the one before T = 2,*
but the drop volume leaks with the time evolution. The reason of this phenomenon is
that there is the spurious velocities in our computational process. The spurious velocities
influence the computation of the new marker positions on the interface .Thus, we only do
our best to reduce the relative error of the drop volume. The relative error of the drop
*volume is smaller than 1% at T = 2.75 in figure 4.1. So, This is an acceptable error range*
for us.

### 4.2 Convergence test of fluid velocity and surfactant concentration

After above subsection, we impose the surfactant on the interface to observe the
drop behavior. Throughout this paper, we choose β = 1 for the surface tension equation
Eq.(3.2). In order to make our research convinced, we should carry out the convergence
*study of the immersed boundary method. We use the different mesh h = ∆r = ∆z = 1/32,*
1/64, 1/128, 1/256 to perform our computation relatively. The Lagrangian mesh is
*cho-sen as ∆s ≈ h and the time step size is ∆t = h/8. The solutions are computed up to time*
*T = 1. Since it is not easy to obtain the analytical solution from our problem, we choose*

*the finest mesh as our reference solution to compute the L*∞ error between the reference

*solution and the solution which is solved from the coarser grid. Here, we choose Re = 1,*

*Table 1: The mesh refinement analysis of the velocity ur, uz*, and the surfactant

concen-tration Γ.

h *kur− URre f*k∞ Rate *kuz− UZ re f*k∞ Rate kΓ − Γ*re f*k∞ Rate

1/32 *8.5478E − 03 -* *1.1470E − 02 -* *1.2303E − 02 *

-1/64 *2.4752E − 03 1.788 4.6027E − 03 1.317 2.4629E − 03 2.321*

*1/128 7.5121E − 04 1.720 1.8433E − 03 1.320 4.9265E − 04 2.322*

*Table 1 shows that the error of kur− URre f*k∞*, kuz− UZ re f*k∞, and kΓ − Γ*re f*k∞converge

*to zero where URre f, UZ re f*, and Γ*re f* are r-axial, z-axial reference velocities, and reference

*surfactant concentration relatively. The convergent rates for ur, uz*, and Γ are about 1.7,

1.3, and 2.3 relatively, but our method is a second order scheme. Note that since the fluid
*velocities are defined at the center of the uniform grid in r − z coordinates, the refined*
mesh will not coincide with the same grid locations. Thus, we approximate the coarser
grid by averaging the velocities obtained from the reference grid near the coarser grid. We
attribute this is part of the reason why the rate of convergence behaves less than
second-order. Similarly, the surfactant concentration is defined at “half-integer” grid. When we
average the reference grid, we should use the weighted average since the concentration is
the mass per unit surface area.

### 4.3 Clean vs. contaminated interface

There is different behavior for clean interface and contaminated interface, so we will study it in this subsection. We compare a clean drop with a contaminated drop in an

*extensional flow, ur* *= −0.5Gr and uz* *= Gz, where G is the principal strain rate [5]. By*

the comparison, we observe that the surface tension reduces when the surfactant is on the
*interface. Here, we use the mesh h = ∆x = ∆y = 1/64, the Lagrangian grid with size*
*∆s ≈ h, and the time step size ∆t = h/8.*

4.3.1 The steady state of the clean interface in the extensional flow

In general, a drop will be stretched along the direction of fluid flow if it is put into the
extensional flow. But there is surface tension on the interface, it may restrict a tendency of
stretching drop. Therefore, before the comparison the clean interface with contaminated
interface, we first observe a property for the clean drop. We fix a small Reynolds number,
*Re = 1, such that our Navier-Stokes equations are similar to the Stokes equations. Then,*
we attempt to change Capillary number and get the numerical result that the drop achieves
a steady state.

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.4 r−axis time =0 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.43685 r−axis time =0.25 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.45546 r−axis time =0.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.46548 r−axis time =0.75 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47121 r−axis time =1 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47455 r−axis time =1.25 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47645 r−axis time =1.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.47748 r−axis time =1.75

Figure 4.2: The drop will form a steady state when the Capillary number is smaller than
*the critical value. Here, the dotted line is the initial interface and Ca = 0.2.*

*If Capillary number is 0.2 such as figure 4.2, then the drop with radius r = 0.4 is*
*stretched slowly at about T = 0.5 ,and the tip points of the interface on the z-axis are*
*about at the positions (r, z) ≈ (0, ±0.47). In figure 4.3, Capillary number is 0.27, and we*
*find that this adjustment can not decelerate the stretching drop after T = 2. The reason*
is that Capillary number has a critical value in this case and it is between 0.2 and 0.27.
Hence, the stretching drop will achieve a steady state quickly as the Capillary number is
smaller than the critical value. But the stretching drop will not stop when the Capillary

*to the viscosity µ, the strain rate G, the drop’s radius r, and the clean surface tension*

σ*c*. Hence, the larger strain rate and smaller surface tension can strengthen the stretching

phenomenon. 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.4 r−axis time =0 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.4949 r−axis time =0.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.54999 r−axis time =1 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.59339 r−axis time =1.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.63507 r−axis time =2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.68175 r−axis time =2.5 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.74252 r−axis time =3 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis , z−right−point = 0.83688 r−axis time =3.5

Figure 4.3: If the Capillary number is larger than the critical value, then the drop continue
*to be stretched after T = 2. Here, the dotted line is the initial interface and Ca = 0.27.*

4.3.2 The effect for the controlled markers on the interface

The initial markers on the interface are defined at the uniform grid, but the markers will
*not be distributed uniformly after the drop is put into the extensional flow with Re = 1*
*and Ca = 0.25. Most of interfacial markers concentrate in the vicinity of the tip points*

*of the drop as the fluid flows along the z-axial direction. For accurate computation, we*
impose an artificial velocity on Eq.(2.19) and Eq.(2.29) to redistribute interfacial marker
positions uniformly.
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time =0, (Va − Ve)/Ve = 0 %
no surfactant
Pes = 0.125
Pes = 12.5
Pes = 1250
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time =0.5, (Va − Ve)/Ve = 0.028247 %
no surfactant
Pes = 0.125
Pes = 12.5
Pes = 1250
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time =1, (Va − Ve)/Ve = 0.047636 %
no surfactant
Pes = 0.125
Pes = 12.5
Pes = 1250
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time =1.5, (Va − Ve)/Ve = 0.10333 %
no surfactant
Pes = 0.125
Pes = 12.5
Pes = 1250
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time =2, (Va − Ve)/Ve = 0.23924 %
no surfactant
Pes = 0.125
Pes = 12.5
Pes = 1250

Figure 4.4: In an extensional flow, compare the clean drop with the contaminated drops whose Peclet numbers are different. Note that the markers on interface are not distributed uniformly. 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 2 Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.5: The time evolution of the surfactant dynamics in a extensional flow. Note that the markers on interface are not distributed uniformly.

*The time and the relative volume error of drop with Pes* = 12.5 is on the top of figure

4.4 and figure 4.6. In figure 4.4, the relative volume error of the drop is 0.23924% at
*T = 2, and it is larger than 0.024377% which is the relative volume error of the drop in*
figure 4.6. So, redistributing markers on the interface uniformly can reduce the volume
error. Both in figure 4.4 and figure 4.6, the clean drop reaches a steady state as the
Capillary number is smaller than the critical value.

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =0, (Va − Ve)/Ve = 0 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =0.5, (Va − Ve)/Ve = 0.014915 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =1, (Va − Ve)/Ve = 0.038418 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =1.5, (Va − Ve)/Ve = 0.036789 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time =2, (Va − Ve)/Ve = 0.024377 % no surfactant Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.6: In an extensional flow, compare the clean drop with the contaminated drops whose Peclet numbers are different. Note that the markers on interface are distributed uniformly. 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 0.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 1.5 Pes = 0.125 Pes = 12.5 Pes = 1250 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 arc−length surfactant concentration time = 2 Pes = 0.125 Pes = 12.5 Pes = 1250

Figure 4.7: The time evolution of the surfactant dynamics in a extensional flow. Note that the markers on interface are distributed uniformly.

We put a drop into a uniaxial extensional flow such that the drop is stretched by the flow. In order to find the differences from the clean drop and the contaminated drop, we impose the surfactant on the clean interface and its surfactant concentration is Γ = 0.5.

*There are three Peclet numbers, Pes* *= 0.125, Pes* *= 12.5, and Pes* = 1250, in figure

4.6 and figure 4.7. The concentration curves has different shapes for different Peclet numbers. Whether we observe figure 4.4 or figure 4.6, we find that the interface move more slowly as the Peclet number is larger. The reason is that the high Peclet number reduce the effect of the surfactant diffusion. According to the direction of the fluid flow, the surfactant concentration in the vicinity of the tip points of the drop is higher than others. The evidence is in figure 4.5 and figure 4.7.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.392 0.3921 0.3922 0.3923 0.3924 0.3925 0.3926 0.3927 0.3928 0.3929 0.393 time

total mass of surfactant

*Figure 4.8: This is the total mass of surfactant for Pes* *= 0.125, Pes* *= 12.5, and Pes* =

1250.

The diffusion term in Eq.(2.29) influences that the surfactant diffuses from the high concentration to the low concentration. Therefore, the surfactant is carried by the fluid flowing, but the diffusion effect is very small if Peclet number is very large. In figure 4.7, it is more obvious that the surfactant concentrate in the vicinity of the interfacial endpoints as Peclet number is larger. Let’s attempt to explain why the interface move more slowly as the Peclet number is larger. It is known for us that the surfactant reduce the surface tension on the interface. Smaller Peclet number makes the surfactant distributed uniformly, so the surface tension is reduced averagely. But larger Peclet number makes the surfactant concentrate in the vicinity of the interfacial endpoints. Most of surface tension on the interface only reduce a little, but the higher concentration near the interfacial endpoints makes the tip points of the drop become sharper. By the comparison, the interface move more slowly as the Peclet number is larger. In figure 4.8, no matter which Peclet number we choose, the total mass of surfactant is conserved very well by our method.

### 4.4 Forming the dumbbell-shaped drop in another fluid velocity field

In order to ensure our scheme and concept, we want to observe whether there is a similar numerical result for the drop in another velocity field. Consider the velocity field

*, ur* *= −0.5Gπr cos(πz) and uz* *= G sin(πz), which satisfies Eq.(2.17) and the boundary*

*condition ur* *= 0, ∂uz/∂r = 0 on r = 0 in the axisymmetric cylindrical coordinates. Put a*

*drop into the velocity field with Re = 1 and Ca = 0.5. Then we impose the surfactant on*
the interface. In the following result, we have imposed the artificial velocity to uniformly
redistribute the marker positions.

0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time = 0, (V_{a} − V_{e})/V_{e} = 0 %
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
z−axis
r−axis
time = 3, (V_{a} − V_{e})/V_{e} = 0.38462 %
0 0.5 1 1.5 2 2.5 3
1.5706
1.5706
1.5706
1.5706
1.5706
1.5706
1.5706
1.5706
1.5706
1.5706
time

total mass of surfactant

*Figure 4.9: The left two figures are the process that the drop is stretched from T = 0 to*
*T = 3. The total mass of surfactant is on the right hand side.*

*Observing the left two subfigures in figure 4.9, the drop with Γ = 0.5 and Pes* = 50

*is stretched along z = −1 and z = 1 in the velocity field. At T = 3, the surfactant*
concentration is higher near the stretched parts than other parts of the interface. The
reason is that the surfactant is carried along the direction of the fluid flowing. We can
compare the drops with different Peclet numbers in figure 4.10. The stretched part of the

*drop with Pes* *= 50 is sharper than the one with Pes* = 0.5 since there is more surfactant

concentrating on the tip parts. In addition, the clean drop is more smooth than others.

0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 0 no surfactant Pes = 0.5 Pes = 50 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 1 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 z−axis r−axis time = 3

0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 0 Pes = 0.5 Pes = 50 0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 1 Pes = 0.5 Pes = 50 0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 2 Pes = 0.5 Pes = 50 0 1 2 3 0 1 2 3 4 5 6 arc−length surfactant concentration time = 3 Pes = 0.5 Pes = 50

Figure 4.11: The time evolution of the surfactant dynamics in this velocity field.

*For the curve of surfactant concentration with Pes* = 0.5 in figure 4.11, it is more

*smooth and more straight than the one with Pes* = 50 since larger Peclet number reduces

the diffusion term in Eq.(2.29). On the right hand side of figure 4.9, the total mass of surfactant is conserved. Hence, this result is similar to the one about the drop in the extensional flow.

### 4.5 Moving contact line on the slip boundary

In this case, we put a clean drop onto a solid, and the initial contact angle is a right
*angle. We restrict the computational domain Ω = [0, 1] × [0, 1] in r − z coordinates. In*
axisymmetric geometries, most of boundary conditions are the same as the previous cases

*except the slip boundary conditions uz* *= 0 and ur* *= (∂ur/∂z) on z = 0, where is a*

coefficient about the slip length. Instead of the solver for the software package Fishpack, we solve Eq.(2.15) by the IMSL subroutine “DLSLXD” since the boundary conditions on

*z = 0 have changed. According to Young’s condition σ cos θs*+σ*sl* = σ*sg*, we get the force

*Fc* = σ*c*(cos θ*s*− cos θ*d) at moving contact line along r-axis where θs*is the equilibrium

contact angle and θ*d*is the dynamic contact angle.

4.5.1 Convergence test of fluid velocity and surfactant concentration

First, we should carry out the convergence study of the immersed boundary method.
*We use the different mesh h = ∆r = ∆z = 1/32, 1/64, 1/128, 1/256 to perform our*
*computation relatively. The Lagrangian mesh is chosen as ∆s ≈ h/2 and the time step size*
*is ∆t = h/8. The solutions are computed up to time T = 1. Since it is not easy to obtain the*

analytical solution from our problem, we choose the finest mesh as our reference solution

*to compute the L*∞ error between the reference solution and the solution which is solved

*from the coarser grid. Here, we choose Re = 1, Ca = 0.5, Pes* = 50, θ*s* = 60◦*, and = h*

for our computation.

*Table 2: The mesh refinement analysis of the velocity ur, uz*, and the surfactant

concen-tration Γ.

h *kur− URre f*k∞ Rate *kuz− UZ re f*k∞ Rate kΓ − Γ*re f*k∞ Rate

1/32 *1.0782E − 03 -* *2.9943E − 03 -* *2.7585E − 03 *

-1/64 *4.3824E − 04 1.299 1.3635E − 03 1.135 1.0808E − 03 1.352*

*1/128 1.5433E − 04 1.506 5.0570E − 04 1.431 3.7413E − 04 1.531*

*Table 2 shows that the error of kur− URre f*k∞*, kuz− UZ re f*k∞, and kΓ − Γ*re f*k∞converge

*to zero where URre f, UZ re f*, and Γ*re f* are r-axial, z-axial reference velocities, and reference

*surfactant concentration relatively. The convergent rates for ur, uz*, and Γ are about 1.5,

1.4, and 1.5 relatively, but our method is a second order scheme. The reason is that it produces some errors as the data location which we want are interpolated by the reference data.

4.5.2 Clean vs. contaminated interface with θs = 60◦

In this subsection, we want to compare the clean interface with contaminated

inter-face as the equilibrium contact angle is 60◦_{. Since the surfactant is insoluble and only}

distributed on the liquid-to-liquid interface, σ*sl* and σ*sg*are not affected by the surfactant.

*Therefore, the force Fc* at the moving contact line must be different from the one on the

*clean interface. That is Fc* = σ*sg*− σ*sl*− σ cos θ*d* where σ*sg*− σ*sl* = σ*c*cos θ*s*.

*In order to achieve the equilibrium contact angle quickly, we set = h = 1/128,*

*∆t = h/8, Re = 100, and Weber number We = ReCa = 0.1. We set Pes* = 500 such that

there is the obvious fluctuation on the surfactant concentration curve. In figure 4.12, we
*find that the tip points of interfaces on z-axis move slowly after about T = 4. The tip point*
*of contaminated interface on z-axis is farther from the original point than The tip point*

*of clean interface on z-axis. It is imaginable since the term σ cos θd* in the force on the

moving contact line is smaller than σ*c*cos θ*s*. Hence, θ*d* must be so small (θ*d* < θ*s*) that it

can achieve the equilibrium and σ*c*cos θ*s* = σ cos θ*d*.

The fluctuation of concentration curve is presented in figure 4.13. We can compare
this figure with figure 4.12. Note that the right endpoint of the curve is defined at the point
which is closest to the moving contact line. In the beginning, the surfactant concentrate
in the vicinity of the moving contact line since the slipping drop produce the similar
*velocity field as the extensional flow out along r-axis. After T = 4, the concentration*
curve becomes straight gradually since the drop becomes a steady state gradually.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 8

Figure 4.12: The drop is on the slip boundary. Clean drop (solid line) and contaminated drop (‘− · − · −·’) are in this figure.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 8

0 1 2 3 4 5 6 7 8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8

Relative Error of Volume : %

Time 0 1 2 3 4 5 6 7 8 0.4707 0.4708 0.4709 0.471 0.4711 0.4712 0.4713 0.4714 0.4715 0.4716 0.4717 Surfactant Mass Time 0 1 2 3 4 5 6 7 8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Cosine Value of the Dynamic Contact Angle

Time 0 1 2 3 4 5 6 7 8 40 45 50 55 60 65 70 75 80 85 90

Dynamic Contact Angle

Time

Figure 4.14: This figure presents the relative volume error, the total mass of surfactant,

the variation of cos θ*d*, and the variation of θ*d* with time.(solid line : clean drop, ‘− · − · −·’

: contaminated drop)

*In figure 4.14, the relative volume error is smaller than 1% from T = 0 to T = 8. The*
total mass of surfactant is conserved by our method. For the clean drop, the contact angle

is about 60◦* _{after T = 4, but the contact angle is about 45}*◦

_{for the contaminated drop since}

the surfactant has changed the original equilibrium contact angle. In fact, our Reynolds number is so large that it influences the smoothness for the curves about concentration and contact angle.

4.5.3 Clean vs. contaminated interface with θs = 120◦

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 12

Figure 4.15: The drop is on the slip boundary. Clean drop (solid line) and contaminated drop (‘− · − · −·’) are in this figure.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 Maker arc−length surfactant concentration time = 12

0 2 4 6 8 10 12 −1 0 1 2 3 4 5 6

Relative Error of Volume : %

Time 0 1 2 3 4 5 6 7 8 0.4707 0.4708 0.4709 0.471 0.4711 0.4712 0.4713 0.4714 0.4715 0.4716 0.4717 Surfactant Mass Time 0 2 4 6 8 10 12 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1

Cosine Value of the Dynamic Contact Angle

Time 0 2 4 6 8 10 12 80 90 100 110 120 130 140

Dynamic Contact Angle

Time

Figure 4.17: This figure presents the relative volume error, the total mass of surfactant,

the variation of cos θ*d*, and the variation of θ*d* with time.(solid line : clean drop, ‘− · − · −·’

: contaminated drop)

*In figure 4.15, we find that the tip points of interfaces on z-axis move slowly after*
*about T = 6. The tip point of contaminated interface on z-axis is closer to the original*

*point than The tip point of clean interface on z-axis. In this case, the term σ cos θd* in the

force on the moving contact line must be equal to σ*c*cos θ*s* to achieve the steady state.

Hence, the new equilibrium contact angle is larger than original contact angle after the interface is contaminated.

The fluctuation of concentration curve is presented in figure 4.16. We can compare this figure with figure 4.15. Note that the right endpoint of the curve is defined at the point which is closest to the moving contact line. In the beginning, the surfactant concentrate in the vicinity of the moving contact line since the slipping drop produce the similar

*velocity field as the extensional flow out along z-axis. After T = 6, the concentration*
curve becomes straight gradually since the drop becomes a steady state gradually.

*In figure 4.17, the relative volume error is larger than 1% after T = 2, and it is over*
*5% at T = 12. This result is a little bad for us. The reason may be that in order to achieve*
the equilibrium contact angle quickly, larger Reynolds number influences the accuracy.
The total mass of surfactant is conserved by our method. For the clean drop, the contact

angle is about 120◦* _{after T = 4, but the contact angle is about 135}*◦

_{for the contaminated}

drop since the surfactant has changed the original equilibrium contact angle. In fact, our Reynolds number is so large that it influences the smoothness for the curves about concentration and contact angle.

4.5.4 Imposing initial velocity on our moving contact line case

*In this subsection, we impose initial velocity, ur* *= Grz*2*and uz= −2Gz*3/3, on our case.

This velocity satisfies the slip boundary condition, axisymmetric boundary condition, and ∇ · u = 0. In order to reduce the effect of the velocity field to make contact line moving,

*we choose G = 0.5 and the equilibrium contact angle θs* = 0◦. Given initial contact angle

120◦_{, Re = 10, We = 0.1 and Pe}

*s*= 625 for the drop and the velocity field. Compare the

clean drop with the contaminated drop as follows.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 12

Figure 4.18: The drop is on the slip boundary. Clean drop (solid line) and contaminated drop (‘− · − · −·’) are in this figure.

In figure 4.18, the clean drop is very close to the contaminated drop since the shear rate is smaller than other previous cases. But the contact angle of the contaminated interface is smaller than the one of clean interface since the surfactant reduces the contact angle more quickly than the clean one. We can observe the phenomenon more clearly in figure 4.19. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z−axis r−axis time = 12

*Figure 4.19: It is an enlargement in figure 4.18 at T = 12.*

In figure 4.20, the surfactant concentration at moving contact line is higher than the
other parts near moving contact line at first. The reason is that the moving contact line
*slips along the direction of positive r−axis, and the fluid flows along the same direction.*
*Therefore, the surfactant concentrate in the vicinity of moving contact line. After T = 9,*
the surfactant diffuses from moving contact line to other parts since the contact angle is
going to achieve the equilibrium contact angle.

*In figure 4.21, the relative volume error is about 1.5% at T = 12. The total mass of*

surfactant is conserved by our method. For the clean drop, the contact angle is about 35◦

*at T = 12, but the contact angle is about 0*◦_{for the contaminated drop since the surfactant}

reduces the surface tension such that the contact angle at moving contact line achieves the equilibrium contact angle more quickly.

In fact, if we choose the larger shear rate, then the effect of the fluid velocity field is greater than the effect of the slip boundary. Hence, it is difficult to see that the contact angle at moving contact line achieves the equilibrium contact angle.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 arc−length surfactant concentration time = 12

Figure 4.20: The time evolution of surfactant concentration depended on figure 4.18.

0 2 4 6 8 10 12 −1.5

−1 −0.5

0 Relative Error of Volume : %

Time 0 2 4 6 8 10 12 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 0.1767 Surfactant Mass Time 0 2 4 6 8 10 12 −0.5 0 0.5 1

Cosine Value of the Dynamic Contact Angle

Time 0 2 4 6 8 10 12 0 20 40 60 80 100 120

Dynamic Contact Angle

Time

Figure 4.21: This figure presents the relative volume error, the total mass of surfactant,

the variation of cos θ*d*, and the variation of θ*d* with time.(solid line : clean drop, ‘− · − · −·’