成果報告
交通大學應用數學系
賴明治 教授
中文摘要：
流體界面問題旨在探討兩種不相混合流體(或者空氣與水)的界面問
題。它在自然界的物理現象或者微流體工業技術中扮演了重要角色。譬如，
薄膜成長的技術，噴墨式印表機的液滴形成、溼化問題，甚至水面上行走
昆蟲的流力問題，皆是界面現象的貢獻。由於它的高度複雜性及應用性，
過去一直深受應用數學家廣泛的研究，這些努力也帶動了對流體界面問題
的更深一層的了解與認識。
在這份報告裡，我們將節錄部分有關界面流體模擬的研究成果，特別
是針對二維溼化問題提出相關的數學模型及其數值方法，其數值模擬結果
顯示與實際物理現象在行為上相當吻合，本結果連同其它成果也陸續發表
在SCI期刊中，在此特別感謝國科會計畫NSC 98-2115-M-009-014-MY3之支
助。
關鍵字：沈浸邊界法, 界面流, Navier-Stokes方程, 界面活性劑, 移動接觸線,
親水液滴, 疏水液滴, 溼化現象
Interfacial flows with surfactant
Ming-Chih Lai
∗Abstract
The study of the incompressible flows with interfaces is of major interests among applied mathematics community. It plays an important role in numerous natural phe-nomena and industrial applications, especially, for the hydrodynamics of micro-fluidic systems. For instance, the thin film flows in coating devices, the dynamics of liquid drops in ink-jet printing, the wetting phenomena on substrate, or even the hydrody-namics of water-walking insect, just to name a few. Since the complexity of effects associated with capillarity gains a more wide applications during the last century, the research effort spent on this topic is overwhelming in applied mathematics community. These effort brought in many remarkable successes in the understanding of the fun-damental physics and the quantitative modeling of capillary phenomena in different problems. In this report, we shall focus our research results on the different issues of mathematical modeling and numerical methods for interfacial flows with surfactant. In particular, we shall report the detailed results in the 2D numerical simulation of moving contact lines with surfactant. Those results were supported by our previous NSC projects, NSC-97-2628-M-009-007-MY3 and NSC-98-2115-M-009-014-MY3.
Keywords: Immersed boundary method; Interfacial flow; Navier-Stokes equations; Sur-factant; Moving contact line; Hydrophilic drop; Hydrophobic drop; Wetting
1
Introduction
In this report, 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 Lagrangian interfacial markers. The immersed boundary force has two components: one from the nonhomogeneous surface tension determined by the distribution of surfactant along the fluid interface, and the other from unbalanced Young’s force at the moving 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 corresponding modified surfactant equation is solved in a way such that the total surfactant 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.
∗Corresponding author. Center of Mathematical Modeling and Scientific Computing & Department of
Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, TAIWAN. E-mail: [email protected]
Surfactant are surface active agents that adhere to the fluid interface and reduce the interface 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 system, surfactant allow small droplets to be formed and used as an emulsion. Surfactant 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 example, 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 surfactant 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 intersection 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 interfaces with insoluble surfac-tant [15]. The mathematical models include the incompressible Navier-Stokes equations and a surfactant concentration equation along the moving interface. 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 tangential 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 biological applications [21, 3, 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 [12, 5, 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 attempts 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 contact 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 remove 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 [32, 19, 7]. 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. Recently, 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 de-veloped 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 Carte-sian 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 Lagrangian 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 section 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.
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, (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. θe is 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 [22, 15], ∂u ∂t + (u· ∇)u + ∇p = 1 Re∇ 2u + 1 ReCaf , in Ω (2) ∇ · u = 0, in Ω (3) f (x, t) = ∫ Σ F (α, t) δ(x− X(α, t)) dα, (4) F (α, t) = ∂ ∂α(σ τ ), τ = ∂X ∂α / ∂X∂α , (5) ∂X(α, t) ∂t = U (α, t) = ∫ Ω u(x, t) δ(x− X(α, t))dx. (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 in-ertial force and the viscous force, and the capillary number (Ca) describing the strength of
the surface tension. Equations (4) and (6) represent the interaction between the immersed interface and the fluids. In particular, Eq. (4) describes the force (f ) acting on the fluid due to the fluid interfacial force F and Eq. (6) states that the interface moves with the fluid velocity. The fluid interfacial force is given by Eq. (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, (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 θ, (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− η Γ)) (9)
where σc is the surface tension of a clean interface, and η satisfying 0 ≤ η < 1 is a di-mensionless number that measures the sensitivity of surface tension to changes in surfactant concentration. The surfactant concentration equation [30] is
DΓ Dt + (∇s· U) Γ = 1 P es ∇2 sΓ, (10)
where ∇s and ∇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α|, (11) ∇sΓ = ∂Γ ∂ττ = ( ∂Γ ∂ατ ) /|Xα|, (12) ∇2 sΓ = ∇s· ∇sΓ = ∂ ∂α ( ∂Γ ∂α/|Xα| ) /|Xα|. (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 redistribution 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 immersed 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 formulation 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. (6) as
∂X(α, t)
∂t = U (α, t) + U
A(α, t) τ , (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(L
b, t) = 0. Here, we need to impose
|Xα|α = 0 so that |Xα| = L1 b
∫Lb
0 |Xα′| dα′. Taking the derivative with respect to t, we have
|Xα|t = L1b ∫Lb 0 |Xα′|tdα′. Since |Xα|t = Xαt· Xα |Xα| = ( ∂U ∂α + ∂UA ∂α τ + U A∂τ ∂α ) · τ = ∂U ∂α · τ + ∂UA ∂α = (∇s· U) |Xα| + ∂UA ∂α (15) we have ∂U ∂α · τ + ∂UA ∂α = 1 Lb ∫ Lb 0 ( ∂U ∂α′ · τ ′+∂UA ∂α′ ) dα′. (16)
Integrating with respect to α and using the fact of UA(0, t) = UA(L
b, t) = 0, we obtain UA(α, t) = α Lb ∫ Lb 0 ∂U ∂α′ · τ ′dα′−∫ α 0 ∂U ∂α′ · τ ′dα′. (17)
We note that an alternative formula for using artificial tangential velocity has been derived by Ceniceros [1] in a similar front-tracking manner. However, unlike the approach in [1], the present formula (17) does not need to compute the curvature of the interface.
2.2
Modified surfactant concentration equation
By taking the artificial velocity UA into account, the material derivative now becomes
DΓ Dt = ∂Γ ∂t − U A τ · ∇sΓ. (18)
Therefore, the surfactant equation (10) becomes
∂Γ ∂t − U Aτ · ∇ sΓ + (∇s· U) Γ = 1 P es∇ 2 sΓ. (19)
Multiplying the stretching factor |Xα| on both sides of the above equation, we obtain
∂Γ ∂t |Xα| − U A∂Γ ∂α + (∇s· U) |Xα| Γ = 1 P es ∇2 sΓ|Xα|. (20)
By writing the surface derivatives in terms of α explicitly and using the identity of Eq. (15), we obtain ∂Γ ∂t|Xα| + ∂|Xα| ∂t Γ− ∂(UAΓ) ∂α = 1 P es ∂ ∂α ( ∂Γ ∂α/|Xα| ) . (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 UA is 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 components 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 Lagrangian 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, respectively. 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 interface 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)), (22) τnk = DαX n k |DαXnk| , (23)
both of which hold for αk = k∆α, and Γnk is simply approximated by the average of Γnk−1/2and Γnk+1/2. Then we compute the interfacial force at the the fluid/fluid markers
Xk, k = 1, . . . , M − 1 by
Fnk = Dα(σknτ n
k) ∆α. (24)
The unbalanced Young’s force at the contact lines (k = 0 and k = M ) can be computed by Fn0 = (σs2− σs1− σ0n cos θ n 0) e1, (25) FnM = −(σs2− σs1− σMn cos θ n M) e1, (26)
where e1 = (1, 0) and cos θkn=−τ n k · e1.
2. Distribute the force from the markers to the fluid by
fn(x) = M∑−1 k=1 Fnkδh(x− Xnk) + ∑ k=0,M Fnkδh(x− Xnk), (27)
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 approx-imated 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(u n· ∇ h)un− (un−1· ∇h)un−1 = −∇hpn+ 1 Re∇ 2 hu∗+ 1 Re Caf n , (28) u∗ = β∂u ∗ ∂y , v ∗ = 0 on y = c, u∗ = 0 otherwise, (29) ∇2 hϕ n+1 = ∇h· u ∗ 2∆t/3 , ∂ϕn+1 ∂n = 0 on ∂Ω, (30) un+1 = u∗− 2∆t 3 ∇hϕ n+1, (31) pn+1 = pn+ ϕn+1− 1 Re∇h· u ∗. (32)
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 equations 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.
Un+1k = ∑ x un+1δh(x− Xnk)h 2, k = 0, 1, . . . M, (33) Xn+1k = Xnk+ ∆t(Un+1k + (UA)n+1k τnk), k = 1, 2, . . . M − 1, (34) Xkn+1 = Xkn+ ∆t Ukn+1, k = 0, M, (35) where (UA)n+1
k is computed by a mid-point integration rule of Eq. (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. (25) and (26). As mentioned before, there is no artificial velocity applied at the moving contact lines; that is, (UA)n+1
0 = (UA)
n+1
M = 0.
5. Update surfactant concentration distribution Γn+1k+1/2. 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.
For the sake of convenience, we use Sα to denote the stretching factor |Xα|. Then we discretize the modified surfactant equation (21) by the Crank-Nicholson scheme in a symmetric way as Γn+1k+1/2− Γn k+1/2 ∆t (Sα)n+1k+1/2+ (Sα)nk+1/2 2 + (Sα)n+1k+1/2− (Sα)nk+1/2 ∆t Γn+1k+1/2+ Γn k+1/2 2 −1 2 ( (UA)n+1 k+1Γ n+1 k+1− (UA) n+1 k Γ n+1 k ∆α + (UA)nk+1Γnk+1− (UA)nkΓnk ∆α ) = 1 2 P es 1 ∆α ( (Γn+1k+3/2− Γn+1k+1/2) ∆α /(Sα) n+1 k+1− (Γn+1k+1/2− Γn+1k−1/2) ∆α /(Sα) n+1 k ) + 1 2 P es 1 ∆α ( (Γn k+3/2− Γ n k+1/2) ∆α /(Sα) n k+1− (Γn k+1/2− Γ n k−1/2) ∆α /(Sα) n k ) , (36)
where (Sα)k+1/2 = ((Sα)k+ (Sα)k+1)/2 and Γk = (Γk−1/2+ Γk+1/2)/2. Since the new interface marker location Xn+1k is 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 Γn+1k+1/2(Sα)n+1k+1/2∆α = ∑ k Γnk+1/2(Sα)nk+1/2∆α. (37)
The above equality can be easily derived by taking the summation of both sides of Eq. (36) and using the no flux boundary condition of Γ.
4
Numerical results
In this section, we present several numerical experiments to test our numerical scheme de-scribed 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 computational 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. (9),
η = 0 implies no surfactant, in which case we do not need to solve the surfactant equation
(21). Throughout this section, we set σc = 1. For the contaminated case, the initial surfac-tant 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 P es = 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 results 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 ratio 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 locations. 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.
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 (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.
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
Table 1: The mesh refinement analysis of the velocity (u, v), and surfactant concentration Γ.
h ∥u − uref∥2 ratio ∥v − vref∥2 ratio ∥Γ − Γref∥2 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 ∥X − Xref∥∞ ratio | cos(θ) − cos(θref)| ratio |A − Aref|/Aref 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
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 (UA̸= 0, solid line), where h = 1/128.
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 (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 substrate. 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. (1). The surfactant is insoluble and affects only the surface tension between fluids. Therefore, both σs2 and σs1 remain 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 result, surface tension near the contact lines is reduced. Eventually, surfactant redistribute due to
−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).
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. 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).
In Fig. 6, we present time evolutionary plots for the comparison of clean and contaminated 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 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 contact 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 evolutions 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 contam-inated 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. (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 surface 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 tension σ is reduced by contamination.
Fig. 8 shows the plots for the drop with surfactant near the left and the right contact 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 corresponding 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 contaminated 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 unbalanced 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 contaminated 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.
References
[1] H. D. Ceniceros, The effects of surfactants on the formation and evolution of capillary waves, Phys. Fluids, 15(1), (2003) 245–256.
[2] S.H. Davis, G.K. Batchelor, H.K. Moffatt, M.G. Worster (Eds.), Perspectives in Fluid
Dynamics, Cambridge University Press, Cambridge, 2002.
[3] P.-G. de Gennes, F. Brochard-Wyart, D. Quere, Capillarity and Wetting Phenomena, Springer, 2003.
[4] E. B. Dussan V., The moving contact line: the slip boundary condition, J. Fluid Mech., 77, (1976), 665–684.
[5] E. B. Dussan V., On the spreading of liquids on a solid surface: static and dynamic contact lines, Ann. Rev. Fluid Mech., 11, (1979), 371–400.
[6] P. G. de Gennes, Wetting: statics and dynamics, Rev. Mod. Phys., 57(3), (1985) 827– 863.
[7] S. Ganesan, L. Tobiska, Modelling and simulation of moving contact line problems with wetting effects, Comput. Visual Sci., DOI 10.1007/s00791-008-0111-3, (2008).
[8] J. L. Guermond, P. Minev, J. Shen, An overview of projection methods for incompress-ible flows, Comput. Methods Appl. Mech. Eng. 195 (2006) 6011-6045.
[9] F. H. Harlow, J. E. Welsh, Numerical calculation of time-dependent viscous incompress-ible flow of fluid with a free surface, Phys. Fluids, 8 (1965) 2181-2189.
[10] T. Y. Hou, J. S. Lowengrub, M. J. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys. 114 (1994) 312-338.
[11] H. Huang, D. Liang, B. Wetton, Computation of a moving drop/bubble on a solid surface using a front-tracking method, Comm. Math. Sci. 2(4) (2004) 535–552.
[12] C. Huh, L. E. Scriven, Hydrodynamic model of steady movement of a solid/liquid/fluid contact line, J. Colloid Interf. Sci. 35 (1971), 85–101.
[13] D. Jacqmin, Contact-line dynamics of a diffuse fluid interface, J. Fluid Mech., 402 (2000), 57–88.
[14] A. J. James, J. S. Lowengrub, A surfactant-conserving volume-of-fluid method for in-terfacial flows with insoluble surfactant, J. Comput. Phys., 201 (2004) 685–722.
[15] M.-C. Lai, Y.-H. Tseng, H. Huang, An immersed boundary method for interfacial flow with insoluble surfactant, J. Comput. Phys. 227 (2008) 7279-7293.
[16] J. Lee, C. Pozrikidis, Effect of surfactants on the deformation of drops and bubbles in Navier-Stokes flow, Comput. Fluids, 35 (2006) 43–60.
[17] V. G. Levich, Physicochemical Hydrodynamics, Prentice Hall, 1962.
[18] Z. Li, K. Ito, The Immersed Interface Method, Frontiers in Applied Mathematics, SIAM (2006), pp 230.
[19] H. Liu, S. Krishnan, S. Marella, H. S. Udaykumar, Sharp interface Cartesian grid method II: A technique for simulating droplet interactions with surfaces of arbitrary shape, J. Comput. Phys., 210 (2005), 32–54.
[20] M. Muradoglu, G. Tryggvason, A front-tracking method for computations of interfacial flows with soluble surfactants, J. Comput. Phys., 227 (2008) 2238–2262.
[21] A. Oron, S. H. Davis, S. G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod.
Phys., 69 (1997), 931–980.
[23] C. S. Peskin, B. F. Printz, Improved volume conservation in the computation of flows with immersed elastic boundaries, J. Comput. Phys. 105 (1993) 33–36.
[24] T. Qian, X.-P. Wang, P. Sheng, Molecular scale contact line hydrodynamics of immis-cible flows, Phys. Rev. E 68, (2003) 016306.
[25] T. Qian, X.-P. Wang, P. Sheng, Power-law slip profile of the moving contact line in two-phase immiscible flows, Phys. Rev. Lett. 93(9), (2004) 094501.
[26] T. Qian, X.-P. Wang, P. Sheng, Molecular hydrodynamics of the moving contact line in two-phase immiscible flows, Commun. Comput. Phys., Vol 1, No 1, pp. 1–52 (2006). [27] W. Ren, W. E, Boundary conditions for moving contact line problem, Phys. Fluids 19,
(2007) 022101.
[28] M. Renardy, Y. Renardy, J. Li, Numerical simulation of moving contact line problems using a volume-of fluid method, J. Comput. Phys., 171 (2001), 243–263.
[29] Y. D. Shikhmurzaev, Capillary Flows with Forming Interfaces, Chapman & Hall/CRC Press, 2008.
[30] H. A. Stone, A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface, Phys. Fluids A 2(1), (1990) 111–112. [31] S. Sikalo, H.-D. Wilhelm, I. V. Roisman, S. Jakirlic, C. Tropea, Dynamic contact angle
of spreading droplets: Experiments and simulations, Phys. Fluids 17, (2005) 062103. [32] P. D.M. Spelt, A level-set approach for simulations of flows with multiple moving contact
lines with hysteresis, J. Comput. Phys. 207 (2005) 389–404.
[33] P. Tabeling, Introduction to Microfluidics, Oxford University Press, 2005.
[34] K. Yokoi, D. Vadillo, J. Hinch, I. Hutchings, Numerical studies of the influence of the dynamic contact angle on a droplet impacting on a dry surface, submitted for publication (2009).
[35] J.-J. Xu, Z. Li, J. S. Lowengrub, H.-K. Zhao, A level-set method for interfacial flows with surfactant, J. Comput. Phys., 212 (2006) 590–616.
