Smoothed profile-lattice Boltzmann method for non-penetration and wetting boundary conditions in two and three dimensions

Computers and Fluids 159 (2017) 64–80

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Smoothed profile-lattice Boltzmann method for non-penetration and wetting boundary conditions in two and three dimensions Takeshi Seta a,∗, Tomomi Uchiyama b, Noboru Takano a a b

Graduate School of Science and Engineering for Research, University of Toyama, Gofuku, Toyama-shi, Toyama 930-8555, Japan Institute of Materials and Systems for Sustainability, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan

a r t i c l e

i n f o

Article history: Received 1 March 2017 Revised 1 July 2017 Accepted 17 September 2017 Available online 20 September 2017 Keywords: Hydrophobicity Hydrophilicity Smoothed profile method Lattice Boltzmann method Non-penetration condition Wetting boundary condition

a b s t r a c t In this study, the smoothed profile-lattice Boltzmann method (SP-LBM) is proposed to determine the contact line dynamics on a hydrophobic or a hydrophilic curved wall. Two types of smoothed indicator functions are introduced, namely a function that identifies the solid domain for non-slip and non-penetration conditions and a function that denotes the boundary layer for no mass-flux and the wetting boundary conditions. In order to prevent fluid penetration into the solid boundary, the fluid-solid interaction force is computed based on the definition of the fluid velocity as proposed by Guo et al. [1]. In order to implement the Neumann boundary conditions for the order parameter and the chemical potential, the fluxes from the solid surfaces are distributed to relevant physical valuables through a smoothed profile. Several two-dimensional and three-dimensional numerical investigations including those determining the Couette flows, flow around a circular cylinder, transition layer on a wetting boundary, and dynamic behavior of a droplet on a flat or curved plate demonstrate the efficiency of the present method in calculating the contact angle of a droplet on curved surfaces with wall impermeability. The present model provides a simple algorithm to compute the surface normal vector and contact line dynamics on an arbitrarily shaped boundary by using a smoothed-profile. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Hydrophobic surfaces have attracted considerable attention because of significant scientific applications such as water repellency, lubricity, and self-cleaning and antifouling properties. Wettability of solid surfaces is governed by the chemical composition and the geometric structure of a surface. For example, the superhydrophobicity of lotus leaves principally results from the presence of binary structures at both micrometer and nanometer scales and the low energy wax-like materials on the surfaces [2]. A previous study developed a phase-field model that couples Cahn–Hilliard and Navier–Stokes equations, and a good agreement was observed between the numerical results and the experimental data with respect to the spreading of a liquid droplet on a smoothed and chemically homogeneous surface [3,4]. Physical properties continuously vary in the interface between two phases in a phase-field model based on free-energy theory [5]. With respect to the study of the contact-line movement, the phase-field model sets the wetting potential that corresponds to a contact angle on the boundary, and satisfies the following three boundary ∗

Corresponding author. E-mail address: [email protected] (T. Seta).

http://dx.doi.org/10.1016/j.compfluid.2017.09.012 0045-7930/© 2017 Elsevier Ltd. All rights reserved.

conditions, namely non-slip condition, non-penetration condition, and no mass-flux condition, ∂ μφ /∂ n = 0. Specifically, μφ denotes the chemical potential, and ∂ /∂ n denotes the normal derivative. The lattice Boltzmann method (LBM) computes fluid dynamics by the evolution of distribution functions on the discrete lattice [6]. The LBM satisfies mass and momentum conservations and is suitable for computing multi-phase and multi-component flows. A color-field model and an inter-particle-potential model were proposed to easily capture irregular topological changes, interface disintegration, and coalescence in a multi-phase flow [7,8]. These models were useful in simulating contact-line dynamics by controlling the concentration on a wall [9–11]. The phase-filed model in combination with LBM was extended to consider the contact angle of a liquid placed on a flat surface by using the relation between the normal derivative of the order parameter and the wetting potential [12–14]. The difficulty in treating large density difference in the lattice Boltzmann scheme was solved such that the spreading phenomena of a droplet on a flat hydrophobic surface could be computed quantitatively [15–19]. The LBMs succeeded in verifying the effect of roughness on the wettability of a wall through the numerical simulation of droplet dynamics on a square pillar microstructure [20–22]. A bounce back on the node or halfway bounce back schemes are used in the fore-mentioned studies,

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

and thus, a smooth curved boundary was approximated by a series of staircases. Filippova proposed an interpolation method to compute the fluid flow around an arbitrarily shaped boundary in a uniform regular Cartesian mesh in the LBM [23]. This was followed by studies that enhanced the accuracy and mass conservation by the modifications of the interpolation approach of the distribution functions between a boundary and fluid nodes [24,25]. In order to obtain the distribution function reflected from a curved wall, it is necessary to compute the distance between the wall surface and the reference nodes along each discrete velocity direction in an interpolation approach. Immersed boundary-lattice Boltzmann methods (IB-LBM) can deal with an arbitrarily complex geometry by imposing a fluidsolid interaction force on Lagrangian points [26–28]. When compared to the sharp interface scheme, the IB-LBM with a diffuseboundary approach simply obtains the fluid-solid interaction by employing a discrete delta-function to transfer quantities between the Eulerian and Lagrangian nodes [29,30]. Feng and Michaelides proposed IB-LBM based on the direct forcing method to obtain the force based on the difference between the fluid velocity and the desired velocity at the boundary without requiring a user-defined stiffness parameter [27]. Although the aforementioned explicit IBLBMs succeeded in simulating the sedimentation of particles in an enclosure, they led to a spurious fluid mass exchange between the interior and exterior of a solid domain due to the inexactness of the non-slip boundary conditions [26–28]. Wu proposed an implicit velocity correction-based IB-LBM to provide force density and velocity correction by solving a system of equations using an inverse of matrix [31]. The IB-LBMs were proposed to accurately satisfy non-slip boundary conditions through an iterative correction of body force on the Lagrangian and Eulerian nodes [32–34]. The non-slip boundary condition was accurately enforced, and therefore, the implicit correction method succeeded in preventing the penetration of a fluid into a solid surface observed in the conventional IB-LBMs and contributed to reducing the velocity slip [35]. Most existing IB-LBMs focused on Dirichlet boundary conditions. In order to investigate hydrodynamics on the complicated binary structured surface, Shao developed an immersed boundaryphase field-lattice Boltzmann method that embodied Neumann boundary conditions [36]. The primary concept of Shao’s method involved interpreting the Neumann boundary conditions as a contribution of flux from the surface to the relevant physical variables in a control volume. The main feature of Shao’s approach corresponded to the IBM that corrected the temperature on the Lagrangian points to implement Neumann (heat flux) conditions [37]. Information on the normal direction on a solid surface is necessary to examine the wettability of a solid boundary, and thus, it is necessary to compute normal directions on the Lagrangian points from the positions of the neighboring points. In the IB-LBM, it is necessary to equidistantly locate the Lagrangian points on the solid surface to form the even distribution of a solid-fluid interaction force. An adequate mesh generator and an algorithm are necessary to calculate the normal vector to the surface to simulate contact line dynamics with a lotus-leaf-like complicated geometry by the IBLBM. With respect to the two-way coupling of an incompressible fluid with rigid bodies of an arbitrary shape, a different approach based on the smoothed-profile method (SPM) was introduced into the LBM field [38,39]. The SPM defines a spatial indicator field to yield the boundary force in Navier–Stokes equations without the Lagrangian points and without the boundary-fitted coordinate system. The indicator profile smoothly transitions between the fluid and solid regions and is a function of the distance from the solid surface [40]. It is expected that the SP-LBM can compute the normal vector to the surface from the Laplacian of the indicator profile in a manner similar to the level set method [41]. In the present study, the SP-LBM is developed to investigate the effect of surface

65

wettability on droplet dynamics under the Dirichlet boundary conditions for velocity and under the Neumann boundary conditions for a chemical potential and for an order parameter (the phase field). The remainder of the paper is organized as follows. Section 2 describes the SP-LBMs for incompressible single-phase and two-phase flows in detail. The section includes the techniques for reducing fluid penetration and implementing Neumann boundary conditions in SP-LBMs in two subsections. The simulation procedures are summarized. Section 3 provides numerical experiments to demonstrate the accuracy and utility of the proposed method. The cylindrical Couette flow, the fluid flow around a circular cylinder, and the axial Couette flow are calculated to validate the non-slip and non-penetration boundary conditions in two and three dimensions. In order to validate the no mass-flux and the wetting boundary conditions, we compute a transition layer on a cylinder and on a sphere, derive the equilibrium state of the contact angle of a droplet on a flat or curved plate, and predict the contact-line motion along a curved surface. Finally, concluding remarks are presented in Section 4. 2. Smoothed profile-lattice Boltzmann method 2.1. Non-slip and non-penetration conditions 2.1.1. Lattice Boltzmann method for single-phase flow The lattice Boltzmann method for the incompressible Navier– Stokes equations uses the following kinetic equations for the distribution function fα :

