II. Numerical Methods 77
7. The geometrical Volume-of-Fluid method 85
7.5. State of the art: Wetting with geometrical VOF methods
Chapter 7. The geometrical Volume-of-Fluid method
cell, which is almost always the case when the slip length is not resolved, the advection method will not be able to recover the correct contact line motion. Moreover, the viscous dissipation in the contact line region is usually underestimated by the numerical method if the flow in the contact line region is not well-resolved. Both effects lead typically to an over-prediction of the contact line velocity (see Chapter 11 for an example).
In fact, the need to resolve the slip length on the scale of nanometers, which is many orders of magnitude smaller than the typical scale of the flow (say mm or cm), poses one of the major challenges in the numerical simulation of dynamic wetting. One way to solve this problem is by massive local mesh refinement (see, e.g., [SS12a]).
However, doing so is computationally expensive and requires a quite sophisticated mesh that is tailor-made for the specific flow configuration at hand. Another approach, which has been introduced for a geometrical VOF method by Afkhami et al. [AZB09], is to make use of the numerical slip on a relatively coarse mesh and to apply a dependent model for the dynamic contact angle based on the analysis of Cox [Cox86] to reduce the mesh-dependence of the overall solution. Satisfactory mesh-independent results have been obtained with that method for a spreading drop and the withdrawing plate problem (see [AZB09]). We introduce a novel adaptation of the discrete Navier slip condition (7.37) called “staggered slip” in Section 11.4. The staggered slip condition shows improved convergence for the capillary rise problem for an under-resolved slip length and might open another route for a “subgrid modeling” of dynamic wetting in the future.
The contact angle boundary condition in geometrical VOF methods: There are several approaches present in the literature to incorporate the contact angle boundary condition in geometrical VOF methods. According to [SDS14], there is a common feature shared by most of the methods:
“Recent developments in VOF methods have allowed the representation of MCLs in droplet impact (...) or spreading (...) and bubble generation in microfluidic T-junctions (...). Despite differences in the contact-line models used, the implementations commonly impose a gradient of the volume fraction at the contact line in terms of a prescribed contact angle.” [SDS14, p.105]
The current implementation in FS3D (see [FB15a, Lip16]) employs the method introduced by Afkhami and Buss-mann [AB08, AB09]. The method prescribes the interface orientation according to the contact angle boundary condition and modifies the surface tension force at the boundary (see Section 7.5.2 for details). It was originally implemented in an early version of the geometrical VOF solverGerrisdeveloped by Popinet [Pop03] and it is still in use in recent works to study the transition to a Landau-Levich-Derjaguin film in forced dewetting [ABG+18]
and the dynamic wetting failure in curtain coating [FZP20]. Interestingly, recent studies compare VOF methods directly with molecular dynamics and/or phase field methods for nano-scale drops [SHGL19, LJF+20].
In fact, the interface orientation is commonly determined directly by the contact angle boundary condition (with a few exceptions, including, e.g., [ ˇSWR+05]). However, since the advection of the volume fraction field is mostly discretized explicitly in time, the contact angle boundary condition is usually not satisfied after an interface advec-tion step. Instead, the contact angle is enforced by an explicit adjustment of the interface orientaadvec-tion at compu-tational cells located at the boundary (see, e.g., [RRL01, AB08, AB09]). Enforcing the contact angle like this is, however, not consistent with the kinematics of moving contact lines as discussed in Chapter 4. We observed that it can also be a source for numerical instabilities in the vicinity of the contact line (see [RBB12] for a discussion of instabilities close to the contact line). Specialized methods to reconstruct the interface close to the boundary are developed in the present work (see Section 8.2). These methods increase the accuracy of the VOF method at the contact line and are consistent with the fundamental kinematics described in Chapter 4.
7.5.1. No-slip and Navier slip boundary conditions
To allow for a motion of the contact line, one can make use of the numerical slip inherent to the method or prescribe the Navier slip condition. Technically, the velocity boundary condition is indirectly enforced within the Finite Volume method using the concept of “ghost cells”.
(a) In the classical approach, the no-slip boundary condition is enforced by extrapolating the tangential velocity field into the layer of ghost cells adjacent to the physical boundary, such that the velocity interpolates to zero
94
7.5. State of the art: Wetting with geometrical VOF methods
v L= 0
vghost
∂Ω
(a) No-slip (“Numerical slip”).
v L >0
vghost
∂Ω
(b) Standard Navier slip.
Figure 7.3.: Ghost-cell based numerical realization of the no-slip and Navier slip boundary conditions in FS3D.
exactly at the boundary (see FIG. 7.3(a)), i.e.
vghost=−v. (7.36)
This velocity gradient is subsequently counteracted by the discrete realization of the viscous forces applied next to the boundary. This approach enforces the no-slip condition only in the limit of the mesh size going to zero. Therefore, the face-centered velocity at the boundary cell layer may be non-zero leading to a motion of the contact line. However, the extent of numerical slip decreases with increasing mesh resolution, leading to a significant mesh dependence of the solution. This numerical effect has been first described in the context of VOF methods by Renardy et al [RRL01].
(b) The Navier slip boundary condition is enforced following the same approach by setting a different velocity in the ghost-cell layer, such that the tangential velocity interpolates to zero at a fixed distanceL>0 from the physical boundary (see FIG. 7.3(b)), i.e.
vghost=v2L−∆x
2L+∆x. (7.37)
It has been demonstrated in the literature that this may lead to mesh-convergent results if the resolution of the computational mesh is well below the slip lengthL(see, e.g., [SS12a, GSA+20a] and Chapter 11). However, physically reasonable slip lengths are on the scale of nanometers [NEB+05]. This length scale is not accessible with the present numerical method without massive computational costs. Moreover, note that for L>0 the magnitude of the “counter velocity”vghostin the ghost cell is always less than or equal to|vghost|in the case of no-slip. This follows from the inequality
2L−∆x 2L+∆x ≤1, which is valid for positiveL.
7.5.2. The height function method for applying contact angles
The height function method described in Section 7.4 was adapted by Afkhami and Bussmann [AB08, AB09] to enforce the contact angle boundary condition in Volume-of-Fluid simulations.
Basic methodology: The purpose of the method is to “impose” the contact angle boundary condition required in the continuum mechanical model (7.25), i.e.
hnΣ,n∂Ωi=−cosθ at Γ(t), (7.38) whereθ satisfies a relation of typeθ = f(VΓ), to a VOF simulation. The boundary condition (7.38) has to two implications for the algorithm, namely “it defines the orientation~nof the VOF reconstruction in the contact line cell, and it influences the calculation ofFst by affecting the calculated curvature.” [AB08, p.459]
Chapter 7. The geometrical Volume-of-Fluid method
fσ θ
H1
H0
H−1
Figure 7.4.: Linear extrapolation of the height function as a numerical driving force for wetting.
(i) The implication of the first statement is that the interface orientation is prescribed at the solid boundary according to the contact angle boundary condition instead of being a result of the interface reconstruction method. As described above, this prescribed orientation might be inconsistent with the interface orientation which results from the interface kinematics.
(ii) The second ingredient in the numerical method to apply the contact angle adapts the numerical approximation of the mean curvature in a contact line cell. A sketch of the method is given in Fig. 7.4. If an interface cell is located at the domain boundary, one cannot apply the central finite differences formulas (7.30) to compute the curvature since, by definition, the height information is missing outside the computational domain (see discrete height H−1 in Fig. 7.4). The suggestion of the method [AB08] is to reconstruct a virtual height in the boundary cell layer using the contact angle boundary condition and to compute the curvature with the formulas (7.30) using the virtual height information. More precisely, the height function is extrapolated linearly into the ghost-cell layer with a slope that is determined by the contact angle boundary condition. For the configuration depicted in Fig. 7.4, the formula to computeH−1is
H−1=H0+∆xcotθ.
The generalization of the method to three spatial dimensions is described in [AB09]. Note that both the two-and three dimensional variant of the method have been implemented in the FS3D solver by Lippert3. It is clear that the linear extrapolation according to the contact angle leads to a kink in the virtual interface if the current contact angle does not coincide with the desired contact angle according to (7.38). This leads to a singularity in the curvature at the contact line in the limit as∆x→0, which might lead to very large numerical values of ˆκeven on coarse meshes. When applied to the numerical approximation of the surface tension force, the singular curvature leads to large velocities at the contact line that drive the interface towards the desired contact angle; see Fig. 7.4 for an example. Hence, the height function method [AB08, AB09] introduces a numerical driving force to indirectly “apply” the contact angle boundary condition. In this sense, the method is similar to the Generalized Navier Boundary condition (GNBC) described in Section 3.4. However, contrary to the GNBC, the driving force is introduced directly in the numerical method with no proper counterpart on the PDE level.
Failure of the method in the equilibrium state: A fundamental drawback of the interface extrapolation method is the missing convergence in stationary state. In fact, thelinearextrapolation of the interface beyond the solid boundary creates a height function representation of classC1only (unless the curvature at the contact line van-ishes). In the case of a two-dimensional spherical droplet, the curvature of the extrapolated height function ˜h admits a jump from−1/Rto zero at the contact line. As a result, the central finite differences scheme doesnot converge to the correct curvature at the contact line, even if the height function is known exactly. In fact, a much higher regularity of the height function (h∈C5) is needed to prove convergence of the method; see Lemma 9.4 in Chapter 9 for details. Example 7.2 illustrates the inconsistency of the method.
Example 7.2(Height Function Extrapolation forθ =π/2). For simplicity, we consider the case of a spherical droplet in 2D with contact angleθ =π/2; see Fig.7.5 for the setup. A height function representation of the interface above the solid boundaryx=0 is given by
{(x,h(x)): 0≤x≤R}
3More details on the implementation in FS3D and some basic validation of the contact line dynamics in two and three dimensions can be found in [Lip16].
96
7.5. State of the art: Wetting with geometrical VOF methods
withh(x) =Rp
1−(x/R)2. The linear extension of the height function reads h(x) =˜ R
(q 1− Rx
2
, 0≤x≤R
1 x≤0
. (7.39)
Clearly, the extended height function is of classC1((−∞,R]). For a finite grid size∆>0, the discrete heightsH±1
H1
H0
H−1 κ= 0
κ=−R1
Figure 7.5.: Linear extrapolation forθ=π/2.
andH0are defined by integration
H0(∆) =1
∆ Z ∆
0
h(x)dx=R 2
s
1− ∆
R 2
+R
∆arcsin ∆
R
,
H1(∆) =1
∆ Z 2∆
∆
h(x)dx=R
s
1− 2∆
R 2
+ R 2∆arcsin
2∆
R
−H0. H−1(∆) =1
∆ Z 0
−∆
h(x)dx=R.
Note that limx→0+1
xarcsin(x) =1 implies lim∆→0+H0=lim∆→0+H1=R. Using l’Hospital’s rule one easily shows that
H1(∆)−H−1(∆)
2∆ →0 as ∆→0+.
Consequently, the numerical approximation of the curvature at the contact line satisfies lim
∆→0+
κˆ= lim
∆→0+
H1(∆)−2H0(∆) +H−1(∆)
∆2
. A short calculation shows that
lim
∆→0+(H1(∆)−2H0(∆) +H−1(∆))0=0 and, by l’Hospital’s rule, it follows that
lim
∆→0+
κˆ = lim
∆→0+
(H1(∆)−2H0(∆) +H−1(∆))00
2 =− 5
6R. (7.40)
So, the curvature at the contact line doesnotconverge to the correct value of−1/R.