An efficient scheme for a phase field model for the moving contact line problem with variable density and viscosity

Accepted Manuscript An efficient Scheme for a Phase Field model for the moving contact line problem with variable density and viscosity Min Gao, Xiao-Ping Wang

PII: DOI: Reference:

S0021-9991(14)00330-1 10.1016/j.jcp.2014.04.054 YJCPH 5244

To appear in:

Journal of Computational Physics

Received date: 24 August 2013 Revised date: 5 March 2014 Accepted date: 29 April 2014

Please cite this article in press as: M. Gao, X.-P. Wang, An efficient Scheme for a Phase Field model for the moving contact line problem with variable density and viscosity, Journal of Computational Physics (2014), http://dx.doi.org/10.1016/j.jcp.2014.04.054

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An efficient Scheme for a Phase Field Model for the Moving Contact Line Problem with Variable Density and Viscosity Min Gaoa , Xiao-Ping Wanga,∗ a Department

of Mathematics, The Hong Kong University of Science and Technology Clear Water Bay, Kowloon, Hong Kong, China

Abstract In this paper, we develop an efficient numerical method for the two phase moving contact line problem with variable density, viscosity, and slip length. The physical model is based on a phase field approach, which consists of a coupled system of the Cahn-Hilliard and Navier-Stokes equations with the generalized Navier boundary condition [1, 2, 5]. To overcome the difficulties due to large density and viscosity ratio, the Navier-Stokes equations are solved by a splitting method based on a pressure Poisson equation [19], while the Cahn-Hilliard equation is solved by a convex splitting method. We show the method is stable under certain conditions. The linearized schemes are easy to implement and introduce only mild CFL time constraint. Numerical tests are carried out to verify the accuracy, stability and efficiency of the schemes. The method allows us to simulate the interface problems with extremely small interface thickness. Three dimensional simulations are included to validate the efficiency of the method. Keywords: Navier-Stokes equations; Cahn-Hilliard equation; Phase field; Moving contact line; Large density ratio; Pressure stabilization scheme

1. Introduction Moving contact line problem, where the fluid-fluid interface intersects the solid wall, is a classical problem that occurs in many physical phenomena. It is well known that classical hydrodynamical models with noslip boundary condition lead to nonphysical singularity in the vicinity of the contact line [3]. A phase field model with generalized Navier boundary condition (GNBC) is proposed in [1, 3] to resolve the issue. It is shown that the numerical results based on the GNBC can reproduce quantitatively the results from the MD simulation. This indicates that the new model can accurately describe the behavior near the contact line. ∗ Corresponding

author Email addresses: [email protected] (Min Gao), [email protected] (Xiao-Ping Wang)

Preprint submitted to J. Comput. Phys.

May 6, 2014

The model involves a coupled system of the Cahn-Hilliard and Navier-Stokes equations. ∂φ + v · ∇φ = M ∇2 μ, ∂t   ∂v + (v · ∇)v = −∇p + ∇ · [ηD(v)] + μ∇φ + ρgext , ρ ∂t ∇ · v = 0,

(1.1) (1.2) (1.3)

here p is the pressure, ηD(v) = η(∇v + ∇vT ) denotes the viscous part of the stress tensor, ρ, η are the fluid mass density and viscosity; ρgext is the external body force density, and M is the phenomenological mobility coefficient; μ = −K∇2 φ − rφ + uφ3 is the chemical potential, and μ∇φ is the capillary force; K, r, u  are the parameters that are related to the interface profile thickness ξ = K/r, the interfacial tension  √ γ = 2 2r2 ξ/3u, and the two homogeneous equilibrium phases φ± = ± r/u (= ±1 in our case). To describe the motion of the contact line, the equation (1.2) is supplemented with the generalized Navier boundary condition at top and bottom boundaries (Fig. 1). βvxslip = −η∂n vx + L(φ)∂x φ,

(1.4)

here L(φ) = K∂n φ + ∂γwf (φ)/∂φ, and γwf (φ) = − 12 γ cos θssurf sin( π2 φ), θssurf is the static contact angle; β is the slip coefficient; vxslip = vx − uw is the slip velocity, uw is the wall speed. The velocity field is denoted by v = (vx , vz ) (see Fig. 1), where vx , vz are velocities along x, z directions, n, τ are unit vectors orthogonal and tangential to the boundaries, vn := v · n, vτ := v · τ . In addition, a relaxation boundary condition is imposed on the phase field variable φ at the top and bottom boundaries: ∂φ + vx ∂x φ = −Γ[L(φ)], ∂t

(1.5)

where Γ is a (positive) phenomenological parameter, together with the following impermeability conditions: vz = 0, ∂n μ = 0.

(1.6)

Here the variable density, viscosity and slip length are taken to be the volume average of those for the two flows, i.e., ρ = ρL (

1+φ 1−φ 1+φ 1−φ ) + ρG ( ), η = ηL ( ) + ηG ( ). 2 2 2 2

Although this may be inconsistent with the continuity equation ρt + ∇ · (ρv) = 0 (when ρ is not a constant), it is shown in [14] that the total mass is still conserved if the boundary conditions (1.6) are imposed. There have been many work on developing efficient numerical methods for two phase flow with general variable density and viscosity [7, 14, 17–19, 22–24]. However, most of the work were on models for problems where the interface does not intersect with the boundary. The main difficulty in those problems comes from the high (fourth) order derivatives and strong nonlinearity in the Cahn-Hilliard equation which introduces 2

uw 1111111111111111111111111 0000000000000000000000000 0000000000000000000000000 1111111111111111111111111

z 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000

−uw

o

x

Figure 1: Two-phase Couette flow with wall speed uw .

a strong stability constraint for the time step. Extra complexity is introduced in the moving contact line model due to the generalized Navier boundary condition (1.4). For two phase system with equal density and viscosity, we have developed an efficient gradient stable scheme [5]. The scheme is unconditionally stable and has the property of total energy decaying. Another efficient scheme based on a least square/finite element method is also proposed in [6]. Solving Navier-Stokes equations with large density ratios using projection method can be very time comsuming since it requires solving a variable-coefficient stiffness matrix (with very large condition number) at each time step. Guermond and Salgado [19] proposed a splitting method for Navier-Stokes equations with variable density based on penalty techniques, which only requires solving pressure Poisson equation with constant coefficient per time step. In this paper, we develop a gradient stable scheme for the system (1.1)-(1.6). The scheme is based on a convex splitting of the bulk free energy functional and the surface energy, and a splitting method based on the pressure Poisson equation for Navier-Stokes equations. We show, under certain conditions, the scheme has the total energy decaying property and is stable. Numerical tests are carried out to verify the stability, accuracy and efficiency of the scheme. We also compared the performance of two types of Navier-Stokes solvers. It is verified that, the splitting scheme based on pressure Poisson equation has similar accuracy as the projection method, but the splitting scheme is significantly more efficient in terms of computational cost, making it more favorable in large scale simulations. The rest of the paper is organized as follows. In section 2, we derive the energy law for the PDE system. The numerical scheme and discrete energy law are derived in section 3. Numerical tests are performed and the results are analyzed in section 4. The paper concludes in section 5 with a few remarks.

