Development of an immersed boundary-phase field-lattice Boltzmann method for Neumann boundary condition to study contact line dynamics

Journal of Computational Physics 234 (2013) 8–32

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Development of an immersed boundary-phase field-lattice Boltzmann method for Neumann boundary condition to study contact line dynamics J.Y. Shao, C. Shu ⇑, Y.T. Chew Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 117576, Singapore

a r t i c l e

i n f o

Article history: Received 5 October 2011 Received in revised form 15 June 2012 Accepted 28 August 2012 Available online 8 September 2012 Keywords: Immersed boundary method Phase field method Lattice Boltzmann method Multiphase flow Contact line dynamics Dirichlet condition Neumann condition

a b s t r a c t The implementation of Neumann boundary condition in the framework of immersed boundary method (IBM) is presented in this paper to simulate contact line dynamics using a phase field-lattice Boltzmann method. Immersed boundary method [10] is known as an efficient algorithm for modelling fluid–solid interaction. Abundance of prominent works have been devoted to refine IBM [1,11,12]. However, they are mainly restricted to problems with Dirichlet boundary condition. Research that implements the Neumann boundary condition in IBM is very limited to the best of our knowledge. This deficiency significantly limits the application of IBM in computational fluid dynamics (CFD) since physical phenomena associated with Neumann boundary conditions are extremely diverse. The difficulty is attributed to the fact that implementation of Neumann boundary condition is much more complex than that of Dirichlet boundary condition. In the present work, we initiate the first endeavour to implement Neumann boundary condition in IBM with assistance of its physical interpretation rather than simple mathematical manipulation. Concretely speaking, rooted from physical conservation law, the Neumann boundary condition is considered as contribution of flux from the boundary to its relevant physical parameter in a control volume. Moreover, the link between the flux and its corresponding flow field variable is directly manipulated through the immersed boundary concept. In this way, the Neumann boundary conditions can be implemented in IBM. The developed method is applied together with phase field-lattice Boltzmann method to study contact line dynamics. The phase field method [27,39], which becomes increasingly popular in multiphase flow simulation, can efficiently capture complex interface topology and naturally resolve the contact line singularity. Meanwhile, the lattice Boltzmann method is known as an alternative to model fluid dynamics and holds good prospect to simulate multiphase flows with complex geometry [38]. In this context, the developed immersed boundary-phase field-LBM is verified in detail for both steady and unsteady contact line problems. Tests show that the proposed method can correctly reproduce both equilibrium results and dynamic process. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction In fluid mechanics, problems involving interactions between fluid and immersed objects are ubiquitous. To simulate such problems, implementation of boundary conditions, including the Dirichlet and Neumann boundary conditions, is an indispensable task. A straightforward way is to use a body-conformal grid. Specifically, a grid is first generated concerning solid ⇑ Corresponding author. E-mail address: [email protected] (C. Shu). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.08.040

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

9

shape as shown in Fig. 1. Subsequently the governing equations are naturally discretized considering solid geometrical information and the boundary conditions can also be readily imposed. In this way, an accurate solution can be yielded [1]. Nevertheless, it is usually difficult to generate a body-fitted grid with good quality, especially when a complex geometry or a moving boundary is presented in the flow field. To overcome this drawback, non-body-conformal methods [2–10] which become increasingly popular, could be employed. Among various non-body-conformal methods, the immersed boundary method (IBM) developed by Peskin [10] has attracted much attention from researchers. Prominent advantages of IBM include easy implementation and robustness to handle complex geometries. In IBM, the governing equations are resolved on a fixed Cartesian (Eulerian) grid. The boundary is presented by a set of Lagrangian points which are independent to the Eulerian grid as shown in Fig. 2. Meanwhile, the influence of boundary on flow field is depicted by adding a force or source term into the governing equations. Hitherto, IBM has undergone continuous development in accurate implementation of boundary conditions. In Peskin’s original work, the boundary is treated as being elastic [10]. The relationship between the restoring force and body deformation is governed by the Hook’s law, in which a user-defined coefficient is involved. To avoid using the arbitrary coefficient, the direct forcing method is introduced by Fadlun et al. [11]. In their method, the Navier–Stokes (NS) equations are employed to compute the force on the boundary. More recently, a boundary conditionenforced IBM is developed by Wu and Shu [12]. Unlike the previous force calculation schemes, in which the restoring force is pre-calculated and there is no mechanism to enforce the no-slip boundary condition. The newly developed method considers the restoring force as unknown and it is determined in a way that the no-slip boundary condition is enforced. An elegant review of kindred algorithms is provided by Mittal and Iaccarino [1]. Tremendous effort has been devoted to refine IBM. However, most of them are restricted to the Dirichlet boundary condition over decades. To the best of our knowledge, the research that aims at extending the IBM to depict the Neumann boundary condition is very limited. This remarkably restricts the application of IBM in computational fluid dynamics (CFD) since physical phenomena associated with Neumann boundary conditions are extremely diverse. One of the instances is the contact line dynamics in multiphase flow problems. This means that there is a practical demanding to develop an efficient IBM for Neumann boundary condition. In multiphase flows, when a solid object is brought into contact with two immiscible fluids, a three-phase contact line will be formed. This class of phenomena [13–15] can be referred as contact line dynamics. The contact line dynamics can be numerically modeled by a slip boundary (Dirichlet) condition [16–23] in sharp interface methods or Robin boundary conditions in the phase field method [24–27]. In the phase field method, the no-slip boundary condition is still utilized. Additionally, two Neumann boundary conditions are used to govern variation of composition on a solid boundary [27]. Although the phase field method is relatively a new method, it becomes increasingly popular because of its efficiency to capture interface topology change, natural regularization of the contact line singularity and a sound physical background [24,27–40]. As an interface capturing scheme, the phase field method can be combined with either Navier–Stokes (NS) solver or lattice Boltzmann method (LBM), which is known as an attractive alternative to simulate fluid dynamics. Advantages of LBM [38] include easy implementation and specifically, the potential to bridge widely dispersed scales in multiphase flows. Therefore, in the present work, a phase field-lattice Boltzmann method is used to capture the interface and study the contact line dynamics. Although many works have been done by both phase field-NS solver and phase field-LBM solver, most of the reported results investigate contact line on perfectly smooth surface or, at most, grooved surface with simple geometries represented by straight lines. This is mainly attributed to the intricacies caused by implementation of the Neumann boundary condition on a complex geometry. To overcome this difficulty, an IBM-based algorithm that can handle Neumann boundary conditions will be instrumental. In fact, some attempts have been made, such as the work done by Francois and Shyy [41]. In their work, the immersed boundary method was used to track fluid–fluid interface through interface markers. On

Fig. 1. Sketch of a body-conformal grid.

10

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 2. Sketch of a grid used in immersed boundary method.

the wall, an ‘‘ad hoc’’ strategy [42] was used. The position of the interface marker near the wall was identified by a prescribed contact angle and no-slip boundary condition was imposed. However, only Dirichlet boundary condition was involved and thus an IBM-based algorithm for Neumann boundary conditions still remains to be developed. To sum up, the contact line dynamics involves both Dirichlet and Neumann boundary conditions. The phase field-lattice Boltzmann method is a good choice to simulate this problem. To apply immersed boundary method to solve this problem with geometric complexity, we need to develop an efficient immersed boundary-phase field-lattice Boltzmann method, which can treat both Dirichlet and Neumann boundary conditions. In this paper, we begin to make the first endeavour to extend the application of IBM with Neumann boundary condition in the context of the contact line dynamics. To achieve this, we start from the physical conservation law, and view the Neumann boundary condition as contribution of flux from the boundary to a dependent variable in a control volume. In fact, the flux contribution from the boundary is directly linked to the correction of dependent variable. We emphasize that the way is rooted from the physical conservation law rather than from mathematical manipulation. The developed algorithms will be applied to study the contact line dynamics on straight and curved boundaries. The rest of paper is organized as follows. The governing equations of phase field model and boundary conditions are presented in Section 2. Section 3 briefly introduces a boundary condition-enforced IBM [12] for the no-slip boundary condition. Section 4 is devoted to elaborate the idea to implement Neumann boundary condition in IBM for contact line dynamics in detail. Section 5 summarizes the overall simulation procedures. Section 6 validates and examines the numerical behaviour of the proposed method. First, the ability of the proposed method to simulate both steady and unsteady contact line problems is verified by simulation of transition layer [43] and droplet dewetting process respectively. Moreover, droplet spreading on a smooth surface in the partial wetting regime is also simulated. Furthermore, the robustness of the current algorithm to handle complex geometries is demonstrated by simulation of the contact line dynamics on a single and two alongside circular cylinders. Conclusions are drawn in Section 7 to end this paper. 2. Phase-field model for incompressible immiscible fluids Before we embark on illustration of IBM for contact line dynamics, a multiphase flow model must be established. In this work, a phase field-lattice Boltzmann method is employed. This section first presents the governing equations of phase field method. Subsequently, the wetting boundary conditions in this framework are described. 2.1. The governing equations for flow field In this work, we primarily follow a model developed by Swift et al. [44] for incompressible immiscible fluids. The governing equations of flow field read

