surfactant by immersed boundary method
Ming-Chih Lai1,∗, Yu-Hau Tseng1, and Huaxiong Huang2
1Center of Mathematical Modeling and Scientific Computing & Department of Ap- plied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, TAIWAN.
2Department of Mathematics and Statistics, York University, Toronto, Ontario, Canada M3J 1P3.
Abstract. In this paper, we present an immersed boundary method for simulating moving contact lines with surfactant. The governing equations are the incompressible Navier-Stokes equations with the usual mixture of Eulerian fluid variables and La- grangian interfacial markers. The immersed boundary force has two components: one from the nonhomogeneous surface tension determined by the distribution of surfac- tant along the fluid interface, and the other from unbalanced Young’s force at the mov- ing contact lines. An artificial tangential velocity has been added to the Lagrangian markers to ensure that the markers are uniformly distributed at all times. The corre- sponding modified surfactant equation is solved in a way such that the total surfac- tant mass is conserved. Numerical experiments including convergence analysis are carefully conducted. The effect of the surfactant on the motion of hydrophilic and hydrophobic drops are investigated in detail.
AMS subject classifications: 65M06, 76D45
Key words: Immersed boundary method; Interfacial flow; Navier-Stokes equations; Surfactant;
Moving contact line; Hydrophilic drop; Hydrophobic drop; Wetting.
1 Introduction
Surfactant are surface active agents that adhere to the fluid interface and reduce the inter- face tension by lowering the surface energy. Surfactant play an important role in many applications in the food, cosmetics and oil industries. For instance, extraction of ore relies on the subtle effects introduced by the presence of surfactant [3]. In a liquid-liquid sys- tem, surfactant allow small droplets to be formed and used as an emulsion. Surfactant
∗Corresponding author. Email addresses: mclai@math.nctu.edu.tw (M.-C. Lai), yhtseng@math.nctu.edu.tw (Y.-H. Tseng), hhuang@yorku.ca (H. Huang)
http://www.global-sci.com/ Global Science Preprint
also play an important role in water purification and other applications where micro- sized bubbles are generated by lowering the surface tension at the liquid-gas interface.
In microsystems with the presence of interfaces, it is extremely important to consider the effect of surfactant since in such cases the capillary effect dominates the inertia of the fluids [33].
When a fluid-fluid interface is in contact with a solid substrate, surfactant can change the wetting properties by altering the value of the contact angle. This simple fact has found may interesting applications in our daily life and industrial processes. For ex- ample, we add detergents (surfactant) in washing machine to clean our clothes more effectively. The detergent helps to remove drops of grease from clothes by increasing the contact angle (measured from inside the drop). An idealization of this problem can be found in a photo shown in [33] or a figure in [3] where it is demonstrated that a drop on clothes can become less “wetting” (from the drop point of view) by increasing sur- factant concentration. By adding surfactant, a drop which sticks to clothes becomes less
“sticky” and the water currents can wash away the drop readily. The physical situation corresponding to this idealized system includes a solid-drop (grease)-water system and the surfactant. Mathematically, as discussed more detail in the main text of this paper, this is a moving contact line problem since the solid-drop (grease)-water triple intersec- tion forms a contact line. The main issue we will try to address in this paper is how surfactant, by changing the contact angle, affects the movement of the contact line.
In the past few years, several numerical methods have been developed for interfacial flows with surfactant, such as the level set method [35], volume-of-fluid method [14], front-tracking method [16, 20], and hybrid methods [1]. Recently, we have proposed an immersed boundary method for simulating the motion of two-dimensional fluid inter- faces with insoluble surfactant [15]. The mathematical models include the incompressible Navier-Stokes equations and a surfactant concentration equation along the moving inter- face. A numerical scheme that conserves the total surfactant mass is introduced and drop deformation under shear flow is studied in detail. In the present work, the Navier-Stokes solver is similar to our previous one in [15] with one modification. Instead of moving with the fluid as in [15], the present Lagrangian markers are uniformly distributed along the interface during the time evolution. This is achieved by using an artificial tangen- tial velocity. Despite this modification, our new scheme for the surfactant equation still conserves the total surfactant mass exactly.
Moving contact line problems arise when one fluid is displaced by another immiscible fluid on a solid surface. They appear in the process of wetting, coating and many biolog- ical applications [3, 21, 29]. It is well-known that the no-slip boundary condition in the vicinity of the moving contact line leads to a non-integrable stress singularity [5, 12, 29].
In other words, if the no-slip boundary condition is imposed on the solid substrate, an infinite force is required to move the contact line. For the past three decades, many at- tempts have been made to resolve this physically unrealistic force singularity. It has been confirmed by the molecular dynamics (MD) simulation that fluid does slip near the con- tact line [24, 25]. In [24], a generalized Navier boundary condition has also been derived
using a diffuse interface approach. For an in-depth discussion, we refer the readers to a recent review [26] and the references therein. In this paper, we avoid force singularity by applying the Navier slip boundary condition in the wall of the moving contact line.
In additional to imposing the Navier slip boundary condition on the solid wall to re- move force singularity, the contact line dynamics must also be prescribed. In the presence of contact angle hysteresis, one can prescribe the contact angle depending on the sign of contact line speed [7, 19, 32]. A more complicated dynamic contact model that takes both low and high capillary number into account has been proposed recently in [34]. In the above models, however, the advancing and receding contact angles must be given. Re- cently, Ren and E derived an effective boundary condition at the contact line from the force balance argument [27]. As mentioned by the authors, the main driving force for the slip is the unbalanced Young’s force which results from the deviation of the contact angle from its static value. Therefore, in this work, we simply apply the unbalanced Young’s force to the contact line directly.
Numerical simulations of two-phase flows with moving contact lines have also been developed in the literature. They include the front tracking method [11], level set method [32], volume-of-fluid method [28, 31, 34], phase field method [13, 26], and sharp interface Cartesian grid method [19]. However, to the best of our knowledge, there is no previous numerical study which takes interface contaminant (surfactant) into account. As a first attempt, we present an immersed boundary method in this paper for simulating moving contact lines with surfactant.
In section 2, the basic mathematical model is introduced. The governing equations are the Navier-Stokes formulation with the usual mixture of Eulerian flow and Lagrangian interfacial variables. The immersed boundary force arises from two parts, namely, the nonhomogeneous surface tension determined by the distribution of surfactant along the fluid interface, and the unbalanced Young’s force applied at the moving contact lines.
An explicit formula for the artificial tangential velocity, which equi-distributes the La- grangian markers, is derived. Numerical method is discussed in section 3, where it is shown that the scheme for the surfactant equation strictly conserves the mass. In sec- tion 4, careful numerical experiments, including convergence analysis are conducted.
The effect of the surfactant on the contact angles for hydrophilic and hydrophobic drops are investigated in detail. We finish the paper by a short conclusion in section 5.
2 Mathematical model
In this paper, we consider two immiscible fluids (fluids 1 and 2) that are placed on a solid substrate in which the interface of the two fluids intersects the substrate at the so called contact lines. (More precisely, it should be called the contact points in two-dimensional flows.) Here, we assume fluid 1 is a drop surrounded by fluid 2 and both fluids rest on the solid substrate; thus, the contact lines are at the bottom of the domain, cf. Figure 1.
Furthermore, we assume that an insoluble surfactant may diffuse along the fluid-fluid
interface and the surface tension varies according surfactant concentration. Throughout this paper, we assume that the two fluids have equal density and viscosity values. The gravity is neglected since the main focus is on the capillary effect, especially, the effect on contact line movements due to surfactant. In the absence of external forces and when the contact lines are static, the surface forces acting at the contact lines follow the well-known Young’s condition [3]
σs2=σs1+σ cosθe, (2.1)
where σs1, σs2, and σ are the corresponding surface tension between the solid and fluid 1, solid and fluid 2 and on the fluid-fluid interface, respectively. θeis the static equilibrium contact angle. When the contact lines are in motion, we need to model both the fluid flow and contact line dynamics.
σs2 σs1
σ
u = 0, v = 0
u = 0
v = 0 u = 0
v = 0
fluid2 fluid1
solid
u = β u
y , v = 0 θ
θ σ σs2 σs1
Figure 1: The diagram of contact lines and the problem set up.
The dynamics of the two-phase fluids over a solid wall y=c in a domain Ω= [a,b]×
[c,d] is governed by the incompressible Navier-Stokes equations using the immersed boundary formulation [15, 22],
∂u
∂t+(u·∇)u+∇p = 1
Re∇2u+ 1
ReCaf , in Ω (2.2)
∇·u = 0, in Ω (2.3)
f(x,t) =
Z
ΣF(α,t)δ(x−X(α,t))dα, (2.4) F(α,t) = ∂
∂α(στ), τ=∂X
∂α/
¯¯
¯¯∂X
∂α
¯¯
¯¯, (2.5)
∂X(α,t)
∂t =U(α,t) =
Z
Ωu(x,t)δ(x−X(α,t))dx. (2.6)
Here, the fluid interface Σ is represented by a parametric form X(α,t) = (X(α,t),Y(α,t)), 0≤α≤Lb, where α is the Lagrangian parameter. As shown in Fig. 1, we can see that X(0,t)and X(Lb,t)are the configurations of right and left contact lines, respectively. 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. Equations (2.4) and (2.6) represent the interaction between the immersed interface and the fluids. In particular, Eq. (2.4) describes the force ( f ) acting on the fluid due to the fluid interfacial force F and Eq. (2.6) states that the interface moves with the fluid velocity. The fluid interfacial force is given by Eq. (2.5), where σ is the fluid/fluid surface tension and τ is the unit tangent vector. The present formulation employs a mixture of Eulerian (x) and Lagrangian (X) variables which are linked by the two-dimensional Dirac delta function δ(x) =δ(x)δ(y).
As mentioned before, in order to avoid a non-integrable singularity on the moving contact line problems, the no-slip boundary condition at the solid wall y=c is replaced by the Navier slip boundary condition [4, 5]
u=β∂u
∂y, v=0, (2.7)
where β is the slip length which can be defined as the distance from fluid-solid interface to where the linear extrapolated tangential velocity vanishes. In addition, the contact line dynamics must be prescribed as well. Instead of using an effective velocity condition at the contact lines as in [27], we simply impose a point force Fcl at the contact line. The force at the moving contact line is only exerted along the tangential direction (that is x direction) of the solid substrate and mainly due to the unbalanced Young’s force [3] and can be written as
Fcl=σs2−σs1−σ cosθ, (2.8) where θ is the dynamic contact angle.
The surface tension σ(α,t)is related to the surfactant concentration Γ(α,t)through the following nonlinear approximation of Langmuir equation of state [17]
σ(Γ) =σc(1+ln(1−η Γ)) (2.9) where σcis the surface tension of a clean interface, and η satisfying 0≤η<1 is a dimen- sionless number that measures the sensitivity of surface tension to changes in surfactant concentration. The surfactant concentration equation [30] is
DΓ
Dt+(∇s·U)Γ= 1
Pes∇2sΓ, (2.10)
where∇sand∇2s are the surface gradient and surface Laplacian operators, respectively.
In this paper, those surface derivatives are given explicitly [15] as
∇s·U = ∂U
∂τ·τ= µ∂U
∂α·τ
¶
/|Xα|, (2.11)
∇sΓ = ∂Γ
∂ττ= µ∂Γ
∂ατ
¶
/|Xα|, (2.12)
∇2sΓ = ∇s·∇sΓ= ∂
∂α µ∂Γ
∂α/|Xα|
¶
/|Xα|. (2.13)
2.1 An equi-distributed technique for Lagrangian markers
In the context of immersed boundary simulations, the interface is tracked in a Lagrangian manner. Once the Lagrangian markers have been chosen initially, the movement of those markers are based on the interpolated local fluid velocity. Very often, as time evolves, those Lagrangian markers will be either clustered together or dispersed so the overall numerical stability or accuracy can be compromised. Therefore, certain grid redistribu- tion technique must be adopted to preserve better resolution. One approach is to add or delete marker points, when they are needed or unwanted, as shown in our previous im- mersed boundary simulation for a drop in a shear flow [15], where the marker points are gradually swept into the tips region. In this paper, we introduce another convenient way to dynamically control the Lagrangian markers by redistributing them equally along the interface as it evolves.
In order to remove the stiffness from the interfacial flows with surface tension, Hou, Lowengrub and Shelly [10] introduced an artificial tangential velocity into their formu- lation of boundary integral methods so that the particles can be uniformly distributed.
Following [10], we propose the following technique to maintain equi-distribution of the Lagrangian markers.
The idea is to introduce an artificial tangential velocity UA(α,t)in Eq. (2.6) as
∂X(α,t)
∂t =U(α,t)+UA(α,t)τ, (2.14)
so that the local stretching factor of the interface |Xα| is independent of α. It is clear that there should be no artificial tangential velocity at moving contact lines (α=0 and α=Lb) since the contact lines movement should not be affected by the marker’s artificial tangential velocity. Thus, we have UA(0,t) =UA(Lb,t) =0. Here, we need to impose
|Xα|α=0 so that|Xα| =L1
b
RLb
0 |Xα0|dα0. Taking the derivative with respect to t, we have
|Xα|t=L1
b
RLb
0 |Xα0|tdα0. Since
|Xα|t = Xαt·Xα
|Xα| = µ∂U
∂α +∂UA
∂α τ+UA∂τ
∂α
¶
·τ
= ∂U
∂α·τ+∂UA
∂α = (∇s·U)|Xα|+∂UA
∂α (2.15)
we have
∂U
∂α·τ+∂UA
∂α = 1 Lb
Z L
b
0
µ∂U
∂α0·τ0+∂UA
∂α0
¶
dα0. (2.16)
Integrating with respect to α and using the fact of UA(0,t) =UA(Lb,t) =0, we obtain UA(α,t) = α
Lb Z L
b
0
∂U
∂α0·τ0dα0−
Z α
0
∂U
∂α0·τ0dα0. (2.17) We note that an alternative formula for using artificial tangential velocity has been de- rived by Ceniceros [1] in a similar front-tracking manner. However, unlike the approach in [1], the present formula (2.17) does not need to compute the curvature of the interface.
2.2 Modified surfactant concentration equation
By taking the artificial velocity UAinto account, the material derivative now becomes DΓ
Dt=∂Γ
∂t−UAτ·∇sΓ. (2.18)
Therefore, the surfactant equation (2.10) becomes
∂Γ
∂t−UAτ·∇sΓ+(∇s·U)Γ= 1
Pes∇2sΓ. (2.19)
Multiplying the stretching factor|Xα|on both sides of the above equation, we obtain
∂Γ
∂t|Xα|−UA∂Γ
∂α+(∇s·U)|Xα|Γ= 1
Pes∇2sΓ|Xα|. (2.20) By writing the surface derivatives in terms of α explicitly and using the identity of Eq. (2.15), we obtain
∂Γ
∂t|Xα|+∂|Xα|
∂t Γ−∂(UAΓ)
∂α = 1 Pes
∂
∂α µ∂Γ
∂α/|Xα|
¶
. (2.21)
The finite difference discretization for the surfactant equation in the next section will be based on the above formulation. Notice that, when the artificial tangential velocity UAis set to be zero, the above equation is recovered to the original surfactant equation derived in [15]. We should also mention that a similar modified surfactant equation (by taking the artificial velocity into account) with a curvature term can be found in [1].
3 Numerical method
In this paper, the fluid flow variables are defined on a staggered marker-and-cell (MAC) mesh introduced by Harlow and Welsh [9]. The pressure is defined on the grid points
labelled as x= (xi,yj) = (a+(i−1/2)h,c+(j−1/2)h)for i,j=1,2...,N, the velocity com- ponents u and v are defined at(xi−1/2,yj) = (a+(i−1)h,c+(j−1/2)h)and(xi,yj−1/2) = (a+(i−1/2)h,c+(j−1)h), respectively, where the spacing h=∆x=∆y. For the immersed interface, we use a collection of discrete points αk=k∆α,k=0,1,... M such that the La- grangian markers are denoted by Xk=X(αk) = (Xk,Yk), where ∆α=Lb/M. Notice that, in present setting, X0 and XM correspond to the right and left contact lines, respec- tively. The surfactant concentration Γk+1/2, is defined at the “half-integer” points given by αk+1/2= (k+1/2)∆α. To maintain the accuracy, we use cubic splines to fit the inter- face such that τ, n, and κ are provided automatically. Without loss of generality, for any function φ(α)defined on the interface, we approximate the partial derivative ∂φ/∂α by Dα(φ(α)), where Dαis the differential operator based on the cubic splines.
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 Xnk=X(αk,n∆t), Γnk+1/2=Γ(αk+1/2,n∆t), un= u(x,n∆t), and pn−1/2=p(x,(n−1/2)∆t)are all given. The details of the numerical time integration are as follows.
1. Compute the surface tension and unit tangent on the interface as
σkn = σc(1+ln(1−ηΓnk)), (3.1) τnk = DαX
nk
|DαXnk|, (3.2)
both of which hold for αk=k∆α, and Γnk is simply approximated by the average of Γnk−1/2 and Γnk+1/2. Then we compute the interfacial force at the the fluid/fluid markers Xk, k=1,..., M−1 by
Fnk=Dα(σknτnk)∆α. (3.3) The unbalanced Young’s force at the contact lines (k=0 and k=M) can be computed by
Fn0 = (σs2−σs1−σ0ncosθ0n)e1, (3.4) FnM = −(σs2−σs1−σMn cosθnM)e1, (3.5) where e1= (1,0)and cosθkn= −τnk·e1.
2. Distribute the force from the markers to the fluid by
fn(x) =M
−1 k
∑
=1Fnkδh(x−Xnk)+
∑
k=0,M
Fnkδh(x−Xnk), (3.6)
where the smooth version of Dirac delta function in [22] is used.
3. Solve the Navier-Stokes equations. This can be done by the following second-order incremental pressure-correction projection method [8]. The time derivative is ap- proximated by the second-order backward difference formula and the nonlinear term is approximated by the extrapolation using previous time step values, so the incremental pressure-correction projection method can be written as follows.
3u∗−4un+un−1
2∆t + 2(un·∇h)un−(un−1·∇h)un−1
= −∇hpn+ 1
Re∇2hu∗+ 1
ReCafn, (3.7) u∗=β∂u∗
∂y ,v∗=0 on y=c, u∗=0 otherwise, (3.8)
∇2hφn+1=∇h·u∗
2∆t/3, ∂φn+1
∂n =0 on ∂Ω, (3.9)
un+1 = u∗−2∆t
3 ∇hφn+1, (3.10) pn+1 = pn+φn+1− 1
Re∇h·u∗. (3.11) Here,∇h is the standard centered difference operator on the staggered grid. One can see that the above Navier-Stokes solver involves solving two Helmholtz equa- tions for velocity u∗= (u∗,v∗)and one Poisson equation for pressure. These elliptic equations are solved by using the geometric multigrid solver.
4. Interpolate the new velocity on the fluid lattice points onto the marker points Uk= (Uk,Vk)and move the marker points Xk= (Xk,Yk)to new positions.
Unk+1 =
∑
x un+1δh(x−Xnk)h2, k=0,1,... M, (3.12) Xnk+1 = Xnk+∆t
³
Unk+1+(UA)nk+1τnk
´
, k=1,2,... M−1, (3.13) Xkn+1 = Xkn+∆tUkn+1, k=0, M, (3.14) where(UA)nk+1is computed by a mid-point integration rule of Eq. (2.17). It is worth noting that at the positions of contact lines k=0 and k=M, we only update the x coordinates by the velocity along the solid wall, and keep their y coordinates fixed as Y0n+1=YMn+1=c. This is consistent with the unbalanced Young’s forces which are applied only along the solid wall (x direction) as shown in Eqs. (3.4) and (3.5). As mentioned before, there is no artificial velocity applied at the moving contact lines;
that is,(UA)n0+1= (UA)nM+1=0.
5. Update surfactant concentration distribution Γnk++1/21 . Since the surfactant is insol- uble, 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.
For the sake of convenience, we use Sαto denote the stretching factor|Xα|. Then we discretize the modified surfactant equation (2.21) by the Crank-Nicholson scheme in a symmetric way as
Γnk++1/21 −Γnk+1/2
∆t
(Sα)nk++1/21 +(Sα)nk+1/2
2 +(Sα)nk++1/21 −(Sα)nk+1/2
∆t
Γnk++1/21 +Γnk+1/2 2
−1 2
Ã(UA)nk++11Γnk++11−(UA)nk+1Γnk+1
∆α +(UA)nk+1Γnk+1−(UA)nkΓnk
∆α
!
= 1 2Pes
1
∆α
Ã(Γnk++3/21 −Γnk++1/21 )
∆α /(Sα)nk++11−(Γnk++1/21 −Γnk−+1/21 )
∆α /(Sα)nk+1
!
+ 1 2Pes
1
∆α
Ã(Γnk+3/2−Γnk+1/2)
∆α /(Sα)nk+1−(Γnk+1/2−Γnk−1/2)
∆α /(Sα)nk
! , (3.15) where(Sα)k+1/2=((Sα)k+(Sα)k+1)/2 and Γk=(Γk−1/2+Γk+1/2)/2. Since the new in- terface marker location Xnk+1is obtained in the previous step, the above discretiza- tion 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,
∑
kΓnk++1/21 (Sα)nk++1/21 ∆α=
∑
k
Γnk+1/2(Sα)nk+1/2∆α. (3.16) The above equality can be easily derived by taking the summation of both sides of Eq. (3.15) and using the no flux boundary condition of Γ.
4 Numerical results
In this section, we present several numerical experiments to test our numerical scheme described in the previous section. The main objective is to investigate how the surfactant affects the drop motion on a solid substrate, by using both clean (without surfactant) and contaminated (with surfactant) drops in these computations. We consider a computa- tional domain Ω= [−1,1]×[0,1]where a half of circular drop with radius 0.5 is initially attached on the bottom of the domain, and both left and right contact angles are set to π/2 initially. The initial fluid velocity is also set to be zero. Using the equation of state given by Eq. (2.9), η=0 implies no surfactant, in which case we do not need to solve the surfactant equation (2.21). Throughout this section, we set σc=1. For the contaminated case, the initial surfactant concentration is uniformly distributed along the interface such that Γ(α,0) =1. Unless otherwise, we set the Reynolds number Re=10, the capillary number Ca=0.1, the surface Peclet number Pes=20. The slip length used in Navier slip boundary condition is chosen as β=h/4, one quarter of the mesh width.
4.1 Convergence test
Before we proceed, we first carry out a convergence study of the present method. Here, we perform four different computations with varying Cartesian mesh h=1/16, 1/32, 1/64, 1/128, 1/256. The Lagrangian mesh is chosen as ∆α≈h and the time step size is
∆t=h/10. The solutions are computed up to time T=6.25. The parameters are set to be σs1=0.5,σs2=1 so the equilibrium contact angle is θe=π/3 for the clean interface.
Since the analytical solution is not available in these simulations, we choose the re- sults obtained from the finest mesh as our reference solution and compute the L2 error between the reference solution and the solution obtained on the coarser grid. Table 1 shows the mesh refinement analysis of the velocity components u, v, and the surfactant concentration Γ. One can see that the error decreases when the mesh is refined, and the ra- tio between two consecutive errors is likely to approach to three indicating the first-order accuracy of the scheme [18]. (Notice that, the fluid variables are defined at the staggered grid and the surfactant concentration is defined at ”half-integer” grid, so when we refine the mesh, the numerical solutions will not coincide with the same grid locations. In these runs, we simply use a linear interpolation to compute the solutions at the desired loca- tions. We attribute this is part of the reason why the ratio shows less than three.) Table 2 shows the L∞ errors of the interface position, the cosine value of the contact angle, and the area loss of the drop for different mesh sizes.
Table 1: The mesh refinement analysis of the velocity (u, v), and surfactant concentration Γ.
h ku−ure fk2 ratio kv−vre fk2 ratio kΓ−Γre fk2 ratio 1/16 5.8079e-03 - 3.3148e-03 - 3.1818e-02 - 1/32 2.9639e-03 1.96 1.9179e-03 1.73 1.7977e-02 1.77 1/64 1.4773e-03 2.00 1.0805e-03 1.77 9.8698e-03 1.82 1/128 5.9179e-04 2.50 4.6628e-04 2.32 4.2087e-03 2.35
Table 2: The mesh refinement analysis of interface positions, the contact angles, and the area of drop.
h kX−Xre fk∞ ratio |cos(θ)−cos(θre f)| ratio |A−Are f|/Are f ratio
1/16 4.7646e-02 - 1.0285e-01 - 1.1476e-01 -
1/32 2.4225e-02 1.97 6.3948e-02 1.61 5.6724e-02 2.02 1/64 1.1426e-02 2.12 3.6915e-02 1.73 2.7473e-02 2.06 1/128 4.1700e-03 2.74 2.1708e-02 1.70 1.3145e-02 2.09
Fig. 2 shows the comparison of the local stretching factor |Xα| =Sα at T=6.25 with (solid line) and without (dotted line) the implementation of equi-distributed technique of Lagrangian markers. Without implementing the artificial tangential velocity in (2.17), the
Lagrangian markers will gradually be swept into the contact lines and become clustered there. Here, we only show the case with surfactant η=0.3. One can see the present equi- distributed technique does preserve the Lagrangian markers more uniformly by keeping the stretching factor less varied along the interface.
0 1
0 0.1 0.2 0.3 0.4 0.5
Sα at t = 6.25
α
Figure 2: Comparison of stretching factor|Xα| =Sαwith (UA=0, dotted line) and (UA6=0, solid line), where h=1/128.
In the following subsections, we investigate the effect of surfactant on the motion of contact lines when the substrate is hydrophilic, hydrophobic and both, from the view point of the drop (fluid 1). We show that in all cases, the presence of surfactant can significantly affect the contact line dynamics.
4.2 Hydrophilic case
To examine the effect of the surfactant on interfacial dynamics, we compare a hydrophilic drop with and without surfactant in a quiescent flow. When the surfactant is present on the interface, surface tension is reduced significantly, cf. equation of state (2.9). The parameters are set to be σs1=0.5,σs2=1 so the equilibrium contact angle is θe=π/3 for the clean interface. The Cartesian mesh size is h=1/128, and the Lagrangian grid size is
∆α≈h. The time step size is set to be ∆t=h/10.
Fig. 3 shows the time evolution of a hydrophilic drop in a quiescent flow field. Here, we distinguish the clean (η=0) and contaminated (η=0.3) drop interfaces by the “dashed”
and “solid” lines, respectively. As expected, both drops start spreading on the solid sub- strate. The clean drop moves gradually to a state with a contact angle approaching to equilibrium θ≈π/3, while the contaminated one spreads further to reach a state with θ≈π/4. This drop behavior can be easily explained by Eq. (2.1). The surfactant is insolu- ble and affects only the surface tension between fluids. Therefore, both σs2and σs1remain the same and so is their difference, but the surface tension σ at the right and left contact lines are smaller than σc. As a result, the dynamic contact angle θ becomes a function of t as θ(t) =arccos((σs2−σs1)/σ) ≤π/3. Thus, the contaminated drop becomes more
hydrophilic than the clean one. It is also interesting to mention that the larger value of η, the more hydrophilic the drop becomes since the surface tension σ decreases more.
−0.5 0 0.5
0 0.2 0.4
T = 0
−0.5 0 0.5
0 0.2 0.4
T = 2.3438
−0.5 0 0.5
0 0.2 0.4
T = 4.6875
−0.5 0 0.5
0 0.2 0.4
T = 7.0313
−0.5 0 0.5
0 0.2 0.4
T = 9.375
−0.5 0 0.5
0 0.2 0.4
T = 12.5
Figure 3: The time evolution of a hydrophilic drop with clean (η=0, dashed line) and contaminated interface (η=0.3, solid line).
Fig. 4 shows the velocity plots for the drop with surfactant near the left and the right contact lines at time T=1.5625. One can see significant contact line velocity pointing outward to the dry region (fluid 2) has been observed during the wetting process.
Plots of the corresponding surfactant concentration (top column) and surface tension (bottom column) versus arc-length for the contaminated case are given in Fig. 5. It can be observed although the initial surfactant is uniformly distributed, concentration at the contact lines becomes larger than that in other locations of the interface. This is because during the wetting process, surfactant has been swept into the contact lines. As a re- sult, surface tension near the contact lines is reduced. Eventually, surfactant redistribute due to diffusion and its concentration reaches a uniform distribution again. The present method preserves the total surfactant mass exactly and the error reaches the machine precision, i.e.|m(t)−m(0)| ≈10−14.
In Fig. 6, we present time evolutionary plots for the comparison of clean and con- taminated drops; namely, (a) left contact line speed; (b) right contact line speed; (c) the contact angles in units of π; and (d) the total length of the drop. Fig. 6-(a) and (b) show the evolution of left and right contact line speed, respectively. Since the initial contact
−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 0
0.1 0.2 0.3 0.4 0.5 0.6
velocity field at T = 1.5625
−0.7 −0.6 −0.5 −0.4 −0.3
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
left contact angle
0.3 0.4 0.5 0.6 0.7
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
right contact angle
Figure 4: The velocity field for the drop with surfactant near the left and right contact lines (η=0.3,T=1.5625).
0 1
0.8 0.85 0.9 0.95 1 1.05
T = 0
0 1
0.64 0.66 0.68 0.7 0.72 0.74
T = 0
0 1
0.8 0.85 0.9 0.95 1 1.05
T = 2.3438
0 1
0.64 0.66 0.68 0.7 0.72 0.74
T = 2.3438
0 1
0.8 0.85 0.9 0.95 1 1.05
T = 4.6875
0 1
0.64 0.66 0.68 0.7 0.72 0.74
T = 4.6875
0 1
0.8 0.85 0.9 0.95 1 1.05
T = 7.0313
0 1
0.64 0.66 0.68 0.7 0.72 0.74
T = 7.0313
0 1
0.8 0.85 0.9 0.95 1 1.05
T = 9.375
0 1
0.64 0.66 0.68 0.7 0.72 0.74
T = 9.375
0 1
0.8 0.85 0.9 0.95 1 1.05
T = 12.5
0 1
0.64 0.66 0.68 0.7 0.72 0.74
T = 12.5
Figure 5: The distribution of the surfactant concentration (top) and the corresponding surface tension (bottom).
angle deviates the equilibrium angle more, the unbalanced Young’s force is larger in the beginning and so is the contact line speed. Notice that, this is a wetting case so the con- tact line speed for left angle is negative while the right one is positive. Furthermore, the magnitude of contact line speed for the contaminated drop is larger than the clean one since the contaminated drop wets more. As time evolves, the contact lines slow down and approach to zero when the drop reaches its steady state. The corresponding evolu- tions for left and right contact angles in units of π are shown in Fig. 6-(c). One can hardly
distinguish the left and right angles since they are matched perfectly well in the plot due to the symmetric setting of the drop. The final dynamic contact angle of the clean drop approaches to π/3 as expected, while the contaminated one approaches to π/4. Fig. 6- (d) confirms that the drop with surfactant wets more than the one without surfactant due to the increase of total length of the interface.
0 5 10 15
1.55 1.6 1.65 1.7 1.75 1.8 1.85
(d)
0 5 10 15
−0.08
−0.04 0 0.04
0.08 (b)
0 5 10 15
−0.08
−0.04 0 0.04
0.08 (a)
0 5 10 15
0.25 0.3 0.35 0.4 0.45 0.5
(c)
Figure 6: The comparison of a hydrophilic drop with clean (η=0, dashed line) and contaminated interface (η=0.3, solid line). (a) left contact line speed; (b) right contact line speed; (c) contact angles in units of π;
(d) total length of the drop.
4.3 Hydrophobic case
In this case, we keep almost the same setup as the previous example but change the surface tension to σs1=1,σs2=0.1557 so that the equilibrium contact angle for a clean interface is around θe=0.82π which corresponds to a hydrophobic drop.
Fig. 7 shows the time evolution plots of a hydrophobic drop in a quiescent flow field.
Again, we distinguish the clean (η=0) and contaminated (η=0.3) drop interfaces by the
“dashed” and “solid” lines, respectively. In this case, both drops start to contract. From Eq. (2.1), the dynamic contact angle θ evolves following θ(t) =arccos((σs2−σs1)/σ) ≥ 0.82π so the contaminated drop becomes more hydrophobic due to the reduction of sur- face tension. Thus, while the clean drop moves gradually to approach to a state with equilibrium contact angle, the contaminated one goes further to become a nearly non- wetting state. This confirms exactly what we describe in the Introduction section. Once again, the larger value of η, the more hydrophobic the drop becomes as the surface ten- sion σ is reduced by contamination.
Fig. 8 shows the plots for the drop with surfactant near the left and the right con- tact lines at time T=2.3438. Significant contact line velocity pointing inward to the wet region (fluid 1) has been observed during the non-wetting process. Plots of the corre- sponding surfactant concentration and surface tension are similar to the hydrophilic case and omitted here.
−0.6 −0.4 −0.20 0 0.2 0.4 0.6 0.2
0.4 0.6
T = 0
−0.6 −0.4 −0.20 0 0.2 0.4 0.6 0.2
0.4 0.6
T = 2.3438
−0.6 −0.4 −0.20 0 0.2 0.4 0.6 0.2
0.4 0.6
T = 4.6875
−0.6 −0.4 −0.20 0 0.2 0.4 0.6 0.2
0.4 0.6
T = 7.0313
−0.6 −0.4 −0.20 0 0.2 0.4 0.6 0.2
0.4 0.6
T = 9.375
−0.6 −0.4 −0.20 0 0.2 0.4 0.6 0.2
0.4 0.6
T = 16.4063
Figure 7: The time evolution of a hydrophobic drop with clean (η=0, dashed line) and contaminated interface (η=0.3, solid line).
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
velocity field at T = 2.3438
−0.6 −0.5 −0.4 −0.3 −0.2
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
left contact angle
0.2 0.3 0.4 0.5
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
right contact angle
Figure 8: The velocity field for the drop with surfactant near the left and right contact lines (η=0.3,T=2.3438).
In Fig. 9, we present time evolutionary plots for the comparison of clean and con- taminated hydrophobic drops; namely, (a) left contact line speed; (b) right contact line speed; (c) the contact angles in units of π; and (d) the total length of the drop. Fig. 9-(a) and (b) show the evolution of left and right contact line speed, respectively. As in the
hydrophilic case, the initial contact angle deviates the equilibrium angle more, the un- balanced Young’s force is larger in the very beginning and so is the contact line speed.
Now, this is a hydrophobic case so the contact line speed for left angle is positive while the right one is negative. Furthermore, the magnitude of contact line speed for the con- taminated drop is larger than the clean one since the contaminated drop becomes more non-wetting. As time evolves, the contact lines slow down and approach to zero when the drop reaches its steady state. The corresponding evolutions for left and right contact angles in units of π are shown in Fig. 9-(c). The final dynamic contact angle of the clean drop approaches to 0.82π as expected, while the contaminated one approaches to π.
0 5 10 15 20 25
1.5 2 2.5
(d)
0 5 10 15 20 25
−0.15
−0.1
−0.05 0
(b)
0 5 10 15 20 25
0 0.05 0.1 0.15
(a)
0 5 10 15 20 25
0.5 0.6 0.7 0.8 0.9 1
(c)
Figure 9: The comparison of a hydrophobic drop with clean (η=0, dashed line) and contaminated interface (η=0.3, solid line). (a) left contact line speed; (b) right contact line speed; (c) contact angles in units of π;
(d) total length of the drop.
4.4 Hydrophobic-hydrophilic case
Now, we consider a substrate that has a spatially inhomogeneous wettability where σs2(x),x∈ [−1,1]is given by
σs2(x) =
0.191, Hydrophobic x∈I1= [−1,−0.65] 5.59(x+0.65)+0.191, x∈I2= (−0.65,−0.45) 1.309, Hydrophilic x∈I3= [−0.45,1].
We choose σs1=1 so that the static contact angles (clean interface) for σs2=0.191 and 1.309 are θe=4π/5 and 2π/5, respectively. This set up mimics the situation where the solid surface has a combining hydrophobic (zone I1) and hydrophilic (zone I3) feature.
Initially, a half circular drop with radius 0.5 and centered at (−0.25,0) is placed in the
domain so that the left and right contact lines are located in zone I1and I3, respectively, cf. Fig. 10.
−10 −0.75 −0.65 −0.45 −0.25 0 0.25 1
0.5 1
Zone I2 Zone I1
Zone I 3
Figure 10: Wettability and the initial drop set up.
Along the circular interface, we choose a nonuniform initial surfactant concentration as
Γ(α) = − 8
5π3α3+ 12
5π2α2+0.5, α∈ (0,π).
Unlike previous two cases, the present set up produces much rich dynamic behavior.
First of all, due to the asymmetric nature of the surface property, the drop behaves dif- ferently near the left and right contact lines. Fig. 11 shows the time evolution of the drop with (η=0.3, solid line) and without (η=0, dashed line) surfactant effect. It can be seen that the motion of the contact lines is mainly driven by unbalanced Young’s force initially. The left contact angle deviates from its equilibrium more than the right con- tact angle, due to a greater unbalanced Young’s force there. Therefore, it moves with a greater speed than the right one, as shown in the top row of Fig. 12, which is also visible in Fig. 14-(a) and (b). Since the static contact angle is greater at the left contact line (less wetting), the center of the mass inside the drop shift towards the left direction, as shown in Fig. 11 and Fig. 12. However, when the left contact line arrives at zone I2(a transition region from hydrophobic to hydrophilic), the unbalanced Young’s force is reduced, and so is the corresponding contact line speed. Meanwhile, the fluid inside is driven by the Marangoni effect and the moving contact lines, it flows from the left contact line to the right one and push right contact angle away from its equilibrium value, which increases the Young’s force and leads an increasing contact line speed, see Fig. 14-(a) and (b). Fi- nally, after the left contact line crosses the middle of I2 (around t=3.6), the entire drop is in the hydrophilic zone, and the drop behaves in a similar way as in the hydrophilic
