An Immersed Boundary Method for Simulating Interfacial Flows with Insoluble Surfactant in Three Dimensions
Yunchang Seol1,∗, Shih-Hsuan Hsu2and Ming-Chih Lai2
1National Center for Theoretical Sciences, No. 1, Sec. 4, Road. Roosevelt, National Taiwan University, Taipei 10617, Taiwan.
2Department of Applied Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan.
Received 12 April 2017; Accepted (in revised version) 29 June 2017
Abstract. In this paper, an immersed boundary (IB) method for simulating the inter- facial flows with insoluble surfactant in three dimensions is developed. We consider a doubly periodic interface separating two fluids where the surfactant exists only along the evolving interface. An equi-arclength parametrization is introduced in order to track the moving interface and maintain good Lagrangian meshes, so stable compu- tations can be performed without remeshing. This surface mesh-control technique is done by adding two artificial tangential velocity components into the Lagrangian marker velocity so that the Lagrangian markers can be equi-arclength distributed dur- ing the time evolution. As a result, the surfactant equation on the interface must be modified based on the new parametrization. A conservative scheme for solving the modified surfactant equation has been developed and proved to satisfy the total sur- factant mass exactly in discrete level. A series of numerical experiments consisting of the validation of Lagrangian mesh control technique, the convergence study, the study of self-healing dynamics, and the simulations of two-layer fluids under Couette flow have been conducted to test our present numerical scheme.
AMS subject classifications: 35Q35, 65M06, 76D45
Key words: Insoluble surfactant, a conservative scheme, mesh control, interfacial flow, Navier- Stokes flow, immersed boundary method.
1 Introduction
Surfactant is a typical compound of amphiphilic molecules with hydrophilic heads and hydrophobic tails. It adheres to the fluid interface and changes the surface tension ac- cordingly. So, in the presence of surfactant, the dynamics of interfacial flows can be
∗Corresponding author. Email addresses: ycseol@ntu.edu.tw (Y. Seol), narisawa.am02g@g2.nctu.edu.tw (S.-H. Hsu), mclai@math.nctu.edu.tw (M.-C. Lai)
http://www.global-sci.com/ 640 2018 Global-Science Pressc
dramatically altered in complex ways due to the variable surface tension. For instance, a recent work of Pan et al. [37] demonstrates that by simply adding a small amount of sur- factant, one can manipulate the impact outcome of droplet collision from coalescence to bouncing. Surfactant also plays a crucial role in numerous industrial applications includ- ing oil recovery, ore extraction, emulsification, pharmaceutical, and food industry [24], etc. Nowadays, we are benefited from scientific success with applications to its practical use; however, the understanding of surfactant effects is far from being complete.
For the past two decades, the interfacial flows with insoluble and soluble surfactant have been extensively studied by numerical computations. Those numerical methods are strongly relevant to how the interface is represented and how the surfactant convection- diffusion equation is solved along the evolving interface. For the soluble surfactant case, another coupling bulk convection-diffusion equation needs to be solved in one of bulk re- gions. In early computational models of surfactant problems, the boundary integral for- mulations in Stokes flow were extensively employed [12, 28, 30, 31, 40, 46], refer to [35] for a review. In order to take the inertia effect into account, various numerical methods were proposed to incorporate with Navier-Stokes equations. These methods are based upon volume-of-fluid (VOF) [4, 11, 21], front-tracking method [9, 32, 33], immersed boundary method [7, 25–27], level-set [47, 48], and arbitrary Lagrangian-Eulerian (ALE) finite ele- ment method [3, 15–17, 19, 36]. Apart from the methods listed above, other numerical studies include diffuse-interface method [43], segmented projection method [22, 23], lat- tice Boltzmann method [13], moving particle semi-implicit method [14], and smoothed particle hydrodynamics (SPH) method [1] etc. Some hybrid approaches combining two different methods were proposed as well [6, 29, 49, 51]. In this paper, we extend the third author’s previous 2D [25, 26] and axisymmetric [27] works and propose a three- dimensional immersed boundary (IB) method for simulating the evolution of a deform- ing interface with insoluble surfactant.
As known, solving the surfactant convection-diffusion equation on an evolving inter- face has already been a numerical challenge, especially in three-dimensional space. From the authors point of view, there are two major numerical issues to be considered carefully.
(1) How to maintain the interfacial mesh quality to avoid Lagrangian markers clustering and surface mesh distortion without re-meshing? (2) How to guarantee that the total sur- factant mass is conserved along the evolving interface without re-scaling? In this paper, we aim to propose a numerical scheme that can handle both issues systematically for an evolving free interface. Our numerical approach is outlined as follows. We first assume a doubly periodic interface separating two fluids where the surfactant exists only along the evolving interface. An equi-arclength parametrization is introduced in order to track the moving interface and maintain good Lagrangian meshes, so stable computations can be performed without re-meshing. This surface mesh-control technique is done by adding two artificial tangential velocity components into the Lagrangian marker velocity so that the Lagrangian markers can be equi-arclength distributed during the time evolution. As a result, the surfactant equation on the interface must be modified based on the new parametrization. A conservative scheme for solving the modified surfactant equation
has been developed and proved to satisfy the total surfactant mass exactly in discrete level. We then apply the developed scheme to simulate the problems of interfacial flows with insoluble surfactant.
The rest of this paper is organized as follows. In the next section, we present the equations of motion under the immersed boundary framework. We also describe the re- lated contents of surface differential geometry needed in our mathematical formulation.
In Section 3, an equi-arclength parametrization for a deforming interface is introduced and the modified surfactant concentration equation is derived accordingly. In Section 4, we first develop a conservative scheme for solving the modified surfactant equation and then outline the time-stepping scheme for solving the whole fluid equations with insol- uble surfactant. A series of numerical tests to validate our present scheme is given in Section 5 which is followed by some conclusions and future work in Section 6.
2 Equations of motion
In this section, we consider the flow of two immiscible viscous incompressible Navier- Stokes fluids in a three-dimensional domain Ω. The sharp interface Σ separating those two fluids is assumed to be doubly periodic in the horizontal x and y directions while it is free in the z direction. Besides the external applying flow, the motion of the inter- face and fluids is mainly driven by the surface tension along the free interface where an insoluble surfactant is distributed and varied on the interface that changes the surface tension in a timely fashion. Using the immersed boundary (or front-tracking) formula- tion, the fluid variables are represented in Eulerian manner while the interface quantities are represented in Lagrangian manner so the dimensionless governing equations for the above two-phase flow with insoluble surfactant can be written as one whole fluid system as follows [9, 25].
∂u
∂t+(u·∇)u= −∇p+ 1
Re∇·µ∇u+(∇u)T+ f
ReCa in Ω, (2.1)
∇·u=0 in Ω, (2.2)
f(x,t) =
Z
ΣF(α,β,t)δ(x−X(α,β,t))dA in Ω, (2.3)
∂X
∂t(α,β,t) =U(α,β,t) =
Z
Ωu(x,t)δ(x−X(α,β,t))dx on Σ, (2.4) F(α,β,t) = ∇sσ−2σHn, σ(α,β,t) =σ0(1−ηΓ(α,β,t)) on Σ, (2.5)
∂Γ
∂t(α,β,t)+(∇s·U)Γ(α,β,t) = 1 Pes
∆sΓ(α,β,t) on Σ, (2.6) Eqs. (2.1) and (2.2) are the Navier-Stokes equations in which u is the velocity and p is the pressure. Here we assume both fluids have equal density but different viscosity de- noting the upper fluid viscosity by µ+and the lower by µ−. So the Reynolds number is defined by Re=ρV L/µ+ where ρ,V and L are the corresponding scales for the density,
velocity and length. We assume µ−≥µ+throughout this paper although this is not a re- striction to our present method. Another dimensionless number in fluid equations is the capillary number Ca=µ+V/σ0where σ0represents the constant surface tension without surfactant. The interface Σ(t)between two fluid layers is a two-dimensional surface in Ω, represented by Σ(t)={X(α,β,t)|0≤α≤2π,0≤β≤2π}, where(α,β)are the Lagrangian coordinates. Throughout the paper, the effect of gravity is neglected so the only force in Eulerian fluid domain arises from the spreading of interfacial force F via the Dirac delta function δ(x) =δ(x)δ(y)δ(z)as shown in Eq. (2.3). Similarly, the interfacial veloc- ity U is also interpolated from the local Eulerian fluid velocity via the delta function in Eq. (2.4). The interfacial force F in Eq. (2.5) consists of tangential and normal components;
namely, the surface gradient of tension∇sσ (Marangoni force) and the mean curvature force 2σHn (capillary force). Here, H is the mean curvature of the free interface and the normal vector n is chosen pointing into the upper fluid direction. The surface tension related to the surfactant concentration is determined by the Langmuir equation of state in Eq. (2.5) where η (0≤η<1) is the dimensionless elasticity number that quantifies the sensitivity of interfacial tension to the changes of surfactant concentration as in [34]. The surfactant is insoluble and its concentration Γ satisfies the convection-diffusion equa- tion [41] as in (2.6), where Pes is the surface Peclet number. Here Γ(α,β,t) is defined in Lagrangian coordinates, so ∂Γ∂t is the material derivative. To be complete, the above governing equations (2.1)-(2.6) must be accompanied with suitable initial and boundary conditions.
In the following, we shall introduce some preliminary surface differential geometry [10] to express the mathematical forms for the surface gradient ∇s, surface divergence
∇s·, and surface Laplace ∆s operators used in above formulation. For a surface patch X(α,β)at some fixed time, the coefficients of the first fundamental form are defined by
E=Xα·Xα, F=Xα·Xβ, G=Xβ·Xβ,
where the subscripts α and β of a function denote its partial derivatives of the function with respect to α and β, respectively. So the local area stretching factor can be written as
|Xα×Xβ| =√
EG−F2. We define the unit normal vector as n=Xα×Xβ/|Xα×Xβ|, where
|·|stands for the usual Euclidean norm. The coefficients of the second fundamental form are defined by
L=Xαα·n= −Xα·nα, M=Xαβ·n= −Xα·nβ= −Xβ·nα, N=Xββ·n= −Xβ·nβ. So the mean curvature of the surface can be written as H=−GL2(−EGEN−+F22FM) .
Using those fundamental forms, the surface gradient ∇sσ of a scalar function σ can be represented by
∇sσ=Gσα−Fσβ
EG−F2 Xα+Eσβ−Fσα
EG−F2 Xβ, (2.7)
while the surface divergence∇s·Uof a vector field U can be represented by
∇s·U=GUα−FUβ
EG−F2 ·Xα+EUβ−FUα
EG−F2 ·Xβ. (2.8)
For a surface vector field Us=PXα+QXβ, its surface divergence can also be alterna- tively written as
∇s·Us= 1
|Xα×Xβ|
∂
∂α |Xα×Xβ|P+ ∂
∂β |Xα×Xβ|Q
. (2.9)
This formula can be derived by using direct substitution of Usinto Eq. (2.8). After com- bining two operators in Eqs. (2.7) and (2.9), the surface Laplacian (or Laplace-Beltrami operator) of Γ has the form
∆sΓ= ∇s·∇sΓ= 1
|Xα×Xβ|
"GΓα−FΓβ
|Xα×Xβ|
α
+
EΓβ−FΓα
|Xα×Xβ|
β
#
. (2.10)
3 An equi-arclength parametrization for a deforming free surface and the modified surfactant equation
As mentioned in the Introduction, three-dimensional interfacial flow simulations pose numerical challenges to track the free surface using the underlying surface meshes. More precisely speaking, it is difficult to maintain a surface parametrization using global co- ordinates that keeps the underlying mesh in good quality (without too much mesh dis- tortion) during the time evolution. In this paper, we propose an approach to maintain such global parametrization by first choosing an equi-arclength parametrization along two principal tangential directions and then try to maintain the property by adding ap- propriate artificial tangential velocities. The validity of such mechanism is theoretically ensured because the interfacial shape remains unchanged by adding any tangential ve- locity, while artificial normal velocity will change the shape accordingly which results in undesired dynamics. This is an extension to the earlier work [26] on a curve interface in two-dimensional interfacial flow by the third author and his coworkers. By exploiting this equi-arclength technique, one can see the Lagrangian markers to track the interface is indeed uniformly distributed [26]. One should notice that, the idea of adding two tan- gential velocities to the surface velocity to preserve some parametrization properties is not new; for instance, a generalized isothermal parametrization was introduced in 3D boundary integral computations for incompressible Darcy’s flow with surface tension in [2]. More relevant works to study different flows using the same idea can be found in the references therein.
To proceed, we modify the surface evolutional equation (2.4) by adding two tangen- tial velocities as
∂X
∂t(α,β,t) =U(α,β,t)+V1(α,β,t)τ1+V2(α,β,t)τ2, (3.1) where τ1=Xα/|Xα|and τ2=Xβ/|Xβ|are the unit tangent vectors satisfying τ1×τ26=0.
The goal is to derive the two evolutional equations for V1and V2so that the equi-arclength
parametrization for the surface can be always preserved if the initial surface is chosen to have such property. Like the free interface X, here, we assume that the velocities U,V1,V2
are all 2π doubly periodic. For a surface X(α,β,t), an equi-arclength parametrization in both principal tangential directions satisfies
∂
∂α|Xα| =0 and ∂β∂ Xβ =0, for all(α,β) ∈ [0,2π]×[0,2π]. Thus, we have
|Xα| = 1 2π
Z 2π
0 |Xα′|dα′ and |Xβ| = 1 2π
Z 2π
0 |Xβ′|dβ′. (3.2) Taking the time derivative in Eq. (3.2), it yields
|Xα|t= 1 2π
Z 2π
0 |Xα′|tdα′ and |Xβ|t= 1 2π
Z 2π
0 |Xβ′|tdβ′, (3.3) where|Xα|tcan be expressed by
|Xα|t=Xαt·Xα
|Xα| =Xαt·τ1= (Xt)α·τ1
=∂U
∂α·τ1+∂V1
∂α +∂V2
∂α
τ2·τ1+V2∂τ2
∂α ·τ1.
Substituting the above equation into Eq. (3.3) and integrating with respect to α, we obtain V1(α,β,t) = α
2π Z 2π
0 Q(α′,β,t)dα′−
Z α
0 Q(α′,β,t)dα′, (3.4) where
Q(α,β,t) =∂U
∂α·τ1+∂V2
∂α
τ2·τ1+V2∂τ2
∂α ·τ1. Similarly, one can derive
V2(α,β,t) = β 2π
Z 2π
0 R(α,β′,t)dβ′−
Z β
0 R(α,β′,t)dβ′, (3.5) where
R(α,β,t) =∂U
∂β·τ2+∂V1
∂β
τ1·τ2+V1
∂τ1
∂β ·τ2
.
In the above derivations, we impose the boundary conditions for V1(0,β,t)=V2(α,0,t)=0, implying that at these boundary points, no additional tangential velocities are needed in their corresponding directions.
Since the surfactant transports and diffuses along the interface, by taking the new surface parametrization Eq. (3.1) into account, the original surfactant equation (2.6) must be modified to
∂Γ
∂t−(V1τ1+V2τ2)·∇sΓ+(∇s·U)Γ= 1 Pes
∆sΓ.
Multiplying the surface stretching factor|Xα×Xβ|on both sides of the above equation and using the following identity derived in [38],
∂|Xα×Xβ|
∂t = |Xα×Xβ|∇s·∂X∂t = |Xα×Xβ|∇s·U+V1τ1+V2τ2
,
we obtain
∂Γ
∂t|Xα×Xβ|−|Xα×Xβ|(V1τ1+V2τ2)·∇sΓ +∂|Xα×Xβ|
∂t Γ−|Xα×Xβ|h∇s·V1τ1+V2τ2
iΓ=|Xα×Xβ| Pes
∆sΓ.
By putting the terms together, one can simplify the above equation as
∂ Γ|Xα×Xβ|
∂t −|Xα×Xβ|∇s·
V1
Xα
|Xα|+V2
Xβ
|Xβ|
Γ
=|Xα×Xβ| Pes
∆sΓ.
Using the formulas for the surface divergence in Eq. (2.9) and the surface Laplacian in Eq. (2.10), we have
∂ Γ|Xα×Xβ|
∂t −
"
|Xα×Xβ|V1Γ
√E
α
+
|Xα×Xβ|V2Γ
√G
β
#
= 1 Pes
"GΓα−FΓβ
|Xα×Xβ|
α
+
EΓβ−FΓα
|Xα×Xβ|
β
#
. (3.6)
One can immediately see that the above modified surfactant equation is written in a form of conservation law which has the advantage of designing a conservative scheme.
Indeed, in next section, Eq. (3.6) will be discretized by a conservative finite difference scheme and numerical evidences in Section 5 demonstrate that the total surfactant mass along the interface is conserved within the machine accuracy.
4 Numerical method
In this section, we present the numerical method and some related implementation de- tails for solving the governing equations introduced in the previous section. To solve the Navier-Stokes equations (2.1)-(2.2) in a computational domain Ω⊂R3, we layout a
Figure 1: The left shows the locations of fluid velocities and pressure on a staggered grid in 3D. In the right, the Lagrangian markers X (filled circles) form a rectangular mesh and the surfactant concentration Γ (open circles) is located at the cell center of the interface.
uniform Cartesian grid with meshwidth h=∆x=∆y=∆z, and allocate the fluid velocity u= (u,v,w)and the pressure p in a staggered grid manner [20] as shown in the left of Fig. 1. The boundary conditions for the x and y directions are 2π periodic while the z direction is Dirichlet. As in [2], we use a Fourier representation to discretize the interface as X(α,β,t) = (α,β,0)T+Y(α,β,t), where the first term represents the flat interface and the second term represents the periodic deviation from the plane written in truncated double Fourier series as
Y(α,β,t) =
N1/2−1
∑
k1=−N1/2
N2/2−1
∑
k2=−N2/2
Yb(k1,k2,t)ei(αk1+βk2). (4.1)
Therefore, the interface is represented by a set of Lagrangian markers Xℓm=X(αℓ,βm), where(αℓ,βm) = (ℓ∆α,m∆β)are the uniform grid points with the corresponding mesh- widths ∆α=2π/N1, ∆β=2π/N2in the(α,β)plane as shown in the right of Fig. 1. The as- sociated geometric quantities listed in previous sections such as the coefficients of the first and second fundamental forms, interface tangents and normals, and the mean curvature are all computed at the Lagrangian markers Xℓm using the Fourier spectral differentia- tions [44]. These spectral derivatives with respect to α and β can be computed efficiently using Fast Fourier Transform (FFT) in a spectral accuracy. For instance, a scalar doubly periodic function f(α,β)can be written as Eq. (4.1) with the truncated Fourier coefficients bf(k1,k2)so the derivative with respect to α only involves the multiplication of the Fourier coefficient bf by ik1. The inversion between the function values and its discrete Fourier co- efficients can be performed efficiently using FFT. The derivative with respect to β can be computed in the same manner. To remove the aliasing error in computations, we adopt the well-known de-aliasing 2/3-rule filter [5].
4.1 A conservative scheme for solving the modified surfactant equation As mentioned earlier, the modified surfactant equation (3.6) is in conservation form so it is quite natural to develop a conservative scheme so that the total surfactant mass along the interface is conserved. Unlike the Lagrangian markers and other geometrical quanti- ties defined at the grid points(αℓ,βm), we define the surfactant concentration at the cell center and denote it by Γℓ+21,m+12 as shown in the right of Fig. 1. In the following, we discretize Eq. (3.6) using forward Euler method in time and centered difference in space as
Γn+1
ℓ+12,m+21(Sαβ)nℓ++11
2,m+12−Γnℓ+21,m+12(Sαβ)nℓ+1 2,m+12
∆t
−∆α1
(Sαβ)nℓ+
1,m+21(V1)nℓ++1
1,m+12Γn
ℓ+1,m+21
qEℓn+1,m+1 2
− (Sαβ)nℓ
,m+12(V1)nℓ+1
,m+12Γn
ℓ,m+21
qEnℓ,m+1 2
−∆β1
(Sαβ)nℓ+1
2,m+1(V2)nℓ++11 2,m+1Γnℓ
+21,m+1
qGnℓ+1 2,m+1
−
(Sαβ)nℓ+1
2,m(V2)nℓ++11 2,mΓnℓ
+12,m
qGnℓ+1 2,m
= 1 Pes∆α
1
Sαβn ℓ+1,m+12
Gℓn+1,m+1 2
Γnℓ
+23,m+12−Γnℓ+12,m+12
∆α −Fℓn+1,m+21
Γnℓ+
1,m+1−Γnℓ+1,m
∆β
!
− 1
Sαβ
n ℓ,m+12
Gℓn,m+1 2
Γn
ℓ+12,m+12−Γnℓ−12,m+12
∆α −Fℓn,m+21
Γnℓ
,m+1−Γnℓm
∆β
!
+ 1 Pes∆β
1
Sαβ
n ℓ+21,m+1
Eℓn+1 2,m+1
Γn
ℓ+12,m+32−Γnℓ+21,m+12
∆β −Fℓn+12,m+1
Γnℓ+
1,m+1−Γnℓ,m+1
∆α
!
− 1
Sαβ
n ℓ+12,m
Eℓn+1 2,m
Γn
ℓ+21,m+12−Γnℓ+21,m−12
∆β −Fℓn+12,m
Γnℓ+1,m−Γnℓm
∆α
!
, (4.2)
where Sαβdenotes|Xα×Xβ|=√
EG−F2, the local stretching factor on the interface. In the above scheme, all the variables except Γ are defined at the Lagrangian marker X shown in Fig. 1. So, when a variable is necessarily shifted by a half-mesh width for relevant computations, the linear interpolation is simply employed. For instance, we define
(Sαβ)nℓ+1
2,m=(Sαβ)nℓm+(Sαβ)nℓ+1,m 2
and
(Sαβ)nℓ+1
2,m+21=(Sαβ)nℓm+(Sαβ)nℓ+1,m+(Sαβ)nℓ,m+1+(Sαβ)nℓ+1,m+1
4 .
In the similar fashion, the surfactant concentration Γnℓmcan be linearly interpolated from the values of neighboring four points.
By taking the double summation with respect to ℓ and m on both sides of Eq. (4.2) and using the doubly periodic boundary condition, one can easily prove that the total surfactant mass in the scheme (4.2) is conserved. That is,
N1−1
∑
ℓ=0 N2−1
∑
m=0
Γn+1
ℓ+21,m+12(Sαβ)nℓ+1
+12,m+12∆α∆β=
N1−1
∑
ℓ=0 N2−1
∑
m=0
Γnℓ
+12,m+21(Sαβ)nℓ+1
2,m+12∆α∆β. (4.3) Notice that, this conservation property in Eq. (4.2) is independent of the numerical scheme for computing the quantities associated with the interfacial geometry and the artificial tangential velocity such as E,F,G,Sαβ and V1,V2, respectively. In the next section, the relative error of total surfactant mass versus time will be presented to verify the mass conservation for various numerical experiments.
4.2 Time-stepping scheme
We here present how to march the Lagrangian markers Xn=X(n∆t)from time level n to obtain Xn+1=X(n∆t+∆t)at time level n+1 with ∆t the time step size. In the present IB method, the surfactant concentration Γn, the fluid velocity un, the pressure pn, and the Lagrangian markers Xnare all given in advance, and from these variables we aim to update Γn+1, un+1, pn+1, and Xn+1. The step-by-step numerical procedure can be done as follows.
1. At the Lagrangian markers Xnℓm, we first compute the tension by σℓm=σ(Xnℓm) = σ0(1−ηΓnℓm), and then compute the interfacial tension force
F(Xnℓm)dA(Xnℓm) = (∇sσℓm−2σℓmHℓmnℓm)dA(Xnℓm),
where the (tangential) Marangoni force density ∇sσℓm and the (normal) capillary force density 2σℓmHℓmnℓm are obtained using the formulas written in earlier sec- tions. The surface area element is computed by dA(Xnℓm) = (Sαβ)nℓm∆α∆β.
2. Distribute the tension force acting on Lagrangian markers into the Eulerian grid by using the smoothed Dirac delta function δhas
fn(x) =
N1−1
∑
ℓ=0 N2−1
∑
m=0
F(Xnℓm)δh(x−Xnℓm)dA(Xnℓm), (4.4)
where x= (x,y,z)is the Eulerian grid point. For δh(x) =h13φ xh φ yh
φ hz
, we em- ploy the 4-point supported function φ developed in [50] to suppress spurious force oscillations in IB method.
3. Solve the Navier-Stokes equations by the second-order incremental pressure-correction projection method in [18] as follows.
3u⋆−4un+un−1
2∆t +2(un·∇h)un−un−1·∇h
un−1
=−∇hpn+ 1 Re
h
λ∆u⋆−λ∆un+∇h·µ∇hun+(∇hun)Ti+ f
n
ReCa,
∆hp⋆= 3
2∆t∇h·u⋆, ∂p
⋆
∂n =0 on ∂ΩD, u⋆=un+1 on ∂ΩD, un+1=u⋆−2∆t3 ∇hp⋆, ∇hpn+1= ∇hp⋆+∇hpn−2λ∆t3Re ∆h(∇hp⋆), where the discrete operators ∇h and ∇h· approximate the gradient and diver- gence operators, respectively, using the second-order finite difference in stag- gered grid. For the nonlinear terms, the skew-symmetric form is employed as (u·∇h)u=12(u·∇h)u+12∇h(uu). The domain boundary is denoted by ∂ΩD. Here, λ=µ−/µ+ is the viscosity contrast while the dimensionless viscosity can be de- fined by µ=1+(λ−1)I(x)after the normalization using µ+. The indicator function I(x) =1 in Ω−and zero elsewhere can be obtained by solving Poisson equation as described in [45].
4. Once un+1 is determined in the Eulerian fluid points, we interpolate the fluid in- terfacial velocity Unℓm+1=∑xun+1(x)δh(x−Xnℓm)h3 in Eq. (3.1). In order to maintain an equi-arclength distribution of updated Lagrangian markers, it requires to solve the artificial velocities(V1)nℓm+1and(V2)nℓm+1in Eqs. (3.4) and (3.5), respectively. This can be simply done using the fixed point iteration by taking V2=0 initially. The error of tolerance is set to be h/100 so that only 2-3 iterations are needed in each time step. (Note that, the magnitude of the error of tolerance affects only the accu- racy of V1and V2here. Since the immersed boundary method is first-order accurate in general, as long as we choose the error of tolerance as the same order of the method accuracy, it has no significant effect on the surfactant and fluid variables as a whole.) For the numerical integration in (3.4) and (3.5), the trapezoidal rule is used. Therefore, the new Lagrangian markers can be updated by
Xnℓm+1=Xnℓm+∆t
Unℓm+1+(V1)nℓm+1
τ1
n
ℓm+(V2)nℓm+1 τ2n ℓm
.
5. After computing (Sαβ)nℓ++11
2,m+12, we update surfactant concentration distribution Γn+1
ℓ+12,m+12 using Eq. (4.2), and then apply the linear interpolation to obtain Γnℓm+1. We summarize the implementing difficulties and differences for the numerical algo- rithm between previous 2D [25,26] and the present 3D version as follows. Besides the computational complexity of the fluid solver used, the significant differences fall into the
handling of the interface (an evolving curve in 2D and an evolving surface in 3D) and solving the surfactant equation on the corresponding interfaces. As discussed in Sec- tion 3, to keep the interface mesh as equi-arclength distributed, two artificial tangential velocity components must be determined in 3D rather than one component in 2D. The modified surfactant equation in 3D involves solving convection-diffusion equation on an evolving surface which is much more difficult to solve than in 2D case as we have emphasized in the third paragraph of the Introduction. Here, we develop an explicit scheme for the modified surfactant equation being able to save the computational cost and to preserve the total surfactant mass exactly.
5 Numerical results
A series of numerical experiments is carried out in this section. Firstly, we examine the effect and necessity of the equi-arclength Lagrangian mesh control developed in Section 3. We then check the rates of convergence for the velocity field, interfacial configuration, and the surfactant concentration. Self-healing dynamics in quiescent flow with non-unity viscosity contrast is investigated in detail and show qualitatively similar results with the ones obtained in literature. Lastly, we study the two-layer fluids under Couette flow with different physical parameters in detail. Throughout the paper, the velocity boundary conditions in both x and y directions are set to be 2π periodic while in the z direction is Dirichlet. Since the interface is parallel to x−y plane initially, the interface is also 2π doubly periodic as described in the beginning of Section 4.
5.1 Effect of Lagrangian mesh control
As mentioned earlier, for a deforming surface with interfacial tension in fluid flows, the underlying Lagrangian mesh will inevitably distort and consequently induce numerical instability of computations. In this subsection, we check the performance and necessity of the enforcement of the equi-arclength parametrization mesh control technique presented in Section 3. The problem setup is the following. We assume the interface is flat initially in the plane z=π with the initial surfactant concentration
Γ(α,β,0) =1−tanh(20(2−r))
2 , (5.1)
where r is the distance of Lagrangian marker from the center(π,π,π). So the centered region is almost surfactant free Γ≈0, while the outside region is Γ≈1, see Fig. 2(a) in detail. The flow is quiescent so the only driving force to the fluid comes from the surface tension. We also choose the same viscosity for the upper and lower fluids; that is, the viscosity contrast λ=µ−/µ+=1. In this test, we use Re=1, Ca=0.01, Pes=1, and η=0.8. We fix σ0=1 for all tests carried out in this paper. The grid size is 1283 in the computational domain [0,2π]3 so the fluid mesh width is h=2π/128. The number of
Figure 2: The underlying Lagrangian mesh and the corresponding surfactant concentration at different times.
Left panel (without mesh control); right panel (with mesh control).
Lagrangian markers on the interface is 2562 so the Lagrangian mesh is ∆α=∆β=h/2.
The time step size is chosen as ∆t=h/512.
Since the interface is flat and 2π doubly periodic, the only mechanism driving the fluid dynamics is the inward Marangoni force due to surface tension gradient (thus sur- factant concentration gradient). Fig. 2 shows the underlying Lagrangian mesh and the corresponding surfactant concentration at different times in which the left panel is the one without imposing equi-arclength mesh control while the right panel is the one with mesh control. One can immediately see from Fig. 2(c) and (e) that the grid lines distort significantly and eventually the computation breaks down at t=0.063, see the erroneous value scale. However, for the results with mesh control in Fig. 2(d) and (f), the Lagrangian mesh remains uniformly and the surfactant undergoes uniform spreading on the flat in- terface.
0 2 4 6 0.01
0.02 0.03 0.04
x
arc length
(a) t=0.031
0 2 4 6
0.01 0.02 0.03 0.04
x
arc length
(c) t=0.063
0 2 4 6
0.01 0.02 0.03 0.04
x (b) t=0.063
0 2 4 6
0.01 0.02 0.03 0.04
x (d) t=2.78
Figure 3: The plots of arclength|Xα|∆α. (a)-(d) correspond to Fig. 2(c)-(f), respectively. The straight red line is the initial uniform arclength.
Fig. 3 shows the corresponding plots of the arclength|Xα|∆α for each α-directed curve at some selected times, in which (a)-(d) correspond to Fig. 2(c)-(f). We can clearly see that the arclength in the vicinity of the circular border of surfactant-free region varies consid- erably and its deviation becomes significant as time evolves as shown in Fig. 3(a) and (c), while the arclength is a nearly constant (overlapped with solid red lines) in Fig. 3(b) and (d).
Fig. 4(a) depicts the fluid velocity field (u,v) on the plane z=π corresponding to Fig. 2(d) which indeed shows the Marangoni force driving the fluid motion inwardly.
Fig. 4(b) depicts the velocity(u,w)on the plane y=π in which one can see the circulatory motion appears. The corresponding streamline of Fig. 4(b) is shown in (c) to demon- strate the interface keeps flat as initially so no penetration flow occurs (normal velocity approximates to zero numerically). Fig. 4(c) also shows that the fluid leakage across the interface is negligible as long as the number of Lagrangian markers and the time step size are chosen appropriately.
5.2 Convergence study
In this subsection, we perform a convergence study of the present numerical method with respect to fluid velocity, Lagrangian marker, and surfactant concentration. The di- mensionless parameters are also chosen by Re=1, Ca=0.01, Pes=1, and η=0.8. The viscosity contrast is set by λ=2, which will consequently lead to interface deformation.
We set the time step size by ∆t=h/256. For the computational domain[0,2π]3, we use
0 2 4 6 0
2 4 6
x
y
(a)
0 2 4 6
0 2 4 6
x
z
(b)
x
z
(c)
0 2 4 6
0 2 4 6
−0.5 0 0.5
Figure 4: At t=0.063 with mesh control corresponding to Fig. 2(d). (a) the velocity field on the plane z=π; (b) the velocity field on the plane y=π; (c) the circulatory streamline corresponding to (b) with the flat interface at z=π. The color represents the stream function magnitude.
four different Cartesian grid sizes N=32,64,128,256 to estimate the convergence rate of the variables. In the interface discretization, as we double the grid size N, we increase the number of Lagrangian markers by a factor of four accordingly. As previous exam- ple, the initial interface is a flat plane and the initial distribution of surfactant is given by Γ(α,β,0) = (1−tanh(20(2−r)))/2, where r is defined as in previous example. The resul- tant motion up to t=2.4 are similar to the inward spreading of surfactant illustrated in Figs. 6 and 7(a).
Table 1 exhibits the rates of convergence of the velocity components u= (u,v,w), the interfacial marker X, and the surfactant concentration Γ at two different times t=1.5 and t=2.4, respectively. In this test, the analytic solution is not known, so we compute the convergence rate by estimating the two successive errors as
Rate=log2(kuN−u2Nk∞/ku2N−u4Nk∞).
The error kuN−u2Nk∞ is the maximum absolute difference of the numerical solutions for the grid N and the grid 2N. The rates for other variables are defined in a similar manner. One can see from Table 1 that the average rate is around one for all variables indicating that our numerical scheme is roughly first-order accurate. This supports that the IB method is generally first-order accurate.