fα (x, t ) − fα(eq ) (x, t ) f˜α (x, t ) = fα (x, t ) − + δt Fα (x, t ),

(1)

fα (x + cα δt , t + δt ) = f˜α (x, t ),

(2)

τn

where f˜α denotes post-collision value, δ t denotes time step, cα denotes discrete velocity, and Fα denotes a forcing term. The fluid density and velocity are conventionally expressed in terms of the distribution function as follows:

ρ=

 α

fα ,

u=

1

ρ

α

f α cα .

(3)

The D2Q9 and D3Q19 models are used in the present study. The D2Q9 model defines the discrete velocities as follows:

⎧ ⎪ ⎨ ( 0, 0 ) c (±1, 0 ), c (0, ±1 ) cα = ⎪c(±1, ±1 ) ⎩

α=0 α =1−4 , α =5−8

(4)

where c denotes lattice velocity magnitude. With respect to the D3Q19 model, the discrete velocity set can be expressed as follows:

⎧ ⎪ ⎨ ( 0, 0, 0 ) c (±1, 0, 0 ), c (0, ±1, 0 ), c (0, 0, ±1 ) cα = ⎪ ⎩c(±1, ±1, 0 ), c(±1, 0, ±1 ), c(0, ±1, ±1 )

α=0 α =1−6 . (5) α = 7 − 18

The equilibrium distribution function is given as follows:



fα(eq ) = ωα ρ 1 +



3u · u 3 cα · u 9 ( cα · u )2 + − , c2 2c 4 2c 2

(6)

where ωα denotes the weight coefficients. With respect to the D2Q9 model, ω0 = 4/9, ω1−4 = 1/9, and ω5−8 = 1/36. With respect to the D3Q19 model, ω0 = 1/3, ω1−6 = 1/18, and ω7−18 = 1/36. Conventionally, the forcing term is defined as follows:

Fα = ωα ρ

3 cα · G , c2

(7)

66

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

where Np denotes the number of the solid regions. The velocity field of the solid region ϕ (1) up is obtained as follows:

ϕ ( 1 ) ( x, t )u p ( x, t ) =

Np 

{Vi (t ) + i (t ) × ri (t )}ϕi(1) (x, t ),

(13)

i=1

where Vi and i denote the translational velocity and the angular velocity of the ith solid particle, respectively. ri denotes distance vector from the rotation point of the ith particle. Jafari defines the fluid-solid interaction force G that acts on the fluid nodes as follows [38]:

u p ( x, t ) − u ( x, t )

G ( x, t ) = ϕ ( 1 ) ( x, t )

Fig. 1. Smoothed profiles for the non-slip Dirichlet boundary conditions, ϕ (1) (solid line), and for the Neumann boundary conditions, ϕ (2) (dashed line).

∇ · u = 0, 1 ∂u + ∇ · (uu ) = − ∇ p0 + ν∇ · [∇ u + (∇ u )T ] + G, ∂t ρ

ϕi ( x, t ) =

⎩

1 2

sin π diξ(x,t ) i

1



+1



di ≤ −ξi /2

| d i | < ξi / 2 ,

(11)

d i ≥ ξi / 2

where di denotes singed distance to the ith solid surface with a negative value in the fluid region and a positive value in the solid region. ξ i denotes interface thickness parameter. With respect to a circular cylinder or a sphere with a radius R, di = R − |x − Ri |. Here, Ri denotes center position of a cylinder or of a sphere. A concentration field ϕ (1) (x, t) is introduced by the summation of individual concentration functions as follows: Np  ϕ ( 1 ) ( x, t ) = ϕi(1) (x, t ), i=1



t+δt

t

Ti (x, t )dt =

 ∀p

 ∀p

ρϕ (1) (x, t )(u(x, t ) − u p (x, t ))dx,

(15)

ρϕ (1) (x, t )ri (t ) × (u(x, t ) − u p (x, t ))dx.

(9)

(10)



t

Fi (x, t )dt =

In order to consider the discrete lattice effect and the contributions of the body force to the momentum flux, Guo redefines the fluid velocity and the forcing term in Eq. (1). This is expressed as follows:

2.1.2. Smooched profile method for the Dirichlet boundary conditions A smoothed profile method is developed to implement nonslip boundary conditions to simulate particulate flows by applying an external boundary force at the solid-fluid interface to the LBM [38]. The boundary force follows the indicator function ϕi(1 ) that changes smoothly within this finite-thickness interface as shown in Fig. 1. The function ϕi(1 ) corresponding to the ith solid region is defined to differentiate between the fluid and solid regions. The following profile is mostly used in the SPM [38,40]: (1 )

t+δt

(16)

The Eulerian mesh size δ x satisfies δx = cδt . In the study, δ x , c, and δ t correspond to unity.

⎧ ⎨0 

(14)

(8)

where p0 denotes pressure√calculated by p0 = cs2 ρ . The speed of sound cs corresponds to c/ 3. The kinematic viscosity ν is determined as follows:

(2τn − 1 )δx2 ν= . 6δt

.

Eq. (14) uses the concept of the direct forcing method in the immersed boundary method proposed by Feng and Michaelides [27]. The force Fi and torque Ti acting on the solid particle are derived from the conservation of momentum between the solid and the fluid. The time integral of the force exerted by the surrounding fluid during the time interval δ t corresponds to the volume integral of the hydrodynamic force over the solid volume ∀p . From the momentum conservation, the following expression is obtained:

 where G denotes a fluid-solid interaction force that is calculated by the smoothed profile method. The Chapman–Enskog expansion is applied to Eqs. (1), (2), (6) and (7) to yield the following macroscopic equations:

δt

(12)

ρ=



u∗ =

α

fα ,

1

ρ

α

(17) u = u∗ + δt

f α cα ,



Fα = ωα ρ 1 −

1 2 τn

 3(c − u ) α c2

+

G , 2

(18)



9 ( cα · u ) cα · G. c4

(19)

In order to ensure that the boundary velocity u is equal to the desired velocity up based on Eq. (18), Kang iteratively corrected the body force and the fluid velocity at the boundary using the following equations [33]:

Gk =

2 ( u p − uk )

δt

,

uk+1 = uk + δt

Gk , 2

(20)

where k denotes the number of iteration steps. Wu proposed the IB-LBM known as the implicit correction method to compute the body force G through simultaneous equations based on Eq. (18) [31]. In the simulation of the dynamic behavior of a droplet on a solid surface in the IB-LBM, Shao used the implicit correction method to satisfy the non-slip and non-penetration boundary conditions. In the numerical calculation of multi-phase fluid flow in the SP-LBM, it is necessary to establish the algorithm to prevent fluid leakage. Following the method proposed by Kang and Hassan [33], u = u p is substituted into Eq. (18) to yield the body force acting on the fluid in the SP-LBM as follows:

G ( x, t ) = ϕ ( 1 ) ( x, t )

2(u p (x, t ) − u∗ (x, t ))

δt

.

(21)

Substituting Eq. (21) into Eq. (18) results in the following expression:

u(x, t ) = ϕ (1) (x, t )u p (x, t ) + (1 − ϕ (1) (x, t ))u∗ (x, t ).

(22)

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

67

Additionally, ϕ (1 ) = 0 is in the fluid region, and thus the body force G is equal to zero. (See Fig. 1). When G = 0 in Eq. (18), u∗ is equal to the fluid velocity uf . Eq. (22) is equivalent to the definition of the total velocity field as follows:

where n, P, φ , M, and μφ correspond to the mean density of two fluids, the dynamic pressure, the order parameter, the mobility, and the chemical potential, respectively [43]. ∇ P is related to the surface tension force φ∇μφ as follows:

u = ϕ ( 1 ) u p + ( 1 − ϕ ( 1 ) )u f ,

∇ P = φ∇ μφ + ∇ p0 ,

(23)

as shown in a previous study [40]. The divergence of both sides of Eq. (22) is considered, and an incompressibility constraint of the fluid velocity uf is used to yield the following equation:

∇ · u = ∇ϕ (1) · (u p − u f ).

(24)

(28)