3

2. Energy law for the system (1.1)-(1.6) We first rewrite (1.1)-(1.6) into the dimensionless form. We scale length by some length scale L, φ by  r/u (=1, in this paper), velocity by the wall speed V , time by L/V , density by liquid density ρL , viscosity by liquid viscosity ηL and pressure by ηL V /L. Then ∂φ + v · ∇φ = Ld ∇2 μ, ∂t ∂v Rρ[ + (v · ∇)v] = −∇p + ∇ · [ηD(v)] + Bμ∇φ + fext , ∂t

(2.1) (2.2)

∇ · v = 0,

(2.3)

here μ = −∇2 φ − φ/ + φ3 / is the chemical potential,  is the ratio between interface thickness ξ and characteristic length L; density ρ =

1+φ 2

+ λρ 1−φ 2 , viscosity η =

1+φ 2

+ λη 1−φ 2 , λρ = ρG /ρL and λη = ηG /ηL

are density and viscosity ratios; external force fext = RGρgext . The boundary conditions at top and bottom boundaries are ∂φ + vx ∂x φ = −Vs [L(φ)], ∂t

(2.4)

vz = 0, ∂n μ = 0,

(2.5)

[Ls ls ]−1 vxslip = BL(φ)∂x φ/η − ∂n vx ,

(2.6)

here L(φ) = ∂n φ + ∂γwf (φ)/∂φ, and γwf (φ) = −

√

2 3

cos θssurf sin( π2 φ); slip length ls =

1+φ 2

+ λls 1−φ 2 , is an

interpolation between two different wall-fluid slip length, and λls = lsG /lsL . Periodic boundary conditions are assumed at the left and right boundaries: φ, v, p are periodic in x direction

(2.7)

√ Eight dimensionless parameters (except G) appear in the above equations. They are (1) Ld = 3M γ/2 2V L2 , √ (2) R = ρL V L/ηL , which is the liquid phase Reynolds number, (3) B = 3γ/2 2ηL V , which is inversely √ proportional to the capillary number Ca = ηL V /γ, (4) Vs = 3γΓL/2 2V , (5) Ls = ηL /βL L, which is the slip length of liquid scaled by L, (6) G = L/V 2 , (7) λls , λρ , λη are slip length, density, viscosity ratios between gas and liquid. We now show that an energy law can be derived from the above phase field model neglecting the external force. Unlike in the equal density case [5], convective term in Navier-Stokes equation cannot be eliminated √ in the derivation of energy law. Guermond and Quatapelle [20] introduced a new variable σ = ρ in place of ρ . Using mass conservation law ρt + ∇ · (ρv) = 0,

(2.8)

one can derive 1 1 σ(σv)t = ρvt + ρt v = ρvt − ∇ · (ρv)v. 2 2

(2.9) 4

Therefore, it is physically consistent to replace ρvt in (2.2) by σ(σv)t + 12 ∇ · (ρv)v, and we can get 1 R[σ(σv)t + ρ(v · ∇)v + ∇ · (ρv)v] = −∇p + ∇ · [ηD(v)] + Bμ∇φ. 2 The main advantage of this form is that the following property holds:  1 [ρ(v · ∇)w + ∇ · (ρv)w] · wdx = 0, 2 Ω

(2.10)

(2.11)

if v · n|T,B = 0 and ρ, v, w are periodic in x. We then introduce the following notations for different energy terms:  1 ρ|v|2 dx, Fk (v) = 2 Ω  1 (φ2 − 1)2 Fb (φ) = [ |∇φ|2 + ]dx, 4 Ω 2  γwf (φ)dA, Fwf (φ) = T,B

Ft = RFk (v) + BFb (φ) + BFwf (φ), here Fk (v) is the bulk kinetic energy, Fb (φ) is the bulk free energy, Fwf (φ) is the surface energy at the fluid-solid boundary and Ft is the total free energy. Theorem 2.1. Suppose φ, v, p are the smooth solutions to the system (2.1), (2.3), (2.10), with boundary conditions (2.4)-(2.7), then the following energy law holds: 1 √ d Ft = −BLd ∇μ2 − BVs L(φ)2T,B −  ηD(v)2 dt 2 1/2 slip 2 slip − L−1 vx T,B − L−1 s (η/ls ) s (η/ls vx , uw )T,B .

(2.12)

Here vxslip = vx − uw , uw is the wall speed and  · , (·, ·) represent L2 norm and inner product on domain Ω, or at the top and bottom boundaries (T , B). Proof. Multiplying equation (2.1) by Bμ = B(−∇2 φ−φ/+φ3 /), dotting equation (2.10) by v, and adding them together, we have, after integrating on the domain Ω, ∂φ 1 ∂ 1 ( Rσv2 ) + R(ρ(v · ∇)v + ∇ · (ρv)v, v) + B( , −∇2 φ − φ/ + φ3 /) + B(v · ∇φ, μ) ∂t 2 2 ∂t ∂  1 ∂φ ∂ 2 2 2 =R Fk (v) + ( B∇φ + Bφ − 1 ) − B( , ∇φ · n)T,B + B(v · ∇φ, μ) ∂t ∂t 2 4 ∂t ∂ ∂φ = [RFk (v) + BFb (φ)] + B(v · ∇φ, μ) − B( , ∂n φ)T,B . ∂t ∂t

LHS =

Note that we have used boundary conditions (2.5), (2.7), and the equality (2.11). Then by using boundary

5

condition (2.4) and GNBC (2.6), we have ∂ ∂φ ∂φ [RFk (v) + BFb (φ)] + B(v · ∇φ, μ) − B( , L(φ))T,B + B( , ∂γwf /∂φ)T,B ∂t ∂t ∂t ∂ = [RFk (v) + BFb (φ) + BFwf (φ)] + B(v · ∇φ, μ) − B(−vx ∂x φ − Vs L(φ), L(φ))T,B ∂t ∂ slip = Ft + B(v · ∇φ, μ) + BVs L(φ)2T,B + (L−1 + η∂n vx , vx )T,B s η/ls vx ∂t ∂ slip slip = Ft + B(v · ∇φ, μ) + BVs L(φ)2T,B + (L−1 + uw )T,B + (η∂n vx , vx )T,B s η/ls vx , vx ∂t ∂ 1/2 slip 2 = Ft + B(v · ∇φ, μ) + BVs L(φ)2T,B + L−1 vx T,B s (η/ls ) ∂t

LHS =

slip + L−1 s (η/ls vx , uw )T,B + (η∂n vx , vx )T,B .

