An augmented method for free boundary problems with moving contact lines

An augmented method for free boundary problems with moving contact lines

Zhilin Li

a,*

, Ming-Chih Lai

b

, Guowei He

c

, Hongkai Zhao

d

a

Center for Research in Scientific Computation & Department of Mathematics, North Carolina State University, Raleigh, NC 27695-8205, United States b

Center of Mathematical Modeling and Scientific Computing & Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 30050, Taiwan c

LNM, Institute of Mechanics, Chinese Academy of Sciences, Beijing 10080, China d

Department of Mathematics, University of California, Irvine, CA 92697, United States

a r t i c l e

i n f o

Article history: Received 11 April 2009

Received in revised form 10 January 2010 Accepted 18 January 2010

Available online 25 January 2010 Keywords:

Moving contact line Free boundary problem Triple junction One-phase flow Navier-Stokes equations Embedding technique Immersed interface method Irregular domain Augmented method

a b s t r a c t

An augmented immersed interface method (IIM) is proposed for simulating one-phase moving contact line problems in which a liquid drop spreads or recoils on a solid substrate. While the present two-dimensional mathematical model is a free boundary problem, in our new numerical method, the fluid domain enclosed by the free boundary is embedded into a rectangular one so that the problem can be solved by a regular Cartesian grid method. We introduce an augmented variable along the free boundary so that the stress balancing boundary condition is satisfied. A hybrid time discretization is used in the projection method for better stability. The resultant Helmholtz/Poisson equations with interfaces then are solved by the IIM in an efficient way. Several numerical tests including an accuracy check, and the spreading and recoiling processes of a liquid drop are presented in detail.

Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction

In this paper, we propose an augmented immersed interface method for free boundary problems with moving contact lines in two space dimensions. This problem arises from the spreading of liquid drop on a solid substrate. We wish to investigate the wetting effects of the drop on the solid surface [9]. The contact line is defined as the intersection of the fluid interface with the solid sur-face. When the contact line is in its equilibrium state (not moving), the contact angle, the angle between the fluid interface and the solid boundary, is called the static contact angle that depends on the surface tension coefficients between the different phases (Laplace–Young relation). In the presence of moving contact lines, it is well-known that the no-slip boundary condition at the moving contact line will lead to a non-integrable force singularity[12,6]. Several different techniques including continuum hydrodynamical or molecular dynamical treatments have been proposed to remove such singularity, see for example, a recent review article by Qian et al.[25]. In this paper, our intention is not to provide a new math-ematical model or analysis, but is to develop an efficient numerical

method for the mathematical model developed by Ren and Weinan

[26]for the moving contact line problem.

The one-phase free boundary problem is a simplified model for a two-phase (liquid–gas) problem with neglecting the influence of one phase (gas), seeFig. 1for an illustration. In this paper, we par-ticularly consider a liquid drop on a solid surface and study its wet-ting behavior. As usual, the governing equations are described by the incompressible Navier-Stokes equations in a time dependent single fluid domainXðtÞ, together with the stress force balancing conditions on the free surfaceCFðtÞ, and the slip boundary condi-tions on the solid surfaceCSðtÞ, that is,

q

@u @t þ ðu 

r

Þu   þ

r

p ¼

l

D

u þ G; in

X

ðtÞ; ð1:1Þ

r

 u ¼ 0; in

X

ðtÞ; ð1:2Þ  p þ nT_

_l

_ð

_r

_{u þ}

_r

_uT_{Þ  n ¼}

_cj

_{ p} air; on

C

FðtÞ; ð1:3Þ sT

l

ð

r

u þ

r

uT Þ  n ¼ 0; on

C

FðtÞ; ð1:4Þ

v

¼ 0; bu ¼

l

@u @y; on

C

SðtÞ: ð1:5Þ

Here, u ¼ ðu;

v

Þ is the fluid velocity, p the pressure, and G is the external force that can include the gravity. The constants

q

and

l

are the fluid density and viscosity, respectively.

0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.01.013

* Corresponding author.

E-mail address:[email protected](Z. Li).

Contents lists available atScienceDirect

Computers & Fluids

Eqs.(1.3) and (1.4)are the stress force balancing conditions on the free boundaryCFðtÞ, where

c

is the surface tension,

j

the cur-vature of the free boundary (negative sign inFig. 1), p_airthe pres-sure outside of the drop, and n and s are the unit normal and tangent vectors at the free boundary. The derivation for above stress balance conditions can be found, for example, in[14]. In additional to the stress force conditions, we should impose the kinematic condition on the free boundary, that is, the normal velocity of the fluid at the free boundary should be equal to the normal velocity of the interface. In this paper, we shall use the level set method[21]to evolve the motion of the free boundary so that the kinematic boundary condition is imposed naturally. The detail for the free boundary representation by level set should be dis-cussed in later section.

Along the fluid–solid boundaryCSðtÞ, as mentioned before, the no-slip boundary condition at the moving contact lines leads to a non-integrable force singularity. Therefore, we allow the solid boundary to be partial slip and impose the Navier slip boundary condition, Eq.(1.5). The first condition in(1.5)is the no penetration boundary condition. The second condition is the slip with friction indicating that the tangential velocity of the fluid is proportional to their tangential stress. Here, as in[26], b is the friction coeffi-cient. Both conditions in Eq.(1.5)can also be easily derived from their general expression as shown in[8]. Notice that, in other liter-ature such as in[8,23], one can define k ¼

