Contact line velocity

Applying a dynamic contact angle poses the challenge of an accurate approximation of the contact line velocity. In the VOF approach the interface is smeared out and no explicit representation of the interface or the contact line exists, which could be taken to calculated the contact line velocity knowing its position over time. Three different approaches to cal-culate the contact line velocity have been implemented.

Simple approach

The simplest approach to estimate the contact line velocity is to take the material point velocityupin the vicinity of the contact line and define the component parallel to the wall as the contact line velocityu_cl (see Figure 3.12a).

u_cl,S= n_Σ−(n_w·n_Σ)n_w

|n_Σ−(n_w·n_Σ)n_w|·u_p. (3.46)

This simple contact line velocity model using the first cell center next to the wall as a reference point is widely used in literature for the simulation of wetting phenomena with interface capturing methods and is also implemented in Open-FOAM.

Roisman et al. (2008) and Linder et al. (2015) argue that this simple model sets the propagation velocity of the contact line equal to the velocity of the material point, which is not physical. This assumption only holds for contact angles of aboutθ≈90^◦. Both propose a new model which we refer to asLinderRoismanmodel.

Projection onto the interface by Roisman et al. and Linder et al.

Based on geometrical considerations shown in Figure 3.12b and assuming a nearly planar interface near the wall, they consider the contact angle rateθ˙and the distancerfrom the pointPat the interface to the contact line and write for the interface velocityu_Σ

u_Σ−rθ˙≈u_clsinθ. (3.47)

If the valuerθ˙is much smaller thanu_cl, the equation can be discretized as

u_cl,LR= u_P·n_Σ

p1−(n_w·n_Σ)². (3.48)

The authors emphasize that for this approximation no creeping flow condition is required since it is only based on geo-metrical considerations. A more in-depth derivation can be found in Linder et al. (2015). The discretized value of the contact line velocity (3.48) has been corrected in this study with a necessary square, which is missing in the derivation and implementation by Linder which was also adapted by Lippert (2016).

Via contact line position

Šikalo et al. (2005) also used a VOF framework to simulate drop impact and they calculated the contact line velocity in an axisymmetric case by differentiation of the drop spreading radius. This application-oriented approach is similar to the more generic approach by Batzdorf (2015) who uses the volume fraction iso-surface reconstruction to estimate the position of the contact line in each time-step. In this approach the contact line velocity can be calculated knowing the position of the current and the last time step (see Figure 3.13a). This study adapts this approach following the steps outlined by Batzdorf.

The iso-surface reconstruction at the wall is a challenge since the volume fraction is a cell-centered field and not properly defined at the wall boundary. To reconstruct the contact line position, the volume fraction at the wall needs a correction.

The gradient of the volume fraction at the wall boundary expressed by the contact angle is used to extrapolate a necessary value ofα^corr_w . This is possible since the gradient of any scalar field in a specific point is orthogonal to its iso-surface,

0.7 0.3

0.5 0.2

(a)

n_w θ

θ

∇αb

∇α^corr_b

c

f

δ

(b)

Figure 3.13:a) The contact line position is reconstructed as an iso-line and by comparing its position in two time steps the contact line velocity is calculated. b) The gradient information of the cell center and the contact angle are used to extrapolate and correct the volume fraction at the boundary face in order to determine the position of the contact line.

which means, the volume fraction extrapolated onto the wall by the contact angle will lead to iso-lines parallel to the contact line. Batzdorf stresses that this is not equivalent to setting the slope of the iso-surface inside the wall adjacent cell.

The extrapolation via the contact angle requires a correction of the volume fraction gradient∇α^corr_w as shown in Figure 3.13b that complies with two conditions. First, the gradient and the wall normal encloses the contact angle

cosθ= ∇α^corr_w

|∇α^corr_w |·n_w, (3.49)

and secondly, the wall tangential component should remain untouched, in particular, be the same as given by the cell-centered volume fraction gradient∇αc assuming a homogeneous Neumann boundary condition. This condition holds fulfilling

∇α^corr_w ×n_w=∇αc×n_w. (3.50)

These equations combine to

∇α^corr_w ·n_w= |∇αc×n_w|cosθ

sinθ , (3.51)

and provide the ingredient to extrapolate the cell-centered volume fractionαc onto the wall faceαf,w

α_f,w=αc+ ∇α^corr_w ·n_w

δf→c. (3.52)

The face-centered volume fractionαf,wis not bounded by zero or one, which is no issue since the flux through the wall is zero and the tangential component to the wall is preserved. To maintain a bounded volume fraction fieldαlimiters are often used on wall boundary conditions. Consequently, the use of limiters should be avoided in combination with this method.

To obtain a smeared out contact line velocity, the distance to the reconstructed contact line is propagated to neighbor cells analogously to the algorithm explained in Section 3.4.2. The current contact line position vectorsxⁿ_cl and the position one from the last time stepxⁿ⁻¹_cl define the contact line velocity

u_cl,pos=− 1

∆t xⁿ_cl−xⁿ⁻¹_cl

· n_Σ−(nw·n_Σ)nw

|nI−(nw·n_Σ)nw|. (3.53)

The contact line position is extrapolated using the current contact angle which itself depends on the contact line velocity and therefore also on the current contact line position itself. This circular dependency is solved by a fix-point iteration in which the value of the contact angle is under-relaxed to promote convergence and the iteration process is stopped if the contact angle change is less then0.01^◦per iteration. This approach is the most accurate implemented here, however, using this method the iteration for almost static contact lines does not converge always due to the contact angle hysteresis and contact line pinning implementation in combination with parasitic currents as will be shown in Section 3.5.5. This model is later referred to ascLVelocityViaCLPosition.