For the right hand side, we have RHS =(−∇p + ∇ · (η(∇v + ∇vT )), v) + B(μ∇φ, v) + B(Ld ∇2 μ, μ) =(p, ∇ · v) + (∂j (η(∂i vj + ∂j vi )), vi ) + B(v · ∇φ, μ) − BLd ∇μ2 + BLd (μ, ∂n μ)T,B = − (η(∂i vj + ∂j vi ), ∂j vi ) + (nj η(∂i vj + ∂j vi ), vi )T,B + B(v · ∇φ, μ) − BLd ∇μ2 1 = − (η(∂i vj + ∂j vi ), ∂j vi + ∂i vj ) + (nj η(∂x vj + ∂j vx ), vx )T,B + B(v · ∇φ, μ) − BLd ∇μ2 2 1 √ = −  η(∇v + ∇vT )2 + (η(∂x vn + ∂n vx ), vx )T,B + B(v · ∇φ, μ) − BLd ∇μ2 2 1 √ = −  ηD(v)2 + (η∂n vx , vx )T,B + B(v · ∇φ, μ) − BLd ∇μ2 . 2 Thus, comparing LHS and RHS gives (2.12). Remark 2.1. Notice that the first four terms on the right hand side in (2.12) are exactly the four dissipations mentioned in [4], and the last term is the work done per unit time by the flow to the wall. Thus the theorem shows exactly the energy conservation of the system.

3. Numerical schemes and discrete energy law 3.1. Schemes based on pressure correction projection method As in [5], our goal is to design numerical schemes satisfying certain discrete energy law. The method in [5] can be extended here as follows: ⎧ ⎨ ⎩

φn+1 −φn + (˜ vn+1 · ∇)φn+1 δt n+1 n φ −φ + v˜xn+1 ∂x φn+1 = δt

= Ld ∇2 μn+1 , −Vs [L(φ

n+1

x ∈ Ω,

)], ∂n μ

n+1

= 0,

x ∈ T, B.

(3.1)