l

=bas the slip length indicating the fictitious distance to the solid surface where the fluid tangential velocity is extrapolated to be zero.

In additional to imposing the Navier slip boundary condition on the solid wall to remove force singularity, the contact line dynam-ics must also be prescribed. In the presence of contact angle hyster-esis, one can prescribe the contact angle depending on the sign of contact line speed[31,22,8,32]. In the above models, however, the advancing and receding contact angles must be given. Recently, Ren and Weinan derived an effective boundary condition at the contact line from the force balance argument[26]. 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 their effective condition to the contact line directly. That is, at the mov-ing contact lines x1ðtÞ and x2ðtÞ, we have

a

u ¼

c

ðcos h cos hÞ; ð1:6Þ

where

a

is the effective friction coefficient, h _{the static} (equilib-rium) contact angle that depends on the surface tensions of the free and solid boundaries, and h is the dynamic contact angle between the moving free boundary and the solid boundary as shown in

Fig. 1.

While a two-phase model, in which the density and viscosity are discontinuous across the interface, is more practical, a direct discretization of the Navier-Stokes equations may be ill-condi-tioned and difficult to solve. The condition number of the discrete system is proportional to the jumps of the physical parameters.

Under certain circumstances, a one-phase model is adequate and can be solved more efficiently than the two-phase model. More discussions about the circumstances will be discussed later.

The solvability and stability of the problem under certain

cir-cumstances are discussed in [10,24,30,29] including simplified

one-dimensional numerical simulations. Most numerical

simula-tions are for the two-phase model. In [17], a volume of fluid

(VOF) method was developed. In [8], an Arbitrary Lagrangian

Eulerian (ALE) finite element method using a body fitted mesh was studied. The Immersed Boundary (IB) method have been proposed in[11,26], in which the boundary condition at the con-tact lines are enforced through an external point force at those lines. These point-sources are then distributed to the nearby grid points via a discrete delta function. In [13], a simplified model using the Laplace equation to obtain the contact angle is studied and a level set method was developed for the numerical simulations.

In this paper, we propose an augmented immersed interface method for the one-phase free boundary problem based on Carte-sian grids. Thus, there is almost no cost in the grid generation even though the free boundary is moving. We first embed the domain into a rectangular one with the solid surface being the bottom side so that we can solve the problem on a rectangular domain using a Cartesian grid. We then use a modified version of the projection method[3,15,16]to solve the Navier-Stokes equations for the mo-tion of the fluid. They key step is to use the jump in the normal derivative of the velocity as an augmented variable so that we can solve the fluid equations efficiently. We will explain our numerical method in detail in the next section; followed by results of numerical experiments and analysis. We draw some conclusions in the last section.

2. The numerical method

Our numerical method for the fluid equations is based on the projection method for solving the incompressible Navier-Stokes equations. There are several versions of the projection method for solving the incompressible Navier-Stokes equations, see for example,[3,15,16]and many others. The projection method that we used in our scheme is the one described in[5]which is based on the pressure increment formulation of[3,15].

As mentioned before, instead of solving the fluid equations in an irregular time dependent domainXðtÞ, we embed the fluid domain

of the droplet into a rectangular computational domain

R ¼ ½a; b  ½c; d which is sufficiently large to enclose the time dependent droplet with the bottom side coinciding with the solid surface. Since the computational domain is a regular rectangle in-stead of an irregular one bounded by the free boundary, the gov-erning equations now can be solved in a Cartesian grid. The apparent advantage of using Cartesian grid is that there are several fast direct solvers for fluid equations available in literature. Here,

the spatial spacings are chosen as hx¼ ðb  aÞ=M and

hy¼ ðd  cÞ=N, where M and N are the number of grid points used in the x and y directions, respectively. For simplicity, we use a stan-dard uniform mesh such that h ¼ hx¼ hy.

As a common practice, from one time level to the next, we use a splitting approach for the free boundary problem. In this approach, we first fix the moving boundary (denoted byCkFbelow) and solve the fluid equations on the irregular domain to get the velocity. We then use the computed velocity to evolve the free boundary using the level set method. We describe the numerical scheme from time tk _{to t}kþ1 _{below. At time step t}k_{, given an approximation of the} velocity uk_{, the pressure p}k_{, the free boundary position C}k

F, and an initial guess of augmented variable qkþ1_{, we carry out the} fol-lowing steps: contact angle n Γ_F(t) Ω(t) Γ_S(t) θ₁ θ₂ τ X₁(t) X₂(t)

Fig. 1. A diagram of a one-phase flow with two moving contact lines X1ðtÞ and X2ðtÞ, and two contact angles h1ðtÞ and h2ðtÞ. The circular portion represents the free boundary CFðtÞ where the stress force balancing conditions are imposed. The bottom line represents the solid boundaryCSðtÞ where the Navier slip boundary conditions are imposed.

Step 1: Prediction step

q

u uk Dt ¼ rpk_

_q

_{ðu ruÞ}kþ1 2_þ

l

2ðDu _þDuk_{Þ þ G}kþ1 2_;x 2 ð_{X\ Z}R_Þ rpk_

