- Email: [email protected]
Oriented ordering parameters for free energy lattice Boltzmann methods using the bounce-back boundary condition
Computers and Mathematics with Applications (
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Oriented ordering parameters for free energy lattice Boltzmann methods using the bounce-back boundary condition Sigvat Stensholt Polytec Research Institute, Sørhauggata 128, 5527 Haugesund, Norway
article
info
Article history: Available online xxxx Keywords: Lattice Boltzmann Boundary conditions Cassie–Baxter Contact angle Simulation stability
abstract We evaluate the Lattice-Boltzmann model that combines the free energy approach for multiphase flow with the bounce-back boundary condition. This method requires a virtual ordering parameter to be assigned to the wall nodes. Existing literature has assigned a single ordering parameter for the entire wall node, but this is too coarse to give accurate simulations since it allows for a spurious interaction of different phases through the wall nodes. The consequences are illustrated by examining rough and hydrophobic surfaces where contact angles can be predicted using the Cassie–Baxter equation. It is found that the single value approach gives a large systemic error in the measured contact angles. Considerable improvement can be obtained by assigning different values of the ordering parameter to each side of the wall node. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction The behavior of fluids on different forms of surfaces is a field of study with multiple scientific and industrial applications. The chemical and structural characteristics of a surface can determine whether water tends to stick to or roll off it. These attributes are significant, for example for applications that require surfaces to be as dry as possible. If a surface attracts water, the droplet will spread out giving a contact angle θ less than 90°. Such a surface is called hydrophilic. In contrast, a water droplet on a surface that repels water will curl up and have a contact angle greater than 90°. Such a surface is called hydrophobic. Very hydrophobic surfaces with contact angles greater than 150° are called superhydrophobic and are of great industrial interest. Several proposals for manufacturing superhydrophobic surfaces exist, many with the aim of achieving characteristics such as ice-phobicity and durability, see for example [1,2]. Many of these publications report on laboratory testing of the results, but fluid dynamics simulations by using mesoscopic methods such as the Lattice Boltzmann method have also been applied [3]. Laboratory tests require the production of a physical test surface, something that is not needed for computer simulations. In this article we examine the feasibility of using the popular bounce-back boundary condition for the solid–fluid interface when using the free energy approach [4,5] for multiphase flow. We aim to improve on the method offered in the previous literature [6,7] and show that this gives simulated contact angles that are more consistent with classical fluid theory. This is applicable to any domain with a complex solid boundary. The classical fluid theory for rough surfaces is covered in Section 2. The Lattice Boltzmann method and the associated boundary conditions that are the focus for this article is covered in Section 3. Simulation results are covered in Section 4, and the conclusions are drawn in Section 5. E-mail address: [email protected]. http://dx.doi.org/10.1016/j.camwa.2015.06.014 0898-1221/© 2015 Elsevier Ltd. All rights reserved.
2
S. Stensholt / Computers and Mathematics with Applications (
)
–
2. Fluid theory of rough surfaces A roughness is required to make a surface superhydrophobic. The contact angle of droplets on heterogeneous and rough surfaces has been the subject of considerable study. The simplest equations for estimating the contact angle on rough surfaces are the Wenzel [8] and Cassie–Baxter [9] equations. The Wenzel equation applies when a droplet fills the cavities of a rough surface. The apparent contact angle θc is then given by cos θc = r cos θ0 ,
(1)
where θ0 is the contact angle the material would give on a smooth surface and the roughness parameter r is the ratio between the actual surface area and the projected surface area. In its original form, the Cassie–Baxter equation applies on a smooth surface with n components with surface fractions ξ1 , ξ2 , . . . , ξn with different contact angle parameters θ1 , θ2 , . . . , θn . For a smooth surface made of different components, the sum over surface fractions will be 1. The resultant contact angle θc on a composite surface is given by cos θc =
n
ξi cos θi .
(2)
i=1
The Cassie–Baxter equation has also been applied to some rough surfaces where the droplet is suspended above the cavities. In this case air, which has a contact angle parameter of 180°, is treated as a separate component. If the surface fraction of the solid component is ξ , one obtains simplified version of Eq. (2) which has been used in a majority of articles that deal with contact angles on a rough surface [10]: cos θc = ξ cos θ − (1 − ξ ).
(3)
Eq. (3) assumes that the solid part of the surface is made up of a single component with an initial contact angle θ that covers a fraction ξ of the surface, and that the droplet does not penetrate into the surface at all. While these assumptions are commonly made, real droplets often exhibit more complex behavior by penetrating somewhat into the surface in which case Eq. (2) should be used without the assumption that the sum over all ξi is equal to 1. As the droplet penetrates deeper into the surface ξsolid approaches the roughness factor and ξair approaches 0. Consequently, Eq. (2) approaches the Wenzel equation (1) [10]. The Cassie–Baxter and Wenzel equations can be derived from thermodynamic analysis [11] and provide a means for analytical predictions of the contact angle for simple cases of stationary droplets that are entirely on top of the surface or entirely penetrated into the surface. However the validity of the equations has been the subject of skepticism [12] and the fluid dynamics for the transitional states or droplets in motion cannot be described fully by the equations. There is therefore an interest in a computational fluid dynamic approach to studying such surfaces and the scales involved invite the use of mesoscopic approaches. 3. Lattice Boltzmann method The Lattice Boltzmann method is a computational fluid dynamics (CFD) model originally introduced by Mcnamara and Zanetti [13] and later interpreted as a discretization of the Boltzmann equation [14,15]. The method is commonly used on a mesoscopic scale and based on the discrete particle distribution function fi (x, t ) which gives the number of particles moving at a given velocity at site x at time t. For two-component flow, two distribution functions are needed. One option is to have a separate distribution function for each component. The other option, used here, is to let the first distribution function fi (x, t ) account for the sum of the components while the second distribution function gi (x, t ) accounts for the difference between them. The distribution functions evolve according to the equations fi (x + ei δ t , t + δ t ) = fi (x, t ) + Ωf ,i + Fi
(4)
gi (x + ei δ t , t + δ t ) = gi (x, t ) + Ωg ,i
(5)
Fi accounts for any external forces such as gravity. Ωi is the collision operator where we have used the single relaxation time operator [16],
Ω f ,i =
fi − fi
τf
eq
eq
,
Ωg ,i =
gi − gi
τg
(6)
where the superscript eq represents an equilibrium distribution and f and g are the relaxation times. The relaxation time for fi gives the kinematic viscosity ν by
ν=
2τf − 1 6
δt .
(7)
S. Stensholt / Computers and Mathematics with Applications (
)
–
3
Space is discretized into a regular cubic lattice. We have opted to use the D3Q19 lattice where particles may move along 19 possible vectors ei . The lattice spacing δ x is set to one. e0 = (0, 0, 0) e1−6 = (±1, 0, 0), e7−18
(0, ±1, 0), (0, 0, ±1) = (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1).
(8)
1 1 , and w7−18 = 36 . The corresponding weight coefficients wi for this lattice are w0 = 13 , w1−6 = 18 The mass, momentum, and ordering parameters at each site can be calculated from the moments
ρ(x, t ) =
f i ( x, t )
(9)
i
ρ(x, t )ui (x, t ) =
fi (x, t )ei
(10)
i
ϕ(x, t ) =
gi (x, t ).
(11)
i
A number of methods to simulate multiphase flow within the lattice Boltzmann framework have been proposed. The oldest method is the color model [17] which redistributes phases to enforce phase separation and imposes a surface tension force. The method provides fairly sharp interfaces and an explicit contact angle parameter to handle wetting phenomena, but it is known to possess a measure of anisotropy and its computational expense is a drawback that becomes serious for three-dimensional simulations [18]. The pseudo-potential method [19] is a simpler method that imposes a force between opposite phases. It is used extensively, but it produces diffuse interfaces and even though the contact angle arising from the model has been studied [20], a priori determination of the contact angle remains difficult. For multiphase flow we have instead employed the free energy approach which was first developed by researchers at Oxford [4,5,21]. This method has previously been used for simulating superhydrophobic surfaces [3] and has the advantage of being thermodynamically consistent. Another significant advantage is the ability to impose an explicit contact angle on a smooth homogenous surface through a single contact angle parameter. Explicit expressions for the free energy model can be found in [22] and for the three-dimensional D3Q19 lattice in [6], which largely uses the method described by [7]. A difference is that Hao and Cheng used the BGK operator (6) while Niu et al. used the multiple relaxation time (MRT) operator [23] which offers more degrees of freedom but is computationally more demanding than the BGK operator. Central to the understanding of the free energy model is the free energy density ψ(ϕ) given by A
B
2
4
ψ(ϕ) = − ϕ 2 + ϕ 4 .
(12)
The parameters A and B are positive and the free energy density function is therefore a quartic function that form a double well potential. The chemical potential µ is given by
µ = −Aϕ + Bϕ 3 + κ∇ 2 ϕ.
(13)
The free energy model uses an explicit thermodynamic pressure tensor Pth which can be written out as the sum of the chemical pressure tensor Pch and hydrostatic pressure P hy giving1 Pth = Pch + P hy
= p0 I + κ(∇ϕ)(∇ϕ)T + ρ I κ = ϕµ − ψ(ϕ) − (∇ϕ)2 + κ(∇ϕ)(∇ϕ)T + ρ I 2
A
=− ϕ + 2
2
3B 4
ϕ − κϕ∇ 2 ϕ − 4
κ 2
(∇ϕ)2 + κ(∇ϕ)(∇ϕ)T + ρ I.
(14)
1 There are typographical errors in some of the equations in Hao and Cheng’s 2009 article. Most notably, the term for p which should be given as 0 p0 = ϕµ − ψ(ϕ) −
κ 2
A
3B
2
4
(∇ϕ)2 = − ϕ 2 +
ϕ 4 − κϕ∇ 2 ϕ −
κ 2
(∇ϕ)2
has substituted the quadratic term 3B ϕ 2 for the correct quartic term, thus destroying the double well potential. Also, µ and κ have been rendered as u and 4 k. The corresponding author of the article has been contacted and indicated that an erratum will be published.
4
S. Stensholt / Computers and Mathematics with Applications (
)
–
The gradient (∇ ) and Laplacian (∇ 2 ) operators are conducted using a 19-point stencil given by [7]. The numerical estimate for the jth component of the gradient is given by 3
18
wi ϕ(x + ei δ t )eij
i =0
∇j ϕ(x) =
δt
18
,
(15)
wi (eij δ t )2
i=0
while the Laplacian is given by 2
∇ ϕ= 2
18
wi (ϕ(x + ei δ t ) − ϕ(x))
i=0 18
.
(16)
wi (eij δ t )2
i =0
The main feature of the free energy model is the expressions for the equilibrium distributions that depend on the pressure tensor and chemical potentials, that are dependent on ϕ and derivatives of ϕ . For the free energy model on a D3Q19 lattice, the explicit expressions for the equilibrium distributions are fi
eq
eq
gi
= Qi + wi = Hi + wi
ei · u cs2
ei · u cs2
+
(ei · u)2
+
(ei · u)2
2cs4 2cs4
− −
u·u
2cs2 u·u 2cs2
,
(17)
,
(18)
√
where cs = c / 3 is the speed of sound and the coefficients Qi and Hi are given by
th th th − Pxx − Pxx )/3c 2 ρ − 2(Pxx th th th Qi = eTi P th ei /4c 4 − (Pxx − Pxx − Pxx )/36c 2 T th th th th ei P ei /8c 4 − (Pxx − Pxx − Pxx )/18c 2 ϕ − 2Γ µ/c 2 i = 0; Hi = Γ µ/6c 2 i = 1, 2, . . . , 6; Γ µ/12c 2 i = 7, 8, . . . , 18.
i = 0; i = 1, 2, . . . , 6; i = 7, 8, . . . , 18,
(19)
(20)
The parameters A, B, and κ that all affect the phase separation equations have already been introduced. The expression for Hi introduced a fourth parameter, Γ , that is related to the mobility. The bulk pressure consists of the ideal component (cs2 ρ ) and the components due to interfaces. It can be written as p = cs2 ρ −
A 2
ϕ2 +
3B 4
ϕ 4 − κϕ∇ 2 ϕ −
κ 2
(∇ϕ)2 .
(21)
3.1. Boundary conditions In our simulations we have employed periodic boundary conditions at the fluid node edges of the domain. For the fluid–solid interface the boundary conditions are required to handle both the physical properties (nonpenetration of the wall and no-slip) and the chemical properties affecting contact angles. The simplest method for handling the solid interface nodes is the bounce-back condition. In this approach any particles that would propagate into the wall node are instead reflected in the opposite direction. This method ensures that no element of the distribution function will remain unknown. The method gives a no-slip boundary condition consistent with a wall placed halfway between the last fluid node and the first wall node [24]. The accuracy of the bounce-back condition has been studied extensively. A particular result to note is that the accuracy is dependent on the relaxation time τ , being more accurate for smaller relaxation times (i.e. less viscous flows), and with an exponential deterioration in accuracy for τ > 1 [25]. The position of the wall may also be affected by τ [26]. Other boundary conditions are possible. The first article to use lattice Boltzmann to simulate superhydrophobic surfaces [27] used a set of velocity boundary conditions were to determine the missing elements of the distribution function near the fluid–wall boundary. The bounce-back condition can be modified so that the reflection to occur on the node instead of between nodes, and this approach has also been applied for studies of superhydrophobic surfaces [28]. Some wall boundary conditions result in artificial parasitic currents at the contact line, but these are eliminated if intermolecular force is handled using the potential form [29].
S. Stensholt / Computers and Mathematics with Applications (
)
–
5
Fig. 1. Wall separating ‘‘red’’ and ‘‘blue’’ phases. In (a) the wall consists of a single row of wall nodes. In (b) the wall consists of a double row of wall nodes. Colors on the wall nodes indicate the value of ϕs assigned to it. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In order to allow the gradient and Laplacian of the ordering parameter to be calculated a virtual ordering parameter for wall nodes, ϕs is needed. There are various approaches to this. The simplest approach is to simply assign a fixed value of ϕs for each wall node. The method has been used by [30] who give the resultant contact angle θ as 3
cos θ =
2
ϕ˜ 1 −
ϕ˜ 2
3
,
(22)
1/2
where ϕ˜ = ϕs / AB . The simplicity of not needing to update the values of ϕs at every time step is an attractive feature of the model. This approach does not account for the ϕ in the surrounding fluid nodes, and will therefore give an artificial gradient and Laplacian around the walls. To avoid the artificial gradients the value of ϕs needs to change with the surrounding fluid. The approach by [7] uses a target contact angle θp . This parameter is used to determine a virtual ϕs for the wall nodes. A wetting potential Θ is given by
Θ = 2 sgn
π
− θp
2
cos
arccos(sin2 θp ) 3
1 − cos
arccos(sin2 θp )
3
(23)
where all angles are given in radians. This wetting potential is used to calculate the parameter λ:
λ = −Θ
κA 2
.
(24)
The method given in [7] approximates the order parameter for the wall nodes ϕs as the average of all the N surrounding fluid nodes, modified by the contact angle parameter and surface tension parameters λ and κ .
ϕs =
(ϕf − δ xλ/κ)
N
N
.
(25)
This value of ϕs can then be used along with the values of ϕf for the fluid nodes so that all nodes are assigned a value of ϕ . One can then proceed with the calculations of ∇ϕ and ∇ 2 ϕ in Eqs. (15) and (16) for the relaxation step of the method. Hao and Cheng [6] write that Eq. (25) ‘‘can easily be applied to complex solid boundary such as porous media’’. However, we observe that averaging over all the surrounding nodes leads to error in cases where two different phases are in contact with the same wall node. In practice this includes any complex wall with corners or edges, as well as any wall node that separates two fluids of different phase. To illustrate the problem, consider a wall separating to different fluid phases. In Fig. 1, a single row of wall nodes, as shown in (a), is in contact with both the red phase with ϕ = 1 and blue phase with ϕ = −1, and so the assigned ordering parameter for the wall nodes is ϕs = 0 − δ xλ/κ , which is incorrect for both sides of the wall. The thin wall allows the unphysical interaction of phases through the wall. The ordering parameter of the wall is correct if a double row of wall nodes is used as seen in (b). In this case each wall node is only in contact with one phase so the wall nodes in contact with the red phase will have an ordering parameter ϕs = 1 − δ xλ/κ while those in contact with the blue phase will have ϕs = −1 − δ xλ/κ , which retains only the desired gradient and Laplacian from the wetting factor δ xλ/κ . A similar problem arises for rough surfaces that feature edges and corners. In Fig. 2 we see that the corner wall node that is in contact with both phases has its value for ϕs affected by both phases and this value is not valid on any side. The problem with incorrect values of ϕs cannot be entirely eliminated in this case, even though increasing the resolution reduces the problem by decreasing the proportion of corner wall nodes.
6
S. Stensholt / Computers and Mathematics with Applications (
)
–
Fig. 2. Details of a fluid interface at a wall edge. x1 and x2 are fluid nodes in the red (ϕ = 1) and blue (ϕ = −1) phases respectively. The ϕs value for each wall node in contact with a fluid phase is noted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In order to avoid introducing an artificial gradient, the edge wall node in Fig. 2 needs to have a different values of ϕs on each side. A significant improvement is to give the ordering parameter of the wall an orientation. Each side and edge of the wall node is assigned an ordering parameter ϕs,i , where subscript i indicates the side or edge of the wall node.
ϕs,i =
(ϕf (x + ei ) − δ xλ/κ)
i∈Eˆ i
,
Nˆ
(26)
where Eˆ i is the set of three (for i = 7, . . . , 18) or five (for i = 1 . . . 6) non-zero indices ˆi so that eˆi − ei ≤ 0. Nˆ is the number of fluid nodes adjacent to the face or edge of the wall node. For example, on the i = 1 surface, which is oriented in the direction (1, 0, 0), the values ϕf (x +(1, 0, 0)), ϕf (x +(1, 1, 0)), ϕf (x + (1, −1, 0)), ϕf (x + (1, 0, 1)), and ϕf (x + (1, 0, −1)) are used in the calculation for ϕs,1 . To assign a value for the i = 7 edge, which is oriented in the direction (1, 1, 0), ϕs,7 is calculated using the fluid node order parameters ϕf (x + (1, 1, 0)), ϕf (x + (1, 0, 0)), and ϕf (x + (0, 1, 0)). The different value of ϕs on different sides of the wall nodes is reflected in a small modification to Eqs. (15) and (16) for the gradient and Laplacian. These are now given by
3
∇j ϕ(x) =
18
wi ϕi* (x + ei δ t )eij
i =0
δt
18
, wi (eij δ t )
(27)
2
i=0
while the Laplacian is given by 2
∇ 2ϕ =
18
wi (ϕi* (x + ei δ t ) − ϕ(x))
i=0 18
(28)
wi (eij δ t )
2
i=0
where i* is the index opposite i, that is ei* = −ei . For fluid nodes, ϕi = ϕ as defined in Eq. (11) is valid for all directions i. In the problem illustrated by Fig. 1, a single row of wall nodes is now sufficient to simulate a separating wall as the top side of the wall will have a ‘‘blue’’ ordering parameter while the bottom side will be ‘‘red’’. For the corner node in Fig. 2 the top side of the node will have a ‘‘red’’ ordering parameter while the vertical side on the right will have a ‘‘blue’’ one. 4. Simulations and results Lattice Boltzmann simulations of sessile droplets on pillared surfaces have been conducted and the resulting contact angles have been measured. Simulations have been conducted with schemes using the original non-oriented ϕs (x, t ) and the proposed oriented ϕs,i (x, t ) for the wall nodes so that the results can be compared.
S. Stensholt / Computers and Mathematics with Applications (
)
–
7
Fig. 3. Top-down view of 2 × 2 pitted (left) and pillared (right) surfaces.
The contact angle parameter θ0 was set to 110°. The aim in this paper is to consider the theoretical effects of modifying the surface rather than simulating real fluids. The relaxation parameters have been set for maximum simulation stability, τf = τg = 1. The interaction parameters are A = B = 0.007, κ = 0.0078, and Γ = 1. For best reliability, the dimensions of the surface details should be larger than the interface width. If this condition is violated, a simulation can still be performed but the results are less reliable. Nonetheless, computational constraints limit the size of the domain and the achievable fineness of the surface details. It is therefore of interest to investigate the reliability of simulations where the details are smaller than recommended. The smallest details (pillars or pits) are of size 2 × 2. Any detail which leaves a fluid node sandwiched between two wall nodes is prone to destabilize the simulation. On the other hand, the chosen surface details are not so far apart that the droplets will penetrate appreciably into the gaps between the pillars. This allows us to compare the results with the simplified Cassie–Baxter equation (3). We have run our simulations on pitted surfaces and pillared surfaces where the pits and pillars cover a quarter of the top surface area. A top-down view of the arrangement of the 2 × 2 surface details is shown in Fig. 3. The arrangement for the 3 × 3 and larger surface details are similar. The interface produced by the free energy model is diffuse, but the ϕ = 0 isosurface can be used to define its location. After several time steps the droplets converge into a final shape. We examined cross sections through the center of the droplets and found, by fitting a circle over the ϕ = 0 contours, that all the droplets assumed a shape of a spherical cap within the degree of precision allowed by the fineness of the lattice. This realization allows us to calculate the contact angle θ by measuring the radius α of the circular contact area and the droplet height β as shown in Fig. 4: 2αβ . (29) β 2 − α2 The measurements α and β are taken on the fitted sphere half a node beyond the last fluid node. This is because the boundary
θ = arctan
is between the last fluid node and the first wall node. 4.1. Fixed value of ordering parameters For the rough pillared and pitted surfaces, simulations using the fixed values of ϕs were unstable with the selected input parameters. Typically, the simulation broke down after about 300 time steps illustrating the problem of using this method on rough surfaces. 4.2. Grid refinement and convergence tests of proposed method A series of tests were run to determine the effects of grid refinement and the evolution of the droplet shape and contact angle when using Eq. (26). The finest domain features a grid 100 nodes long and wide, and 70 nodes tall. The floor consists of a base and pillars that have a 6 node × 6 node cross section and a height of 4 nodes. The solid surface fraction ξ = 0.25. The initial contact angle θ0 is 110°. According to Eq. (3), the theoretical contact angle of a droplet in the Cassie–Baxter state is then 146.7°. A droplet of radius r is initiated to rest atop the pillars and the simulation is started. Tests have also been run on three coarser grids where the pillar widths, and droplet radii are scaled down proportionally. The pillar height was not scaled down however. If a pillar is only one or two nodes tall, the interface width will cause the droplet to make contact with the surface base, causing the simulated droplet to assume a Wenzel state. Therefore, the pillar height is kept at four nodes for all domain sizes. The results are given in Table 1.
8
S. Stensholt / Computers and Mathematics with Applications (
)
–
Fig. 4. Measurements used to determine the contact angle θ .
Fig. 5. Diagram of a simulated droplet on a pillared surface. The droplet surface indicates the ϕ/ρ = −0.8 isosurface instead of the ϕ = 0 isosurface which is inside this surface and entirely over the top of the pillars. Table 1 Contact angles of a sessile droplet at different time steps and varying degrees of resolution. Domain size
Pillar width
Droplet radius r
Time steps 10 000
20 000
30 000
40 000
50 000
100 × 100 × 70 84 × 84 × 60 67 × 67 × 50 50 × 50 × 42
6×6 5×5 4×4 3×3
35.0 29.4 23.5 17.5
123.1° 127.3° 133.2° 139.7°
142.1° 140.6° 144.3° 143.7°
146.2° 146.5° 146.1° 144.7°
147.4° 147.3° 147.1° 145.5°
147.8° 147.1° 146.2° 146.3°
The maintenance of the spherical shape of the droplet is apparent already at 10 000 time steps, but about 40 000 time steps are required before the contact angle converges to a figure where the change between 10 000 time steps is less than 1°. Since the contact angle is greater than 110°, the contact area between the droplet and the surface is small, and only
S. Stensholt / Computers and Mathematics with Applications (
)
–
9
covers the tops of seven pillars. Despite this limited contact area, the resultant contact angles are in close agreement with the theory of Cassie and Baxter. 4.3. 2 × 2 and 3 × 3 pitted surfaces These surfaces contain square pits carved into an otherwise flat surface. The pits have a size of 2 × 2 or 3 × 3 and spaced so that they cover 25% of the surface area. Using the simplified Cassie–Baxter Equation (3), ξ = 0.75. The results are given in Table 2. Table 2 Theoretical and measured contact angles of droplets on pitted surfaces. Pit size
Initial contact angle
Theoretical contact angle
Measured contact angle using Eq. (25)
Measured contact angle using Eq. (26)
2×2 3×3
110° 110°
120.4° 120.4°
140.4° 132.9°
126.8° 122.0°
4.4. 2 × 2 and 3 × 3 pillared surfaces These surfaces contain square, flat-topped pillars imposed over the surrounding surface. The pillars have a size of 2 × 2 or 3 × 3 and spaced so that they cover 25% of the surface area. All simulations were conducted in Matlab on a cubic D3Q19 lattice with dimensions 80 × 80 × 80 for 50 000 time steps. In all simulations the fluid mass of both phases remained constant throughout the simulation. The methods are all mass conservative. Using the simplified Cassie–Baxter equation (3), ξ = 0.25. The results are given in Table 3. Table 3 Theoretical and measured contact angles of droplets on pillared surfaces. Pillar size
Initial contact angle
Theoretical contact angle
Measured contact angle using Eq. (25)
Measured contact angle using Eq. (26)
2×2 3×3
110° 110°
146.7° 146.7°
≈180°, see Fig. 5 159.1°
147.9° 149.0°
Using the Niu et al. method on hydrophobic 2 × 2 pillars resulted in a droplet that was barely in contact with the surface at all, and where it was impossible to make measurements of the contact angle. In Fig. 4 the ϕ/ρ = −0.8 isosurface is shown, the base of the spherical ϕ = 0 isosurface was above the top of the pillars. 5. Conclusions When the free energy model is used together with the bounce-back boundary condition, a value for the virtual ordering parameter ϕs in the wall node needs to be determined. The primary aim is to avoid or minimize the artificial gradients of ϕ . Such gradients may destabilize a simulation and cause serious and systemic error in the fluid flow calculations. The simplest method of assigning each wall node a fixed value of the virtual ordering parameter ϕs . This method has the non-trivial benefit of not needing to continuously update ϕs which equates to a large saving of computational effort and simulation time. The drawback is that the approach does nothing to prevent the large artificial gradients and Laplacians that can cause stability issues. A step to avoid the artificial gradients is to calculate a new value of ϕs at each time step where the value varies depending on which phase is near the wall nodes. The method provided by Niu et al. averages the ordering parameter of all surrounding fluid nodes to provide a baseline for ϕs . However, a severe limitation to this approach is that it is only valid where the wall node is not affected by different phases on different sides. If the wall node is in contact with different phases on different sides, calculating ϕs as an average of different phases gives a value that is valid nowhere and allows phases to interact through the walls. This method generated large errors of the contact angles on surfaces that had many edges and corners. This was most pronounced in the simulations where the surface details were small so that there were many edges and corners on the surface. Even when the surface details were larger (3 × 3 lattice sites) the contact angle error was around 10°. In this article we proposed a further refinement to the Niu et al. method by assigning each side and edge of the wall node a separate ordering parameter ϕs,i . Computationally this is a more costly method since many more values of ϕs need to be calculated. However, the problem with artificial gradients is considerably mitigated since the fluids on one side of the wall node no longer affects φs on another side of the wall node. Through the grid refinement study on a pillared surface, we found that the measured contact angles remained accurate for smaller grids, despite the surface details being smaller than the interface width. The method uses an average of phases at the edges of the wall nodes. This causes a small error in the measured contact angles but for structural details on the order of 3 × 3 nodes, the difference from the Cassie–Baxter prediction are only on
10
S. Stensholt / Computers and Mathematics with Applications (
)
–
the order of 3°. For the surfaces where the details (pillars or pits) were of size 2 × 2 lattice nodes, the error in contact angle was reduced to something of the order of 5°. The method makes it feasible to use finer surface details than what is possible using the Niu et al. method. There are no issues with mass conservation in the method. Cassie–Baxter type-surfaces were considered in this article but the applications extend to any system with complex boundaries, such as porous media problems. The method that was implemented still assigns a single target contact angle θp that as valid over the whole of each wall node. If more flexibility is needed, it is possible to assign different values of θp on each side of each wall node similar to how different values of ϕs were assigned to each side. Acknowledgments The research has been conducted as part of the Polytec Research Foundation’s ANTIS project, with funding from Gassco, Statoil, and the Research Council of Norway. References [1] G. Polizos, K. Winter, M.J. Lance, H.M. Meyer, B.L. Armstrong, D.A. Schaeffer, J.T. Simpson, S.R. Hunter, P.G. Datskos, Scalable superhydrophobic coatings based on fluorinated diatomaceous earth: Abrasion resistance versus particle geometry, Appl. Surf. Sci. 292 (2014) 563–569. [2] X. Zhu, Z. Zhang, X. Men, J. Yang, K. Wang, X. Xu, X. Zhou, Q. Xue, Robust superhydrophobic surfaces with mechanical durability and easy repairability, J. Math. Chem. 21 (2011) 15793–15797. [3] A. Dupuis, J.M. Yeomans, Modeling droplets on superhydrophobic surfaces: Equilibrium states and transitions, Langmuir 21 (6) (2005) 2624–2629. [4] E. Orlandini, M. Swift, J. Yeomans, A lattice Boltzmann model of binary-fluid mixtures, Europhys. Lett. 32 (6) (1995) 463–468. [5] M.R. Swift, W. Osborn, J. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75 (5) (1995) 830–834. [6] L. Hao, P. Cheng, Lattice Boltzmann simulations of liquid droplet dynamic behavior on a hydrophobic surface of a gas flow channel, J. Power Sources 190 (2) (2009) 435–446. [7] X.-D. Niu, T. Menukata, S.-A. Hyodo, K. Suga, An investigation of water–gas transport processes in the gas-diffusion-layer of a PEM fuel cell by a multiphase multiple-relaxation-time lattice Boltzmann model, J. Power Sources 172 (2007) 542–552. [8] R.N. Wenzel, Resistance of solid surfaces to wetting by water, Ind. Eng. Chem. 28 (8) (1936) 988–994. [9] A.B. Cassie, S. Baxter, Wettability of porous surfaces, Trans. Faraday Soc. 40 (1944) 546–551. [10] A.J. Milne, A. Amirfazli, The Cassie equation: How it is meant to be used, Adv. Colloid Interface Sci. 170 (2012) 48–55. [11] G. Whyman, E. Bormashenko, T. Stein, The rigorous derivation of Young, Cassie–Baxter and Wenzel equations and the analysis of contact angle hysterisis phenomenon, Chem. Phys. Lett. 450 (2008) 355–359. [12] L. Gao, T.J. McCarthy, How Wenzel and Cassie were wrong, Langmuir 23 (2007) 3762–3765. [13] G.R. McNamara, G. Zanetti, Use of the Boltzmann equation to simulate lattice–gas automata, Phys. Rev. Lett. 61 (1988) 2332–2335. [14] X. He, L.-S. Luo, A priori derivation of the lattice Boltzmann equation, Phys. Rev. E 55 (1997) R6333–R6336. [15] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997) 6811–6817. [16] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [17] A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (8) (1991) 4320–4327. [18] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, 2001. [19] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819. [20] H. Huang, D.T. Thorne, M.G. Schaap, M.C. Sukop, Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models, Phys. Rev. E 76 (2007) 066701. [21] M.R. Swift, E. Orlandini, W. Osborn, J. Yeomans, Lattice Boltzmann simulations of liquid–gas and binary fluids, Phys. Rev. E 5 (1996) 5041–5052. [22] R. van der Smaan, S. van der Graaf, Emulsion droplet deformation and breakup with lattice Boltzmann model, Comput. Phys. Comm. 178 (2008) 492–504. [23] D. D’Humieres, I. Gizburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. A 360 (2002) 437–451. [24] Q. Zou, S. Hou, G.D. Doolen, Analytical solutions of the lattice Boltzmann BGK model, J. Stat. Phys. 81 (1–2) (1995) 319–334. [25] M.A. Gallivan, D.R. Noble, J.G. Georgiadis, R.O. Buckius, An evaluation of the bounce-back boundary condition for lattice Boltzmann simulations, Internat. J. Numer. Methods Fluids 25 (1997) 249–263. [26] I. Ginzburg, D. D’Humiéres, Multireflection boundary conditions for lattice Boltzmann models, Phys. Rev. E 68 (2003) 066614. [27] A. Dupuis, J. Yeomans, Mesoscopic modelling of droplets on topologically patterned substrates, in: Computational Science-ICCS, 2004, pp. 556–563. [28] K. Connington, T. Lee, Lattice Boltzmann simulations of forced wetting transitions of drops on superhydrophobic surfaces, J. Comput. Phys. 250 (2013) 601–615. [29] T. Lee, L. Liu, Wall boundary conditions in the lattice Boltzmann equation for nonideal gases, Phys. Rev. E 78 (2008) 017702. [30] S. van der Graaf, T. Nisisako, C.G. Schroën, R.G. van der Sman, R.M. Boom, Lattice Boltzmann simulations of droplet formation in a T-Shaped microchannel, Langmuir 22 (2006) 4144–4152.
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Oriented ordering parameters for free energy lattice Boltzmann methods using the bounce-back boundary condition Sigvat Stensholt Polytec Research Institute, Sørhauggata 128, 5527 Haugesund, Norway
article
info
Article history: Available online xxxx Keywords: Lattice Boltzmann Boundary conditions Cassie–Baxter Contact angle Simulation stability
abstract We evaluate the Lattice-Boltzmann model that combines the free energy approach for multiphase flow with the bounce-back boundary condition. This method requires a virtual ordering parameter to be assigned to the wall nodes. Existing literature has assigned a single ordering parameter for the entire wall node, but this is too coarse to give accurate simulations since it allows for a spurious interaction of different phases through the wall nodes. The consequences are illustrated by examining rough and hydrophobic surfaces where contact angles can be predicted using the Cassie–Baxter equation. It is found that the single value approach gives a large systemic error in the measured contact angles. Considerable improvement can be obtained by assigning different values of the ordering parameter to each side of the wall node. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction The behavior of fluids on different forms of surfaces is a field of study with multiple scientific and industrial applications. The chemical and structural characteristics of a surface can determine whether water tends to stick to or roll off it. These attributes are significant, for example for applications that require surfaces to be as dry as possible. If a surface attracts water, the droplet will spread out giving a contact angle θ less than 90°. Such a surface is called hydrophilic. In contrast, a water droplet on a surface that repels water will curl up and have a contact angle greater than 90°. Such a surface is called hydrophobic. Very hydrophobic surfaces with contact angles greater than 150° are called superhydrophobic and are of great industrial interest. Several proposals for manufacturing superhydrophobic surfaces exist, many with the aim of achieving characteristics such as ice-phobicity and durability, see for example [1,2]. Many of these publications report on laboratory testing of the results, but fluid dynamics simulations by using mesoscopic methods such as the Lattice Boltzmann method have also been applied [3]. Laboratory tests require the production of a physical test surface, something that is not needed for computer simulations. In this article we examine the feasibility of using the popular bounce-back boundary condition for the solid–fluid interface when using the free energy approach [4,5] for multiphase flow. We aim to improve on the method offered in the previous literature [6,7] and show that this gives simulated contact angles that are more consistent with classical fluid theory. This is applicable to any domain with a complex solid boundary. The classical fluid theory for rough surfaces is covered in Section 2. The Lattice Boltzmann method and the associated boundary conditions that are the focus for this article is covered in Section 3. Simulation results are covered in Section 4, and the conclusions are drawn in Section 5. E-mail address: [email protected]. http://dx.doi.org/10.1016/j.camwa.2015.06.014 0898-1221/© 2015 Elsevier Ltd. All rights reserved.
2
S. Stensholt / Computers and Mathematics with Applications (
)
–
2. Fluid theory of rough surfaces A roughness is required to make a surface superhydrophobic. The contact angle of droplets on heterogeneous and rough surfaces has been the subject of considerable study. The simplest equations for estimating the contact angle on rough surfaces are the Wenzel [8] and Cassie–Baxter [9] equations. The Wenzel equation applies when a droplet fills the cavities of a rough surface. The apparent contact angle θc is then given by cos θc = r cos θ0 ,
(1)
where θ0 is the contact angle the material would give on a smooth surface and the roughness parameter r is the ratio between the actual surface area and the projected surface area. In its original form, the Cassie–Baxter equation applies on a smooth surface with n components with surface fractions ξ1 , ξ2 , . . . , ξn with different contact angle parameters θ1 , θ2 , . . . , θn . For a smooth surface made of different components, the sum over surface fractions will be 1. The resultant contact angle θc on a composite surface is given by cos θc =
n
ξi cos θi .
(2)
i=1
The Cassie–Baxter equation has also been applied to some rough surfaces where the droplet is suspended above the cavities. In this case air, which has a contact angle parameter of 180°, is treated as a separate component. If the surface fraction of the solid component is ξ , one obtains simplified version of Eq. (2) which has been used in a majority of articles that deal with contact angles on a rough surface [10]: cos θc = ξ cos θ − (1 − ξ ).
(3)
Eq. (3) assumes that the solid part of the surface is made up of a single component with an initial contact angle θ that covers a fraction ξ of the surface, and that the droplet does not penetrate into the surface at all. While these assumptions are commonly made, real droplets often exhibit more complex behavior by penetrating somewhat into the surface in which case Eq. (2) should be used without the assumption that the sum over all ξi is equal to 1. As the droplet penetrates deeper into the surface ξsolid approaches the roughness factor and ξair approaches 0. Consequently, Eq. (2) approaches the Wenzel equation (1) [10]. The Cassie–Baxter and Wenzel equations can be derived from thermodynamic analysis [11] and provide a means for analytical predictions of the contact angle for simple cases of stationary droplets that are entirely on top of the surface or entirely penetrated into the surface. However the validity of the equations has been the subject of skepticism [12] and the fluid dynamics for the transitional states or droplets in motion cannot be described fully by the equations. There is therefore an interest in a computational fluid dynamic approach to studying such surfaces and the scales involved invite the use of mesoscopic approaches. 3. Lattice Boltzmann method The Lattice Boltzmann method is a computational fluid dynamics (CFD) model originally introduced by Mcnamara and Zanetti [13] and later interpreted as a discretization of the Boltzmann equation [14,15]. The method is commonly used on a mesoscopic scale and based on the discrete particle distribution function fi (x, t ) which gives the number of particles moving at a given velocity at site x at time t. For two-component flow, two distribution functions are needed. One option is to have a separate distribution function for each component. The other option, used here, is to let the first distribution function fi (x, t ) account for the sum of the components while the second distribution function gi (x, t ) accounts for the difference between them. The distribution functions evolve according to the equations fi (x + ei δ t , t + δ t ) = fi (x, t ) + Ωf ,i + Fi
(4)
gi (x + ei δ t , t + δ t ) = gi (x, t ) + Ωg ,i
(5)
Fi accounts for any external forces such as gravity. Ωi is the collision operator where we have used the single relaxation time operator [16],
Ω f ,i =
fi − fi
τf
eq
eq
,
Ωg ,i =
gi − gi
τg
(6)
where the superscript eq represents an equilibrium distribution and f and g are the relaxation times. The relaxation time for fi gives the kinematic viscosity ν by
ν=
2τf − 1 6
δt .
(7)
S. Stensholt / Computers and Mathematics with Applications (
)
–
3
Space is discretized into a regular cubic lattice. We have opted to use the D3Q19 lattice where particles may move along 19 possible vectors ei . The lattice spacing δ x is set to one. e0 = (0, 0, 0) e1−6 = (±1, 0, 0), e7−18
(0, ±1, 0), (0, 0, ±1) = (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1).
(8)
1 1 , and w7−18 = 36 . The corresponding weight coefficients wi for this lattice are w0 = 13 , w1−6 = 18 The mass, momentum, and ordering parameters at each site can be calculated from the moments
ρ(x, t ) =
f i ( x, t )
(9)
i
ρ(x, t )ui (x, t ) =
fi (x, t )ei
(10)
i
ϕ(x, t ) =
gi (x, t ).
(11)
i
A number of methods to simulate multiphase flow within the lattice Boltzmann framework have been proposed. The oldest method is the color model [17] which redistributes phases to enforce phase separation and imposes a surface tension force. The method provides fairly sharp interfaces and an explicit contact angle parameter to handle wetting phenomena, but it is known to possess a measure of anisotropy and its computational expense is a drawback that becomes serious for three-dimensional simulations [18]. The pseudo-potential method [19] is a simpler method that imposes a force between opposite phases. It is used extensively, but it produces diffuse interfaces and even though the contact angle arising from the model has been studied [20], a priori determination of the contact angle remains difficult. For multiphase flow we have instead employed the free energy approach which was first developed by researchers at Oxford [4,5,21]. This method has previously been used for simulating superhydrophobic surfaces [3] and has the advantage of being thermodynamically consistent. Another significant advantage is the ability to impose an explicit contact angle on a smooth homogenous surface through a single contact angle parameter. Explicit expressions for the free energy model can be found in [22] and for the three-dimensional D3Q19 lattice in [6], which largely uses the method described by [7]. A difference is that Hao and Cheng used the BGK operator (6) while Niu et al. used the multiple relaxation time (MRT) operator [23] which offers more degrees of freedom but is computationally more demanding than the BGK operator. Central to the understanding of the free energy model is the free energy density ψ(ϕ) given by A
B
2
4
ψ(ϕ) = − ϕ 2 + ϕ 4 .
(12)
The parameters A and B are positive and the free energy density function is therefore a quartic function that form a double well potential. The chemical potential µ is given by
µ = −Aϕ + Bϕ 3 + κ∇ 2 ϕ.
(13)
The free energy model uses an explicit thermodynamic pressure tensor Pth which can be written out as the sum of the chemical pressure tensor Pch and hydrostatic pressure P hy giving1 Pth = Pch + P hy
= p0 I + κ(∇ϕ)(∇ϕ)T + ρ I κ = ϕµ − ψ(ϕ) − (∇ϕ)2 + κ(∇ϕ)(∇ϕ)T + ρ I 2
A
=− ϕ + 2
2
3B 4
ϕ − κϕ∇ 2 ϕ − 4
κ 2
(∇ϕ)2 + κ(∇ϕ)(∇ϕ)T + ρ I.
(14)
1 There are typographical errors in some of the equations in Hao and Cheng’s 2009 article. Most notably, the term for p which should be given as 0 p0 = ϕµ − ψ(ϕ) −
κ 2
A
3B
2
4
(∇ϕ)2 = − ϕ 2 +
ϕ 4 − κϕ∇ 2 ϕ −
κ 2
(∇ϕ)2
has substituted the quadratic term 3B ϕ 2 for the correct quartic term, thus destroying the double well potential. Also, µ and κ have been rendered as u and 4 k. The corresponding author of the article has been contacted and indicated that an erratum will be published.
4
S. Stensholt / Computers and Mathematics with Applications (
)
–
The gradient (∇ ) and Laplacian (∇ 2 ) operators are conducted using a 19-point stencil given by [7]. The numerical estimate for the jth component of the gradient is given by 3
18
wi ϕ(x + ei δ t )eij
i =0
∇j ϕ(x) =
δt
18
,
(15)
wi (eij δ t )2
i=0
while the Laplacian is given by 2
∇ ϕ= 2
18
wi (ϕ(x + ei δ t ) − ϕ(x))
i=0 18
.
(16)
wi (eij δ t )2
i =0
The main feature of the free energy model is the expressions for the equilibrium distributions that depend on the pressure tensor and chemical potentials, that are dependent on ϕ and derivatives of ϕ . For the free energy model on a D3Q19 lattice, the explicit expressions for the equilibrium distributions are fi
eq
eq
gi
= Qi + wi = Hi + wi
ei · u cs2
ei · u cs2
+
(ei · u)2
+
(ei · u)2
2cs4 2cs4
− −
u·u
2cs2 u·u 2cs2
,
(17)
,
(18)
√
where cs = c / 3 is the speed of sound and the coefficients Qi and Hi are given by
th th th − Pxx − Pxx )/3c 2 ρ − 2(Pxx th th th Qi = eTi P th ei /4c 4 − (Pxx − Pxx − Pxx )/36c 2 T th th th th ei P ei /8c 4 − (Pxx − Pxx − Pxx )/18c 2 ϕ − 2Γ µ/c 2 i = 0; Hi = Γ µ/6c 2 i = 1, 2, . . . , 6; Γ µ/12c 2 i = 7, 8, . . . , 18.
i = 0; i = 1, 2, . . . , 6; i = 7, 8, . . . , 18,
(19)
(20)
The parameters A, B, and κ that all affect the phase separation equations have already been introduced. The expression for Hi introduced a fourth parameter, Γ , that is related to the mobility. The bulk pressure consists of the ideal component (cs2 ρ ) and the components due to interfaces. It can be written as p = cs2 ρ −
A 2
ϕ2 +
3B 4
ϕ 4 − κϕ∇ 2 ϕ −
κ 2
(∇ϕ)2 .
(21)
3.1. Boundary conditions In our simulations we have employed periodic boundary conditions at the fluid node edges of the domain. For the fluid–solid interface the boundary conditions are required to handle both the physical properties (nonpenetration of the wall and no-slip) and the chemical properties affecting contact angles. The simplest method for handling the solid interface nodes is the bounce-back condition. In this approach any particles that would propagate into the wall node are instead reflected in the opposite direction. This method ensures that no element of the distribution function will remain unknown. The method gives a no-slip boundary condition consistent with a wall placed halfway between the last fluid node and the first wall node [24]. The accuracy of the bounce-back condition has been studied extensively. A particular result to note is that the accuracy is dependent on the relaxation time τ , being more accurate for smaller relaxation times (i.e. less viscous flows), and with an exponential deterioration in accuracy for τ > 1 [25]. The position of the wall may also be affected by τ [26]. Other boundary conditions are possible. The first article to use lattice Boltzmann to simulate superhydrophobic surfaces [27] used a set of velocity boundary conditions were to determine the missing elements of the distribution function near the fluid–wall boundary. The bounce-back condition can be modified so that the reflection to occur on the node instead of between nodes, and this approach has also been applied for studies of superhydrophobic surfaces [28]. Some wall boundary conditions result in artificial parasitic currents at the contact line, but these are eliminated if intermolecular force is handled using the potential form [29].
S. Stensholt / Computers and Mathematics with Applications (
)
–
5
Fig. 1. Wall separating ‘‘red’’ and ‘‘blue’’ phases. In (a) the wall consists of a single row of wall nodes. In (b) the wall consists of a double row of wall nodes. Colors on the wall nodes indicate the value of ϕs assigned to it. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In order to allow the gradient and Laplacian of the ordering parameter to be calculated a virtual ordering parameter for wall nodes, ϕs is needed. There are various approaches to this. The simplest approach is to simply assign a fixed value of ϕs for each wall node. The method has been used by [30] who give the resultant contact angle θ as 3
cos θ =
2
ϕ˜ 1 −
ϕ˜ 2
3
,
(22)
1/2
where ϕ˜ = ϕs / AB . The simplicity of not needing to update the values of ϕs at every time step is an attractive feature of the model. This approach does not account for the ϕ in the surrounding fluid nodes, and will therefore give an artificial gradient and Laplacian around the walls. To avoid the artificial gradients the value of ϕs needs to change with the surrounding fluid. The approach by [7] uses a target contact angle θp . This parameter is used to determine a virtual ϕs for the wall nodes. A wetting potential Θ is given by
Θ = 2 sgn
π
− θp
2
cos
arccos(sin2 θp ) 3
1 − cos
arccos(sin2 θp )
3
(23)
where all angles are given in radians. This wetting potential is used to calculate the parameter λ:
λ = −Θ
κA 2
.
(24)
The method given in [7] approximates the order parameter for the wall nodes ϕs as the average of all the N surrounding fluid nodes, modified by the contact angle parameter and surface tension parameters λ and κ .
ϕs =
(ϕf − δ xλ/κ)
N
N
.
(25)
This value of ϕs can then be used along with the values of ϕf for the fluid nodes so that all nodes are assigned a value of ϕ . One can then proceed with the calculations of ∇ϕ and ∇ 2 ϕ in Eqs. (15) and (16) for the relaxation step of the method. Hao and Cheng [6] write that Eq. (25) ‘‘can easily be applied to complex solid boundary such as porous media’’. However, we observe that averaging over all the surrounding nodes leads to error in cases where two different phases are in contact with the same wall node. In practice this includes any complex wall with corners or edges, as well as any wall node that separates two fluids of different phase. To illustrate the problem, consider a wall separating to different fluid phases. In Fig. 1, a single row of wall nodes, as shown in (a), is in contact with both the red phase with ϕ = 1 and blue phase with ϕ = −1, and so the assigned ordering parameter for the wall nodes is ϕs = 0 − δ xλ/κ , which is incorrect for both sides of the wall. The thin wall allows the unphysical interaction of phases through the wall. The ordering parameter of the wall is correct if a double row of wall nodes is used as seen in (b). In this case each wall node is only in contact with one phase so the wall nodes in contact with the red phase will have an ordering parameter ϕs = 1 − δ xλ/κ while those in contact with the blue phase will have ϕs = −1 − δ xλ/κ , which retains only the desired gradient and Laplacian from the wetting factor δ xλ/κ . A similar problem arises for rough surfaces that feature edges and corners. In Fig. 2 we see that the corner wall node that is in contact with both phases has its value for ϕs affected by both phases and this value is not valid on any side. The problem with incorrect values of ϕs cannot be entirely eliminated in this case, even though increasing the resolution reduces the problem by decreasing the proportion of corner wall nodes.
6
S. Stensholt / Computers and Mathematics with Applications (
)
–
Fig. 2. Details of a fluid interface at a wall edge. x1 and x2 are fluid nodes in the red (ϕ = 1) and blue (ϕ = −1) phases respectively. The ϕs value for each wall node in contact with a fluid phase is noted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In order to avoid introducing an artificial gradient, the edge wall node in Fig. 2 needs to have a different values of ϕs on each side. A significant improvement is to give the ordering parameter of the wall an orientation. Each side and edge of the wall node is assigned an ordering parameter ϕs,i , where subscript i indicates the side or edge of the wall node.
ϕs,i =
(ϕf (x + ei ) − δ xλ/κ)
i∈Eˆ i
,
Nˆ
(26)
where Eˆ i is the set of three (for i = 7, . . . , 18) or five (for i = 1 . . . 6) non-zero indices ˆi so that eˆi − ei ≤ 0. Nˆ is the number of fluid nodes adjacent to the face or edge of the wall node. For example, on the i = 1 surface, which is oriented in the direction (1, 0, 0), the values ϕf (x +(1, 0, 0)), ϕf (x +(1, 1, 0)), ϕf (x + (1, −1, 0)), ϕf (x + (1, 0, 1)), and ϕf (x + (1, 0, −1)) are used in the calculation for ϕs,1 . To assign a value for the i = 7 edge, which is oriented in the direction (1, 1, 0), ϕs,7 is calculated using the fluid node order parameters ϕf (x + (1, 1, 0)), ϕf (x + (1, 0, 0)), and ϕf (x + (0, 1, 0)). The different value of ϕs on different sides of the wall nodes is reflected in a small modification to Eqs. (15) and (16) for the gradient and Laplacian. These are now given by
3
∇j ϕ(x) =
18
wi ϕi* (x + ei δ t )eij
i =0
δt
18
, wi (eij δ t )
(27)
2
i=0
while the Laplacian is given by 2
∇ 2ϕ =
18
wi (ϕi* (x + ei δ t ) − ϕ(x))
i=0 18
(28)
wi (eij δ t )
2
i=0
where i* is the index opposite i, that is ei* = −ei . For fluid nodes, ϕi = ϕ as defined in Eq. (11) is valid for all directions i. In the problem illustrated by Fig. 1, a single row of wall nodes is now sufficient to simulate a separating wall as the top side of the wall will have a ‘‘blue’’ ordering parameter while the bottom side will be ‘‘red’’. For the corner node in Fig. 2 the top side of the node will have a ‘‘red’’ ordering parameter while the vertical side on the right will have a ‘‘blue’’ one. 4. Simulations and results Lattice Boltzmann simulations of sessile droplets on pillared surfaces have been conducted and the resulting contact angles have been measured. Simulations have been conducted with schemes using the original non-oriented ϕs (x, t ) and the proposed oriented ϕs,i (x, t ) for the wall nodes so that the results can be compared.
S. Stensholt / Computers and Mathematics with Applications (
)
–
7
Fig. 3. Top-down view of 2 × 2 pitted (left) and pillared (right) surfaces.
The contact angle parameter θ0 was set to 110°. The aim in this paper is to consider the theoretical effects of modifying the surface rather than simulating real fluids. The relaxation parameters have been set for maximum simulation stability, τf = τg = 1. The interaction parameters are A = B = 0.007, κ = 0.0078, and Γ = 1. For best reliability, the dimensions of the surface details should be larger than the interface width. If this condition is violated, a simulation can still be performed but the results are less reliable. Nonetheless, computational constraints limit the size of the domain and the achievable fineness of the surface details. It is therefore of interest to investigate the reliability of simulations where the details are smaller than recommended. The smallest details (pillars or pits) are of size 2 × 2. Any detail which leaves a fluid node sandwiched between two wall nodes is prone to destabilize the simulation. On the other hand, the chosen surface details are not so far apart that the droplets will penetrate appreciably into the gaps between the pillars. This allows us to compare the results with the simplified Cassie–Baxter equation (3). We have run our simulations on pitted surfaces and pillared surfaces where the pits and pillars cover a quarter of the top surface area. A top-down view of the arrangement of the 2 × 2 surface details is shown in Fig. 3. The arrangement for the 3 × 3 and larger surface details are similar. The interface produced by the free energy model is diffuse, but the ϕ = 0 isosurface can be used to define its location. After several time steps the droplets converge into a final shape. We examined cross sections through the center of the droplets and found, by fitting a circle over the ϕ = 0 contours, that all the droplets assumed a shape of a spherical cap within the degree of precision allowed by the fineness of the lattice. This realization allows us to calculate the contact angle θ by measuring the radius α of the circular contact area and the droplet height β as shown in Fig. 4: 2αβ . (29) β 2 − α2 The measurements α and β are taken on the fitted sphere half a node beyond the last fluid node. This is because the boundary
θ = arctan
is between the last fluid node and the first wall node. 4.1. Fixed value of ordering parameters For the rough pillared and pitted surfaces, simulations using the fixed values of ϕs were unstable with the selected input parameters. Typically, the simulation broke down after about 300 time steps illustrating the problem of using this method on rough surfaces. 4.2. Grid refinement and convergence tests of proposed method A series of tests were run to determine the effects of grid refinement and the evolution of the droplet shape and contact angle when using Eq. (26). The finest domain features a grid 100 nodes long and wide, and 70 nodes tall. The floor consists of a base and pillars that have a 6 node × 6 node cross section and a height of 4 nodes. The solid surface fraction ξ = 0.25. The initial contact angle θ0 is 110°. According to Eq. (3), the theoretical contact angle of a droplet in the Cassie–Baxter state is then 146.7°. A droplet of radius r is initiated to rest atop the pillars and the simulation is started. Tests have also been run on three coarser grids where the pillar widths, and droplet radii are scaled down proportionally. The pillar height was not scaled down however. If a pillar is only one or two nodes tall, the interface width will cause the droplet to make contact with the surface base, causing the simulated droplet to assume a Wenzel state. Therefore, the pillar height is kept at four nodes for all domain sizes. The results are given in Table 1.
8
S. Stensholt / Computers and Mathematics with Applications (
)
–
Fig. 4. Measurements used to determine the contact angle θ .
Fig. 5. Diagram of a simulated droplet on a pillared surface. The droplet surface indicates the ϕ/ρ = −0.8 isosurface instead of the ϕ = 0 isosurface which is inside this surface and entirely over the top of the pillars. Table 1 Contact angles of a sessile droplet at different time steps and varying degrees of resolution. Domain size
Pillar width
Droplet radius r
Time steps 10 000
20 000
30 000
40 000
50 000
100 × 100 × 70 84 × 84 × 60 67 × 67 × 50 50 × 50 × 42
6×6 5×5 4×4 3×3
35.0 29.4 23.5 17.5
123.1° 127.3° 133.2° 139.7°
142.1° 140.6° 144.3° 143.7°
146.2° 146.5° 146.1° 144.7°
147.4° 147.3° 147.1° 145.5°
147.8° 147.1° 146.2° 146.3°
The maintenance of the spherical shape of the droplet is apparent already at 10 000 time steps, but about 40 000 time steps are required before the contact angle converges to a figure where the change between 10 000 time steps is less than 1°. Since the contact angle is greater than 110°, the contact area between the droplet and the surface is small, and only
S. Stensholt / Computers and Mathematics with Applications (
)
–
9
covers the tops of seven pillars. Despite this limited contact area, the resultant contact angles are in close agreement with the theory of Cassie and Baxter. 4.3. 2 × 2 and 3 × 3 pitted surfaces These surfaces contain square pits carved into an otherwise flat surface. The pits have a size of 2 × 2 or 3 × 3 and spaced so that they cover 25% of the surface area. Using the simplified Cassie–Baxter Equation (3), ξ = 0.75. The results are given in Table 2. Table 2 Theoretical and measured contact angles of droplets on pitted surfaces. Pit size
Initial contact angle
Theoretical contact angle
Measured contact angle using Eq. (25)
Measured contact angle using Eq. (26)
2×2 3×3
110° 110°
120.4° 120.4°
140.4° 132.9°
126.8° 122.0°
4.4. 2 × 2 and 3 × 3 pillared surfaces These surfaces contain square, flat-topped pillars imposed over the surrounding surface. The pillars have a size of 2 × 2 or 3 × 3 and spaced so that they cover 25% of the surface area. All simulations were conducted in Matlab on a cubic D3Q19 lattice with dimensions 80 × 80 × 80 for 50 000 time steps. In all simulations the fluid mass of both phases remained constant throughout the simulation. The methods are all mass conservative. Using the simplified Cassie–Baxter equation (3), ξ = 0.25. The results are given in Table 3. Table 3 Theoretical and measured contact angles of droplets on pillared surfaces. Pillar size
Initial contact angle
Theoretical contact angle
Measured contact angle using Eq. (25)
Measured contact angle using Eq. (26)
2×2 3×3
110° 110°
146.7° 146.7°
≈180°, see Fig. 5 159.1°
147.9° 149.0°
Using the Niu et al. method on hydrophobic 2 × 2 pillars resulted in a droplet that was barely in contact with the surface at all, and where it was impossible to make measurements of the contact angle. In Fig. 4 the ϕ/ρ = −0.8 isosurface is shown, the base of the spherical ϕ = 0 isosurface was above the top of the pillars. 5. Conclusions When the free energy model is used together with the bounce-back boundary condition, a value for the virtual ordering parameter ϕs in the wall node needs to be determined. The primary aim is to avoid or minimize the artificial gradients of ϕ . Such gradients may destabilize a simulation and cause serious and systemic error in the fluid flow calculations. The simplest method of assigning each wall node a fixed value of the virtual ordering parameter ϕs . This method has the non-trivial benefit of not needing to continuously update ϕs which equates to a large saving of computational effort and simulation time. The drawback is that the approach does nothing to prevent the large artificial gradients and Laplacians that can cause stability issues. A step to avoid the artificial gradients is to calculate a new value of ϕs at each time step where the value varies depending on which phase is near the wall nodes. The method provided by Niu et al. averages the ordering parameter of all surrounding fluid nodes to provide a baseline for ϕs . However, a severe limitation to this approach is that it is only valid where the wall node is not affected by different phases on different sides. If the wall node is in contact with different phases on different sides, calculating ϕs as an average of different phases gives a value that is valid nowhere and allows phases to interact through the walls. This method generated large errors of the contact angles on surfaces that had many edges and corners. This was most pronounced in the simulations where the surface details were small so that there were many edges and corners on the surface. Even when the surface details were larger (3 × 3 lattice sites) the contact angle error was around 10°. In this article we proposed a further refinement to the Niu et al. method by assigning each side and edge of the wall node a separate ordering parameter ϕs,i . Computationally this is a more costly method since many more values of ϕs need to be calculated. However, the problem with artificial gradients is considerably mitigated since the fluids on one side of the wall node no longer affects φs on another side of the wall node. Through the grid refinement study on a pillared surface, we found that the measured contact angles remained accurate for smaller grids, despite the surface details being smaller than the interface width. The method uses an average of phases at the edges of the wall nodes. This causes a small error in the measured contact angles but for structural details on the order of 3 × 3 nodes, the difference from the Cassie–Baxter prediction are only on
10
S. Stensholt / Computers and Mathematics with Applications (
)
–
the order of 3°. For the surfaces where the details (pillars or pits) were of size 2 × 2 lattice nodes, the error in contact angle was reduced to something of the order of 5°. The method makes it feasible to use finer surface details than what is possible using the Niu et al. method. There are no issues with mass conservation in the method. Cassie–Baxter type-surfaces were considered in this article but the applications extend to any system with complex boundaries, such as porous media problems. The method that was implemented still assigns a single target contact angle θp that as valid over the whole of each wall node. If more flexibility is needed, it is possible to assign different values of θp on each side of each wall node similar to how different values of ϕs were assigned to each side. Acknowledgments The research has been conducted as part of the Polytec Research Foundation’s ANTIS project, with funding from Gassco, Statoil, and the Research Council of Norway. References [1] G. Polizos, K. Winter, M.J. Lance, H.M. Meyer, B.L. Armstrong, D.A. Schaeffer, J.T. Simpson, S.R. Hunter, P.G. Datskos, Scalable superhydrophobic coatings based on fluorinated diatomaceous earth: Abrasion resistance versus particle geometry, Appl. Surf. Sci. 292 (2014) 563–569. [2] X. Zhu, Z. Zhang, X. Men, J. Yang, K. Wang, X. Xu, X. Zhou, Q. Xue, Robust superhydrophobic surfaces with mechanical durability and easy repairability, J. Math. Chem. 21 (2011) 15793–15797. [3] A. Dupuis, J.M. Yeomans, Modeling droplets on superhydrophobic surfaces: Equilibrium states and transitions, Langmuir 21 (6) (2005) 2624–2629. [4] E. Orlandini, M. Swift, J. Yeomans, A lattice Boltzmann model of binary-fluid mixtures, Europhys. Lett. 32 (6) (1995) 463–468. [5] M.R. Swift, W. Osborn, J. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75 (5) (1995) 830–834. [6] L. Hao, P. Cheng, Lattice Boltzmann simulations of liquid droplet dynamic behavior on a hydrophobic surface of a gas flow channel, J. Power Sources 190 (2) (2009) 435–446. [7] X.-D. Niu, T. Menukata, S.-A. Hyodo, K. Suga, An investigation of water–gas transport processes in the gas-diffusion-layer of a PEM fuel cell by a multiphase multiple-relaxation-time lattice Boltzmann model, J. Power Sources 172 (2007) 542–552. [8] R.N. Wenzel, Resistance of solid surfaces to wetting by water, Ind. Eng. Chem. 28 (8) (1936) 988–994. [9] A.B. Cassie, S. Baxter, Wettability of porous surfaces, Trans. Faraday Soc. 40 (1944) 546–551. [10] A.J. Milne, A. Amirfazli, The Cassie equation: How it is meant to be used, Adv. Colloid Interface Sci. 170 (2012) 48–55. [11] G. Whyman, E. Bormashenko, T. Stein, The rigorous derivation of Young, Cassie–Baxter and Wenzel equations and the analysis of contact angle hysterisis phenomenon, Chem. Phys. Lett. 450 (2008) 355–359. [12] L. Gao, T.J. McCarthy, How Wenzel and Cassie were wrong, Langmuir 23 (2007) 3762–3765. [13] G.R. McNamara, G. Zanetti, Use of the Boltzmann equation to simulate lattice–gas automata, Phys. Rev. Lett. 61 (1988) 2332–2335. [14] X. He, L.-S. Luo, A priori derivation of the lattice Boltzmann equation, Phys. Rev. E 55 (1997) R6333–R6336. [15] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997) 6811–6817. [16] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [17] A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (8) (1991) 4320–4327. [18] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, 2001. [19] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819. [20] H. Huang, D.T. Thorne, M.G. Schaap, M.C. Sukop, Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models, Phys. Rev. E 76 (2007) 066701. [21] M.R. Swift, E. Orlandini, W. Osborn, J. Yeomans, Lattice Boltzmann simulations of liquid–gas and binary fluids, Phys. Rev. E 5 (1996) 5041–5052. [22] R. van der Smaan, S. van der Graaf, Emulsion droplet deformation and breakup with lattice Boltzmann model, Comput. Phys. Comm. 178 (2008) 492–504. [23] D. D’Humieres, I. Gizburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. A 360 (2002) 437–451. [24] Q. Zou, S. Hou, G.D. Doolen, Analytical solutions of the lattice Boltzmann BGK model, J. Stat. Phys. 81 (1–2) (1995) 319–334. [25] M.A. Gallivan, D.R. Noble, J.G. Georgiadis, R.O. Buckius, An evaluation of the bounce-back boundary condition for lattice Boltzmann simulations, Internat. J. Numer. Methods Fluids 25 (1997) 249–263. [26] I. Ginzburg, D. D’Humiéres, Multireflection boundary conditions for lattice Boltzmann models, Phys. Rev. E 68 (2003) 066614. [27] A. Dupuis, J. Yeomans, Mesoscopic modelling of droplets on topologically patterned substrates, in: Computational Science-ICCS, 2004, pp. 556–563. [28] K. Connington, T. Lee, Lattice Boltzmann simulations of forced wetting transitions of drops on superhydrophobic surfaces, J. Comput. Phys. 250 (2013) 601–615. [29] T. Lee, L. Liu, Wall boundary conditions in the lattice Boltzmann equation for nonideal gases, Phys. Rev. E 78 (2008) 017702. [30] S. van der Graaf, T. Nisisako, C.G. Schroën, R.G. van der Sman, R.M. Boom, Lattice Boltzmann simulations of droplet formation in a T-Shaped microchannel, Langmuir 22 (2006) 4144–4152.