@n þ r  ðnuÞ ¼ 0 @t @nu þ r  ðn uuÞ ¼ rP þ lr2 u þ F b @t

ð1Þ ð2Þ

where n is the mean density of two fluids. It must be stressed that this model [44] could only be applied when the density contrast is small. Moreover, t, u, l and Fb take their conventional meaning as time, velocity vector, dynamic viscosity and body force respectively. Additionally, P is the dynamic pressure and rP is related to the surface tension force

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

rP ¼ /rl/ þ rp0

11

ð3Þ

In Eq. (3), / is the order parameter, p0 is the ideal gas pressure. Substituting Eq. (3) into Eq. (2) gives

@ðnuÞ þ r  ðn uuÞ ¼ rð/l/ þ p0 Þ þ l/ r/ þ lr2 u þ F b @t

ð4Þ

where l/ is the chemical potential derived from free energy density functional, l/r/ is surface tension force. The freeenergy functional [31] for an immiscible two-phase fluid system without a solid boundary takes the form of

F¼

Z

½AWð/Þ þ 0:5jðr/Þ2 dV

ð5Þ

V

Eq. (5) shows that the free energy includes two parts: a bulk energy density W(/) weighted by a coefficient A and an interface 2 gradient energy ðr2/Þ scaled by a coefficient j. A commonly applied double-well form bulk energy density is adopted

Wð/Þ ¼ ð/ þ /# Þ2 ð/  /# Þ2 #

ð6Þ

#

where / and / indicate two bulk regions of order parameter field. According to the definition of chemical potential, which is the rate of change of the free energy with respect to the phase concentration, it could be derived that

l/ ¼ 4Að/3  /#2 /Þ  jr2 /

ð7Þ

Additionally, in the phase-field model, the order parameter profile along the normal direction of an equilibrium surface could be written as

/ ¼ /# tanhð21=wÞ

ð8Þ

where 1 is the coordinate perpendicular to the interface and w is interface width which can be evaluated through

w¼

pffiffiffiffiffiffiffiffiffiffiffiffi 2j=A /#

ð9Þ

Meanwhile, the surface tension could be calculated by an equilibrium flat interface through

r¼j

 2 @/ d1 @1 1

Z

1

ð10Þ

It could further be derived by

r¼

pffiffiffiffiffiffiffiffiffiffiffiffi 32jA #3 / 3

ð11Þ

In the numerical simulation, the values of order parameter in the bulk region /#, interface thickness w and surface tension r are given. Eqs. (9) and (11) are then used to determine j and A. Eqs. (1) and (4) are the governing equations of multiphase flows, normally called Navier–Stokes (NS) equations. In the realm of the lattice Boltzmann equation (LBE), they could be presented by

fa ðx þ ea dt; t þ dtÞ ¼ fa ðx; tÞ þ Xa

ð12Þ

In the above equation, the hydrodynamic process proceeds through streaming and collision of density distribution function fa. In Eq. (12), a is the discrete lattice direction and ea represents lattice velocity in each lattice direction. In addition, Xa is the collision operator including additional force terms [45]

Xa ¼

faeq ðx; tÞ  fa ðx; tÞ

sn

þ

   wa 1 ðea  uÞ ðe 1   uÞ þ e a a  ðl/ r/ þ F b Þdt 2sn c2s c2s

ð13Þ

with wa being coefficient that depends on the lattice model used. Moreover, the relationship between the relaxation time sn and the kinematic viscosity t is expressed as



t ¼ sn 

 1 2 c dt 2 s

ð14Þ

Using conservation laws, the macroscopic variables could be evaluated through

n¼

X faeq

ð15Þ

a

X 1 faeq ea þ ðl/ r/ þ F b Þdt 2 a   X 2 faeq eal eak /l/ þ cs n dlk þ nul uk ¼ un ¼

a

ð16Þ ð17Þ

12

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

and the equilibrium distribution function faeq that satisfies the conservation laws reads

  3 9 faeq ¼ wa Aa þ wa n 3eal ul  u2 þ ul uk eal eak 2 2

ð18Þ

If D2Q9 lattice model, as shown in Fig. 3, is used, the coefficients in Eq. (18) are