_q

_{ðu ruÞ}k þ

l

Du_{þ G}k_{;x 2 ð} X\ ZIÞ

l

2ðDu _þDuk_{Þ; x 2 ðX}c \ ZRÞ

l

Du_{;x 2 ðX}c \ ZIÞ 8 > > > > > > < > > > > > > : ð2:7Þ ½u_ Ck F¼ 0; @u @n   Ck F ¼ qkþ1_{; ð2:8Þ} @u @n¼ 0 at x ¼ a;x ¼ b; and y ¼ d;

v

_{¼ 0 and bu}_¼

_l

@u @y at y ¼ c: ð2:9Þ

Step 2: Projection step

D

/kþ1¼

r

 u 

D

t ; x 2 R; ½/kþ1Ck F¼ 0; @/kþ1 @n " # Ck F ¼ 0; @/ kþ1 @n      @R ¼ 0; 8 > > > > < > > > > : ð2:10Þ ukþ1_{¼ u}_

_D

_t

_r

/kþ1; x 2 R; ð2:11Þ

r

pkþ1_¼

_r

_pk_þ

_r

_/kþ1_{; x 2}

_X

_ð2:12Þ

Step 3: Computing the residual of the boundary condition(1.3) and (1.4). Ekþ1 1 ¼ pkþ1þ nT

l

ðrukþ1þrðukþ1Þ T Þ  n  ð

cj

 pairÞ; onCkF ð2:13Þ Ekþ1 2 ¼ sT

l

ðrukþ1þrðukþ1Þ T Þ  n ¼ 0; on CkF: ð2:14Þ

The first two steps consist of the traditional steps of the projec-tion method for the Navier-Stokes solver. The augmented variable qkþ1_{defined only along the free boundaryC}k

Fis introduced in the prediction step and should be determined suitably such that the free boundary conditions(1.3) and (1.4)are satisfied at time level tkþ1_{. If the boundary conditions are satisfied, that is, E}kþ1

1 ¼ 0 and

Ekþ1

2 ¼ 0, then we have computed the velocity ukþ1, otherwise, we need to find better qkþ1 _{until the residual kE}kþ1

1 k and kE

kþ1

2 k are

smaller enough. This can be done iteratively, or by solving a linear system of equations. This is explained in detail below and in Sec-tion2.3. Notice that, the normal vector n, the tangent vector s and the curvature

j

are all evaluated using the free boundary posi-tion at the time level tk_.

In the algorithm above, X\ ZR and Xc\ ZR are regular grid points inside and outside of the domain X, respectively, while

X\ ZIandXc\ ZIare irregular grid points in the corresponding do-mains. The classification of the regular or irregular grid points can be found in Section2.2. The non-linear term ðu ruÞkþ12_is

approx-imated by ðu 

r

uÞkþ12_¼3 2ðu k 

r

Þuk1 2ðu k1_

_r

_Þuk1_{; ð2:15Þ}

at regular grid points. We use a first order discretization for ðu ruÞkat irregular grid points as well to avoid cross differentia-tion for the explicit terms.

The pressure term in the prediction step is treated explicitly and we are only interested in the solution in the domainX, so we do not need to computerp outside the domain. As we can see, the resultant equation in the prediction step inXis an elliptic equa-tion, so it is natural to embed the domainXinto a computational rectangular domain R by setting another elliptic equation (heat equation here) inXc¼ R X. In such way, we can discretize the whole equations and form a linear system based on Cartesian grids.

Given the augmented variable q, we can see that in Step 1 and 2, there involve solving elliptic interface problem with singular sources on a Cartesian domain R with different boundary conditions. The detail on how to solve those equations using immersed interface method can be found in Section2.4. To solve the unknown qkþ1_{, we} can choose q ¼ eito get coefficient matrix for qkþ1, see Section2.3. In other words, to get the coefficient matrix, we need to run the algo-rithm Nbtimes at one time level, where Nbis the dimension of Qkþ1, the discrete approximation of qkþ1_{. Alternatively, we can use the} GMRES iterative method that only requires the matrix–vector mul-tiplication. Each matrix–vector multiplication requires to solve the Navier-Stokes equation once, and to interpolate the boundary condi-tions once. More details will be followed later.

In this paper, we use the backward Euler scheme near/on the free boundary while use the Crank-Nicholson scheme away from the boundary in the prediction step. The reason is for the stability consideration. The scheme is at least second order accurate in space and first order in time. When Dt  Oðhx;hyÞ, the scheme seems to be still second order accurate since we can use one order lower scheme along the boundary, see for example,[2,4,19]. The stability reason is to avoid using the explicit Laplacian term Duk for the moving boundary. Even if uk_{is second order accurate, the} approximation of the explicit termDuk _{may be inaccurate near/} on the boundary because the errors are often not smooth. Such a treatment is very useful if the boundary is moving. The disadvan-tage is that a fast Helmholtz solver based on FFT cannot be applied directly. Here, we use the structured multi-grid solver DMGD9V

[7]. Note that, it is possible to use a second-order time-stepping scheme such as the BDF scheme. The difficulty is the extra book-keeping effort needed to deal with the irregular grid points in three different time levels rather than the two levels.