where p0 denotes ideal gas pressure that is defined as p0 = cs2 n. A double-well potential, (φ ) = β (φ 2 − φ #2 )2 , with two equilibrium states φ # and −φ # yields the chemical potential as follows:

Nakayama proved that the divergence-free condition, ∇ · u = 0, yields the solid-fluid impermeability condition, u p = u f , and vice versa in Eq. (24). The conventional SP-LBM proposed by Jafari does not follow a smooth combination between the particle velocity field, up , and the fluid velocity field, uf , as shown in Eq. (23) [38]. In the definition of the force and the velocity in the present study as shown in Eqs. (21) and (22), the total velocity changes smoothly from the particle velocity, up , to the fluid velocity, uf , in the fluidsolid interface (0 < ϕ (1) < 1). With respect to Eq. (24), the present SP-LBM is assumed to satisfy the impermeability condition within the same degree of the compressible error as that of the LBM.

μφ = 4β (φ 3 − φ #2 φ ) − κ∇ 2 φ ,

2.1.3. Simulation procedure The solution procedure of the present smoothed profile-lattice Boltzmann method for the single phase flow is performed as follows:

φ = φ # tanh

(1) Give the initial values including ρ , u, ϕi(1 ) , Vi , and i . Initially, the body force G is zero. (2) Compute the density field ϕ (1) and velocity field up using Eqs. (12) and (13), respectively. (eq ) (3) Compute fα that corresponds to the initial values by using (eq ) Eq. (6), and obtain the initial distribution function as fα = fα . (4) Compute the forcing term Fα using Eq. (19). Add the SRT collision term and the effect of the body force G to the distribution function fα by using Eq. (1). (5) f˜α streams following Eq. (2). (6) Obtain ρ , u∗ using Eqs. (17) and (18). (7) Evaluate the force density G using Eq. (21), and correct the fluid velocity u using u∗ and G in Eq. (18). (eq ) (8) Compute the equilibrium distribution function fα using Eq. (6), and return to step (4).

(29)

where β denotes a parameter that modulates interaction energy between the two phases, and κ denotes a coefficient related to the surface tension. In the study, φ # equals to unity. The phase-field method considers the fluid-fluid interface width w that is given by the following expression:



w=

2κ /β

φ#

.

(30)

The order parameter profile along the normal direction is given as follows:

2η

w

,

(31)

where η denotes a coordinate that is perpendicular to the fluidfluid interface. Zheng evaluated the surface tension, σ , for a flat plate as follows [43]:



4 2κβ

σ=

3

φ #3 .

(32)

Eqs. (30) and (32) yield the equations for the coefficients β and κ that are used in the numerical experiments discussed in a later section.

3σ φ #4 , 4w

β=

κ=

3σ w . 8φ #2

(33)

Eqs. (25) and (26) are solved by using the lattice Boltzmann equations (1) and (2) for fα with the modified equilibrium distribution function and forcing term. The equilibrium distribution function for the two-phase flows is given as follows:





3u · u 3 cα · u 9 ( cα · u )2 + − + ωα Aα , c2 2c 4 2c 2

2.2. No mass-flux and the wetting boundary conditions

fα(eq ) = ωα n 1 +

2.2.1. Lattice Boltzmann method for gas-liquid two-phase flows Several novel phase-field lattice Boltzmann models were proposed to compute the viscosity contrast and high density difference without numerical stability and without artificial oscillations at a fluid-fluid interface [15,16,18,19,42]. These phase-field lattice Boltzmann methods successfully solved representative problems including head-on collisions of binary micro-droplets and droplet splashing on a thin film for various Weber and Reynolds numbers. In the present study, a phase-field model developed by Swift et al. [12] is used for the small density contrast fluids to compare the accuracy of present method with that of the IB-LBM proposed by Shao et al. [36]. The following governing equations are computed by using the lattice Boltzmann equations.

where A0 = − 15 4 φμφ and A1−8 = 3φμφ for the D2Q9 model. With respect to the D3Q19 model, A0 = −6φμφ and A1−19 = 3φμφ . In order to derive the dynamic pressure, P, defined by Eq. (28) in the Navier-Stokes Eq. (26), the forcing term includes μφ ∇φ as follows:

∂n + ∇ · ( nu ) = 0, ∂t ∂ nu + ∇ · (nuu ) = −∇ P + μ∇ 2 u + nG, ∂t ∂φ + u · ∇ φ = M ∇ 2 μφ , ∂t

(25)

Fα =

(34)



 3(c − u ) 9(c · u )  1 α α ωα 1 − + cα · (μφ ∇φ + nG ), 2 τn c2 c4

(35) where G denotes the interaction force calculated by Eq. (21). The macroscopic valuables are evaluated as follows:

n=

 α

fα ,

nu =

 α

fα cα + δt

μφ ∇φ + nG 2

.

(36)

(26)

In the phase field model, the fluid-fluid interface is governed by the Cahn–Hilliard Eq. (27). The lattice Boltzmann equation for Eq. (27) assumed the following form:

(27)

gα (x + cα δt , t + δt ) = gα (x, t ) −

gα (x, t ) − g(αeq ) (x, t )

τφ

,

(37)

68

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

where τ φ denotes the relaxation time that is selected as τφ = τn by following the method in a previous study [36]. The order parameter, φ , is calculated as the local sum with respect to the discrete-velocity directions according to the following expression:

φ=

 α

gα .

(38)

The mobility, M, is determined as follows:

(2τφ − 1 )δt

M=

,

4

(39)

where  denotes a coefficient to alter the mobility. With respect to D2Q9, the equilibrium distribution function is given as follows:

5  μφ , 3c 2

g(0eq ) = φ −

(40)





3ωα  μφ + φ c α · u . c2

) g(1eq = −8

(41)

In order to compute the three-dimensional flow fields, the D3Q7 model is used for simplicity and high computational efficiency as follows:

3  μφ , c2

g(0eq ) = φ −

(42)





1  μφ + φ c α · u . 2c 2

) g(1eq = −6

(43)

In the D3Q7 model, the discrete velocity directions are defined as follows:

cα =

( 0, 0, 0 ) c (±1, 0, 0 ), c (0, ±1, 0 ), c (0, 0, ±1 )

α=0 . α =1−6

(44)

2.2.2. Smoothed profile method for the Neumann boundary conditions In order to control the wetting behavior of a solid surface, Briant proposed Neumann boundary conditions enforced on a boundary [14]. This is expressed as follows:

 ˜ ∂φ   =− , ∂ n s κ

∂μφ  = 0, ∂ n s

(45)

˜ denotes the wetting powhere s denotes the solid surface, and  ∂μ

tential. Additionally, ∂ nφ |s = 0 corresponds to a no mass-flux condition of the chemical potential. The following non-dimensional wetting potential  determines the equilibrium contact angle θ as follows:

= cos θ

(φ )

˜ 

# 2

2κβ

,

(46)





1 2 2 = (1 + ) 3 − (1 − ) 3 . 2

(47)

In order to satisfy the Neumann boundary conditions shown in Eq. (45) in the IB-LBM scheme, Shao proposed a correction method for the order parameter and for the chemical potential by the surface fluxes that are defined as follows:

 ∂μ  φ  δφ = −2M ∂n 

∂μ∗φ  − , ∂n Specified

δμφ = 2κ

 ∂φ   ∂n 

Specified

−

∂φ ∗  , ∂n

(48)

(49)

where φ ∗ and μ∗φ denote pre-collected values computed by

Eqs. (29) and (38). The IB-LBM interpolates the Eulerian values of

φ ∗ and μ∗φ to the Lagrangian values, and computes the fluxes δμφ , and δφ on the Lagrangian points using Eqs. (48) and (49). After spreading δφ and μφ to the Eulerian points, the IB-LBM corrects φ and μφ as given by the following equations:

φ (x, t ) = φ ∗ (x, t ) + δφ (x, t ),

(50)

μφ (x, t ) = μ∗φ (x, t ) + δμφ (x, t ).

(51)

The IB-LBM succeeds in adequately controlling the surface wettability of the curved surface by treating the Neumann boundary conditions as the surface fluxes. In the present study, the procedure shown in Eqs. (48) through (51) is introduced into the SP-LBM to impose the Neumann boundary conditions. In order to distribute the effect of the fluxes δφ and μφ to the relevant physical values in the scheme of the SP-LBM, the following density profile is defined for the fluid-solid interface (see Fig. 1) as follows:

ϕ (2) (x, t ) = 1 − tanh2

2d ( x, t )

i

ξi

.

(52)

The same values are used for the distance di from the boundary surface, and the width ξ i of the boundary layer in Eqs. (11) and (52). The present SP-LBM does not follow the interpolation procedure for Eqs. (50) and (51) and directly corrects φ and μφ on the Eulerian nodes as follows:

φ (x, t ) = φ ∗ (x, t ) + ϕ (2) (x, t )δφ (x, t ),

(53)

μφ (x, t ) = μ∗φ (x, t ) + ϕ (2) (x, t )δμφ (x, t ).

(54)

As shown in Fig. 1, the contributions of the fluxes, ϕ (2) δφ and ϕ (2) δμφ , are the strongest in the middle of the boundary layer, and decrease in a manner proportional to the distance from the solid surface in Eqs. (53) and (54). In the smoothed profile method, Eqs. (48) and (49) are calculated on the Eulerian nodes. The gradients of φ ∗ are obtained on the Eulerian nodes near the solid surface by using the second-order central difference approximation as follows: ∂ φ ∗ /∂ x ≈ (φi∗+1, j − φi∗−1, j )/2δx . The directional derivative in the direction of the normal vector is given by the inner prod∗ ∂φ ∗ ∂φ ∗ uct as follows: ∂φ ∂ n = nx ∂ x + ny ∂ y . Although the normal vector

(nx , ny ) can be computed by ∇ ϕ (1) /|∇ ϕ (1) | in the smoothed profile method, this study targets only a flat plate, a cylinder, and a sphere in which the normal direction is easily calculated by the boundary geometry without considering the influence of the computation accuracy of ∇ϕ (1) . For example, x − Ri denotes surface normal vector of a sphere with a center position Ri . In the following computations, the second-order central difference approximation is utilized to compute ∇μφ , ∇φ , and ∇ 2 φ . 2.2.3. Simulation procedure The solution procedure of the present SP-LBM for the twophase flows with Neumann boundary conditions is performed as follows: (1) Give the initial values including n, u, φ , μφ , σ , ϕi(1 ) , ϕi(2 ) , Vi ,

i ,

∂μφ ∂φ ∂ n |Specified , and ∂ n |Specified . Initially, the body force G is

zero. (2) Compute the density field ϕ (1) and velocity field up using Eqs. (12) and (13), respectively. (eq ) (eq ) (3) Compute fα and gα corresponding to the initial values using Eqs. (34), and (40) through (43), respectively. Give the initial (eq ) (eq ) distribution functions as fα = fα and gα = gα . (4) Compute the forcing term Fα using Eq. (35). fα and gα evolve following Eqs. (1), (2), and (37), respectively. (5) Obtain n, u∗ using Eq. (36).

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

(6) Evaluate the force density G using Eq. (21), and update the velocity u using Eq. (36). (7) Obtain φ ∗ using Eq. (38), and compute μ∗φ using Eq. (29).

(8) Evaluate the fluxes δφ and δμφ using Eqs. (48) and (49). Correct φ and μφ using φ ∗ , μ∗φ , δφ , and δμφ in Eqs. (53) and (54), respectively. (eq ) (eq ) (9) Compute the distribution function fα and gα using Eqs. (34), and (40) through (43), respectively and return to step (4). The fluid-solid interaction force and the velocity definition are given by Eqs. (14) and (18) in the SP-LBM proposed by Hu et al. [39]. The combination of Eq. (14) with Eq. (18) does not satisfy the non-penetration conditions in Eq. (24). In the present work, we theoretically and numerically demonstrate that the fluid-solid interaction force should be given by Eq. (21) to completely implement the non-penetration condition in the framework of the SP-LBM. In Section 3.2, the computation of flow around a circular cylinder numerically reveals the necessity of Eq. (21) for the elimination of fluid leakage into the solid region. The novelty of the present method is the introduction of an indicator profile ϕ (2) in the SP-LBM, which can impose the Neumann boundary condition on a curved surface without using Lagrangian points. Through the implementation of the Neumann boundary condition, ϕ (2) allows the SP-LBM to satisfy no mass-flux and wetting boundary conditions on surfaces of arbitrary shape. 3. Results and discussion 3.1. Cylindrical Couette flow This is followed by examining the accuracy of the present SPLBM, that can enforce the non-slip boundary conditions in the calculation of the two-dimensional cylindrical Couette flow. The exact solution for the tangential velocity of the cylindrical Couette flow is given as follows:

uˆθ (R ) = udθ

R/Ro − Ro/R , Ri /Ro − Ro/Ri

(55)

where R denotes the radial coordinate. The two rings are placed at the centre of the simulation domain covered by 400 × 400 grid points. The radii of the outer ring, Ro , and the inner ring, Ri , correspond to 140δ x and 90δ x , respectively. The outer ring is at rest, and the tangential velocity of the inner cylinder, udθ , corresponds to 0.01c. The exact solution for the torque Tˆ acting on the inner cylinder is given as follows:

Tˆ =

4π ν udθ Ro

Ri /Ro − Ro/Ri

.

(56)

Fig. 2(a) shows the contour lines of the tangential velocity that are calculated using the SP-LBMs and the IB-LBMs at τn = 0.75. The solid-fluid interface width ξ i is equal to 1.0δ x to obtain the best accuracy for the non-moving solid boundary [40]. In order to estimate the accuracy of numerical methods, the cylindrical Couette flow is computed by two IB-LBMs, the direct forcing method proposed by Feng and Michaelides [26], and the implicit correction method proposed by Wu and Shu [31]. Additionally, 565 Lagrangian points are set on the inner cylinder and 879 points are set on the outer cylinder, which yields a discrete length of s ≈ δ x in the IB-LBMs. Fig. 2(a-1) shows the numerical results calculated by the direct forcing method and the implicit correction model. The left-hand and right-hand sides of Fig. 2(a-2) show the velocity contours calculated by the conventional SP-LBM proposed by Jafari et al. [38] and by the SP-LBM used in the present study, respectively. Minute fluctuations are observed around the surface of the cylinder in the numerical results by the conventional SP-LBM in

69

Figs. 2(a-2). The present SP-LBM eliminates this noise in the solid regions. Fig. 2(b) shows the profiles of the tangential velocity at the horizontal centerline of the computational domain. The velocity profile for all the methods agree well with the analytical solution in the fluid domain in Fig. 2(b). The relative errors are defined as follows:





 ErrorT =

(uθ (x ) − uˆθ (x ))2 , 2 ˆ θ (x ) x∈fluid u

x∈fluid

Erroru =

(T − Tˆ )2 Tˆ 2

,

(57)

(58)

where the fluid represents the region between the outer and inner rings as shown in Fig. 2. The exact solutions uˆθ and Tˆ are given by Eqs. (55) and (56), respectively. Fig. 3 shows the relation between the accuracy of the numerical methods and the number of grid points on the side of the square computational domain in the simulation of the cylindrical Couette flows. The present SP-LBM exhibits approximately first-order accuracy in space for both the velocity and the torque in a manner similar to the IB-LBMs. We investigate the computational effort involved in running the proposed SP-LBM. The calculations were performed on a 3.1 GHz Linux workstation. When the discrete length for each Lagrangian point on the cylinder, s, is approximately equal to 2δ x , and the number of grid points is 300 × 300, the implicit velocity correction-based IB-LBM requires 315.643 [s] and 15,820 iterations to satisfy the convergence criteria, max ||uk+1 − uk || < 10−6 . It takes 43.648 [s] and 2037 iterations for the proposed SP-LBM to get the convergent solution with 300 × 300 grid points. Fig. 4(a) shows the computational time under different grid size required to obtain the convergent solution. The SP-LBM is about eight times faster than the IB-LBM for any grid size to converge a calculation. Fig. 4(b) shows the temporal evolution of the norm given by max ||uk+1 − uk || with 300 × 300 grid points. Since the IB-LBM takes a long time to obtain a convergent solution due to the slight oscillation near the convergence state in Fig. 4(b), we measure the computational time after 20,0 0 0 iterations. Table 1 indicates the computational times needed by the proposed SP-LBM and the implicit velocity correction-based IB-LBM, with an underline indicating when the IB-LBM is faster than the SP-LBM. Since the circumference of the cylinder is given by 2π R, the number of Lagrangian points increases with an increase in radius of the cylinders, R. The radius of the cylinder, R, is proportional to the number of the grid points, N. The computational cost of the IBM procedures increases with order O(N). While the IB-LBM obtains the fluid-solid interaction force on the Lagrangian points, the SP-LBM computes the force on the Eulerian points inside the solid region. Since the area of the region inside the cylinder is given by π R2 , the computational time of the SPM procedures is proportional to N2 . The SPM procedures constitute a large amount of all computational processes in the SPLBM for the large number of Eulerian points. Although the SP-LBM does not always produce the computational effect for any grid size, the diffusive nature of its fluid-solid interface profile allows a significant reduction in the computer time required to obtain convergent solutions. Since the elements of matrix A are based on the nearby interaction between two Lagrangian points in the square area 2δ x × 2δ x , the coefficient matrix A is sparse in the implicit velocity correction-based IBM. We must implement a high-performance solution procedure and optimization for the zero and nonzero elements in the sparse matrix of the IB-LBM [44]. To compute the IBM code in parallel, a given processor usually deals with the Lagrangian points. This strategy raises some issues including loadbalance due to the unequal distribution of Lagrangian points, and communication overhead produced by shared processing of the

70

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

Fig. 2. Tangential velocity profiles in the cylindrical Couette flows. Table 1 Comparison of computational time needed by SP-LBM and IB-LBM. An underline indicates that the IB-LBM is faster than the SP-LBM. Grid size N × N

502

1002

2002

3002

4002

SP-LBM (Present)[s] IB-LBM (s ≈ 2δ x )[s] IB-LBM (s ≈ δ x )[s] IB-LBM (s ≈ 2δ x /3)[s]

7.863 10.743 13.960 17.716

32.493 37.955 45.324 54.586

131.062 140.641 163.058 196.775

444.201 366.486 471.412 514.245

947.654 835.173 915.169 1327.930

points among processes [45]. Since the SP-LBM is a simple and local computational algorithm on the Cartesian grid, it is more suitable than the IB-LBM for parallelization using the domain decomposition approach. The advantages of the proposed SP-LBM over the IB-LBM are fast computation without numerical oscillation near the convergence state, the simplicity of the algorithm without the cost of sparse matrix computation, and efficient and easy parallel implementation.

3.2. Flow past a circular cylinder The second test corresponds to a calculation of the fluid flow past a circular cylinder at Re = 20 and 40. Numerous solutions reported by extant theoretical, experimental, and numerical approaches are compared with the results of the present study. The relaxation time, τ n , corresponds to 0.65. The diameter of the cylinder is considered as D = 100δx . The characteristic free-stream velocity u0 is calculated by u0 = Reν /D. The cylinder is located at

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

Fig. 3. Relative error with respect to the number of grid points in the calculation of the cylindrical Couette flows.

(16D, 20D) in the computational domain with 40D × 40D. The fluid density is given by ρ = 1.0. All the boundary values of fα are (eq ) given by the equilibrium distribution functions fα for ρ = 1.0 and u0 = Reν /D. Fig. 5 shows the streamlines in the half computational domain. The direct forcing IB-LBM allows for fluid leakage into the solid region in Fig. 5(a). Fig. 5(b) shows a small penetration into the cylinder in the numerical results as computed by the conventional SP-LBMs proposed by Jafari et al. [38] and Hu et al. [39]. A previous study [35] demonstrated that the fluid velocity, u, should not equal the solid velocity, up , both in the direct forcing IB-LBM and in the conventional SP-LBMs in which the fluid-solid interactions are given by the difference between u and up , as shown in Eq. (14). In order to accurately satisfy the non-slip boundary conditions, u = u p , the implicit correction IB-LBM and the present SPLBM define the body force based on the definition of the fluid velocity as shown in Eq. (18). Although the SP-LBM proposed by Hu is based on the fluid velocity and the forcing term defined by Guo as shown in Eqs. (18) and (19), it does not satisfy the non-

71

Fig. 4. Computational time under different grid size and the temporal evolution of the norm to obtain the convergent solutions in the calculation of the cylindrical Couette flows.

penetration condition due to use of the fluid-solid interaction force given by Eq. (14). Using Eq. (21), the method in the present study successfully reduces the fluid leakage in Fig. 5(c). Since the flow velocity is completely equal to zero in the solid region in these numerical results as computed by the proposed SP-LBM, the streamline is constant inside the circular cylinder in Fig. 5(c). The results in Fig. 5 demonstrate that the combination of Eqs. (18), (19), and (21) are necessary to satisfy the non-slip and non-penetration condition given by Eqs. (23) and (24). The drag and pressure coefficients are defined as follows:

Cd =

FD , (1/2 )ρ u20 D

Cp =

ps − p0 , (1/2 )ρ u20

(59)

where ps denotes pressure on the surface of a circular cylinder. The free stream fluid pressure p0 is given by p0 = c2 ρ /3. The drag force Fd is calculated as follows:



Fd = −

∀p

Gx (x, t )dx,

(60)

72

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

Fig. 5. Streamlines for a flow past a circular cylinder at Re = 20 and 40. The thick line indicates the zero streamline.

where Gx denotes x component of the force density G. Since the SP-LBM uses only Eulerian points, we utilize a spatial interpolation technique to obtain a computed value on the surface of a wall. The pressure at the surface of a cylinder can be obtained through the interpolation of the pressure at the neighboring fluid nodes x using a discrete delta function D(x) as follows:

ps ( xs ) =



p(x )D(x − xs )δx2 ,

(61)

f

dinate system and Nieuwstadt and Keller used a semi-analytical method with a finite Fourier series [46,47]. He and Doolen solved the Navier–Stokes equations by extending the lattice Boltzmann method to a cylindrical coordinate system [48]. The findings reveal good agreement between the results of the present method and the results of previous studies [46–48]. At a low Reynolds number, the present SP-LBM appropriately computes the drag force and pressure for the steady fluid flows around a cylindrical cylinder by satisfying the impenetration condition.

where xs denotes the Lagrangian points on the circular cylinder surface. With respect to a two-dimensional coordinate system x = (x, y ), D(x) is expressed as follows:

δ (r ) =

 1 4

0

  1 + cos π|2r|

|r | ≤ 2 , |r | > 2

3.3. Axial circular Couette flow

(62)

and

D ( x − xs ) = δ ( x − xs )δ ( y − ys ).

(63)

Table 2 shows the drag coefficient Cd , the wake length L normalized with respect to the radius of the cylinder, the separation angle θ s , and the pressure coefficient Cp in the calculation shown in Fig. 5. In order to solve stream-function vorticity equations, Dennis and Chang used a finite difference method in the polar coor-

The axial circular Couette flow is calculated to examine the accuracy of the present SP-LBM for three-dimensional flows. With respect to the axial circular Couette flow, an inner circular cylinder moves to the z-direction with a constant velocity, udz . The outer circular cylinder remains at rest. The periodic boundary conditions are applied in all directions. The simulation domain is covered by 200 × 200 × 5 grid points. The radius of the outer ring, Ro = 70δx , and the radius of the inner ring, Ri = 45δx . udz correspond to 0.01c. The relaxation time is set as τn = 0.75. The exact solution for the

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

73

Table 2 Comparison of drag coefficient Cd , wake length 2L/D, separation angle θ s , pressure coefficient at the front stagnation point Cp (π ), and pressure coefficient at the rear stagnation point Cp (π ) as calculated by the present SP-LBM and in previous studies. Re

References

Cd

2L/D

θs

Cp (π )

−C p (0 )

20

Dennis and Chang [46] Nieuwstadt and Keller [47] He and Doolen [48] IB-LBM (Implicit Correction) [31] SP-LBM (Jafari et al.) [38] SP-LBM (Hu et al.) [39] SP-LBM (Present) Dennis and Chang [46] Nieuwstadt and Keller [47] He and Doolen [48] IB-LBM (Implicit Correction) [31] SP-LBM (Jafari et al.) [38] SP-LBM (Hu et al.) [39] SP-LBM (Present)

2.045 2.053 2.152 2.103 2.106 2.106 2.110 1.522 1.550 1.499 1.566 1.569 1.569 1.572

1.88 1.786 1.842 1.900 1.900 1.900 1.920 4.69 4.357 4.490 4.600 4.640 4.640 4.660

43.7 43.37 42.96 40.89 41.85 41.78 41.83 53.8 53.34 52.84 50.70 51.43 51.36 51.43

1.269 1.274 1.233 1.220 1.266 1.266 1.266 1.144 1.117 1.133 1.105 1.158 1.158 1.167

0.589 0.582 0.563 0.580 0.587 0.587 0.585 0.509 0.554 0.487 0.505 0.509 0.509 0.508

40

steady velocity is given as follows:

uˆz (R ) = udz

ln(Ro/R ) . ln(Ro/Ri )

(64)

The exact solution for the force Gˆ z acting on the inner cylinder is given as follows:

Gˆ z = −

2π ν udz . ln(Ro/Ri )

(65)

Fig. 6(a) shows the contour lines of the z component of the velocity that is calculated by the conventional SP-LBM and by the method proposed in the present study. A significant amount of noise is observed in the solid region in the numerical results of the conventional method as shown in the left-hand side of Fig. 6(a) in a manner similar to the results in Fig. 2. In the right-hand side of Fig. 6(a), the present method reduces these types of unwanted noises by Eq. (22) in which the velocity field is given by the solid velocity when ϕ (1 ) = 1. Fig. 6(b) shows the velocity profile on the horizontal centerline for the following different interface widths: ξi = δx and ξi = 0. When ξi = 0, the indicator function, ϕi(1) , discontinuously changes between the fluid region and the solid region. As shown in Fig. 6(b), the velocity profiles for both values of ξ i agree well with the exact solution (64). The relative error in the velocity and in the force acting on the inner cylinder is calculated by the following equations:





 ErrorG =

(uz (x ) − uˆz (x ))2 , 2 ˆ z (x ) x∈fluid u

x∈fluid

Erroru =

(Gz − Gˆ z )2 Gˆ 2z

.

(66)

(67)

Fig. 7 shows the dependence of the relative error in the velocity and in the force on the number of grid points. Both the conventional and the present methods yield the first order accuracy space convergence as measured by the relative errors. When ξi = δx , the error of the present method exceeds that of the conventional method in Fig. 7. Conversely, the accuracy of the method in the present study exceeds that of the conventional method for ξi = 0. Although the proposed method exactly enforces the nonslip boundary conditions, the SP-LBM considers the solid velocity in the boundary layer for ξ i > 0. The solid velocity is considered in the fluid region through ϕ (1) , and thus the solid-fluid interface width affects the accuracy of the SP-LBM [40]. Jafari demonstrated that the nonzero interface width of the profile was necessary to suppress the velocity oscillation around a wall in a moving boundary problem [38]. Although the present study does not involve

Fig. 6. Velocity profiles as calculated by the SP-LBMs in the axial circular Couette flows. The dashed line denotes the position of the circular cylinder in (a). The solid line indicates the exact solution in (b).

74

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

where φ # = 1 in the following calculations. In the two-dimensional case, the ratio between the length of the computational domain and the cylinder radius is maintained at 3.5. In order to use the same values as that used in a previous study [36], 251 × 251 grid points and a circular cylinder of radius corresponding to 500/7δ x are used. The transition layer thickness w is equal to 9δ x , and the fluid-solid width ξ i is equal to 2δ x . In the three-dimensional case, the radius of a sphere and the length of the cubic domain are equal to 40δ x and 140δ x , respectively. The natural boundary condition, ∇φ = 0, is applied at the outer boundaries of the domain. Substituting φ = 1, φ # = 1 and ∇ 2 φ = 0 into Eq. (29) gives μφ = 0. When ∇φ = 0 in Eq. (35), Fα is equal to zero. That is, ∇φ = 0 and ∇ 2 φ = 0 are imposed on the boundaries via Eqs. (29) and (35), respectively. The distribution functions fα and gα at the outer boundaries are given by the equilibrium distribution functions over all (eq ) (eq ) the lattice directions as fα = fα and gα = gα for n = 1, u = 0, φ = 1, and μφ = 0. In all the following computations, it is assumed that σ = 0.001, τn = τφ = 0.75, and  = 10. The contact angle θ is given in the 5° to 175° range. The wetting potential  on the cylinder or the sphere surface is given by the following inverse function of Eq. (47) for 0 < θ < π [13] as follows:

π λ 

 λ  1/2  = 2sgn − θ cos 1 − cos , 2

3

3

(69)

where λ = arccos(sin2 θ ), and sgn denotes the signam function. The convergence criterion is set as follows:

max

Fig. 7. Relative error with respect to the number of grid points in the calculation of the axial circular Couette flows.

solving the moving boundary problem, the nonzero interface width ξ i is necessary to control the surface wettability by the proposed indicator profile ϕ (2) in Eq. (52). The present method can calculate the fluid velocity and the force acting on the solid walls with an acceptable precision for three-dimensional flow simulations. 3.4. Transition layers around a circular cylinder and a sphere The previous section demonstrated that the proposed SP-LBM satisfies the non-slip and non-penetration conditions. In this section, the implementation of the Neumann boundary conditions for φ and μφ in the SP-LBM are examined using the calculation of the transition layer around a circular cylinder or a sphere. A circular cylinder is located in the center of the two-dimensional square domain and is filled with a quiescent fluid with a mean density of n = 1. With respect to the three-dimensional problem, a sphere is centered in the cubic computational domain. The exact solution for the order parameter on a solid surface is given as follows [36]:

φˆ = φ # 1 ± ||,

φˆ = −φ # 1 ± ||,

(68)

|φ k+1 − φ k | < 10−6 , |φ k |

(70)

where the superscript k denotes the time level. Fig. 8 shows the order parameter distribution around a hydrophilic curved wall as calculated by the SP-LBM used in the present study. The method in the present study applies a gradient of the order parameter on a wall to force the contact angle to correspond to 60°, which corresponds to  = 0.334933. The increase in the order parameter within the transition layer is observed both in two-dimensional and three-dimensional calculations in Fig. 8. The accuracy of the numerical results is verified in the transition layer via a comparison with the theoretical solution as shown in Eq. (68). The values of the non-dimensional wetting potential,  = 0.678525 and  = −0.678525, correspond to the contact angles, θ = 5◦ and θ = 175◦ , respectively. The numerical values of the order parameter depend on the non-dimensional wetting potential as shown in Fig. 9. Fig. 9 demonstrates that the present results approximately agree with the theoretical values predicted by Eq. (68), while achieving an accuracy that is almost identical to that of the IB-LBM proposed by Shao et al. [36]. The difference between the numerical and theoretical solutions increases with increases in the deviation of  from zero. Shao observed a similar problem in the IB-LBM calculation and assumed that this problem could be attributed to the high sensitivity of the  to the contact angle θ as shown in Eq. (47) [36]. 3.5. A contact angle of a static droplet on a flat plate In the study, the contact angle of a static droplet is verified on a hydrophilic or hydrophobic plate in two and three dimensions. A droplet of diameter 75δ x is centered in a computational domain with 150m grid points in m-dimensions. With respect to the comparative study, the contact angle is computed using the IB-LBM proposed by Shao in two dimensions. In order to eliminate the influence of fluid dynamics in the solid region on the contact angle, a set of boundary points is located along a straight line at y = 75δx which corresponds to the position of the center of the droplet (Fig. 10). The boundary points of the IB-LBM do not possess

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

75

1.15

1.00

1.05

1.00

0.95

(a) two dimensions

Fig. 9. Comparison between numerical and theoretical values of φ on a boundary. The solid line denotes the theoretical solution (68). The symbols (◦) show the numerical values as computed by the IB-LBM proposed by Shao et al. [36].

1.15

1.10

1.05

1.00

0.95

(b) three dimensions Fig. 8. The order parameter profiles of the convergent solution in the transition layers under a hydrophilic boundary condition (θ = 60◦ ). With respect to the threedimensional simulations, the order parameter distribution is shown at the midplane at z = 80δx .

thickness, and thus a droplet evolves based on the given wettability and forms symmetrical shapes on both sides of the boundary. As shown in Fig. 1, the boundary of the SP-LBM possesses thickness. When the fluid-slid interface width ξi = 2δx in the SP-LBM, the minimum thickness of the plate corresponds to 4δ x . In order to implement the same boundary conditions as shown in Fig. 10 in the calculation by the SP-LBM, a side of the plate is located at the position corresponding to y = 75δx , and the other side of the plate is located at the position corresponding to y = 71δx . It is assumed that the fluid-fluid interface thickness w corresponds to 5δ x . The Neumann boundary conditions corresponding to the equilibrium contact angle are implemented from 30° to 150°. The initial values are set as n = 1, u = 0,  = 0.1, σ = 0.04, and τn = τφ = 0.9. In order to compute the droplet contact angle, the following equation is used:

θ = 2 arctan

2h

y

hx

,

(71)

where hy denotes height of a droplet from the plate, and hx denotes the width of the droplet along the plate. Fig. 11 shows an

Fig. 10. A schematic diagram of a droplet on a flat plate. The open circles indicate the boundary nodes along the flat plate in the IB-LBM.

equilibrium shape of the droplets above the flat plate. As shown in Fig. 11, the contact angles of the droplet are adequately controlled by the present SP-LBM in two and three dimensions. Fig. 12 shows the relation between the non-dimensional wetting potential and the equilibrium contact angle. The numerical results computed by the present method are almost identical to the theoretical values in the range from 30° to 150°. Fig. 12 shows that the numerical values of the SP-LBM agree well with those of the IB-LBM based on the correction for the order parameter and for the chemical potential in Eqs. (50) and (51). The method in the present study treats arbitrarily shaped boundary conditions on the fixed Cartesian coordinate system in a more simple manner when compared to that of the IB-LBM. Additionally, the present method can implement the Neumann boundary conditions with an accuracy almost equal to that of the IB-LBM.

76

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

Fig. 11. The droplet shape on a hydrophilic or hydrophobic flat plate. The fluid-fluid interfaces are indicated by the contour line of φ = 0.

3.6. Static contact line on a circular cylinder or on a sphere

Fig. 12. Comparison between the numerical results and the theoretical solution for the droplet contact angle on a flat plate. The solid lines denote the theoretical solution (47).

This section demonstrates the capacity of the present scheme to compute contact line dynamics on a single circular cylinder or a single sphere. A circular cylinder of radius of 40δ x is located at the center of square domain covered with 160 × 160 grid points. The radius of a sphere corresponds to 40δ x , and the length of a side of the cubic domain corresponds to 160δ x . The upper-half and lower-half regions of the computational domain are filled with fluids for φ = −1, and for φ = 1, respectively. Initially, the flat fluidfluid interface with w = 5δx is set at the center of the computational domain. As shown in Section 3.4, neutral wetting boundary conditions (natural boundary conditions) for φ = −1 and for φ = 1 are imposed on the top and bottom walls, respectively. The periodic boundary conditions are applied at the other boundaries. The following parameters are given with respect to the computation: n = 1, σ = 0.001, τn = 0.75, τφ = 0.75,  = 10, and ξi = 2δx . Fig. 13(a) shows the fluid-fluid interface around a sphere at the steady state. As shown in Fig. 13(a-1), the fluid-fluid interface increases in the vicinity of the solid surface when the contact angle equals 60°. When the contact angle is set as 120°, the interface decreases due to the hydrophobic property as shown in Fig. 13(a-2). In order to verify the accuracy of the present method on the curved walls, the results of the present study are compared with the two-dimensional numerical results obtained by

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

77

Fig. 13. Contact lines on curved surfaces. The dashed lines represent the surfaces of a circular cylinder or those of a sphere in (b). The dotted lines denote the numerical results as calculated by the IB-LBM in two dimensions in (b).

the IB-LBM that employs adequate contact angles as used in a previous study [36]. It is not possible to equidistantly locate Lagrangian points on the surface of a sphere, and the numerical accuracy depends on the discrete length s in the implicit correction method [49]. Therefore, the two-dimensional results as calculated by the IB-LBM method are used. The three-dimensional and twodimensional results of the present study are shown in the figures on the left-hand and right-hand side of Fig. 13(b-1) and (b-2), respectively. The dotted lines represent the fluid-fluid interfaces as calculated by the IB-LBM in the D2Q9 model. In the right hand side figures, the interfaces calculated by the present SP-LBM agree well with the those in the IB-LBM as shown in Fig. 13(b-1) and (b-2).

In the left-hand side figures, the three-dimensional present results indicate a slight deviation from the two-dimensional IB-LBM results due to the three-dimensional effects. The interface profile on the surface of the sphere affects the shape of the order parameter on the opposite surface through a periodic boundary. The distance between the surface of the sphere and the opposite surface varies as a function of direction in the cubic domain through a periodic boundary. The shape of the order parameter around the surface depends on the direction in three dimensions. These results demonstrate that the SP-LBM in the present study can be used to easily calculate the static contact angle on the curved boundary with the same accuracy as that of the IB-LBM.

78

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

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

3.7. Contact line dynamics on a plate or on a circular cylinder We calculate the spreading of a water droplet on a solid surface in the computational domain covered by 220 × 150 grid points with periodic boundary conditions. A flat plate with a width 4δ x is located at the position corresponding to y = 18δx . The initial contact angle of a droplet on the plate is 160°. A droplet with radius 45δ x spreads on a flat plate in which the equilibrium contact angle θ is varied between 20° and 150° with the following parameters: n = 1,  = 1.0, σ = 0.05, τn = τφ = 0.9, and w = 5δx . The non-dimensional time tˆ is given by δ x σ t/Dnν [36]. Here, t is an iteration number of the LBM calculation, and D = 90δx is the initial diameter of the droplet. The sequential snapshots of the spreading droplets on the plates with the equilibrium contact angle of θ = 60◦ and θ = 120◦ are shown in Fig. 14(a) and (b), respectively. The droplet on the hydrophilic wall elongates along the horizontal direction at a higher speed for θ = 60◦ than for θ = 120◦ . The temporal evolutions of the droplet height hy and width hx normalized by the initial droplet diameter D are shown in Fig. 15(a). The present results agree well with a previous investigation using a body-conformal grid [36]. Fig. 15(b) shows the effect of the equilibrium contact angle on the spreading behavior of a droplet. The numerical solutions resulting from a phase-field model were shown to substantiate the following equation: hx /D ∝ tˆm for 1 < tˆ < 10 [4]. The curve-fitting results give m = 0.385, 0.380, 0.369, 0.352, 0.335, 0.313, and 0.237 for θ = 20°, 30°, 45°, 60°, 75°, 90°, and 120° in Fig. 15(b). These values agree well with the previous results: m = 0.38, 0.37, 0.35, 0.34, 0.333, 0.33, and 0.215, for θ = 20°, 30°, 45°, 60°, 75°, 90°, and 120° in Reference [4]. The results shown in Fig. 15 demonstrate that the proposed SP-LBM can correctly predict the contact line dynamics in drop spreading on a flat plate. Finally, we compute the motion of a fluid-fluid interface along a curved surface in two dimensions, and compare the results obtained using the proposed SP-LBM and the IB-LBM. The initial conditions of the fluid-fluid interface are given by the convergent solutions for θ = 60◦ or for θ = 120◦ , as shown in Fig. 13(b). In Fig. 16(a), the fluid-fluid interface initially makes contact with the solid surface with θ = 120◦ , and forms the shape corresponding to the equilibrium contact angle of θ = 60◦ . Fig. 16(b) shows the temporal evolution of the contact angle of the interface from 60° to 120°. Although the numerical results obtained by the SP-LBM are almost identical with those by the IB-LBM for wetting boundary condition, the interface motion of the proposed SP-LBM is slightly faster than that of the IB-LBM in Fig. 16. The diffusive nature of the solid-fluid interface of the SPM affects the dynamic contact angle induced by the chemical diffusion in the phase field model on the curved surface. We can conclude that the proposed SP-LBM computes the contact line dynamics on the curved surface with acceptable accuracy.

Fig. 15. Normalized droplet height hy /D and width hx /D plotted against nondimensional time tˆ in the calculation of the droplet spreading on a flat plate.

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

79

Fig. 16. Temporal evolution of the fluid-fluid interface along the curved surface. Dashed lines represent the surfaces of a circular cylinder. The computational domain in these figures corresponds to the enlarged figures of Fig. 13(b).

4. Conclusions In this study, a smoothed profile-lattice Boltzmann method is proposed for the non-slip, non-penetration, no mass-flux, and wetting boundary conditions in the phase field method by introducing two types of indicator functions. The fluid-solid interaction is computed from the definition of fluid velocity including the effects of the body force as proposed by Guo et al. [1] to successfully prevent fluid penetration into the solid wall and reduce the velocity fluctuation in the solid region as observed in the conventional model. The present method achieved first order accuracy space convergence in the computation of the fluid velocity and the force acting on the solid wall in the numerical experiments including those determining the cylindrical Couette flow, the flow around a circular cylinder, and the axial Couette flow. The Neumann boundary conditions (i.e., the no mass-flux and the wetting boundary conditions) are implemented by correcting the order parameter and the chemical potential by the fluxes from the surfaces [36]. An indicator function is formulated such that it considers the contribution of fluxes only from the surfaces within the boundary layer. The fluxes involved in the method used in the present study directly corrected the relevant physical variables at Eulerian points without requiring an interpolation procedure as used in the IBLBM. The numerical values on flat or curved walls agreed well with the theoretical ones in the transition layer around a circular cylinder or a sphere. The proposed SP-LBM can reasonably predict spreading behavior of a drop on a flat plate and the motion of a fluid-fluid interface along a curved surface. The numerical experiments on contact line dynamics demonstrated the accuracy and efficiency of the proposed method to handle hydrophilic and hydrophobic surface conditions. To improve the accuracy of the prediction, we must investigate the effective correction of the dynamic contact angle by the introduction of a slip length on the order of the mesh size in the SP-LBM [50].

Although the present SP-LBM is applied in a study involving contact line dynamics on a arbitrarily shaped boundary, the method is unable to solve the moving boundary problem due to the pressure oscillation around the solid wall. A previous study confirmed the efficiency of the convectional SP-LBM to simulate a neutrally buoyant cylinder in a simple shear flow, or sedimentation of circular cylinders in a viscous fluid [38]. The direct forcing approach employed in the fore-mentioned study is a better choice to apply the SP-LBM to the moving boundary with the wetting boundary condition. With respect to the comparative study with the experimental data, a future study will compute the droplet behavior with a high density ratio and high viscosity contrast by using LBMs [15–19,22]. The SP-LBM for simulating high density ratios can investigate the complicated phenomena such as droplet impact dynamics on micro-pillared hydrophobic surfaces at high Weber and Reynolds numbers. Additionally, the advantages of the present model including the computation of the normal direction to the surface from the Laplacian of the indicator profile ϕ (1) , will be used to explore the wettability of a lotus-leaf-like complicated structured wall in numerical experiments. Acknowledgements This work was supported by JSPS KAKENHI Grant Number JP16K06070. The authors would like to thank Prof. R. Takahashi, Prof. E. Takegosi, Prof. K. Okui, Prof. Y. Abe, and Prof. A. Tomiyama for their helpful suggestions. References [1] Guo Z, Zheng C, Shi B. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys Rev E 2002;65:046308. [2] Xue CH, Jia ST, Zhang J, Ma JZ. Large-area fabrication of superhydrophobic surfaces for practical applications: an overview. Sci Technol Adv Mater 2010;11:033002.

80

T. Seta et al. / Computers and Fluids 159 (2017) 64–80

[3] Jacqmin D. Calculation of two-phase Navier–Stokes flows using phase-field modeling. J Comput Phys 1999;155:96–127. [4] Khatavkar VV, Anderson PD, Meijer HEH. Capillary spreading of a droplet in the partially wetting regime using a diffuse-interface model. J Fluid Mech 2007;572:367–87. [5] Cahn JW, Hilliard JE. Free energy of a nonuniform system. I. Interfacial free energy. J Chem Phys 1958;28:258–67. [6] Chen S, Doolen GD. Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 1998;30:329–62. [7] Gunstensen AK, Rothman DH, Zaleski S, Zanetti G. Lattice Boltzmann model of immiscible fluids. Phys Rev A 1991;43:4320. [8] Shan X, Chen H. Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation. Phys Rev E 1994;49:2941. [9] Lishchuk SV, Care CM, Halliday I. Lattice Boltzmann algorithm for surface tension with greatly reduced microcurrents. Phys Rev E 2003;67:036701. [10] Sergi D, Grossi L, Leidi T, Ortona A. Surface growth effects on reactive capillary-driven flow: lattice Boltzmann investigation. Eng Appl Comput Fluid Mech 2014;8:549–61. [11] Liu H, Ju Y, Wang N, Xi G, Zhang Y. Lattice Boltzmann modeling of contact angle and its hysteresis in two-phase flow with large viscosity difference. Phys Rev E 2015;92:033306. [12] Swift MR, Orlandini E, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of liquid-gas and binary fluid systems. Phys Rev E 1996;54:5041–52. [13] Briant AJ, Wagner AJ, Yeomans JM. Lattice Boltzmann simulations of contact line motion. I. Liquid-gas systems. Phys Rev E 2004;69:031602. [14] Briant AJ, Yeomans JM. Lattice Boltzmann simulations of contact line motion. II. Binary fluids. Phys Rev E 2004;69:031603. [15] Lee T, Lin C-L. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. J Comput Phys 2005;206:16–47. [16] Inamuro T, Ogata T, Tajima S, Konishi N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J Comput Phys 2004;198:628–44. [17] Mazloomi A, Shyam M, Chikatamarla SS, Karlin IV. Entropic lattice Boltzmann method for multiphase flows: fluid-solid interfaces. Phys Rev E 2015;92:023308. [18] Fakhari A, Geier M, Lee T. A mass-conserving lattice Boltzmann method with dynamic grid refinement for immiscible two-phase flows. J Comput Phys 2016;315:434–57. [19] Wang Y, Shu C, Yang LM. An improved multiphase lattice Boltzmann flux solver for three-dimensional flows with large density ratio and high Reynolds number. J Comput Phys 2015;302:41–58. [20] Moradi N, Varnik F, Steinbach I. Roughness-gradient induced spontaneous motion of droplet on hydrophobic surfaces: a lattice Boltzmann study. Europhys Lett 2010;89:26006. [21] Connington K, Lee T. Lattice Boltzmann simulations of forced wetting transitions of drops on superhydrophobic surfaces. J Comput Phys 2013;250:601–15. [22] Takada N, Matsumoto J, Matsumoto S, Kurihara K. Phase-field model-based simulation of two-phase fluid motion on partially wetted and textured solid surface. J Comput Sci 2016;17:315–24. [23] Filippova O, Hanel D. Grid refinement for lattice-BGK models. J Comput Phys 1998;147:219–28. [24] Mei R, Luo L-S, Shyy W. An accurate curved boundary treatment in the lattice Boltzmann method. J Comput Phys 1999;155:307–30. [25] Bao J, Peng Y, Schaefer L. A mass conserving boundary condition for the lattice Boltzmann equation method. J Comput Phys 2008;227:8472–87. [26] Feng ZG, Michaelides EE. The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems. J Comput Phys 2004;195:602–28.

[27] Feng ZG, Michaelides EE. Proteus: a direct forcing method in the simulations of particulate flows. J Comput Phys 2005;202:20–51. [28] Niu XD, Shu C, Chew YT, Peng Y. A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows. Phys Lett A 2006;354:173–82. [29] Chen L, Yu Y, Hou G. Sharp-interface immersed boundary lattice Boltzmann method with reduced spurious-pressure oscillations for moving boundaries. Phys Rev E 2013;87:053306. [30] Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys 1972;10:252–71. [31] Wu J, Shu C. Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications. J Comput Phys 2009;228:1963–79. [32] Suzuki K, Inamuro T. Effect of internal mass in the simulation of a moving body by the immersed boundary method. Comput Fluids 2011;49:173–87. [33] Kang SK, Hassan YA. A comparative study of direct-forcing immersed boundary-lattice Boltzmann methods for stationary complex boundaries. Int J Numer Methods Fluids 2011;66:1132–58. [34] Zhang C, Cheng Y, Zhu L, Wu J. Accuracy improvement of the immersed boundary-lattice Boltzmann coupling scheme by iterative force correction. Comput Fluids 2016;124:246–60. [35] Seta T, Rojas R, Hayashi K, Tomiyama A. Implicit-correction-based immersed boundary-lattice Boltzmann method with two relaxation times. Phys Rev E 2014;89:023307. [36] Shao JY, Shu C, Chew YT. Development of an immersed boundary-phase fieldlattice Boltzmann method for Neumann boundary condition to study contact line dynamics. J Comput Phys 2013;234:8–32. [37] Ren W, Shu C, Yang W. An efficient immersed boundary method for thermal flow problems with heat flux boundary conditions.. Int J Heat Mass Transf 2013;64:694–705. [38] Jafari S, Yamamoto R, Rahnama M. Lattice-Boltzmann method combined with smoothed-profile method for particulate suspensions. Phys Rev E 2011;83:026702. [39] Hu Y, Li D, Shu S, Niu X. An efficient smoothed profile-lattice Boltzmann method for the simulation of forced and natural convection flows in complex geometries. Int Commun Heat Mass Transfer 2015;68:188–99. [40] Nakayama Y, Yamamoto R. Simulation method to resolve hydrodynamic interactions in colloidal dispersions. Phys Rev E 2005;71:036707. [41] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59. [42] Liang H, Shi BC, Guo ZL, Chai ZH. Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows. Phys Rev E 2014;89:053320. [43] Zheng HW, Shu C, Chew YT. A lattice Boltzmann model for multiphase flows with large density ratio. J Comput Phys 2006;218:353–71. [44] Bunch JR, Rose DJ. Sparse matrix computations. 2014. Academic Press. [45] Wang S, He G, Zhang X. Parallel computing strategy for a flow solver based on immersed boundary method and discrete stream-function formulation. Comput Fluids 2013;88:210–24. [46] Dennis SCR, Chang GZ. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J Fluid Mech 1970;42:471–89. [47] Nieuwstadt F, Keller HB. Viscous flow past circular cylinders. Comput Fluids 1973;1:59–71. [48] He X, Doolen GD. Lattice Boltzmann method on curvilinear coordinates system: flow around a circular cylinder. J Comput Phys 1997;134:306–15. [49] Seta T. Implicit temperature correction-based immersed boundary-thermal lattice Boltzmann method for the simulation of natural convection. Phys Rev E 2013;87:063304. [50] Jacqmin D. Contact-line dynamics of a diffuse fluid interface. J Fluid Mech 20 0 0;402:57–88.