Article pubs.acs.org/Langmuir
Adaptive Phase-Field Modeling of Anisotropic Wetting with Line Tension at the Triple Junction S.Y. Yeh and C.W. Lan*
Downloaded by RUTGERS UNIV on August 29, 2015 | http://pubs.acs.org Publication Date (Web): August 21, 2015 | doi: 10.1021/acs.langmuir.5b02175
Department of Chemical Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan ABSTRACT: Line tension could affect the contact angle at triple junction, especially in micro- to nanoscale wetting. We have developed an adaptive phase-field model to consider the line tension quantitatively. This model is coupled to the smoothed boundary method for treating the contact line with the solid phase, while the volume constraint is imposed. Our calculated contact angles are in good agreement with the modified Young’s equation. Further examples are illustrated for the anisotropic wetting on hydrophilic/hydrophobic stripes and rectangular grooves.
1. INTRODUCTION The line tension plays an important role on heterogeneous nucleation,1 nanodrop wetting,2 grain boundary grooving,3,4 anisotropic wetting on grooved surface,5,6 and so on. When the wetting dimension is small, the line tension could affect the contact angle ξ,7,8 as described by the generalized Young’s equation as cos ξ = cos θc −
γslv 1 1 dγslv − γlv rb γlv drb
The control of such anisotropic wetting is becoming important in many applications, such as microfluidics and nanofluidic devices.12 We take the anisotropic wetting on stripes in Figure 1a as an example. The contact angle in the parallel view (Figure 1b) is very different from the perpendicular one (Figure 1c) due to different contact angels. The aspect ratio (L/W) departs from unity due to the anisotropic wetting. Because the liquid drop or the strip spacing is near the line tension length, the wetting phenomena would be affected. A few approaches13−17 have been reported to simulate the equilibrium shape of a liquid drop during anisotropic wetting. Among them, the phase-field model is the most versatile one. Unlike the thermodynamics approaches,17 in addition to the simulation of equilibrium shapes, the phase-field model could be easily coupled to Navior−Stokes equation,18 as well as heat and mass transfer equations,19 for the simulation of dynamic drops. Furthermore, even drop splitting and merging that have complex morphological changes could also be considered by the phase-field model without the limit of mesh deformation in the front tracking techniques. Recently, Huang et al.20 compared different types of boundary conditions for the contact angles. Warren et al.21 also proposed a method to model the contact angle by using a three-phase model. To deal with the complex contact surface, the smoothed boundary method22 was also proposed by introducing a stationary phase, where the contact angle with the phase field was taken into account. Yu et al.23 further extended the smoothed boundary method to complex geometries. Badillo18 proposed a new phase-field model by changing the form of chemical potentials, and the model could capture the detailed experimental
(1)
where θc is the Young’s contact angle, γslv is the line energy at the triple junction line, γlv is the liquid−vapor interfacial energy, and rb is the base radius. The Young’s contact angle θc satisfies cos θc = (γsv − γlv)/γlv, where γsv and γsl are the solid−vapor and solid−liquid interfacial energies, respectively. When the drop radius is larger than the thickness of contact region except for small contact angles, the curvature effect of line tension might be neglected.8 Moreover, such data are still not available. Therefore, the modified Young’s equation is often recognized as cos ξ = cos θc −
γslv 1 γlv rb
(2)
1,5
As a result, most people determined the line tension by eq 2, and the magnitude of line tension varied in a wide range, from 10−5 to 10−11 J/m.1,5 The sign of line tension could be positive or negative9 based on the measurement, but the physical meaning of the negative line tension is not yet clear. When a drop wets a surface with grooves5,6 or patterns having hydrophobic and hydrophilic regions,10,11 as illustrated in Figure 1a, where the gray region is hydrophobic and the white region is hydrophilic, for example, the wetting phenomena are very different from the homogeneous one. © XXXX American Chemical Society
Received: June 15, 2015 Revised: August 9, 2015
A
DOI: 10.1021/acs.langmuir.5b02175 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
model considering the line tension would be necessary to model a realistic anisotropic wetting on a complex surface; however, so far no such a model has been proposed. In addition to the models, efficient computations would also be necessary. To have good numerical resolutions at the interface and the contact line, the adaptive mesh refinement (AMR)27,28 would also be necessary. We propose an adaptive phase-field model with line tension at the triple junction line based on AMR for the anisotropic wetting problems. In addition to the validation of the calculated wetting angle by the modified Young’s equation, some results are compared with the experimental observations.5 The phase-field model and the adaptive numerical scheme are described briefly in Section 2. The validation of the present phase-field model with analytical solutions is presented in Section 3. Section 4 is devoted to the simulation of anisotropic wetting on single and multiple stripes, as well as grooves, before short conclusions are given in Section 5. Downloaded by RUTGERS UNIV on August 29, 2015 | http://pubs.acs.org Publication Date (Web): August 21, 2015 | doi: 10.1021/acs.langmuir.5b02175
2. ADAPTIVE PHASE-FIELD MODEL The present model is based on the smoothed boundary method23 proposed by Yu et al. As previously mentioned, the smoothed boundary method is suitable to model the drop wetting on a complex substrate. The governing equation of the smoothed boundary method for the phase-field model proposed by Yu et al.23 could be rewritten as the following τ0
⎛ W2 ⎞ W2 ∂ϕ = ⎜ 0 ∇·(ψ ∇ϕ) + fϕ ⎟ + 0 |∇ψ ||∇ϕ|cos θc ∂t ψ ⎝ ψ ⎠ (3)
where ϕ is the phase field variable; the liquid phase is located at ϕ = 1, the vapor phase at ϕ = −1, and the liquid−vapor interface at ϕ = 0. In addition, τ0 is the time constant, t is the time, and W0 is the liquid−vapor interface thickness in the phase field. In addition, ψ is the stationary phase variable used in the smoothed boundary method; the solid phase is at ψ = 0, the liquid phase or the vapor phase at ψ = 1, and the solid− liquid or solid−vapor interface at ψ = 0.5. The solid−liquid− vapor triple junction is located at ϕ = 0 and ψ = 0.5. In addition, fϕ is the derivative of double-well function f(ϕ) with respect to ϕ, that is, fϕ = ϕ(1 − ϕ2). The right-hand side could be viewed as the driving force for the wetting, that is, the gradient of chemical potential. As a steady state is reached, the equilibrium shape is obtained. The time constant τ0 on the left is only chosen for numerical stability. Moreover, the stationary field variable ψ is designed for the smooth boundary. It is similar to the phase-field variable and is given by ψ = [1 + tanh(n/√2ζ)]/2, where ζ is the interface thickness of the smooth boundary; n is the distance normal to the boundary surface. The volume of the drop is conserved to model the drop wetting. The volume constraint term u(t) could be added into the model as previously described29 as
Figure 1. Schematic of an anisotropic wetting on multiple stripes (gray stripes are hydrophobic): (a) top view; (b) parallel (the drop width is W); and (c) perpendicular view (the drop length is L).
phenomena with good agreements; however, these phase-field models did not consider the line tension effects. Recently, Nestler et al.24 proposed a multiple phase-field model to simulate the nucleation during solidification. They tried to consider the line tension effects by using an isotropic interfacial energy; however, their results were not consistent with the modified Young’s equation. Their calculated contact angle remained constant at different line tension strengths. Recently, Vedantam and Panchagnula25 proposed a 2D phase-field model to simulate the hysteresis wetting of a sessile drop on a substrate. Their model could simulate the contact-angle hysteresis on a chemically heterogeneous surface, and the results were consistent with the experimental observation.26 Unfortunately, their model was not able to simulate the case having negative line tension, and the model was only 2D. A 3D
τ0
⎛ W2 ⎞ W2 ∂ϕ = ⎜ 0 ∇·(ψ ∇ϕ) + fϕ ⎟ − u(t )gϕ + 0 |∇ψ ||∇ϕ|cos θc ∂t ψ ⎝ ψ ⎠ (4)
where u(t) is given by29
u(t ) =
B
∫V fϕ dV̅ ∫V gϕ dV̅
(5) DOI: 10.1021/acs.langmuir.5b02175 Langmuir XXXX, XXX, XXX−XXX
Article
Downloaded by RUTGERS UNIV on August 29, 2015 | http://pubs.acs.org Publication Date (Web): August 21, 2015 | doi: 10.1021/acs.langmuir.5b02175
Langmuir
Figure 2. Field variables in the computational domain and sample meshes: (a) t* = 0; (b) t* = 200.
where gϕ = (1 − ϕ2) and V is the volume of the overall computational domain. To consider the line tension effect, we could replace cos(c) by cos() according to the modified Young’s equation from eq 2. Then, the model with line tension could be obtained τ0
where the 1/rb could be replaced by (1 −ϕ2)/rb for better numerical stability during computation; (1 − ϕ2) has no effect on the equilibrium shape. In eq 6, the scale relation between the interfacial energies γlv and the interface thickness W0 have been derived by Yu et al.23 by using the sharp interface limit;30 however, for the phase-field variable considered here, which ranges from −1 to 1,31 the relation of the liquid−vapor surface tension with the interface thickness W0 could be obtained as γlv = 2√2/3 HW0, where H is a constant chosen for the free-energy function having a unit of energy per volume. Most people chose 1 for convenience;
⎛ W2 ⎞ W2 ∂ϕ = ⎜ 0 ∇·(ψ ∇ϕ) + fϕ ⎟ − u(t )gϕ + 0 |∇ψ ||∇ϕ| ∂t ψ ⎝ ψ ⎠ ⎛ γ 1⎞ ⎜⎜cos θc − slv ⎟⎟ γlv rb ⎠ ⎝
(6) C
DOI: 10.1021/acs.langmuir.5b02175 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir the state with the minimum energy is not affected by H. Furthermore, the computation of the base radius of the curvature is straightforward by taking the surface divergence of the normal vector of the liquid−vapor phase nlv projected on the substrate as 1 = ∇s ·nlv,s rb
numerical mesh is also shown in Figure 2a. As shown, the mesh is refined near the interfaces for −0.99 < ϕ < 0.99 and 0.05 < ψ < 0.95. In addition, in Figure 2, the initial shape of the liquid phase is a sphere (a radius of 20) in contact with the solid phase having a contact angle of 150°. The minimum dimensionless mesh size Δx/W0 = 0.4, and the domain lengths in the x, y, and z directions are 32, 32, and 48, respectively. The dimensionless time step Δt/τ0 = 0.1, and it is chosen simply for the numerical stability. A sample steady-state mesh and the field variables are illustrated in Figure 2b. Figure 3 illustrates a sample simulation of an anisotropic wetting; however, such an evolution is only a numerical
(7)
where the surface gradient ∇s is defined by ∇s = (I − n sv ⊗ n sv) ·∇
(8)
And nlv,s is given by nlv,s = (I − n sv ⊗ n sv) ·nlv
(9)
Downloaded by RUTGERS UNIV on August 29, 2015 | http://pubs.acs.org Publication Date (Web): August 21, 2015 | doi: 10.1021/acs.langmuir.5b02175
where I is the unit tensor and the notation ⊗ is dyadic product. Again, with the stationary field variable ψ, nsv can be calculated easily by n sv =
∇ψ |∇ψ |
(10)
Similarly, the surface normal of the liquid−vapor interface could be obtained by using the phase-field variable as nlv =
∇ϕ |∇ϕ|
Figure 3. Time evolution of an anisotropic wetting on multiple strips from t* = 0 to 200.
(11)
Then, the last term in the modified Young’s equation is rewritten by γslv 1 3γslv = ∇*s ·nlv,s = εij∇*s ·nlv,s γlv rb 2 2 W02
convergence process. Nevertheless, as the fluid flow is considered,20 the dynamic wetting could be simulated. The equilibrium shape is obtained when a steady state is reached. Moreover, the volume of the liquid phase Vl is calculated by
(12)
where ∇*s is the dimensionless surface gradient. The parameter εtj is related to the line tension strength, that is, l = γslv/γlv = εtjW0. In other words, the interfacial energies could be obtained by given εtj and W0. If we rescale the length by W0 and time by τ0, the present model in dimensionless form (with asterisk superscript) becomes
Vl =
⎛ εij ⎞ 1 ⎟ |∇*ψ ||∇*ϕ|⎜cos θc − ψ rb* ⎠ ⎝
⎞
(14)
The volume difference in the whole computation is