Note that a fully implicit approach (backward Euler)

q

u _{ u}k

D

t ¼ 

r

pk_

_q

_{ðu }

_r

_uÞk þ

l

D

u_{þ G}kþ1_{; x 2 ð}

_X

_Þ

l

D

u_{; x 2 ðR n}

_X

_Þ ( ð2:16Þ

is simpler and may be appropriate if we are only interested in the equilibrium solution. The method is first order in time and second order in space and has better stability. Another advantage is that, fast Helmholtz/Poisson solvers based FFT can be applied. This ap-proach has also been tested and worked well.

2.1. The implementation of the contact line condition

As discussed in[26], on the solid boundary and away from the contact line, the Navier slip boundary condition(1.5)is imposed. However, on the solid boundary but at the contact line, the effective condition Eq. (1.6) must be imposed. There are several ways to implement the moving contact line condition (1.6) numerically. The simplest approach is to enforce the condition smoothly by using

uðxÞ ¼

c

a

ðcos h



 cos hÞeðxxiÞ2=w2_{; ð2:17Þ}

where xiis the position of the contact line on the solid boundary, and w  OðhxÞ is the effective width of the transition region, see [26]. In summary, we have Navier slip boundary condition at almost all grid points along the solid boundary except for the two neigh-boring grid points next to the contact line where the condition

(2.17)is imposed. Notice that, the above numerical boundary con-dition is implemented for u_{on the solid boundary in the projection} step in which the contact angle h and the contact line position xiare both obtained explicitly (that is, from previous time step).

The second approach is to add an external point force fidðx  xiÞdðyÞ at the contact line, say ðxi;0Þ, i ¼ 1; 2, for a pre-deter-mined force strength fiarising from the unbalanced Young’s force,

then use a discrete delta function to distribute the force term to the nearby grid points of the contact lines as in the Peskin’s Immersed Boundary method, see[11,26].

2.2. The level set representation and the orthogonal projections We use the level set method[23,28]to evolve the free boundary for convenience since some of the codes used for the results in[13]

can be handily applied. We believe that the front tracking method would work equally well in two space dimensions. One difficulty is to deal with the boundary condition of the free boundary at the contact lines. In the level set method, the free boundary is repre-sented by the zero level set of a two-dimensional Lipschitz contin-uous function

u

ðx; y; tÞ, for example, the signed distance function. The level set function then is evolved according to the Hamilton– Jacobi equation

u

tþ u 

r

u

¼ 0: ð2:18Þ

In discretization, the level set function is defined at grid points. One advantage of the level set representation is that we can clas-sify grid points easily. In reference to the standard central five-point stencil, a grid five-point xij¼ ðxi;yjÞ is regular if

u

minij

u

maxij >0, where

u

max

ij ¼ max

u

i1;j;

u

ij;

u

iþ1;j;

u

i;j1;

u

i;jþ1

n o

; ð2:19Þ

u

min

ij ¼ min

u

i1;j;

u

ij;

u

iþ1;j;

u

i;j1;

u

i;jþ1

n o

: ð2:20Þ

Otherwise it is called an irregular grid point which means that the free boundary cuts through the finite difference stencil.

In discretization, the augmented variable is defined at the orthogonal projection of the irregular grid points from the domain

Xþ, that is the outside of the fluid domain. Let be xijbe such a grid point, the orthogonal projection can be approximated by

x

ij¼ xijþ

a

r

u

; ð2:21Þ

where

a

is an approximate signed distance between xijand the free boundary, which is obtained by solving the quadratic equation,

u

ðxÞ þ ð

r

u

ðxÞ  pÞ

a

þ1 2ðp T_Heð

_u

_{ðxÞÞ pÞ}

_a

2_{¼ 0; ð2:22Þ} where p ¼r

u

and pT_Heð

_u

Þ p ¼

u

2 x

u

xxþ 2

u

x

u

y

u

xyþ

u

2 y

u

yy; ð2:23Þ

in Cartesian coordinates. The partial derivativesr

u

ðxÞ, the Hessian matrix Heð

u

Þ are computed at the grid point xij. The computed projec-tions have third order of accuracy ifr

u

ðxÞ and Heð

u

Þ are computed using the standard centered five-point finite difference formulas.

The contact angle between the interface and the x-axis can be easily obtained from n ¼r

u

=jr

u

j. It is well-known that we should keep the level set function as a good approximation to the signed distance function through a re-initialization process. A simple way is to solve the following Hamilton–Jacobi equation

u

tþ sgnð

u

Þ jð

r

u

j  1Þ ¼ 0: ð2:24Þ

through an artificial time t, where sgnð

u

Þ is the sign function of

u

. Both Eqs.(2.18) and (2.24) are solved by the third order WENO scheme.

Note that since the free boundary cuts the bottom of the bound-ary, a simple treatment of the boundary for the level set function as used for closed interface cannot be used here. When we solve the level set function(2.18), we apply the boundary condition

v

¼ 0 and thus the level set equation on y ¼ c becomes a one-dimen-sional Hamilton–Jacobi equation that is solved by the WENO scheme. For the re-initialization process along y ¼ c, after we up-dated the level set function for the interior points, we use the

WENO scheme again in the x-direction, and approximate the par-tial derivative in y-direction using

@

u

@yð:; 0Þ ¼

u

:;2

u

:;1

hy

ð2:25Þ

in which

u

:;2and

u

:;1have already been updated as they are defined at interior grid points. InFig. 2, we show a contour plot of the level set function after the proposed re-initialization process. The con-tour lines have almost the same distance, or jr

u

j  1 as we desired. 2.3. The augmented method and the Schur complement system

Since it is sufficient to consider the solution from one time level tk _{level to the next t}kþ1_{, we will ignore the time index for} conve-nience if there is no confusion. Below we describe how to get the linear system of equations for the augmented variable Q. Let U be the vector whose components are Uijand Pij, the approximate solution to the problem at one particular time step. Thus U’s dimension is Oð3MNÞ, where M and N are the number of grid lines in the x and y direction respectively. Let Q be the vector of the dis-crete values of q at the orthogonal projections of the irregular grid

points from the Xþ _{side. Thus Q’s dimension is OðNÞ assuming}

M  N. Then the discrete solution of(2.9)–(2.11) and (2.12)given Q can be written as

AU þ BQ ¼ F1 ð2:26Þ

for some vector F1 and sparse matrices A and B. It requires solving three Helmholtz1_{/Poisson equations with different source terms and}

jump conditions to get U.

Once we know the solution U given Q, we can interpolate U and P to getrU_{and P}_{to approximate the one-phase boundary} con-ditions at those points where the discrete values Q are defined. We denoterU_{and P}_{as the limiting values ofrU and P from the}

X

side. The interpolation scheme depends on U; Q linearly. There-fore we can write

EU þ TQ  F2¼ 0; ð2:27Þ

where E and T are two sparse matrices, and F2is a vector. We need to choose such a vector Q that the free boundary conditions(1.3) and (1.4)are satisfied along the free boundary @X. If we put the two matrix–vector Eqs.(2.26) and (2.27)together we get

A B E T   U Q   ¼ F1 F2   : ð2:28Þ

Note that Q is defined only on a set of points fXlg on the free bound-ary while U is defined at grid points. The Schur complement for Q is

−10 −0.5 0 0.5 1

0.2 0.4 0.6 0.8

Fig. 2. A snap shot of a free boundary and the level set function after the proposed re-initialization process. The solid red line in the middle is the free boundary. We can see that the level set function is a good approximation of the signed distance function after the re-initialization. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

1

They are actually generalized Helmholtz equations Du  k2_{u ¼ f . For simplicity,} we simply call them the Helmholtz equations in this paper if there is no confusion.

ðT  EA1BÞQ ¼ F2 EA1F1¼ F: ð2:29Þ

If we can solve the system above to get Q, then we can get U easily. Because the dimension of Q is much smaller than that of U, we expect to get a reasonably fast algorithm if we can solve

(2.29)efficiently. We refer the readers to[19]for other details of an augmented immersed interface method.

Since we know how to get the matrix–vector multiplication, we can get the Schur complement matrix by taking Q ¼ ei, where eiis the ith base vector. Then we can use the Gaussian elimination meth-od to solve the Schur complement system. Or alternatively, we can use the GMRES iterative method[27]to solve the Schur complement system. The GMRES method only requires the matrix–vector multi-plication. Since the free boundary is moving, the Schur complement matrix is changing with time, in most of our tests, it is more efficient to use the GMRES iterative method than that of forming the matrix, then applying the LU decomposition, seeFig. 4b. In our tests, we do not apply preconditioning techniques because we do not form the matrix. It is still an interesting numerical issue on how to develop efficient preconditioning techniques with only the knowledge of the matrix–vector multiplication. More details about the augmented method can be found in[19,21].

2.3.1. The least squares interpolation for the free boundary condition One crucial step of the algorithm is to interpolate the free boundary conditions(1.3) and (1.4)given an approximate solution Uij and Pij. At an orthogonal projection xij corresponding to an irregular grid point xij, we need to interpolate Uijto get @u=@n and @u_=@

_s

_{. This is computed using a least squares interpolation,} for example, @u @n ¼ X ns1 k¼0

c

kUiþik;jþjk Cij ð2:30Þ

where ikand jkare integers taken from 0; 1; 2; . . . ; nsis the num-ber of neighboring grid points involved, often taken between 9 and 16; Cijis a correction term to offset the jump in the normal derivative of u. The coefficients

c

_kare chosen so that the interpola-tion scheme is second order accurate. The linear system of equa-tions

c

k is obtained by expanding uðxijÞ at xijand then matching up to second order partial derivative terms. If uðxijÞ involved is from outside ofX, then it is replaced by the extended value

uðxijÞ ¼ uþþ uþ_nnþ uþg

g

þ12 u þ nnn 2 þ 2uþ ngn

g

þ uþgg

g

2   þ Oðh3Þ ¼ u_{þ u} nnþ ug

g

þ12 unnn 2 þ 2u ngn

g

þ ugg

g

2  

þ½unn þ ½ug

g

þ12 ½unnn2þ 2½ungn

g

þ ½ugg

g

2þ Oðh3Þ

where ðn;

g

Þ is the local coordinate system in the normal and tan-gential directions, and all the values are defined at x

ij. Since we know the jump conditions ½u and ½un ¼ ½un, we can obtain all other jump conditions, see[19]. Note that we need to take nsP6 to get a consistent linear system. Thus the linear system of equations for the coefficients

c

kis under-determined, which can be solved by the sin-gular value decomposition (SVD).

2.4. Solving an elliptic interface problem with singular sources Our proposed algorithm consists of solving two generalized Helmholtz equations for the intermediate velocity u_{in(2.7), and} one Poisson equation for the pressure increment /kþ1_in(2.10)with a given jump condition in the normal derivative of the velocity. The details for solving Helmholtz/Poisson equations with jump condi-tions in the solution and its normal derivative (or the elliptic PDEs with singular sources) can be found in[18–20]. Here, we just give a brief sketch on how we solve such problems in an efficient and accurate way.

Without loss of generality, we consider the following general-ized Helmholtz equation

wxxþ wyy kw ¼ f ; ðx; yÞ 2 R; ½w@X¼ 0; @w @n   @X ¼ q: ð2:31Þ

In our application, we have k ¼ 2=ð

l

DtÞ for the prediction step, and k ¼ 0 for /kþ1_{in(2.10). Since k is a constant, a fast Poisson} sol-ver can be used.

The finite difference discretization using the immersed inter-face method can be simply written as

w_iþ1;jþ wi1;j 2wi;j

h2x

þwi;jþ1þ wi;j1 2wi;j h2y

 kwij¼ f ðxi;yjÞ þ Cij;

ð2:32Þ

where the correction term Cijis zero at regular grid points where the boundary @Xdoes not cut through the standard centered 5-point stencil. The correction term Cijat those irregular grid points can be determined in a dimension by dimension fashion, see[20]

for the formula of Cij.

We use the structured multi-grid solver DMGD9V[7]to solve the discrete system for u_{. While using the fast Poisson solver from} [1]for /. This is because of the hybrid discretization for u _that changes the coefficient in the Helmholtz equations for u_{, and} dif-ferent boundary conditions for u at the bottom side.

3. Numerical experiments

In this section, we present some numerical results for the free boundary problem with moving contact lines. All the computations were performed at the North Carolina State University using either notebook or desktop computers. Most simulations are done within minutes to a couple of hours depending on the mesh, surface ten-sion, the initial geometry, and the static contact angle. We use a le-vel set method to evolve the free boundary for convenience since some of the codes used for the results in[13]can be used. We be-lieve that the front tracking method would work equally well in two space dimensions. In most of simulations, we take the density

q

¼ 1, the viscosity

l

¼ 2, and the surface tension

c

¼ 0:5. The other parameters are the slip coefficient b ¼ 2, the effective friction coefficient

a

¼ 0:5, and the smoothing length w ¼ h, unless stated otherwise. The time step size is take as

D

tk_{¼ min} h 2; h 2Uk max ( ) ; ð3:33Þ

where h ¼ maxfhx;hyg and Ukmax¼ maxijjUkijj unless specified otherwise.

Example 1 (Validation of the method against an exact solution). As a first numerical test for our scheme, we consider an example in a stationary irregular domain in which the exact solution is known analytically. This is an example without the effect of moving contact lines which we simply use it as an accuracy check. The analytic solution is uðx; y; tÞ ¼ sinðtÞ _y r 2y  ; if r P 1=2 sinðtÞ r2_1 4 y; otherwise; 8 > < > : ð3:34Þ

v

ðx; y; tÞ ¼ sinðtÞ  x rþ 2x  ; if r P 1=2  sinðtÞ r2_1 4 x; otherwise; 8 > < > : ð3:35Þ

pðx; y; tÞ ¼ sinðtÞ sin x þ sin yð Þ; ð3:36Þ

where r ¼pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffix2_{þ y}2_{. The domain of the interest is the domain} bounded by r 6 1=2 and y P 0. We extend the solution to the

rect-angular domain R ¼ ½1; 1  ½0; 1 so that we have a consistent initial condition and a problem with the exact jump condition in the normal derivative of the velocity. The source term G is derived directly from the exact solution.

InTable 1, we show the grid refinement analysis to check the order of the accuracy of our method. We use the exact non-homo-geneous normal derivative boundary condition at @R. Since we are interested in the computed solutions in the domainX, we set

kEuk1¼_r max ij61=2; yP0 jUkij uðxi;yj;TÞj n o þ max rij61=2; yP0 jVkij

v

ðxi;yj;TÞj n o kEpk1¼_r max ij61=2; yP0 jPkij pðxi;yj;TÞj n o

to be the error in the velocity at time T. The numbers order are the approximated order of accuracy from the two consecutive errors for the velocity. Second order accuracy is clearly seen for the velocity. The pressure is at least first order accurate. The number cond is the condition number of the coefficient matrix for the unknown jump in the normal derivative of the velocity at the first step. Example 2 (Drop spreading and recoiling). We started with the same example that was used in[26]. The initial drop is a semi-circle centered at the origin with radius r ¼ 0:5, and

l

¼ 2. Thus

the initial contact angle between the boundary and the x-axis is

p

=2. The surface tension is taken as

c

¼ 0:5. InFig. 3a, we show the initial droplet and its shape at t ¼ 3:9060. The static contact angle is h

¼

p

=4 so the drop would spread as we can see from the

simulation. The simulation was computed using a 128 by 64 grid. We use homogeneous Neumann boundary condition on @R except for the bottom side where we use the Navier BC for u and

v

¼ 0. We set u ¼ 0;

v

¼ 0, and p ¼ 0, initially. We keep the boundary fixed and solve the problem to get consistent initial data before letting the free boundary to move. The condition number of the Schur complement matrix is well under 200. InFig. 3b, we show a recoiling case by plotting the initial free boundary, r ¼ 0:4 and y P 0, and its shape at t ¼ 0:9586 without the gravity force. The static angle now is h

¼ 3

p

=4. In this case, the drop contracts at the solid surface. Also for this case, the condition number of the Schur complement matrix is larger than that of the spreading case.

InFig. 4a, we show the time evolution of the contact angle and the contact line speed in the spreading case which the static con-tact angle is

p

=4. The contact angle and line speed do not decrease exactly in a monotonic way to the equilibrium state due the inter-action of the fluid and the surface tension.

InFig. 4b, we show the number of iterations of the GMRES iter-ation at each time level for the spreading case. The mesh size is 128 by 64. The size of the Schur complement system is about 178 by 178 for the augmented variable. The convergence tolerance is 105_{. The number of iterations of the GMRES fluctuates but seldom} exceeds its full dimensions. It is faster than that when we form the matrix of the Schur complement and then use the LU decomposi-tion to solve the system of equadecomposi-tions. However, while the condidecomposi-tion number is small to modest, the number of iterations of the GMRES is not close to a constant and can be quite large if the tolerance is very small. This is because the matrix is far from a normal matrix ðAAT¼ ATAÞ. It is a challenging problem to develop an efficient lin-ear solver for the Schur complement system of equations knowing only the matrix–vector multiplications.

Table 1

A grid refinement analysis against the exact solution at a final time T ¼ 0:5:kEuk1is the sum of the maximal error in the velocity component u and v, order is the approximated convergence order computed from the two consecutive errors, cond is the condition number of the coefficient matrix for the unknown jump in the normal derivative of the velocity at the first step.

M  N kEuk1 orderu kEpk1 orderp cond 32  16 9:4533  103 _{7:2203  10}3 _11.90 64  32 2:2300  103 _{2.0838 1:8997  10}3 _{1.9278 27.42} 128  64 6:0532  104 _{1.8813 6:0351  10}4 _{1.6528 29.66} 256  128 9:7791  105 _{2.6290 2:8304  10}4 _{1.0923 190.0} −10 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 t=0 t=3.9060 −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

(a)

(b)

Fig. 3. (a) A drop spreads along the contact line with h

¼p=4. The initial free boundary is r ¼ 0:5 above the x-axis. The red-line plot is the initial free boundary. (b) A drop recoils along the contact line with h

¼ 3p=4 with t ¼ 0:9586. The initial free boundary is r ¼ 0:4 above the x-axis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Example 3 (A perturbed surface). In Fig. 5a, we show the free boundary at t ¼ 0, and at a time t ¼ 10. The initial boundary was a perturbed semi-circle r ¼ 0:5 þ 0:05 sinð8hÞ; 0 6 h 6

p

. In this case, the surface tension smooths the free boundary while the contact line is moving. The static contact angle is h

¼

p

=4.

Example 4 (The effect of the gravity). InFig. 6, we show the numer-ical experimental results of the simulation with and without grav-ity. The initial free boundary is the half circle r ¼ 0:4 and y P 0 centered at the origin. Thus, the initial angle h ¼

p

=2. We assume that the static angle is h

¼

p

=4. Without the gravity, the motion is slower, see the plot of the green line . With the gravity, the motion is faster and flatter, see the plot of the blue line.

Example 5. This example is adapted from [17]. The initial free boundary is r ¼ 0:3 and y P 0 centered at (0, 0.2). The initial angle h >

p

=2. We assume again that the static angle is h_¼

_p

_{=4.Fig. 7} shows some snap shots of the free boundary at different time. The GMRES method does take more number of iterations.

4. Conclusions

In this paper, we propose an augmented immersed interface method for solving one-phase free boundary problems with moving contact lines. The one-phase model is a simplification of two-phase model under some assumptions. The one-phase model

−0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

π/4

tip speed contact angle 0 100 200 300 400 500 0 20 40 60 80 100 120 140 160 180 200

178

(a)

(b)

Fig. 4. (a) A history of the contact angle and the contact line speed as functions of time. (b) The number of GMRES iterations for the Schur complement system whose size is about 178 by 178 for a 128 by 64 grid.

−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 t=0 t=10

Fig. 5. A drop with perturbation spreads along the contact 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 t=0 t=3.4493, with g t=6.4191, no g

is valid if, the boundary conditions at the contact lines are nearly consistent. The method can track the interactions of the fluid equa-tions and the moving contact lines. By introducing the augmented variable along the boundary, we can solve the problem in a rectan-gular domain. We use a hybrid time discretization with backward Euler at irregular grid points and Crank-Nicholson at regular grid points for the prediction step to have better stability. The linear system of equations is well-conditioned too.

Acknowledgement

We would like to thank Dr. Weinan E (Princeton) for bringing us this interesting problem. We would also like to thank Drs. Weiqing Ren (NYU), Kazufumi Ito and Michael Shearer (NCSU) for valuable discussions.

The author would like to acknowledge the following supports. The first author is partially supported by US-ARO Grants 56349-MA and 49308-56349-MA, US-AFSOR Grant FA9550-09-1-0520, and NSF Grant 0911434. The second author is partially supported by Na-tional Science Council of Taiwan under Grant NSC-97-2628-M-009-007-MY3 and MoE-ATU project.

References

[1] Adams J, Swarztrauber P, Sweet R. Fishpack: efficient Fortran subprograms for the solution of separable elliptic partial differential equations. %3chttp:// www.netlib.org/fishpack/%3e.

[2] Beale JT, Layton AT. On the accuracy of finite difference methods for elliptic problems with interfaces. Math Comput Sci 2006;1:91–119.

[3] Bell JB, Colella P, Glaz HM. A second-order projection method for the incompressible Navier-Stokes equations. J Comput Phys 1989;85:257–83. [4] Beyer RP, LeVeque RJ. Analysis of a one-dimensional model for the immersed

boundary method. SIAM J Numer Anal 1992;29:332–64.

[5] Brown DL, Cortez R, Minion ML. Accurate projection methods for the incompressible Navier-Stokes equations. J Comput Phys 2001;168:464. [6] Dussan EBV, Davis SH. On the motion of a fluid–fluid interface along a solid

surface. J Fluid Mech 1974;65:71.

[7] De Zeeuw D. Matrix-dependent prolongations and restrictions in a blackbox multigrid solver. J Comput Appl Math 1990;33:1–27.

[8] Ganesan S, Tobiska L. Modeling and simulation of moving contact line problems with wetting effects. Comput Visual Sci 2008. doi:10.1007/ s00791008-0111-3.

[9] de Gennes P-G, Brochard-Wyart F, Quere D. Capillarity and wetting phenomena. Springer; 2003.

[10] Haley PJ, Miksis MJ. The effect of the contact line on droplet spreading. J Fluid Mech 1991;223:57–81.

[11] Huang H, Liang D, Wetton B. Computation of a moving drop/bubble on a solid surface using a front-tracking method. Commun Math Sci 2004;2:535–52. [12] Huh C, Scriven LE. Hydrodynamics model of steady movement of a solid/

liquid/fluid contact line. J Colloid Interface Sci 1971;35:85.

[13] Hunter J, Li Z, Zhao H. Autophobic spreading of drops. J Comput Phys 2002;183:335–66.

[14] Joseph DD, Renardy YY. Fundamentals of two-fluid dynamics. Springer-Verlag; 1993.

[15] Van Kan J. A second-order accurate pressure-correction scheme for viscous incompressible flow. SIAM J Sci Comput 1986;7:870–91.

[16] Kim J, Moin P. Application of a fractional-step method to incompressible Navier-Stokes equations. J Comput Phys 1985;59:308–23.

[17] Li J, Renardy M. Numerical simulation of moving contact line problems using a volume-of-fluid method. J Comput Phys 2001;171:243–63.

[18] Li Z. A fast iterative algorithm for elliptic interface problems. SIAM J Numer Anal 1998;35:230–54.

[19] Li Z, Ito K. The immersed interface method – numerical solutions of PDEs involving interfaces and irregular domains. SIAM Front Ser Appl Math 2006;FR33.

[20] Li Z, Lai M-C. The immersed interface method for the Navier-Stokes equations with singular forces. J Comput Phys 2001;171:822–42.

[21] Li Z, Wan X, Ito K, Lubkin SR. An augmented approach for the pressure boundary condition in a Stokes flow. Commun Comput Phys 2006;1:874–85. [22] Liu H, Krishnan S, Marella S, Udaykumar HS. Sharp interface Cartesian grid method II: a technique for simulating droplet interactions with surfaces of arbitrary shape. J Comput Phys 2005;210:32–54.

[23] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. New York: Springer; 2002.

[24] Pismen LM, Rubinstein BY, Bazhlekov I. Spreading of a wetting film under the action of van der waals forces. Phys Fluids 2000;12:480–3.

[25] Qian T, Wang X-P, Sheng P. Molecular hydrodynamics of the moving contact line in two-phase immiscible flows. Commun Comput Phys 2006;1(1):1–52. [26] Ren W, WE. Boundary conditions for the moving contact line problem. Phys

Fluids 2007;19.

[27] Saad Y. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986;7:56–869. [28] Sethian JA. Level set methods and fast marching methods. 2nd ed. Cambridge

University Press; 1999.

[29] Socolowsky J. The solvability of a free boundary problem for the stationary Navier-Stokes equations with a dynamic contact line source. Nonlinear Anal 1993;21:763–84.

[30] Spaid MA, Homsy GM. Stability of Newtonian and viscoelastic dynamic contact lines. Phys Fluids 1996;8:460–78.

[31] Spelt PDM. A level-set approach for simulations of flows with multiple moving contact lines with hysteresis. J Comput Phys 2005;207:389–404.

[32] Wu C, Lei S, Qian T, Wang X. Stick-slip motion of moving contact line on chemically patterned surfaces. Commun Comput Phys 2010;7:403–22.

−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 t=0 t=0.4343 t=1.1525 t=2.2859