(

15ð/l þn=3Þ

/ A0 ¼ 94 n  4 A18 ¼ 3ð/l/ þ n=3Þ

ð19Þ

8 > < w0 ¼ 4=9 w1;3;5;7 ¼ 1=9 > : w2;4;6;8 ¼ 1=36

ð20Þ

2.2. The interface capturing equation In the phase field model, the deformation and aggregation of interface is governed by the convective-diffusive Cahn– Hilliard equation. It is expressed as

@/ þ ur  / ¼ M r2 l/ @t

ð21Þ

where M represents diffusion rate, and it is also named as mobility. The corresponding lattice Boltzmann equation takes the form of [46]

g b ðx þ eb dt; t þ dtÞ ¼ g b ðx; tÞ þ Xb

ð22Þ

with collision term

Xb ¼

g eq b ðx; tÞ  g b ðx; tÞ

ð23Þ

s/

In Eq. (23), b is lattice direction in D2Q9 model and s/ is the relaxation time which can be chosen as

s/ ¼ sn

ð24Þ

Moreover, the distribution function satisfies the following conservation laws

/¼

X eq gb

ð25Þ

b

/ul ¼

X b

Cl/ dlk ¼

g eq b ebl X

ð26Þ

g eq b ebl ebk

ð27Þ

b

Fig. 3. Sketch of D2Q9 model in the lattice Boltzmann method for flow field.

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

13

In Eq. (27), C is a coefficient relevant to mobility M by

M¼



s/ 

 1 Cdt 2

ð28Þ

In this context, the equilibrium distribution function is

  /elb ul g eq ¼ w B þ b b b c2s

ð29Þ

with relevant coefficients set as

h i 8 < B0 ¼ w1 /  ð1  w0 Þ Cl2/ c 0 s

:B

18

¼

ð30Þ

Cl/ c2s

where the coefficients wb take the same values as those in Eq. (20). 2.3. Partial wetting boundary conditions A general phase field framework to model multiphase flows without considering a solid boundary has been laid previously. When a solid surface is present in multiphase flows, and the interactions among fluids and solid surface are assumed to be short-range [47], the surface characteristics is represented by a wall free energy in the total free energy functional

F¼

Z

½ AWð/Þ þ 0:5kðr/Þ2 dV þ

V

Z

Ws ð/s ÞdS

ð31Þ

S

where S denotes the solid surface, /s is the order parameter on the wall and Ws(/s) is the wall free energy density which could take the form of [32]

~ /s Ws ð/s Þ ¼ x

ð32Þ

~ describes the essential feature of a solid surface and it is also called as wetting potential. Minimizing the total free where x energy functional subjects to the natural boundary condition (the system is assumed to be in equilibrium) gives the boundary condition

~ jn  ðr/Þs ¼ W0s ð/s Þ ¼ x

ð33Þ

where n indicates the local unit outward normal on a solid surface. From a mathematical perspective of view, this is a Neumann boundary condition for the order parameter. Moreover, four equilibrium solutions of the order parameter on the wall are (a systematic derivation could be found in reference [43])

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ j xj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 1  j xj /s2 ¼ /

/s1 ¼ /#

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  jxj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 1 þ jxj /s4 ¼ / /s3 ¼ /#

ð34aÞ ð34bÞ ð35aÞ ð35bÞ

with the non-dimensional wetting potential x defined as

x¼

~ x pffiffiffiffiffiffiffiffiffi ð/# Þ2 2jA

ð36Þ

The relation between x and the equilibrium contact angle reads

cos heq ¼

i 3 3 1h ð1 þ xÞ2  ð1  xÞ2 2

ð37Þ

Based on the monotonous relationship (as shown in Fig. 4) established by Eq. (37), the value of x could be uniquely deter~ could also be calculated in a straightforward way. mined with a given equilibrium contact angle. Thereafter x In addition to the boundary condition expressed by Eq. (33), a no mass flux condition of the chemical potential is applied

n  ðrl/ Þs ¼ 0

ð38Þ

It is a Neumann boundary condition for the chemical potential. Conclusively, the boundary conditions that should be enforced on a solid surface are

14

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 4. The non-dimensional wetting potential versus the equilibrium contact angle.

us ¼ 0   ~ @/ x ¼ @n s j   @ l/ ¼0 @n s

ð39Þ ð40Þ ð41Þ

3. Immersed boundary method for Dirichlet boundary condition In the present study, both the Dirichlet and Neumann boundary conditions are consistently implemented through IBM. For the sake of self-completeness, a brief introduction on implementation of no-slip (Dirichlet) boundary condition is provided in the first place. To enforce the no-slip boundary condition, a recently developed boundary condition-enforced IB-LBM [12] is adopted in this work. In this method, the velocity correction term is determined in a way that guarantees the no-slip boundary condition to be enforced. To illustrate the idea, the set of LBEs for the flow field are rearranged as

1 eq fa ðx þ ea dt; t þ dtÞ  fa ðx; tÞ ¼  fa ðx; tÞ  fa ðx; tÞ þ F a dt s n    wa 1 ðea  uÞ ðea  uÞ þ  ðl/ r/ þ F b Þ e Fa ¼ 2 1  a 2sn c2s cs X 1 ea fa þ ðl/ r/ þ F b Þdt nu ¼ 2 a

ð42Þ ð43Þ ð44Þ

As can be seen from Eq. (44), the velocity consists of three parts. One arises from the density distribution function. The other two are attributed to the force density F = (l/r/ + Fb). It should be stressed that two distinct types of forces contribute to F, namely, the interfacial force l/r/[39] and the forcing term Fb that represents the effect of immersed boundary. Hence we could define three velocity components as below

un ¼

1X ea fa n a

1 l r/dt 2n / 1 F b dt du ¼ 2n u/ ¼

ð45Þ ð46Þ ð47Þ

with this definition, the overall velocity is

u ¼ un þ u/ þ du

ð48Þ

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

15

where u is the overall velocity, un is from the density distribution function without considering the presence of solid boundary, u/ is caused by interfacial force and du is the velocity correction term caused only by no-slip boundary condition. Based on this analysis, an intermediate velocity u⁄ is defined as

u  ¼ un þ u /

ð49Þ

Therefore, the overall velocity can be rearranged as

u ¼ u þ du

ð50Þ

In the boundary condition-enforced IB-LBM, the velocity correction term du is set as unknown and evaluated implicitly. To solve du on Eulerian point, dulB is first defined to represent velocity correction on Lagrangian point. In this context, du can be obtained by Dirac delta function

Z

duðx; tÞ ¼

C

    dulB X lB ; t d x  X lB ds

ð51Þ

  where d x  X lB is smoothly approximated by continuous kernel distribution Dij[48]

   1  D xij  X lB ¼ 2 d xij  xlB d yij  ylB h

dðrÞ ¼

ð52Þ

ð1=4Þ½1 þ cosðpjrj=2Þ; jrj 6 2 0;

ð53Þ

jrj > 2

Using the continuous delta function, Eq. (51) could be written as

X

duðx; tÞ ¼

    dulB X lB ; t Dij xij  X lB Ds

ð54Þ

l¼1...m

In the equation above, Ds is the arc length between two neighbouring boundary points. Recalling Eq. (50), the corrected velocity in Eulerian domain can also be written as

uðxij ; tÞ ¼ u ðxij ; tÞ þ duðxij ; tÞ

ð55Þ

To satisfy  the no-slip boundary condition, the interpolated fluid velocity on a boundary point must equal to the wall velocity U lB X lB ; t at the same place. The mathematical description is

  X   XX U lB X lB ; t ¼ u ðxij ; tÞDij DxDy þ dulB X lB ; t Dij Dsl Dij DxDy i;j

i;j

ð56Þ

l

or in a matrix form

AX ¼ B

ð57Þ

with A, B and X expressed as

0

d11 Bd B 21 A¼B B .. @ .

d12 d22 .. .

  .. .

dm1

dm2

   dmn

0

1

0 d11 C B C B d21 C B CB . C @ .. A m dm1 UB U 1B

B 2 B UB B B¼B . B . @ .

d1n d2n .. .

10 B d11 CB dB21 CB CB B CB .. A@ . dBn1

dB12 dB22 .. . dBn2

   dB1m

C    dB2m C C .. .. C C . . A    dBnm

d12



d1n

10

d22 .. .

 .. .

d2n .. .

C B u C CB 2 C CB C CB .. C A@ . A

dm2

   dmn

 T X ¼ du1B ; du2B ; . . . ; dum B

1

u1

ð58Þ

1 ð59Þ

un ð60Þ

where subscript m is the number of Lagrangian points; and n is the number of adjacent Eulerian points around the boundary. With dulB determined in this system, it can then be distributed back to Eulerian points and the resulting force density can be evaluated through

F b ¼ 2ndu=dt

ð61Þ

With the velocity corrected in this manner, the no-slip boundary condition can be enforced and streamline penetration to the boundary, which is commonly observed in the conventional IBM, could be eliminated [12,49–51].

16

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

4. Immersed boundary method for Neumann boundary condition The Dirichlet boundary condition is implemented through a boundary condition-enforced IBM in the previous section. Nonetheless, in a phase field description of contact line dynamics, both Dirichlet and Neumann boundary conditions are present. As introduced formerly, although numerous works have been devoted to refine the implementation of Dirichlet boundary condition in IBM, an IBM application with Neumann boundary condition is still unavailable at present to the best of our knowledge. In this work, the Neumann boundary condition is implemented based on physical conservation law. For simplicity, we will first consider a linear diffusion equation with Neumann boundary condition to illustrate the implementation process. Then we demonstrate how to apply the approach to implement wetting condition in contact line dynamics.

4.1. Flux contribution at the control surface to dependent variable in a control volume To elaborate the implementation process, we consider the following diffusion equation

@/ ¼ ar2 / @t

ð62Þ

with Neumann boundary condition vector q as

@/ @n

specified. In Eq. (62), a is a constant, r2 is the Laplacian operator. If we define a flux

q ¼ ðqx ; qy ; qz Þ ¼ ar/

ð63Þ

then, Eq. (62) can be rearranged as

@/ ¼ r  q @t

ð64Þ

As shown in Fig. 5, Eq. (64) can be obtained by applying the physical conservation law to a control volume. There are 6 control surfaces for the control volume.  The flux at control surface of x = 0 is qxdydz, which is into the control volume, while the x flux at control surface of x = dx is qx þ @q dx dydz, which is out the control volume. The net flux into the control volume by @x x these two surfaces is  @q dxdydz. Overall, the net flux into the control volume by the 6 control surfaces is (r  q)dxdydz. @x From physical conservation law, this net flux must be equal to the rate of change of / within the control volume, @/  dxdydz, that is @/ dxdydz ¼ ðr  qÞdxdydz. As a result, Eq. (64) is obtained. From the above process, for a general control @t @t surface, which has an out normal direction n = (nx, ny, nz), its flux contribution to the control volume is

ðqx nx þ qy ny þ qz nz ÞdS ¼ qn dS

ð65Þ

where dS is the area of the control surface, and qn ¼ a @/ . This means that the flux qndS directly contributes to @/ when the @n @t single control surface is considered. We will use this feature to correct / when the Neumann boundary condition is implemented.

Fig. 5. Sketch of a control volume with flux along three directions respectively.

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

17

4.2. Implementation of Neumann boundary condition in the context of IBM In the context of IBM, the solution of Eq. (62) can be obtained by the following steps. In the first step, we solve Eq. (62) in the whole domain including interior and exterior of the immersed object. For simplicity, we note its solution as /⁄. That is, /⁄ satisfies

@/ ¼ ar2 / @t

ð66Þ

with /⁄, we can calculate its normal derivative at the boundary point by the following way. Using interpolation, the first order derivatives at the boundary point can be calculated by

  @/  i  X @/ 2 XB ; t ¼ ðxj ; tÞD xj  X iB h @x @x j   @/  i  X @/ 2 XB ; t ¼ ðxj ; tÞD xj  X iB h @y @y j where

@/ @x



 X iB ; t and



@/ @y



ð67Þ ð68Þ

 X iB ; t represent the first order derivatives of /⁄ with respect to x and y at the boundary point X iB ,



ðxj ; tÞ and @/ ðxj ; tÞ are the first order derivatives of /⁄ with respect to x and y at Eulerian point xj. Note that the derivwhile @/ @x @y atives at Eulerian points are obtained by the second order central difference schemes. Finally, the normal derivative at the boundary point is calculated by

@/  i  @/  i  @/  i  XB; t ¼ X B ; t nxi þ X B ; t nyi @n @x @y

ð69Þ

For a general case, the computed @/⁄/@n is not equal to the given @//@n at the boundary point. As discussed in Section 4.1, their difference will contribute as a surface flux to correct / value at surrounding Eulerian points. In the context of IBM, the whole domain including interior and exterior of the immersed object is used as the computational domain. Thus, at a boundary point, there are two normal directions. One is to point to the flow domain while the other is to direct into the inside of immersed object. The surface fluxes from two directions will affect / field at surrounding Eulerian points. As shown in Fig. 6, due to the feature of Dirac delta function interpolation, the surface flux on a small area dS will only affect its surrounding Eulerian points in the box of jx  XBj 6 2h. In fact, for any Eulerian point in the box, its control volume must enclose the small surface area dS. For this case, the surface fluxes from two opposite directions of dS will both contribute positively to / in the control volume. Thus, due to unsatisfying of Neumann boundary condition (offset of normal derivative), the surface flux from dS can be written as



dqn X iB ; t



 @/ ¼ 2a @n 

@/  @n Specified

! ð70Þ

The flux dqn in Eq. (70) will be used to correct / at Eulerian points in the box. Suppose that the correction is noted as d/. Following the concept of IBM for distributing the surface force (flux) at a boundary segment to its surrounding Eulerian points by the delta function (the effect of the boundary to the surrounding Eulerian points is decreased as the distance of a Eulerian point from the boundary is increased) [1,10], we can have

Fig. 6. Illustration of flow domain, immersed boundary points and influence region of boundary points to surrounding fluids.

18

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

 d/ðxj ; tÞ X  i   ¼ dqn X B ; t Dij xj  X iB Dsi dt i

ð71Þ

The above process is similar to the velocity correction in Eq. (54). Once d/(xj, t) is obtained by Eq. (71), we can correct / field by

/ðxj ; tÞ ¼ / ðxj ; tÞ þ d/ðxj ; tÞ

ð72Þ

Eq. (71) is applied for the time-dependent diffusion Eq. (62). For this case, IBM takes an explicit way to update the solution in time due to the effect of boundary. For a time-independent problem, an iterative process must be taken first before IBM is applied. To illustrate the iterative process, we consider the following time-independent equation

/ ¼ ar2 /

ð73Þ

with Neumann boundary condition. After numerical discretization, the resultant equation system of Eq. (73) can be solved by a direct method. It can also be solved by an iterative method. One of iterative processes can be written as

/mþ1 ¼ ar2 /m

ð74Þ

where m is the iteration number. To start the iteration, we need to give an initial guess of /. Then from Eq. (74), we can get /⁄. Like time-dependent case, the computed @ /⁄/@n is in general not equal to the given @//@n at the boundary point. Then we can follow the same procedure as for the time-dependent problem to compute dqn by using Eq. (70). With dqn, the correction d/ can be calculated by

d/ðxj Þ ¼

 X  i  dqn X B Dij xj  X iB Dsi

ð75Þ

i

Eq. (75) is used to determine d/(xj) and it is the same as Eq. (71) if one takes dt as 1 in Eq. (71). Next, / is updated through Eq. (72) at Eulerian points around the immersed boundary. After that, we move to the next iteration. With updated / value given by Eq. (72), we can get a new /⁄ from Eq. (74) (updated / value is applied to the right side of Eq. (74)). We carry on this iteration process until convergence criterion is satisfied. 4.3. Application to wetting boundary condition in contact line dynamics The two Neumann boundary conditions (40) and (41) are related to the solution of Eq. (21) and Eq. (7) respectively. For demonstration, we write Eq. (21) and Eq. (7) to the following two forms

@/ þ u  r/ ¼ r  q/ @t l/ ¼ 4ð/3  /#2 /Þ  r  ql

ð76Þ ð77Þ

where

q/ ¼ Mrl/ ;

ql ¼ jr/

ð78Þ

Following the procedure in Section 4.2, we first solve Eqs. (76) and (77) in the whole domain without consideration of wetting condition to obtain /⁄ and l/ . In other words, /⁄ and l/ satisfy the following equations

@/ þ u  r/ ¼ r  q/ @t l/ ¼ 4ð/3  /#2 / Þ  r  ql

ð79Þ ð80Þ

Following Eqs. (67)–(69), with obtained /⁄ and l/ , we can easily compute @/⁄/@n and @ l/ =@n. In general, the calculated @/⁄/@n and @ l/ =@n are not equal to the given Neumann boundary conditions. Their differences will generate the surface fluxes dq/n and dqln. According to Eqs. (70) and (78), dq/n and dqln can be expressed as

   @ l/  @ l/  dq/n X iB ; t ¼ 2M   @n Specified @n !     @/ @/ dqln X iB ; t ¼ 2j  @n Specified @n

!

ð81Þ ð82Þ

Note that dq/n is to correct @ //@t in Eq. (76), while dqln is to correct l/ in Eq. (77). If we set / correction as d/ and l/ correction as dl/, following Eqs. (71) and (75), we have

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

    d/ðxj ; tÞ X ¼ dq/n X iB ; t Dij xj  X iB Dsi dt i     X dl/ ðxj ; tÞ ¼ dqln X iB ; t Dij xj  X iB Dsi

19

ð83Þ ð84Þ

i

After obtaining d/ and dl/, the corrected / and l/ at Eulerian points can be computed by

/ðxj ; tÞ ¼ / ðxj ; tÞ þ d/ðxj ; tÞ

ð85Þ

lðxj ; tÞ ¼ l/ ðxj ; tÞ þ dl/ ðxj ; tÞ

ð86Þ

At present, the updated l/ and / have included the influence of the solid boundary. Thereafter, they are readily adopted to evaluate the resulting forcing and velocity terms in Eqs. (43) and (44). In this manner, the Neumann boundary conditions are implemented in the framework of immersed boundary method. 5. Simulation procedures In summary, the simulation of contact line problems through IBM involves: (1) the Dirichlet boundary condition is treated by a boundary condition-enforced IBM and (2) the Neumann boundary condition is embodied analogous to the way that a flux affects its relevant physical parameter in a control volume as elaborated in Section 4. To provide an outline of the algorithm, the simulation procedures are presented below, (1) Set the initial flow fields, compute the coefficient matrix A in Eq. (57) and evaluate A1; (2) Using Eqs. (12) and (22) to obtain the density distribution function at time level t = tn (initial values of du, dl/ and d/ are set to zero) and compute the macroscopic variables; (3) Solve equation system (57) to determine the velocity correction term du at all boundary points and distribute them to Eulerian points; (4) Apply Eq. (82) to evaluate dqln on the boundary and Eq. (84) to compute dl/ on Eulerian points. Update l/ using Eq. (86); (5) Obtain dq/n through Eq. (81), compute d//dt on the Eulerian points using Eq. (83) and update the overall / according to Eq. (85); (6) Evaluate force density in Eq. (43), update velocity in Eq. (44) and compute the equilibrium distribution function; (7) Repeat steps (2) to (6) until convergence criteria are reached. 6. Results and Discussion In this section, the proposed IBM for Neumann boundary conditions is first verified via simulation of transition layer around a circular cylinder, which has analytical solutions to be compared with. Subsequently, the numerical behavior of simulating dynamic problems is examined in detail through droplet dewetting. In addition, the method is further validated through simulation of droplet spreading on a flat surface over a wide range of partial wetting regimes. Furthermore, the ability of the present method to handle Neumann boundary condition on complex geometries, probably the most desirable feature of proposed IBM, is demonstrated through simulation of the contact line dynamics on circular cylinder (s). In this study, we restrict ourselves to investigate two density matched fluids’ interaction with solid boundary. The present computation is conducted in the lattice Boltzmann system. To make the computation in the lattice Boltzmann system be consistent with the computation in the physical system, the corresponding non-dimensional variables have to be kept the same. 6.1. Transition layers on hydrophilic and hydrophobic walls  When a solid wall is not neutral, that is, @/ has non-zero value, the value of order parameter on the wall will deviate @n s from that in the bulk region [43] and a transition layer can be formed along the solid surface. In this case, analytical solutions exist for the order parameter on the solid surface. Hence it is adopted to test the performance of the proposed algorithm. The problem considered in this section is a square domain with a circular cylinder located in the centre of the domain. In the physical system, the cylinder radius is chosen as the reference length, that is, the non-dimensional radius is 1, and the non-dimensional length of computational domain is taken as 3.5. As mentioned above, to make the computation be consistent between the physical system and the lattice Boltzmann system, the non-dimensional radius and the length of computational domain in the lattice Boltzmann system must also be kept as 1 and 3.5 respectively. For the present LBM calculation, dx is set as 1 for simplicity. Thus in the lattice Boltzmann system, the number of mesh points for the radius, Nradius, is the radius of the cylinder, and the number of mesh points used in each direction, N mesh size is N  N), is actually the length of domain. Since the cylinder radius Nradius is taken as the reference length, in the lattice Boltzmann system, the non-dimensional radius of the cylinder is Nradius/Nradius = 1, and the non-dimensional length of the square computational domain should

20

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

be N/Nradius = 3.5. This general rule must be obeyed when different mesh size and radius of cylinder are used in the lattice Boltzmann system. Hereafter, all the computations are conducted in the lattice Boltzmann system. 6.1.1. Effect of the transition layer thickness A solid circular cylinder with radius of 40 represented by 148 Lagrangian points is centered at (70, 70) in a 140  140 computational domain. The flow domain is initialized as: q = 1 for both fluids, /# is set as 1 in the bulk region, r = 0.001, sn and s/ take the same value as 0.75. (These physical parameters will be used for the cases hereafter unless otherwise stated). In this case, the gradient of the order parameter is fixed as x = 0.334933 which corresponds to static contact angle of 60°, and the transition layer thickness (w in Eq. (9)) is varied from 4.0 to 6.0 with increment of 0.5. Moreover, neutral wetting boundary conditions are applied on the outer boundaries. To be concrete, for density distribution functions:  fa jBoundary ¼ faeq Boundary and gajBoundary = gaeqjBoundary with a running over 0 to 8 lattice directions. For the macroscopic variables, natural boundary conditions (or neutral wetting boundary conditions) are applied. To be concrete: ujBoundary = 0, qjBoundary = 1, j/nþ1 /n j /jBoundary = 1, r/jBoundary = 0 and rl/jBoundary = 0. In addition, the convergence criterion is set as s /n s < 105 for this case. s

Fig. 7 displays the initial flow field where / takes the same value in the whole flow domain. Fig. 8 demonstrates the flow domain after applying the boundary condition through IBM when the transition layer thickness is set as 4.0. The transition layer attached to the solid boundary could be seen clearly in this figure. As introduced previously, the analytical solutions exist for this problem. With a given non-dimensional wetting potential x, the order parameter on the boundary could be evaluated through Eq. (35). Table 1 compares the / value on the boundary obtained through numerical simulation with the analytical solution when the surface thickness varies from 4.0 to 6.0. The relative error is defined as

ERelative ¼

j/Theoretical  /numerical j /Theoretical

ð87Þ

It can be found in Table 1 that the present results agree well with the theoretical prediction. The present results also demonstrate a common feature of diffuse interface method. Concretely speaking, it could be observed that the relative error reduces as transition layer thickness increases. This is due to the fact that, in the diffuse interface approach, the transition layer should be thick enough to resolve the interfacial dynamics accurately. Moreover, on account of the circular cylinder being a central-symmetric geometry, the solution of this problem should also be isotropic. Hence, it is natural to raise a question, that is, whether the / obtained on the solid boundary demonstrates isotropy and if there is any oscillation of the numerical solution. To examine this issue, one could first look at Fig. 8. No oscillation could be observed in the flow domain. To be more precise, quantitative comparison is made in Table 2. It lists the maximum error between the local and the average / value. It shows that the maximum errors are less than 6  104 and the errors monotonously decrease as the thickness increases. The above comparisons demonstrate good isotropy of the obtained numerical results. In addition, the influence of grid size on numerical solution is examined. The interface width is chosen as 7 grid length and other computational parameters are set the same. Three sets of grids (140  140,210  210 and 280  280) are tested and listed together with results in Table 3. It can be seen that the numerical results steadily approach theoretical solution

Fig. 7. Initial flow field with solid particle located at the centre.

21

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 8. Transition layer generated along the solid surface due to implementation of wetting boundary conditions through immersed boundary method.

Table 1 Comparison of / on cylinder surface with theoretical prediction. /theoretical = 1.155 Interface thickness

/numerical

Relative error (%)

4.0 4.5 5.0 5.5 6.0

1.091 1.095 1.098 1.101 1.104

5.5 5.1 4.9 4.6 4.4

Table 2 Maximum error on the boundary. Interface thickness

4.0

4.5

5.0

5.5

6.0

j/  /averagejmax

5.634  104

5.573  104

5.543  104

5.403  104

5.336  104

Table 3 Influence of mesh size on / value on the boundary. /theoretical = 1.155 Mesh size

/numerical

Relative error (%)

140  140 210  210 280  280

1.114 1.122 1.127

3.54 2.86 2.42

as grid size is refined from 140  140 to 280  280. Note that for the three sets of grids, the non-dimensional length of the domain is kept as 3.5.

6.1.2. Effect of (@//@n)s Previous subsection simulates the case with a fixed Neumann boundary condition that is equivalent to contact angle of 60°. This subsection demonstrates capacity of the present method to handle different (@//@n)s (corresponding to 5° < heq < 175°). The surface thickness is taken as 9.0. The computation is carried on a 251  251 mesh size and the radius of the cylinder is represented by 65 grid points. Fig. 9 compares the numerical results with the theoretical values corresponding to different wetting potential predicted by Eq. (35). The data are listed in Table 4. It can be seen that the numerical results compare well with theoretical values when contact angle varies from 5° to 135°. It could also be found that both the absolute

22

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 9. Theoretical and numerical / values on the boundary versus the non-dimensional wetting potential x.

Table 4 Comparison of / value on the boundary. heq(°)

x

/Theoretical

/IBM

EAbsolute

ERelative(%)

5 30 45 60 75 90 105 120 135 150 175

0.679 0.586 0.476 0.335 0.173 0 0.173 0.335 0.476 0.586 0.679

1.295 1.259 1.214 1.155 1.083 1.0 0.91 0.816 0.724 0.643 0.567

1.23 1.202 1.168 1.122 1.065 1.0 0.928 0.853 0.779 0.716 0.652

0.065 0.057 0.046 0.033 0.018 0.0 0.018 0.037 0.055 0.073 0.085

5.062 4.566 3.789 2.857 1.662 0.0 1.978 4.534 7.597 11.331 14.996

and relative errors increase from zero when heq deviates from 90°. This is may be attributed to the nature of this problem. According to the relationship provided by Eq. (37) which is also shown in Fig. 4, it could be seen that when the equilibrium contact angle is close to 0° or 180°, the curve becomes very sharp. This indicates that in the region nearby 0° or 180°, the equilibrium contact angle is very sensitive to the x value. That is, when x is changed a little bit, the equilibrium contact angle will change a lot. So, in the region nearby 0° or 180°, a small numerical error for x would cause a large numerical error for the equilibrium contact angle. To further study this issue, the same problem is simulated by direct implementation of boundary conditions. It was found that the numerical results given by direct implementation of boundary conditions also show larger errors when the contact angle is close to 0° or 180°, but are more accurate than the IBM results. This means that the present solver may not be efficient to resolve the problem when the contact angle is close to 0° or 180°. It would need a very fine mesh to solve the problem in order to get an accurate solution.

6.2. Droplet dewetting Another contact line problem, droplet dewetting, is used to test numerical behaviour of current method for dynamic problems. Both accuracy and grid-independency are examined in the first place. Thereafter, comparison between results of direct implementation of boundary conditions and those of immersed boundary method are carried out at different surface wettability. In dewetting problems, a droplet is initialized with a relatively small contact angle (60° for instance), while the plate where the droplet is placed is set as neutral wetting or hydrophobic (with large equilibrium contact angle equal or larger than 90°). Consequently, there is an initial contact angle difference D h = heqlibirum  hinitial. In this case, the droplet will move upward due to this initial difference. This situation is similar to the electrowetting experiment in which a voltage is suddenly applied. The speed of droplet motion depends on parameters such as the initial difference. To quantitatively characterize the droplet motion, two quantities are defined in the first place. They are averaged Y-position and Y-velocity component of the droplet:

23

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32 Table 5 Three sets of Lagrangian grid and ratio of Lagrangian grid spacing over Eulerian grid spacing DL/DE. Number of Lagrangian grid points

168

126

112

DL/DE

1.2

1.6

1.8

Fig. 10. Time evolution of YDrop (Eulerian grid: 300  210; DIBC: Direct implementation of boundary conditions; DL/DE: Mesh spacing ratio between Lagrangian grid and Eulerian grid).

Fig. 11. Time evolution of YDrop (Three different Eulerian grids: 300  210, 360  250 and 400  290, are used; DIBC: Direct implementation of boundary conditions).

P yij /ij Y Drop ¼ P /ij P v ij /ij V Drop ¼ P /ij

ð88aÞ ð88bÞ

First, the accuracy and grid-independency for this dynamic case are verified by changing both Eulerian and Lagrangian grids. Then dewetting on solid with different wettability is also examined. In these cases, the time evolutions of droplet height and velocity are compared with results by direct implementation of boundary conditions.

24

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 12. Time evolution of droplet VDrop (DIBC: Direct implementation of boundary conditions; IBM: immersed boundary method).

Fig. 13. Equilibrium state of the spreading droplet on a flat plate.

6.2.1. Grid-independency test for dynamic cases Grid-independency study is first performed, and the numerical results are compared with those obtained by direct implementation of boundary conditions.

6.2.1.1. Lagrangian grid-independency test. To verify the grid-independency of Lagrangian grid, three different sets of Lagrangian grid (enumerated in Table 5) are tested. The computational domain is set as 300  210 with a smooth plate (length equals to 200 lattice units) centered at (150, 30). Concurrently, a body-conformal grid of 201  151 with a plate centered at (100, 0) is used for direct implementation of boundary conditions. Moreover, a droplet of radius 60 is centered on the plate with initial contact angle of 60°. Meanwhile, the equilibrium contact angle of the plate is set as 120°. Neutral wetting boundary conditions are applied on upper and lower walls, while periodic boundary conditions are used for the left and right boundaries. Furthermore, sf is fixed as 0.65 (corresponding to kinematic viscosity of 0.05) hereafter unless mentioned otherwise. The time evolution of YDrop obtained by immersed boundary method on three sets of Lagrangian grid are compared with that by direct implementation of boundary conditions in Fig. 10. It can be seen from Fig. 10 that the dynamic processes provided by IBM are almost identical when Lagrangian grid size varies in a wide range (normally, the value of DL/DE is larger than 1.0 and less than 2.0 in IBM). Furthermore, all results of IBM are compared well with those given by direct implementation of boundary conditions. This shows that the Lagrangian resolution in this study is adequate.

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

25

Fig. 14. Two ways to calculate equilibrium contact angle for droplet spreading on smooth surface. (a) Local / level contours when the droplet approaches heq = 60°. (b) Sketch of arc shape of an equilibrium droplet.

Fig. 15. The non-dimensional wetting potential versus the theoretical and numerical equilibrium contact angles.

26

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32 Table 6 Comparison of equilibrium contact angle on a flat plate. Theoretical value

IBM

Direct implementation of wetting BCs

60° 90° 120°

60.6° 90° 120.5°

60.4° 90° 120°

6.2.1.2. Eulerian grid-independency test. The grid-independency of Eulerian grid is also examined. In this part, numerical simulation is performed on three Eulerian grids. The evolution of YDrop for three cases is compared with that by direct implementation of boundary conditions in Fig. 11. It can be seen that the change in Eulerian grid size also has little influence on the dynamic process obtained and all numerical results by IBM show good agreement with those by direct implementation of boundary conditions. This study also shows that the Eulerian grid size used in this work is fine enough to get accurate numerical results. 6.2.2. Influence of surface wettability Besides examining the accuracy and grid-independency of the proposed algorithm to simulate dynamic process, the surface wettability is varied in this part and droplet velocity obtained by IBM is compared with that by direct implementation of boundary conditions. As demonstrated previously, Eulerian grid of 301  211 with DL/D E varying from 1.2 to 1.8 can provide stable and accurate solution. The surface wettability is set as contact angle of 90°, 120° and 150° respectively. Fig. 12 compares IBM results of these three different cases with those by direct implementation of boundary conditions. It can be observed that the time evolution of VDrop provided by IBM basically compares well with that on the body conformal grid even when the surface demonstrates super-hydrophobicity (i.e. heq = 120° and heq = 150°). 6.3. Droplet spreading on a plate in the partial wetting regime We now apply the solver to simulate another practical problem: droplet spreading in the partial wetting regime. In this case, the numerical contact angles are compared with theoretical predictions. Moreover, the time evolution of resolved droplet height and base diameter are compared with that obtained by direct implementation of the wetting boundary conditions. For this problem, a flat plate with length of 400 is centered at (250, 50) in a computational domain of 500  230. A droplet with radius of 45 is initially located above a plate with contact angle of 160°. The boundary conditions corresponding to the equilibrium contact angle heq from 30° to 150° are tested. The interface thickness is set as 4.5. The non-dimensional time is defined as t = Nhr/Dnt[52], in which N indicates time steps, h means the mesh spacing and D is the diameter of the droplet. Figs. 13(a)–(e) display the equilibrium status of the spreading droplet with different heq. To identify droplet contact angle, two ways are available. One way is to use a protractor and directly measure it on a phase field contour (only one contour line / = 0 that represents the interface will be used) after the droplet evolves to its equilibrium state, as shown in Fig. 14(a), in which scl is tangent line of level curve / = 0 on the wall. The other way is based on the fact that the droplet takes an arc shape to minimize the free energy after it approaches equilibrium on a smooth surface. It is the way adopted

Fig. 16. The time evolution of non-dimensional droplet height and diameter (normalized by the initial droplet diameter) when heq = 60°.

27

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 17. Spreading process of a droplet with three different equilibrium contact angles.

to calculate heq in this work. As sketched in Fig. 14 (b), when heq < 90°, we denote radius of droplet, height of droplet and  2 r with radius of the circle as RDrop, HDrop and R respectively. It could be derived that the contact angle is h ¼ p  cos1 1k 1þk2 r

kr = RDrop/HDrop. To be specific, it can be found that cos heq ¼ H2Drop þR2Drop

RH

RHDrop . R

2 2 2 Moreover, ðR  HDrop  Þ þ RDrop  ¼ R which gives H2

R2

Drop . Substituting this relation into cos heq ¼ RDrop , we can have heq ¼ p  cos1 HDrop . If we denote kr = 2 2 Drop þRDrop  2 1 1kr . Following this procedure, the same equation could be derived RDrop/HDrop, it could be simplified as heq ¼ p  cos 1þk2

R¼

2HDrop

r

when heq P 90°. (In this work, the values of RDrop and HDrop are tracked in the code based on the contour line of / = 0). Fig. 15 compares numerical results of current method with analytical solution given by Eq. (37). It demonstrates that the numerical results agree well with theoretical values in a wide range, especially when the contact angle is less than 135°. This may be due to the reason that for a droplet with a fixed volume, the base radius decreases for larger contact angle and hence the accuracy is negatively affected. This tendency is also observed when the boundary condition is directly implemented [53]. In addition, quantitative comparisons are made in Table 6 between theoretical values and those obtained by direct

28

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 18. Level curves of order parameter together with velocity vector field during spreading process. The left panel is the result by direct implementation of boundary conditions and the right panel is result by IBM (heq = 60°).

Fig. 19. Contact line on a single circular cylinder.

implementation of wetting boundary conditions. It could be seen that the numerical results obtained by IBM are almost identical with those by direct implementation of wetting boundary conditions. Beside comparison of the equilibrium contact angle, the time evolution of the droplet height and base diameter (normalized by initial droplet diameter) by IBM is compared with that by direct implementation of wetting boundary conditions on a body-conformal grid in Fig. 16. It could be found that dynamics provided by IBM is also in good agreement with that obtained by direct implementation of boundary conditions. Fig. 16 also demonstrates overall behavior of a spreading droplet. Initially, the base diameter undergoes dramatic changes while the droplet height stays almost the same when t < 1. Subsequently, the development of both the droplet height and diameter can be observed until the droplet relaxes to its equilibrium status [52]. This process can be seen more directly in Fig. 17. In addition, the level curves for the phase variable with velocity vector field at two time stages are plotted together with those by direct implementation of boundary conditions in Fig. 18. It can be seen that the flow field provided by IBM matches very well with that by direct implementation of boundary conditions. To further validate the numerical results, the spreading rate of this process is evaluated. Fitting the curves in Fig. 16 between 100 < t < 101 gives r / tn where n = 0.351. This value compares reasonably well with that reported in [52] which is n = 0.34. These results demonstrate that the present algorithm could not only reproduce the equilibrium results, but also can correctly predict dynamic process of moving contact line phenomenon.

29

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 20. Schematic depiction of contact angle definition on a circular cylinder.

Table 7 Comparison of equilibrium contact angle on circular cylinder. Theoretical value

60°

90°

120°

IBM

61° ± 1°

90° ± 1°

120° ± 1°

6.4. Contact line dynamics on a single and two alongside circular cylinders Most of the cases examined previously only consider the simple smooth surface in order to compare the results with those by direct implementation of boundary conditions. This subsection is devoted to demonstrate robustness of the proposed algorithm to handle problems with curved boundary. First, the contact line dynamics on a single circular cylinder is simulated. Then, contact line dynamics around two alongside cylinders are studied with different surface wettability. 6.4.1. Single cylinder In this case, a circular cylinder with radius of 40 is fixed at the center of flow domain of 200  200. Periodic boundary condition is applied at left and right boundaries and neutral wetting boundary condition is used for upper and lower wall. The fluid–fluid interface is initialized as a flat surface located at the middle of computational domain. The lower region is set as liquid phase with / = 1 and the ambient fluid is set as / = 1 (the order parameter in the cylinder is set the same as ambient fluid). Figs. 19(a)–(f) show the equilibrium status when theoretical contact angles are set as 60°, 90°and 120°. It could be seen that the fluid–fluid interface evolves along the solid boundary. For hydrophilic surface, the fluid–fluid interface will raise above the initial horizontal line. On the contrary, the fluid–fluid interface is lowered below the initial horizontal line. In this case, the numerical contact angle is measured directly based on the obtained flow field (contour line of / = 0). The definition of the equilibrium contact angle [54] is defined in Fig. 20. In this figure, nb is the out normal of the circular cylinder and sb is the tangential at the same position. These two directions are plotted at the three-phase contact line (contact point in this figure). Table 7 compares the numerical equilibrium contact angles with theoretical values. It is validated that the current numerical algorithm could accurately produce equilibrium contact angle when curved boundary is involved in the multiphase flow-solid interactions. 6.4.2. Two alongside cylinders To further unfold the robustness of the present algorithm to manipulate arbitrary number and shape of complex solid boundaries, the cases with two alongside circular cylinders holding the same as well as different wettability are simulated. The two circular cylinders with radius of 40 are fixed at (150, 50) and (250, 50) respectively in an expanded computational domain of 400  200. Fig. 21 shows the equilibrium flow field. The two cylinders in Fig. 21(a) have the same wettability with equilibrium contact angle of 60°. Meanwhile, the two cylinders in Fig. 21(b) have the same equilibrium contact angle of 120°. It could be

30

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

Fig. 21. Contact line on two alongside cylinders with the same as well as different surface wettability. (a) Both cylinders with the same surface wettability of heq = 60°. (b) Both cylinders with the same surface wettability of heq = 120°. (c) The left cylinder with surface wettability of heq = 60° and the right one with heq = 120°.

found in these figures that, for each cylinder, the fluid–fluid interface deformation is asymmetric. At two ends far away from the interval between two cylinders (the left and right free ends), the fluid–fluid interface relaxes to the prescribed heq. However, the fluid–fluid interface in the interval is raised (heq = 60°) or lowered (heq = 120°) as compared with initial horizontal line. This is due to the fact that when two cylinders are located close enough, a capillary interaction in the interval is generated in response to the overlap of perturbations in the meniscus shape [55]. In addition, another case with two cylinders having different wettability was also simulated. As shown in Fig. 21(c), the left one has heq = 60° and the right one has heq = 120°. In this case, the fluid–fluid interface also relaxes to the prescribed heq at two free ends. However, different from two cylinders with the same wettability, the interface is raised near the hydrophilic cylinder and lowered in the vicinity of the hydrophobic cylinder. In this case, the two cylinders, if they are allowed to move freely, will be pushed away from each other due to repulsive long-range force generated. In-depth investigation [54–57] of these phenomena is beyond the scope of present work and interesting readers could refer to the literature listed above and references therein. In summary, numerical examples in this subsection demonstrate that, the present algorithm could easily be adapted to contact line problems with curved boundaries as well as surfaces with chemical inhomogeneous characteristics. Hence, it could serve as an efficient approach to dig into the multiphase fluid–solid interaction problems. 7. Conclusions An immersed boundary-phase field-lattice Boltzmann method is developed to embody Neumann boundary condition for the first time in this study. The primary concept of current method is to utilize physical mechanism and interpret Neumann

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

31

boundary condition as contribution of flux from the surface to its relevant physical variables in a control volume. Using the concept of IBM, the flux is directly related to the correction of flow variables at Eulerian points. Several numerical experiments are carried out to examine the performance of the developed method for both steady and dynamic problems. Numerical results show that the method can accurately reproduce equilibrium parameters and can also correctly predict dynamic processes as compared with direct implementation of the same boundary conditions. Furthermore, its capacity to be adapted to geometrical and/or chemical patterned surface is also demonstrated. This work releases IBM from the long existed restriction that it can only handle Dirichlet boundary condition and sheds light on the implementation of IBM to ubiquitous fluid– solid interactions associated with Neumann boundary conditions. References [1] R. Mittal, G. Iaccarino, Immersed boundary methods, Annual Review of Fluid Mechanics 37 (2005) 239–261. [2] R.J. Leveque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM Journal on Numerical Analysis 31 (4) (1994) 1019–1044. [3] Z.L. Li, M.C. Lai, The immersed interface method for the Navier–Stokes equations with singular forces, Journal of Computational Physics 171 (2) (2001) 822–842. [4] L. Lee, R.J. Leveque, An immersed interface method for incompressible Navier–Stokes equations, SIAM Journal on Scientific Computing 25 (3) (2003) 832–856. [5] D.V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, Journal of Computational Physics 220 (1) (2006) 109–138. [6] R.P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method), Journal of Computational Physics 152 (2) (1999) 457–492. [7] X.D. Liu, R.P. Fedkiw, M.J. Kang, A boundary condition capturing method for Poisson’s equation on irregular domains, Journal of Computational Physics 160 (1) (2000) 151–178. [8] D. Dezeeuw, K.G. Powell, An adaptively refined Cartesian mesh solver for the Euler equations, Journal of Computational Physics 104 (1) (1993) 56–68. [9] H.S. Udaykumar, H.C. Kan, W. Shyy, R. TranSonTay, Multiphase dynamics in arbitrary geometries on fixed Cartesian grids, Journal of Computational Physics 137 (2) (1997) 366–405. [10] C.S. Peskin, flow patterns around heart valves-numerical method, Journal of Computational Physics 10 (2) (1972) 252–&. [11] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, Journal of Computational Physics 161 (1) (2000) 35–60. [12] J. Wu, C. Shu, Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications, Journal of Computational Physics 228 (6) (2009) 1963–1979. [13] A.M. Worthington, A Study of Splashes, Longmans, Green, and Co., London, New York, Bombay, Calcutta, 1908. [14] A.L. Yarin, Drop impact dynamics: splashing, spreading, receding, bouncing, Annual Review of Fluid Mechanics 38 (2006) 159–192. [15] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, E. Rolley, Wetting and spreading, Reviews of Modern Physics 81 (2) (2009) 739–805. [16] E.B. Dussan, S.H. Davis, Motion of a fluid–fluid interface along a solid-surface, Journal of Fluid Mechanics 65 (AUG12) (1974) 71–95. [17] E.B.V. Dussan, Spreading of liquids on solid surfaces: static and dynamic contact lines, Annual Review of Fluid Mechanics 11 (1979) 371–400. [18] M.Y. Zhou, P. Sheng, Dynamics of immicible-fluid displacement in a capillary-tube, Physical Review Letters 64 (8) (1990) 882–885. [19] P.J. Haley, M.J. Miksis, The effect of the contact line on droplet spreading, Journal of Fluid Mechanics 223 (1991) 57–81. [20] M. Renardy, Y. Renardy, J. Li, Numerical simulation of moving contact line problems using a volume-of-fluid method, Journal of Computational Physics 171 (2001) 243–263. [21] M. Bussmann, J. Mostaghimi, S. Chandra, On a three-dimensional volume tracking model of droplet impact, Physics of Fluids 11 (1999) 1406–1417. [22] M. Bussmann, J. Mostaghimi, S. Chandra, Modeling the splash of a droplet impacting a solid surface, Physics of Fluids 12 (2000) 3121–3132. [23] P.D.M. Spelt, A level-set approach for simulations of flows with multiple moving contact lines with hysteresis, Journal of Computational Physics 207 (2005) 389–404. [24] P.T. Yue, C.F. Zhou, J.J. Feng, Sharp-interface limit of the Cahn–Hilliard model for moving contact lines, Journal of Fluid Mechanics 645 (2010) 279–294. [25] J. S. Rowlinson, B. Widom, Molecular Theory of Capillarity. 1982: Oxford, Oxfordshire: Clarendon Press, 1982. [26] P. Seppecher, Moving contact lines in the Cahn–Hilliard theory, International Journal of Engineering Science 34 (9) (1996) 977–992. [27] D. Jacqmin, Contact-line dynamics of a diffuse fluid interface, Journal of Fluid Mechanics 402 (2000) 57–88. [28] L.K. Antanovskii, A phase field model of capillarity, Physics of Fluids 7 (4) (1995) 747–753. [29] R. Chella, J. Vinals, Mixing of a two-phase fluid by cavity flow, Physical Review E 53 (4) (1996) 3832–3840. [30] D.M. Anderson, G.B. McFadden, A diffuse-interface description of internal waves in a near-critical fluid, Physics of Fluids 9 (7) (1997) 1870–1879. [31] A.J. Briant, A.J. Wagner, J.M. Yeomans, Lattice Boltzmann simulations of contact line motion. I. Liquid–gas systems, Physical Review E 69 (3) (2004). [32] A.J. Briant, J.M. Yeomans, Lattice Boltzmann simulations of contact line motion. II. Binary fluids, Physical Review E 69 (3) (2004). [33] V.V. Khatavkar, P.D. Anderson, P.C. Duineveld, H.E.H. Meijer, Diffuse-interface modelling of droplet impact, Journal of Fluid Mechanics 581 (2007) 97– 127. [34] H. Ding, P.D.M. Spelt, Wetting condition in diffuse interface simulations of contact line motion, Physical Review E 75 (4) (2007). [35] J.J. Huang, C. Shu, Y.T. Chew, Lattice Boltzmann study of droplet motion inside a grooved channel, Physics of Fluids 21 (2) (2009). [36] P.T. Yue, J.J. Feng, Wall energy relaxation in the Cahn–Hilliard model for moving contact lines, Physics of Fluids 23 (1) (2011). [37] J.J. Huang, C. Shu, J.J. Feng, Y.T. Chew, A phase-field based hybrid Lattice–Boltzmann finite-volume method and its application to droplet motion under electrowetting control, Journal of Adhesion Science and Technology 26 (2012) 1825–1851. [38] R.R. Nourgaliev, T.N. Dinh, T.G. Theofanous, D. Joseph, The lattice Boltzmann equation method: theoretical interpretation, numerics and implications, International Journal of Multiphase Flow 29 (1) (2003) 117–169. [39] D. Jacqmin, Calculation of two-phase Navier–Stokes flows using phase-field modeling, Journal of Computational Physics 155 (1) (1999) 96–127. [40] H. Ding, T.G. Theofanous, The inertial regime of drop impact on an anisotropic porous substrate, Journal of Fluid Mechanics 691 (2012) 546–567. [41] M. Francois, W. Shyy, Computations of drop dynamics with the immersed boundary method, Part 2: Drop impact and heat transfer, Numerical Heat Transfer Part B-Fundamentals 44 (2) (2003) 119–143. [42] S. Manservisi, R. Scardovelli, A variational approach to the contact angle dynamics of spreading droplets, Computers and Fluids 38 (2) (2009) 406–424. [43] P. Papatzacos, Macroscopic two-phase flow in porous media assuming the diffuse-interface model at pore level, Transport in Porous Media 49 (2) (2002) 139–174. [44] M.R. Swift, E. Orlandini, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann simulations of liquid–gas and binary fluid systems, Physical Review E 54 (5) (1996) 5041–5052. [45] Z.L. Guo, C.G. Zheng, B.C. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Physical Review E 65 (4) (2002). [46] J.J. Huang, C. Shu, Y.T. Chew, Numerical investigation of transporting droplets by spatiotemporally controlling substrate wettability, Journal of Colloid and Interface Science 328 (2008) 124–133. [47] P.G. de Gennes, Wetting: statics and dynamics, Reviews of Modern Physics 57 (3) (1985) 827–863.

32

J.Y. Shao et al. / Journal of Computational Physics 234 (2013) 8–32

[48] Z.G. Feng, E.E. Michaelides, Proteus: a direct forcing method in the simulation of particulate flows, Journal of Computational Physics 202 (1) (2005) 20– 51. [49] J. Wu, C. Shu, Particulate flow simulation via a boundary condition-enforced immersed boundary-lattice Boltzmann scheme, Communications in Computational Physics 7 (4) (2010) 793–812. [50] J. Wu, C. Shu, Y.H. Zhang, Simulation of incompressible viscous flows around moving objects by a variant of immersed boundary-lattice Boltzmann method, International Journal for Numerical Methods in Fluids 62 (3) (2010) 327–354. [51] J. Wu, C. Shu, An improved immersed boundary-lattice Boltzmann method for simulating three-dimensional incompressible flows, Journal of Computational Physics 229 (13) (2010) 5022–5042. [52] V.V. Khatavkar, P.D. Anderson, H.E.H. Meijer, Capillary spreading of a droplet in the partially wetting regime using a diffuse-interface model, Journal of Fluid Mechanics 572 (2007) 367–387. [53] T. Lee, L. Liu, Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces, Journal of Computational Physics 229 (20) (2010) 8045–8062. [54] P. Singh, D.D. Joseph, Fluid dynamics of floating particles, Journal of Fluid Mechanics 530 (2005) 31–80. [55] P.A. Kralchevsky, K. Nagayama, Capillary forces between colloidal particles, Langmuir 10 (1) (1994) 23–36. [56] P.A. Kralchevsky, K. Nagayama, Capillary interactions between particles bound to interfaces, liquid films and biomembranes, Advances in Colloid and Interface Science 85 (2–3) (2000) 145–192. [57] P.C. Millett, Y.U. Wang, Diffuse-interface field approach to modeling arbitrarily-shaped particles at fluid–fluid interfaces, Journal of Colloid and Interface Science 353 (1) (2011) 46–51.