Here μn+1 = −∇2 φn+1 +(sφn+1 −(1+s)φn +(φn )3 )/ or μn+1 = −∇2 φn+1 −φn /+(φn+1 )3 /; L(φn+1 ) = ˜ n+1 − φn ) and γwf (φ) = − ∂n φn+1 + ∂γwf (φn )/∂φ + α(φ

√

2 3

cos θssurf sin( π2 φ). Notice that a convex splitting

for the boundary term L(φn+1 ) is also needed in order to have the energy decaying property. 6

(ρn+1 , η n+1 , lsn+1 ) =

 1 + φn+1 1 − φn+1 + (λρ , λη , λls ) , σ n+1 = ρn+1 . 2 2

⎧ n+1 n+1 n n ⎪ ⎪ vn+1 + 12 (∇ · (ρn vn ))˜ vn+1 ], R[σ n+1 σ v˜ δt −σ v + ρn (vn · ∇)˜ ⎪ ⎨ vn+1 )] + Bμn+1 ∇φn+1 , = ∇ · [η n+1 D(˜ ⎪ ⎪ ⎪ ⎩ v˜n+1 = 0, [L ln+1 ]−1 (˜ v slip )n+1 = BL(φn+1 )∂ φn+1 /η n+1 − ∂ v˜n+1 , n

s s

x

x

⎧ ⎨ Rρn+1 ( vn+1 −˜vn+1 ) + ∇pn+1 = 0, ∇ · vn+1 = 0, δt ⎩ vn+1 · n = 0,

n x

(3.2)

x ∈ Ω,

(3.3)

x ∈ T, B.

x ∈ Ω,

(3.4)

x ∈ T, B.

In the following theorem, we give the discrete energy law for the total free energy Ft = RFk (v)+BFb (φ)+ BFwf (φ), defined as in Section 3, which also provides the stability of the above scheme. Theorem 3.1. Let vn , φn , pn , n = 0, · · · , k, satisfy equations (3.1)-(3.4) and assume N = | Also if linear splitting is applied,

sup x∈Ω, n=0,1,··· ,k

√

2 π 2 3 (2)

cos θssurf |.

|φn (x)| ≤ M is needed. Then if s ≥ (3M 2 − 1)/2 (only necessary

in linear splitting case), α ˜ ≥ N/2, we have n+1 n+1 slip n+1 /ls (˜ vx ) , uw )T,B Ftn+1 − Ftn + δtL−1 s (η

 1 ≤ − δt( η n+1 D(˜ vn+1 )2 +BLd ∇μn+1 2 +BVs L(φn+1 )2T,B +L−1  η n+1 /lsn+1 (˜ vxslip )n+1 2T,B ), s 2 (3.5)

for n = 0, 1, · · · , k − 1. Moreover, taking summation of (3.5) gives the following bound Ftk + δt

k   1 √ 1 n n v slip )n 2 (  η n D(˜ vn )2 + BLd ∇μn 2 + BVs L(φn )2T,B + L−1 s  η /ls (˜ x T,B ) ≤ C, 2 2 n=1

for any δt and k. Here C is a constant depending only on initial conditions and T = kδt. From the theorem, we see that the total energy Ft decreases when the wall speed uw = 0. Proof. Using the idea of convex spitting in [5], we can decompose bulk free energy Fb (φ) and surface energy Fwf (φ) into contractive and expansive parts: Fb (φ) = Fbc (φ) − Fbe (φ),    s 2 s 1 c 2 e Fb (φ) = [ (∇φ) + φ ]dx, Fb (φ) = [ φ2 − (φ2 − 1)2 ]dx, 2 2 2 4 Ω Ω c e (φ) − Fwf (φ), Fwf (φ) = Fwf   1 2 1 c e Fwf α ˜ φ2 dA, Fwf ˜ φ − γwf (φ)]dA, (φ) = (φ) = [ α T,B 2 T,B 2

7

where s ≥ (3M 2 − 1)/2, α ˜ ≥ N/2, N = |

√

2 π 2 3 (2)

cos θssurf | as in [5]. Then

Fb (φn+1 ) − Fb (φn ) ≤(μn+1 , φn+1 − φn ) + (∂n φn+1 , φn+1 − φn )T,B ,

(3.6)

Fwf (φn+1 ) − Fwf (φn ) ≤(∂γwf (φn )/∂φ + α ˜ (φn+1 − φn ), φn+1 − φn )T,B =(L(φn+1 ) − ∂n φn+1 , φn+1 − φn )T,B .

(3.7)

Adding above two inequalities (3.6), (3.7) and using boundary conditions (3.1), (3.3), we have Fb (φn+1 )+Fwf (φn+1 ) − Fb (φn ) − Fwf (φn ) ≤(μn+1 , φn+1 − φn ) + (L(φn+1 ), φn+1 − φn )T,B =(μn+1 , φn+1 − φn ) − δt(L(φn+1 ), v˜xn+1 ∂x φn+1 + Vs L(φn+1 ))T,B =(μn+1 , φn+1 − φn ) − δtVs L(φn+1 )2T,B n+1 n+1 slip n+1 − δt/B(L−1 /ls (˜ vx ) + η n+1 ∂n v˜xn+1 , v˜xn+1 )T,B s η

n+1 /ln+1 (˜ =(μn+1 , φn+1 − φn ) − δtVs L(φn+1 )2T,B − δt/BL−1 vxslip )n+1 2T,B s s  η n+1 n+1 slip n+1 − δt/BL−1 /ls (˜ vx ) , uw )T,B − δt/B(η n+1 ∂n v˜xn+1 , v˜xn+1 )T,B . s (η

Multiplying equation (3.1) by Bδtμn+1 and using the above inequality, we have B[Fb (φn+1 )+Fwf (φn+1 ) − Fb (φn ) − Fwf (φn )]

 η n+1 /lsn+1 (˜ vxslip )n+1 2T,B ≤ − δtBLd ∇μn+1 2 − δtBVs L(φn+1 )2T,B − δtL−1 s n+1 n+1 slip n+1 − δtL−1 /ls (˜ vx ) , uw )T,B − δt(η n+1 ∂n v˜xn+1 , v˜xn+1 )T,B s (η

− δtB(˜ vn+1 · ∇φn+1 , μn+1 ).

(3.8)

˜ n+1 δt, we obtain Multiplying equation (3.3) by v 1 ˜ n+1 2 − σ n vn 2 + σ n+1 v ˜ n+1 − σ n vn 2 ) R(σ n+1 v 2 1  ˜ n+1 ). vn+1 )2 + δt(η n+1 ∂n v˜xn+1 , v˜xn+1 )T,B + δt(Bμn+1 ∇φn+1 , v = − δt η n+1 D(˜ 2 Dotting equation (3.4) by δtvn+1 gives ˜ n+1 )2 = 0. ˜ n+1 2 + σ n+1 (vn+1 − v σ n+1 vn+1 2 − σ n+1 v Combing above two equations, we have 1  1 R(σ n+1 vn+1 2 − σ n vn 2 ) ≤ − δt η n+1 D(˜ vn+1 )2 + δt(η n+1 ∂n v˜xn+1 , v˜xn+1 )T,B 2 2 ˜ n+1 ) + δtB(μn+1 ∇φn+1 , v

8

(3.9)

Adding (3.8) and (3.9) together gives inequality (3.5). 1  Ftn+1 − Ftn + δt(  η n+1 D(˜ vn+1 )2 + BLd ∇μn+1 2 + BVs L(φn+1 )2T,B 2

n+1 n+1 slip n+1 n+1 /ln+1 (˜ + L−1 vxslip )n+1 2T,B ) ≤ −δtL−1 /ls (˜ vx ) , uw )T,B s s  η s (η



1 1 n+1 /ln+1 (˜ n+1 /ln+1 u 2 ≤ δtL−1 vxslip )n+1 2T,B + δtL−1 s s w T,B . s  η s  η 2 2

Thus, 1  Ftn+1 + δt(  η n+1 D(˜ vn+1 )2 + BLd ∇μn+1 2 + BVs L(φn+1 )2T,B 2



1 1 n+1 /ln+1 (˜ n+1 /ln+1 u 2 + L−1 vxslip )n+1 2T,B ) ≤ Ftn + δtL−1 s s w T,B . s  η s  η 2 2 Taking summation of the above inequality from n = 0 to n = k − 1 gives Ftk + δt

k  1 √ (  η n D(˜ vn )2 + BLd ∇μn 2 + BVs L(φn )2T,B 2 n=1 k    1 1 0 −1 n n v slip )n 2 + L−1  η n /lsn uw 2T,B ≤ C1 + C2 T ≡ C, s  η /ls (˜ x T,B ) ≤ Ft + δtLs 2 2 n=1

where slip length ls and viscosity η are bounded by properties of boundaries and flows and wall speed uw is also assumed to be bounded. Remark 3.1. The inequality (3.5) shows that the total energy decreases in time if wall speed is zero and indicates that the numerical scheme inherits discrete energy decreasing property from the PDE system, and thus is gradient stable. Remark 3.2. This scheme requires to solve an elliptic equation as follows to update pn+1 in equation (3.4):

∇·(

1

ρ

∇pn+1 ) = n+1

˜ n+1 R∇ · v , δt

(3.10)

which is very time-consuming, especially when density ratio is large. 3.2. Schemes based on pressure stabilization As mentioned in Remark 3.2, the above scheme is computationally expensive when density ratio is large. To improve the efficiency, Guermond and Salgado [19] proposed a splitting method for Navier-Stokes equations involving only pressure Poisson equation with constant coefficient, which is based on pressure stabilization [17], namely, the divergence free condition is replaced by ∇ · v − τ ∇2 pt = 0, where τ is a small parameter. Based on these ideas, we propose the following pressure stabilization scheme: 9

Given are initial conditions φ0 , p−1 = p0 = 0, v0 , and set ρ¯ = min(1, λρ ). We compute (φn+1 , vn+1 , pn+1 ) for n ≥ 0 by ⎧ ⎨ φn+1 −φn + vn+1 · ∇φn+1 = L ∇2 μn+1 , d δt ⎩ φn+1 −φn + v n+1 ∂ φn+1 = −V [L(φn+1 )], ∂ lμn+1 = 0, x s n x δt

x ∈ Ω,

(3.11)

x ∈ T, B,

where μn+1 = −∇2 φn+1 + (sφn+1 − (1 + s)φn + (φn )3 )/ or μn+1 = −∇2 φn+1 − φn / + (φn+1 )3 /; ˜ (φn+1 − φn ) and γwf (φ) = − L(φn+1 ) = ∂n φn+1 + ∂γwf (φn )/∂φ + α

(ρn+1 , η n+1 , lsn+1 ) =

√

2 3

cos θssurf sin( π2 φ).

1 + φn+1 1 − φn+1 + (λρ , λη , λls ) . 2 2

(3.12)

⎧ 1 n+1 n n+1 n n

(ρ +ρ )v −ρ v ⎪ ⎪ R 2 + ρn+1 (vn · ∇)vn+1 + 12 (∇ · (ρn+1 vn ))vn+1 ⎪ δt ⎨ = ∇ · [η n+1 D(vn+1 )] − ∇(2pn − pn−1 ) + Bμn+1 ∇φn+1 , ⎪ ⎪ ⎪ ⎩ v˜n+1 = 0, (L ln+1 )−1 (˜ vxslip )n+1 = BL(φn+1 )∂x φn+1 /η n+1 − ∂n v˜xn+1 , s s n ⎧ ⎨ ∇2 (pn+1 − pn ) = ⎩ ∂ pn+1 = 0,

ρ¯ δt R∇

· vn+1 ,

x ∈ Ω,

(3.13)

x ∈ T, B.

x ∈ Ω,

(3.14)

x ∈ T, B.

n

For this scheme, we can prove the following stability result. Theorem 3.2. Let vn , φn , pn , n = 0, · · · , k, satisfy equations (3.11)-(3.14) and assume N = | Also if linear splitting is applied,

sup x∈Ω, n=0,1,··· ,k

√

2 π 2 3 (2)

cos θssurf |.

|φn (x)| ≤ M is needed. Then if s ≥ (3M 2 − 1)/2 (only necessary

in linear splitting case), α ˜ ≥ N/2, we have n+1 n+1 slip n+1 Ftn+1 − Ftn + δtL−1 /ls (vx ) , uw )T,B + s (η

δt2 (∇(pn − pn−1 )2 + ∇pn+1 2 − ∇pn 2 ) 2¯ ρR

1  ≤ −δt(  η n+1 D(vn+1 )2 + BLd ∇μn+1 2 + BVs L(φn+1 )2T,B 2

n+1 /ln+1 (v slip )n+1 2 + L−1 s s  η x T,B )

(3.15)

Proof. Multiplying equation (3.13) by vn+1 δt, we obtain 1 R( (ρn+1 + ρn )vn+1 − ρn vn , vn+1 ) 2 1 = R(σ n+1 vn+1 2 + (σ n (vn+1 − 2vn ), σ n vn+1 )) 2 1 = R(σ n+1 vn+1 2 − σ n vn 2 + σ n (vn+1 − vn )2 ) 2 1  = − δt η n+1 D(vn+1 )2 − δt(∇(2pn − pn−1 ), vn+1 ) 2 + δt(η n+1 ∂n vxn+1 , vxn+1 )T,B + δt(Bμn+1 ∇φn+1 , vn+1 ). 10

(3.16)

2

Taking inner product of equation (3.14) with − δtρ¯ (pn+1 − 2pn + pn−1 ), we get δt2 (∇(pn+1 − pn ), ∇(pn+1 − 2pn + pn−1 )) ρ¯ δt2 (∇(pn+1 − pn )2 − ∇(pn − pn−1 )2 + ∇(pn+1 − 2pn + pn−1 )2 = 2¯ ρ = −δtR(∇ · vn+1 , pn+1 − 2pn + pn−1 ) = δtR(vn+1 , ∇(pn+1 − 2pn + pn−1 )). Taking inner product of equation (3.14) with −

(3.17) δt2 n+1 , ρ¯ p

we have

δt2 (∇pn+1 2 − ∇pn 2 + ∇(pn+1 − pn )2 ) = −δtR(vn+1 , ∇pn+1 ). 2¯ ρ

(3.18)

Taking the difference of equation (3.14) at time step n + 1 and n, we obtain ∇2 (pn+1 − 2pn + pn−1 ) =

ρ¯ R∇ · (vn+1 − vn ) δt

then, δt2 ∇(pn+1 − 2pn + pn−1 )2 ≤ ρ¯R2 vn+1 − vn 2 ≤ R2 σ n (vn+1 − vn )2 , ρ¯

(3.19)

where we have used ρ¯ = min ρ. Combing above three equations together, we derive δt(vn+1 , ∇(−2pn + pn−1 ) ≤

δt2 (−∇(pn − pn−1 )2 − ∇pn+1 2 + ∇pn 2 ) 2¯ ρR 1 + Rσ n (vn+1 − vn )2 , 2

(3.20)

˜ ), we have Combing above two equations (3.16), (3.20), and (3.8) in theorem 3.1 (v instead of v δt2 (∇(pn − pn−1 )2 + ∇pn+1 2 − ∇pn 2 ) 2¯ ρR

 η n+1 /lsn+1 (vxslip )n+1 2T,B ≤ − δtBLd ∇μn+1 2 − δtBVs L(φn+1 )2T,B − δtL−1 s 1  n+1 n+1 slip n+1 − δtL−1 /ls (vx ) , uw )T,B − δt η n+1 D(vn+1 )2 . s (η 2

Ftn+1 −Ftn +

(3.21)

3.3. The second order in time scheme Both the pressure correction method and the pressure stabilization method above can be improved to second order in time. Although the stability and energy law for the second order scheme are not easy to prove analytically numerical results have shown that it is quite stable with good energy decaying property.

11

The pressure correction method can be extended to the following two step incremental pressure correction method [20] ⎧ ⎨ ⎩

3φn+1 −4φn +φn−1 + (˜ vn+1 · ∇)φ∗,n+1 = Ld ∇2 μ ˜n+1 , 2δt n+1 n φ −φ ˜ n+1 )], ∂n μ + v˜xn+1 ∂x φ∗,n+1 = −Vs [L(φ ˜n+1 δt

x ∈ Ω, = 0,

(3.22)

x ∈ T, B.

˜ n+1 ) = ∂n φn+1 +∂γwf (φ∗,n+1 )/∂φ+ where μ ˜n+1 = −∇2 φn+1 +[sφn+1 −(1+s)φ∗,n+1 +(φ3 )∗,n+1 ]/, and L(φ α ˜ (φn+1 − φ∗,n+1 ) and γwf (φ) = −

(ρn+1 , η n+1 , lsn+1 ) =

√

2 3

cos θssurf sin( π2 φ), (·)∗,n+1 := 2(·)n − (·)n−1 .

1 + φn+1 1 − φn+1 + (λρ , λη , λls ) . 2 2

⎧ n+1 n +vn−1 ⎪ ⎪ R[ρn+1 3˜v −4v + ρn+1 (v∗,n+1 · ∇)˜ vn+1 ) + 12 (∇ · (ρn+1 v∗,n+1 ))˜ vn+1 ] ⎪ 2δt ⎨ vn+1 )], = −∇pn + Bμn+1 ∇φn+1 + ∇ · [η n+1 D(˜ ⎪ ⎪ ⎪ ⎩ v˜n+1 = 0, [L ln+1 ]−1 (˜ vxslip )n+1 = BL(φn+1 )∂x φn+1 /η n+1 − ∂n v˜xn+1 , s s z

(3.23)

x ∈ Ω,

(3.24)

x ∈ T, B,

where μn+1 = −∇2 φn+1 − φn+1 / + (φn+1 )3 /. ⎧ 1 ⎪ ⎪ ∇ · ( ρn+1 ∇ψ n+1 ) = ⎪ ⎨

3 2δt R∇

˜ n+1 , ·v

∂n ψ n+1 = 0, ⎪ ⎪ ⎪ ⎩ pn+1 = pn + ψ n+1 − η n+1 ∇ · v ˜ n+1 ,

x ∈ Ω, x ∈ T, B.

(3.25)

x ∈ Ω.

For the pressure stabilization scheme, the second order in time is achieved as following ⎧ ⎨ 3φn+1 −4φn +φn−1 + (vn+1 · ∇)φ∗,n+1 = L ∇2 μ ˜n+1 , x ∈ Ω, d 2δt n+1 n+1 ⎩ φn+1 −φn + v n+1 ∂ φ∗,n+1 = −V [L(φ ˜ )], ∂ μ ˜ = 0, x ∈ T, B, δt

x

x

s

(3.26)

n

˜ n+1 ) = ∂n φn+1 + ∂γwf (φ∗,n+1 )/∂φ + where μ ˜n+1 = −∇2 φn+1 + [sφn+1 − (1 + s)φ∗,n+1 + (φ3 )∗,n+1 ]/, L(φ α ˜ (φn+1 − φ∗,n+1 ), γwf (φ) = −

(ρn+1 , η n+1 , lsn+1 ) =

√

2 3

cos θssurf sin( π2 φ), and (·)∗,n+1 := 2(·)n − (·)n−1 .

1 + φn+1 1 − φn+1 + (λρ , λη , λls ) . 2 2

⎧ n+1 n +vn−1 ⎪ ⎪ R[ρn+1 3v −4v + ρn+1 (v∗,n+1 · ∇)vn+1 )] + ∇(pn + 43 ψ n − 13 ψ n−1 ) ⎪ 2δt ⎨ = Bμn+1 ∇φn+1 + ∇ · [η n+1 D(vn+1 )], ⎪ ⎪ ⎪ ⎩ (L ln+1 )−1 (v slip )n+1 = B L(φ ˜ n+1 )∂x φn+1 /η n+1 − ∂n vxn+1 , vnn+1 = 0, s s x 12

(3.27)

x ∈ Ω, x ∈ T, B,

where μn+1 = −∇2 φn+1 − φn+1 / + (φn+1 )3 /. ⎧ 3ρ¯ ⎪ ⎪ R∇ · vn+1 , ∇2 ψ n+1 = 2δt ⎪ ⎨ ∂n ψ n+1 = 0, ⎪ ⎪ ⎪ ⎩ pn+1 = pn + ψ n+1 − η n+1 ∇ · vn+1 ,

x ∈ Ω, x ∈ T, B,

(3.28)

x ∈ Ω.

3.4. Numerical implementation To implement the above numerical schemes, we use the finite difference method to discretize space, and vn+1 + 12 (∇ · treat the nonlinear terms explicitly, i.e, (˜ vn+1 · ∇)φn+1 as (vn · ∇)φn in (3.1), ρn (vn · ∇)˜ (ρn vn ))˜ vn+1 as ρn (vn · ∇)vn + 12 (∇ · (ρn vn ))vn in (3.3), which will decouple CH equation (3.1) from NS equations (3.3). However, if the density ratio is large, we need to pay special attention to the discretization of the convection terms. We will use third order weighted essentially non-oscillatory (WENO) schemes [25] to reduce the undershoot and overshoot around the interface. Numerical truncation of φ is also needed. Here we redistribute the phase variable φ as follows [32]: 1. Truncate the undershoot/overshoot values: ⎧ ⎨ 1.0, φ > 1.0 + 10−7 , φˆ = ⎩ −1.0, φ < −(1.0 + 10−7 ). 2. Compute the total mass M , M0 and the difference G:  φdx, M0 = Ω ˆ φdx, M= Ω

G = M0 − M. 3. Uniformly distribute G at NG points as G/NG , where NG is the total number of cells that fall into the interface region (−0.999 < φ < 0.999). By using this redistribution, the undershoot/overshoot are eliminated and the total mass is conserved. To solve (3.1), we decompose the original problem into two equations − ∇2 φn+1 + sφn+1 / − μn+1 = RHS1

(3.29)

φn+1 − δtLd ∇2 μn+1 = RHS2

(3.30)

here RHS1 = [(1 + s)φn − (φn )3 ]/, RHS2 = φn − δtF (vn , φn ), where F (vn , φn ) represents WENO scheme for convection term. In this form, the boundary conditions (3.1) can be easily incorporated. If periodic boundary conditions are applied in one direction (2D) or two directions (3D) where fast Fourier transform 13

are applicable, the resulting linear systems (3.29), (3.30) are band matrices which can be solved by direct elimination. Otherwise, they can be solved by iterative methods, e.g., generalized minimal residual method (GMRES) with preconditioner using the periodic boundary conditions. Once φn+1 is solved, it can be inserted into NS equations to solve velocity and pressure. When solving velocity, it is better to decouple different components of velocity through treating the first term of the right hand side of the following equality implicitly and second term explicitly: ∇ · (ηD(v)) = η∇2 v + ∇η · D(v). This will greatly remove stiffness of viscous term and the second term introduce some mild CFL constraint. The resulting system has only variable coefficients in diagonal positions, which is solved efficiently by preconditioned conjugate gradient method.

4. Numerical results In this section, we present several numerical examples to verify the accuracy, stability and the energy decay property of the above numerical methods and also to show the efficiency and capability of our method. In the simulations, the second order semi-implicit schemes: incremental pressure correction scheme (IPCS) (3.22)-(3.25) and pressure stabilization scheme (PSS) (3.26)-(3.28) will be used. 4.1. Two phase Couette flow In the first test, we simulate the Couette flow with different density, viscosity and slip length as in Fig. 2, where the top and bottom boundaries are moving with uw = 1.0 × (1.0 − e−t/0.1 ) in opposite directions. GNBC is used for velocity and relaxation boundary condition is used for phase variable at top and bottom boundaries while periodic boundary conditions are used at left and right boundaries. The parameters used are: Ld = 5.0 × 10−4 , R = 10, B = 40, Vs = 500, ls = 0.0038,  = 0.01, θssurf = 77.6◦ , s = 1.5, α = 0.125,

(4.31)

and viscosity, density and slip length ratios are λη = 0.2, λρ = 0.1, λls = 10. The domain is (x, z) ∈ [0, 1.2] × [0, 0.4], and grids are (nx, nz) = (3N, N ), N = 40, 80, 160, 320. The time step is δt = 0.1h and final time is T = 1.0. In this test, third order WENO is applied but phase variable redistribution is not used due to low density ratio. Accuracy check The error in L2 norm is calculated by 3N/2,N/2

EN (·) = [



[(·)N (2i, 2j) − (·)N/2 (i, j)]2 dxN/2 dzN/2 ]1/2 ,

i=1,j=1

14

Figure 2: Initial condition of φ for the two phase Couette flow.

Figure 3: Density plot of (φ) at T = 0.05, 0.15, 1.0.

and the convergence rate is calculated by log2 (EN /E2N ). The results calculated by both PSS (Table 1) and IPCS (Table 2) show second order accuracy for velocity and phase variable. However, IPCS is computationally much more expensive than PSS due to the cost of solving elliptic equation (3.10). The slip velocity along lower boundary is plotted in Fig. 4, which shows large slipping around the contact line. Slip is nonsymmetric due to different slip length at two sides. Profiles of interfaces at various time are plotted in Fig. 3 and the interface reaches steady state at time T = 1.0. N

EN (φ)

rate

EN (u)

rate

80

3.19e-002

160

2.47e-003

3.69

5.83e-004

3.66

5.32e-004

3.55

3.74e-001

3.20

320

5.96e-004

2.05

1.34e-004

2.12

1.32e-004

2.01

8.61e-002

2.12

640

1.49e-004

2.00

3.20e-005

2.07

3.33e-005

1.99

2.20e-002

1.97

7.38e-003

EN (v)

rate

6.23e-003

EN (p)

rate

3.44e+000

Table 1: L2 norm of the error and convergence rate for phase function φ, velocity v, pressure p at T = 1.0 with different grids. The results are obtained with the pressure stabilization scheme for shear flow with static contact angle θssurf = 77.6◦ . The results show second order accuracy for all quantities φ, v, p. .

Discrete energy law We now verify the discrete energy law of our scheme proved in theorem 3.2. The physical configuration is as in Fig. 2 and parameters are listed in (4.31) with larger B = 40. Two cases are simulated: one is without shearing on the flow and the other is with shearing on the flow with shear velocity uw = 1.0. We first look at the case when the wall speed is zero. The behavior of the various energy terms is shown in Fig 15

1 0.8

v

s

0.6 0.4 0.2 0 0

0.5

1 x

Figure 4: Slip velocity ux along lower boundary at T = 1.0.

N

EN (φ)

rate

EN (u)

rate

80

7.01e-02

160

3.81e-03

4.20

2.36e-03

4.12

3.71e-03

4.47

1.72e-02

3.59

320

1.09e-03

1.80

6.41e-04

1.88

1.12e-03

1.73

4.15e-03

2.05

640

3.01e-04

1.86

1.69e-04

1.92

3.17e-04

1.82

1.15e-03

1.85

4.11e-02

EN (v)

rate

8.20e-02

EN (p)

rate

2.08e-01

Table 2: L2 norm of the error and convergence rate for phase function φ, velocity v, pressure p at T = 1.0 with different grids. The results are obtained with the incremental pressure correction scheme for shear flow with static contact angle θssurf = 77.6◦ . The results show second order accuracy for all quantities φ, v, p.

16

5, where Fb (φ), Fwf (φ), Fk (v), Ft are those defined in section 2 multiplied by B, B, R, 1, respectively, and n  k k slip k (Ft + W )n = Ftn + δtL−1 s (η /ls (vx ) , uw )T,B is the total free energy at time t = nδt plus work done k=0

to the wall. As indicated by inequality (3.15), when the wall speed is zero, the total energy decreases in time (here pressure gradient term is neglected because it is lower order in time). This is shown in Fig. 5(a). Notice that surface energy at the fluid-solid boundary decays very fast at the beginning, and is transferred to other energies and dissipations. The bulk free energy then decreases slowly while the surface energy increases as the interface reaches steady state with the prescribed static contact angle (see the density plot of φ at time step n = 50 and n = 200). When the wall speed is not zero (uw = 1.0), the total free energy Ft may increase, but Ft + W should decrease as indicated by (3.15). These are shown in Fig. 5(b). The density plots of φ are also shown for time step n = 50 and n = 200, which shows steady state structure of the interface.

4.2. Droplet spreading on a flat boundary To study the contact line dynamics, we simulate the droplet spreading on a flat surface and also study the effect of the interface thickness . Again, the parameters used are as follows: Ld = 5.0 × 10−4 , R = 1000, B = 12, Vs = 500, ls = 0.038,  = 0.01, θssurf = 50◦ , s = 1.5, α = 0.374,

(4.32)

and viscosity, density and slip length ratios are λη = 0.1, λρ = 0.001, λls = 1. The domain is (x, z) ∈ [0, 4.0] × [0, 1.0]. The initial condition is a spherical drop centered at (2.0, 0.1) with radius 0.4:  √ φ = − tanh(( (x − 2.0)2 + (y − 0.1)2 − 0.4)/( 2)).

(4.33)

The interface profiles and contact point positions at different time between t = 0 and t = 10 are shown in Fig.7. The droplet initially spreads quickly due to flow inertial and the contact point reaches a maximum distance, then it starts to retract due to surface tension. The process then repeats and contact point oscillate around x = 2.703 and eventually reaches the steady state position with a contact angle equal to the static contact angle θssurf = 50◦ . We then study the effect of the interface thickness with decreasing  = 0.02, 0.01, 0.005, 0.0025, and corresponding grids (nx, nz) = (4N, N ), N = 50, 100, 200, 400. The time step is δt = 0.1h. The droplet interface positions at the same time t = 0.5 but for different  are shown in Fig. 6. The results show that the interface thickness  does not seem to affect much the dynamics and structure of the interface. 4.3. Droplet spreading on a patterned surface To demonstrate the efficiency of our method, we also carry out a three dimensional simulation of a droplet spreading on a chemically patterned surface. As shown in Fig. 8, two intersecting hydrophobic strips of width 0.1 with contact angle θssurf = 2π/3 are placed at the center of otherwise hydrophilic 17

(a) Without shear

(b) With shear

Figure 5: Four energy terms as functions of time with and without shear applied on the top and bottom boundaries. The shear velocity is uw = 1.0 and contour plots of phase variable φ are embedded in the figure to show interface structure at time step n = 50, 200. In the figure, the values are adjusted to zero at n = 0 by subtracting their initial values.

18

1 ε=0.02

0.495 0.8

ε=0.01 0.49

ε=0.005

0.6 y

1.99

ε=0.0025

2.01

x 10

0.4

0.2

2

−3

2 0

0 0

1.5 0.5

1.51 1

1.5

2 x

2.5

3

3.5

4

Figure 6: Interface positions at time T = 0.5 for four different  = 0.02, 0.01, 0.005, 0.0025. The embedded figures show the detailed structure around moving contact line and at the center of the droplet.

0.5

3

t=0.0

0.4

2.8

t=2.0

x

z

0.3 2.6

0.2 2.4

t=5.0

0.1 t=10.0

0 2

2.2

2.4

2.6

2.8

2.2 0

3

x

10

20

30

40

50

t

(a) Interface profile

(b) Position of contact point

Figure 7: Interface position at different times from t = 0 to t = 10.0 with time interval 1.0 (left) and the contact point position as a function of time (right). The equilibrium contact point position is x = 2.703. The result is calculated for  = 0.01.

19

background with θssurf = π/4. We use the same parameters as in (4.31). The computational domain is (x, y, z) ∈ [0, 0.8] × [0, 0.8] × [0, 0.6] with grid (nx, ny, nz) = (160, 160, 80) and the time δt = 0.05h. Initially, a half sphere with radius r = 0.2 is located at the center and the shape of droplet at time t = 0.0, 0.04, 0.12, 0.40 is shown in Fig. 8 . The contact line positions for a sequence of time from t = 0.0 to t = 0.16 are shown in Fig. 9. As expected, the results show that the droplet quickly spreads outwards on hydrophilic surface and contracts inwards on hydrophobic strips and eventually reaches an equilibrium state.

(a) t=0.00

(b) t=0.04

(c) t=0.12

(d) t=0.40

Figure 8: droplet surface at time t = 0.0, 0.04, 0.12, 0.40.

5. Concluding remarks In this paper, a splitting method based on a pressure Poisson equation [19] is proposed for the moving contact line problems with larger density and viscosity ratio. The method requires only one elliptic solver with constant coefficient per time step. We show energy decaying and unconditional stability property for the method under certain conditions. The linearized schemes are easy to implement and introduce only mild CFL time constraint. Numerical tests are carried out to verify the accuracy, stability and efficiency of the schemes. The behavior of the solution near the contact line is examined and three dimensional simulations are included to demonstrate the efficiency of the numerical schemes. 20

0.8

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

Figure 9: Positions of contact line between time t = 0 to t = 0.16 with time interval 0.02.

Acknowledgments This publication was based on work supported in part by Award No. SA-C0040/UK-C0016, made by King Abdullah University of Science and Technology (KAUST), the Hong Kong RGC-GRF Grants 605311, 605513 and NNSF of China grant 91230102

References [1] T. Qian, X.P. Wang, P. Sheng, Molecular scale contact line hydrodynamics of immiscible flows, Phys. Rev. E 68 (2003) 016306. [2] X.P. Wang, T. Qian, P. Sheng, Moving contact line on chemically patterned surfaces, J. Fluid Mech. 605 (2008) 59-78. [3] T. Qian, X.P. Wang, P. Sheng,

Molecular hydrodynamics of the moving contact line in two-phase immiscible flows,

Commun. Comput. Phys. 1 (2006) 1-52. [4] T. Qian, X.P. Wang, P. Sheng, A variational approach to the moving contact line hydrodynamics, J. Fluid Mech. 564 (2006) 333-306. [5] M. Gao, X.P. Wang, A gradient stable scheme for a phase field model for the moving contact line problem, J. Comput. Phys. 231 (2012) 1372-1386. [6] Q. He, R. Glowinski, X.P. Wang, A least square/finite element method for the numerical solution of the Navier-StokesCahn-Hilliard system modeling the motion of the contact line, J. Comput. Phys. 230 (2011) 4991-5009. [7] J. Shen, X. Yang and Q. Wang, Mass and volume conservation in phase field models for binary fluids, Commun. Comput. Phys. 13 (2013) 1045-1065. [8] D.J. Eyre, An unconditionally stable one-step scheme for gradient systems, Unpublished article, June 1998. [9] D. Anderson, G.B. McFadden, A.A. Wheeler, Diffuse interface methods in fluid mechanics, Ann. Rev. Fluid Mech. 30 (1998) 139. [10] J. Shen, X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations, DCDS-A 28 (4) (2010) 1669-1691. [11] H. Johnston, J.G. Liu, Accureate, stable and efficient Navier-Stokes solvers based on explicit treatment of the pressure term, J. Comput. Phys. 199 (2004) 221-259.

21

[12] J.L. Guermond, P. Minev, J. Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Eng. 195 (2006) 6011-6045. [13] D.L. Brown, R. Cortez, M.L. Minion, Accurate projection methods for the incompressible Navier-Stokes equations, J. Comput. Phys. 168 (2) (2001) 464-499. [14] H. Ding, P.D.M. Spelt, C. Shu, Diffuse interface model for incompressible two-phase flows with large density ratios, J. Comput. Phys. 266 (2007) 2078-2095. [15] H. Ding, P.D.M. Spelt, Onset of motion of a three-dimensional droplet on a wall in shear flow at moderate Reynolds numbers, J. Fluid Mech. 599 (2008) 341-362. [16] P. Yue, C. Zhou, J. J. Feng, Sharp interface limit of the Cahn-Hilliard model for moving contact lines, J. Fluid Mech. 645 (2010) 279-294. [17] J. Shen, X. Yang, A phase-field model and its numerical approximation for two-phase incompressible flows with different densities and viscosities, SIAM J. Sci. Comput. 32 (3) (2010) 1159-1179. [18] J. Pyo and J. Shen, Gauge-Uzawa methods for incompressible flows with variable density, J. Comput. Phys. 221 (2007) 181-197. [19] J.-L. Guermond and A. Salgado, A splitting method for incompressible flows with variable density based on a pressure Poisson equation, J. Comput. Phys. 228 (2009) 2834-2846. [20] J.-L. Guermond and L. Quartapelle, A projection FEM for variable density incompressible flows, J. Comput. Phys. 165 (2000) 167-188. [21] W. E and J.-G. Liu, Gauge method for viscous incompressible flows, Commun. Math. Sci. 1 (2003) 317-332. [22] C.M. Pooley, Olga Kuksenok, Anna C. Balazs, Convection-driven pattern formation in phase-separating binary fluids, Phys. Rev. E 71 (2005) 030501. [23] L.K. Antanovskii, A phase field model of capillary, Phys. Fluids 7 (4) (1995) 747. [24] J. Lowengrub, I. Truskinovsky,

Cahn-Hilliard fluids and topological transitions, Proc. R. Soc. London A 454 (1998)

2617-2654. [25] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200-212. [26] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes II, J. Comput. Phys. 83 (1989) 32-78. [27] S. Wise, J. Kim, J. Lowengrub,

Solving the regularized, strongly anisotropic Cahn-Hilliard equation by an adaptive

nonlinear multigrid method, J. Comput. Phys. 226 (1) (2007) 414-446. [28] S.M. Wise, J.S. Lowengrub, V. Cristini, An adaptive multigrid algorithm for simulating solid tumor growth using mixture models, Mathematical and Computer Modelling, (2010). [29] G. Tryggvason, Numerical simulations of the Rayleigh-Taylor instability, J. Comput. Phys. 75 (1988) 253-282. [30] B.J. Daly, Numerical study of two fluid Rayleigh-Taylor instabilities, Phys. Fluids 10 (1967) 297. [31] X.Y. He, S.Y. Chen, R.Y. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, J. Comput. Phys. 152 (1999) 642-663. [32] P.H. Chiu, Ya.T. Lin, A conservative phase field method for solving incompressible two-phase flows, J. Comput. Phys. 230 (2011) 185-204. [33] T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for incompressible two-phase flows with large density differences, J. Comput. Phys. 198 (2) (2004) 628-644. [34] Y.Y. Yan, Y.Q. Zu, A lattice Boltzmann method for incompressible two-phase flows on partial wetting surface with large density ratio, J. Comput. Phys. 227 (1) (2007) 763-775. [35] Y. Tanaka, Y. Washio, M. Yoshino, T. Hirata, Numerical simulation of dynamic behavior of droplet on solid surface by the two-phase lattice Boltzmann method, Computers & Fluids 40 (1) (2011) 68-78.

22