國
立
交
通
大
學
應用數學系
碩
士
論
文
模擬帶有表面活性劑之軸對稱界面流問題
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 HuangAdvisors：Dr.
Ming-Chih LaiDepartment﹙Institute﹚of
Applied MathematicsNational 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 ∂ui ∂t + (ui· ∇) ui ! = ∇ · Ti+ ρig, in Ωi, (2.1) ∇ · ui = 0, in Ωi, (2.2) u = ub, in ∂Ω, (2.3)
where for i = 1, 2 in each fluid domain, Ti = −piI + µi(∇ui + ∇uTi) is the stress tensor,
pi is the pressure, ui 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 = r0x∗, u = U∞u∗, t = (r0/U∞)t∗, p = ρU∞2 p∗,
f = f0f∗, and Γ = Γ∞Γ∗. Then, the non-dimensional equations are presented as follows.
∂u∗ ∂t∗ + (u ∗· ∇∗) u∗ = −∇∗p∗+ 1 Re∇ ∗· (2E∗) + r0f0 ρU2 ∞ f∗, (2.9) ∇∗· u∗ = 0, (2.10) ∂Γ∗ ∂t∗ + ∇ ∗ s· u∗ Γ∗ = 1 Pes ∇∗2 s Γ∗, (2.11)
where Re = ρU∞r0/µ and Pes = r0U∞/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 r2 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∞ = Gr0where G is the principle strain rate and 2r0is
the equivalent diameter of the bubble. Hence, we have Re = ρ(Gr0)r0/µ, Ca = µGr0/σ0,
and Pes = r0U∞/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 ∂2u r ∂r2 + 1 r ∂ur ∂r + ∂2u r ∂z2 − ur r2 ! + 1 ReCafr, (2.15) ∂uz ∂t + ur ∂uz ∂r + uz ∂uz ∂z + ∂p ∂z = 1 Re ∂2u z ∂r2 + 1 r ∂uz ∂r + ∂2u z ∂z2 ! + 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 ( uris the radial velocity
and uzis 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, σcis 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∂s ∂X ∂s = ur r r ∂X∂s + ∂u∂τ · τ ! 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 ( 5s· u ) Γ = 1 Pes ∂ ∂s " r∂Γ ∂s ! , ∂X∂s # (2.29)
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) − UA(0, t) = s 2π Z 2π 0 ∂U ∂s0 · τ 0ds0− Z s 0 ∂U ∂s0 · τ 0ds0, (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 ( 5s· u ) Γ − ∂ ∂s rUAΓ= 1 Pes ∂ ∂s " r∂Γ ∂s ! , ∂X∂s # . (2.32)
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 |DsXk| can approximate the interface stretching factor by using this finite difference
convention, the unit tangent vector τkare 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) , Γnk = Γ sk−1/2, n∆t
, un =
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 σnk = σc 1 + ln 1 − βΓnk , (3.2) τn k = DsXnk DsXnk , (3.3)
both of which hold for sk−1/2 = (k − 1/2)∆s. Then we define the interfacial force as
Fn
k = Ds σnkτnk
, (3.4)
at point Xk.
2. Distribute the force from the markers to the fluid by
fn(x) =X
k
Fnkδh x − Xnk
∆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 ˜un+1− 4un+ un−1 2∆t + 2un· ˜∇ h un−un−1· ˜∇ h un−1+ ˜∇pn = 1 Re˜∆˜u n+1+ 1 ReCaf n, (3.7) ˜u = ub, on ∂Ω, (3.8) ˜∇h· un+1= 0, (3.9) ˜∇2 hφn+1= 3 ˜∇h· ˜u 2∆t , ∂φ ∂n = 0, on ∂Ω (3.10) ˜∇hφn+1= 3˜un+1− un+1 2∆t , (3.11)
pn+1= pn+ φ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· = (1r∂r∂r,∂z∂)· and
the Laplace operator ˜ ∆h˜un+1= ∂2un+1 r ∂r2 + 1 r ∂un+1 r ∂r + ∂2un+1 r ∂z2 − un+1 r r2 , ∂2un+1 z ∂r2 + 1 r ∂un+1 z ∂r + ∂2un+1 z ∂z2 ! ˜∇2 hφn+1 = ∂2φn+1 ∂r2 + 1 r ∂φn+1 ∂r + ∂2φn+1 ∂z2
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 ˜un+1and 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.
Un+1 k = X x un+1δ h(x − Xnk)h2, (3.13) Xn+1k = Xnk + ∆tUn+1k . (3.14)
In order to redistribute the interfacial markers more uniformly, we should add the artificial velocity, UAn+1 k = k∆s 2π M X k0=0 DsUn+1k · τn+1k ∆s − k X k0=0 DsUn+1k · τn+1k ∆s. (3.15)
From the above equations, we have the modified marker positions,
Xn+1
k = Xnk + ∆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 + ∂t∂ ∂X∂s r ! Γ = 1 Pes ∂ ∂s " r∂Γ ∂s ! , ∂X∂s # (3.17)
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 DsXn+1k + rnk DsXnk 2 + rn+1 k DsXn+1k − rnk DsXnk ∆t Γn+1 k + Γnk 2 = 1 Pes 1 ∆s r n+1 k+1/2 |DsXn+1k+1| + |DsXn+1k | Γn+1 k+1− Γn+1k ∆s − rn+1 k−1/2 |DsXn+1k | + |DsXn+1k−1| Γn+1 k − Γn+1k−1 ∆s + r n k+1/2 DsXnk+1 + DsXnk Γn k+1− Γnk ∆s − rn k−1/2 DsXnk + DsXnk−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 Xn+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 |DsXn+1k |rn+1k−1/2∆s = X k Γn k|DsXnk|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 DsXn+1k + rnk DsXnk 2 + rn+1 k DsXn+1k − rnk DsXnk ∆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 |DsXn+1k+1| + |DsXn+1k | Γn+1 k+1− Γn+1k ∆s − rn+1 k−1/2 |DsXn+1k | + |DsXn+1k−1| Γn+1 k − Γn+1k−1 ∆s + r n k+1/2 DsXnk+1 + DsXnk Γn k+1− Γnk ∆s − rn k−1/2 DsXnk + DsXnk−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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 , (Va − Ve)/Ve = 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 fk∞ Rate kuz− UZ re fk∞ Rate kΓ − Γre fk∞ 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 fk∞, kuz− UZ re fk∞, and kΓ − Γre fk∞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, (Va − 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 z−axis r−axis time = 3, (Va − Ve)/Ve = 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 θsis the equilibrium
contact angle and θdis 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 fk∞ Rate kuz− UZ re fk∞ Rate kΓ − Γre fk∞ 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 fk∞, kuz− UZ re fk∞, and kΓ − Γre fk∞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 σsgare 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 = σccos θ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 σccos θs. Hence, θd must be so small (θd < θs) that it
can achieve the equilibrium and σccos θ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 σccos θ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 = Grz2and uz= −2Gz3/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, ‘− · − · −·’
