Contact line advection using the geometrical Volume-of-Fluid method

TL;DR

This paper develops a geometrical Volume-of-Fluid method for accurately simulating contact line advection in dynamic wetting, ensuring correct contact angle transport and analyzing convergence properties of different interface reconstruction techniques.

Contribution

It introduces a boundary-adapted VOF method for contact line advection that achieves linear convergence of the contact angle, improving upon existing methods.

Findings

Boundary ELVIRA method achieves linear convergence.

Boundary Youngs method exhibits oscillations and does not converge.

The proposed method correctly transports the contact angle with prescribed velocity.

Abstract

We consider the interface advection problem by a prescribed velocity field in the special case when the interface intersects the domain boundary, i.e. in the presence of a contact line. This problem emerges from the discretization of continuum models for dynamic wetting. The kinematic evolution equation for the dynamic contact angle (Fricke et al., 2019) expresses the fundamental relationship between the rate of change of the contact angle and the structure of the transporting velocity field. The goal of the present work is to develop an interface advection method that is consistent with the fundamental kinematics and transports the contact angle correctly with respect to a prescribed velocity field. In order to verify the advection method, the kinematic evolution equation is solved numerically and analytically (for special cases). We employ the geometrical Volume-of-Fluid (VOF) method on a structured Cartesian grid to solve the hyperbolic transport equation for the interface in two spatial dimensions. We introduce generalizations of the Youngs and ELVIRA methods to reconstruct the interface close to the domain boundary. Both methods deliver first-order convergent results for the motion of the contact line. However, the Boundary Youngs method shows strong oscillations in the numerical contact angle that do not converge with mesh refinement. In contrast to that, the Boundary ELVIRA method provides linear convergence of the numerical contact angle transport.

Equations

$$\dot{\theta }=\frac{\partial v}{\partial \tau }\cdot n_{\Sigma }$$

$$\Gamma (t)=\partial \Sigma (t)\subset \partial \Omega .$$

$$\cos \theta (t,x)=-⟨n_{\Sigma }(t,x),n_{\partial \Omega }(x)⟩,$$

$$\gamma (t)\in \overline{\Sigma }(t)\text{and}\gamma (t_0)=x_0\in \overline{\Sigma }(t_0),$$

$$V_{\Sigma }(t_0,x_0)=⟨{\gamma }^{\prime }(t_0),n_{\Sigma }(t_0,x_0)⟩.$$

$$n_{\Gamma }=\frac{{\mathcal{P}}_{\partial \Omega }{n}_{\Sigma }}{\Vert {\mathcal{P}}_{\partial \Omega }{n}_{\Sigma }\Vert },$$

$$t_{\Gamma }=n_{\Gamma }\times n_{\partial \Omega }.$$

$$V_{\Sigma }=⟨v,n_{\Sigma }⟩\text{on}\mathrm{gr}\overline{\Sigma },$$

$$\nabla \cdot v=0$$

$$v\cdot n_{\partial \Omega }=0$$

$$\left\llbracket v\right\rrbracket=0\ \ \text{on}\ \operatorname{gr}\overline{\Sigma},\quad\nabla\cdot v=0\ \ \text{in}\ \ \Omega\setminus\Sigma(t),\quad v\cdot n_{\partial\Omega}=0\ \ \text{on}\ \ \partial\Omega$$

$$\dot{x}(t)=v(t,x(t)),x(t_0)=x_0\in \overline{\Sigma }(t_0).$$

$$\dot{\psi }(t_0,x_0)=\frac{d}{dt}\psi (t,x(t){)}_{|t=t_0}.$$

$${\dot{n}}_{\Sigma }(t,x)$$

$$\dot{\theta }(t,x)=⟨{\partial }_{\tau }v(t,x),n_{\Sigma }(t,x)⟩,$$

$$\tau (t,x)=-\cos [\theta (t,x)]n_{\Gamma }(t,x)-\sin [\theta (t,x)]n_{\partial \Omega }(t,x).$$

$$v(x_1,x_2)=(v_0+c_1{x}_1+c_2{x}_2,-c_1{x}_2).$$

$$\frac{v(x_1,x_2)}{L/T}=({\hat{v}}_0+{\hat{c}}_1{\hat{x}}_1+{\hat{c}}_2{\hat{x}}_2,-{\hat{c}}_1{\hat{x}}_2)=:\hat{v}({\hat{x}}_1,{\hat{x}}_2)$$

$$\nabla v=\left(\begin{array}{cc}{c}_1& c_2\\ 0& -c_1\end{array}\right).$$

$${\dot{x}}_1(t)=v_1(x_1(t),0)=v_0+c_1{x}_1(t),x_1(0)=x_1^0.$$

$$x_1(t)=x_1^0{e}^{c_{1}t}+\frac{v_0}{c_1}\left(e^{c_{1}t}-1\right)\text{for}{c}_1\ne 0$$

$$n_{\Sigma }^l=\left(\begin{array}{c}-\sin \theta \\ \cos \theta \end{array}\right)\text{and}{n}_{\Sigma }^r=\left(\begin{array}{c}\sin \theta \\ \cos \theta \end{array}\right).$$

$${\tau }^l=\left(\begin{array}{c}\cos \theta \\ \sin \theta \end{array}\right)\text{and}{\tau }^r=\left(\begin{array}{c}-\cos \theta \\ \sin \theta \end{array}\right).$$

$$\dot{\theta }(t)=\pm c_2{\sin }^2\theta -2c_1\sin \theta \cos \theta ,\theta (0)={\theta }_0$$

$$\theta (t)=\frac{\pi }{2}+\arctan \left(-\cot {\theta }_0{e}^{2c_{1}t}\pm c_2\frac{e^{2c_{1}t}-1}{2c_1}\right).$$

$$\underset{c_1\to 0}{\lim }\left(-\cot {\theta }_0{e}^{2c_{1}t}\pm c_2\frac{e^{2c_{1}t}-1}{2c_1}\right)=-\cot {\theta }_0\pm c_{2}t.$$

$$\chi (t,x)=\{\begin{array}{ll}1& \text{if}x\in {\Omega }^{-}(t)\\ 0& \text{if}x\notin {\Omega }^{-}(t)\end{array},$$

$${\partial }_t\chi +v\cdot \nabla \chi =0$$

$${\partial }_t\chi +\nabla \cdot (v\chi )=\chi (\nabla \cdot v).$$

$$\frac{d}{dt}{\alpha }_{ij}(t)=-\frac{1}{|V_{ij}|}{\int }_{\partial V_{ij}}\chi v\cdot ndA+\frac{1}{|V_{ij}|}{\int }_{V_{ij}}\chi (\nabla \cdot v)dV,$$

$${\alpha }_{ij}(t)=\frac{1}{|V_{ij}|}{\int }_{V_{ij}}\chi dV.$$

Full text

We employ the geometrical Volume-of-Fluid (VOF) method on a structured Cartesian grid to solve the hyperbolic transport equation for the interface in two spatial dimensions. We introduce generalizations of the Youngs and ELVIRA methods to reconstruct the interface close to the domain boundary. Both methods deliver first-order convergent results for the motion of the contact line. However, the Boundary Youngs method shows strong oscillations in the numerical contact angle that do not converge with mesh refinement. In contrast to that, the Boundary ELVIRA method provides linear convergence of the numerical contact angle transport.

This preprint was submitted and accepted for publication in the Journal of Computational Physics.

When citing this work, please refer to the journal article: DOI: 10.1016/j.jcp.2019.109221.

Keywords: Volume-of-Fluid, Interface Reconstruction, Dynamic Contact Angle, Kinematics

1 Introduction

The present article deals with the passive transport of an interface by a prescribed velocity field, a well-known problem in the numerical description of multiphase flows called interface advection (see, e.g., [1, 2]). In particular, we consider the special case when the interface intersects the domain boundary and a so-called contact line is formed. This problem arises, for example, as part of the mathematical description of wetting processes, where a liquid displaces another liquid (or a gas) in order to wet the surface of a third solid phase. In the latter case, the velocity field is a solution of the two-phase Navier Stokes equations with appropriate boundary and transmission conditions modeling the physics of wetting.

In the present work, we focus on the kinematics of wetting, hence treating the velocity field as given, while keeping in mind that it is a solution to some continuum mechanical model. The continuum mechanical modeling of wetting is a separate issue that is not addressed here (see, e.g., [3, 4]). The kinematics of wetting is studied analytically in [5, 6] where it is shown that the contact line advection problem is a well-posed initial value problem if the velocity field is sufficiently regular and tangential to the domain boundary. Hence the full dynamics of the interface can be inferred from the (time-dependent) velocity field and the initial interface configuration. Note that this is true even though in the full continuum mechanical model the interface shape is typically strongly coupled to the flow. The reason is that the velocity field contains enough information to reconstruct the evolution of the interface from the initial configuration (see Figure 2 for a sketch of the idea). In particular, the motion of the contact line as well as the evolution of the contact angle, i.e. the angle of intersection between the fluid-fluid interface and the solid boundary (see Figure 2), can be computed from the knowledge of the velocity field and the initial geometry. In fact, the rate of change of the contact angle $\dot{\theta}$ satisfies the kinematic evolution equation

[TABLE]

derived in [5] (see Figure 2). The goal of the present work is to develop an advection scheme for the geometrical VOF method that allows to transport the contact angle consistently with the fundamental kinematic relation (1), which will serve as a reference to validate the method (see Section 2).

On the other hand, mathematical models for dynamic wetting usually prescribe a boundary condition for the contact angle. The contact angle is an important physical parameter characterizing the wettability of the solid surface. Mathematically, it is typically prescribed as a fixed value or as a function of the speed of the contact line (see, e.g., [7]). Note that in view of (1), a boundary condition for the contact angle leads to a consistency condition for the velocity field. In order to avoid the moving contact line paradox [8], tangential slip is usually introduced on the solid boundary, at least close to the contact line.

The contact angle boundary condition also plays an important role in the numerical simulation of moving contact lines (see [9] for a recent review). In the present study, we employ the Volume-of-Fluid (VOF) method (see [10, 1]) to numerically track the location of the interface. The VOF method has been successfully adapted for the simulation of moving contact lines by several authors (see, e.g., [11, 12, 13, 14, 15, 16, 17, 18, 19]). Since the advection of the volume fraction field is mostly discretized explicitly in time, the contact angle boundary condition is usually not satisfied after a transport step. Instead, the contact angle is enforced by an explicit adjustment of the interface orientation at computational cells located at the boundary. Enforcing the contact angle as described above is, however, not consistent with the kinematics of moving contact lines as discussed in [5],[6]. We observed that it can also be a source for numerical instabilities in the vicinity of the contact line (see [17] for a discussion of instabilities close to the contact line).

Notation:

We consider a moving $\mathcal{C}^{1,2}$-hypersurface $\{\overline{\Sigma}(t)\}_{t\in I}$ with boundary (as defined in the appendix) embedded in two or three-dimensional Euclidean space. Each instantaneous interface $\overline{\Sigma}(t)$ is assumed to be an orientable $\mathcal{C}^{2}$-hypersurface with unit normal vector field $n_{\Sigma}(t,\cdot)$ and boundary $\partial\Sigma(t)$ that is contained in the planar domain boundary denoted by $\partial\Omega$. The contact line at time $t$ is denoted by

[TABLE]

Given a point $x\in\Gamma(t)$, the contact angle $\theta\in(0,\pi)$ is defined by the relation

[TABLE]

where $\left\langle\cdot,\cdot\right\rangle$ is the Euclidean scalar product and $n_{\partial\Omega}$ is the outer unit normal field to $\partial\Omega$. In order to describe the motion of the interface, the notion of the speed of normal displacement is required. If $\gamma$ is a continuously differentiable curve such that

[TABLE]

then the speed of normal displacement at the point $x_{0}$ at time $t_{0}$ is given as111It can be shown that this definition is independent of the choice of the curve $\gamma$, see e.g. [20].

[TABLE]

It is convenient to introduce a contact line normal vector (see Figure 2) via projection as

[TABLE]

where $\mathcal{P}_{\partial\Omega}=\mathds{1}-n_{\partial\Omega}\otimes n_{\partial\Omega}$ is the local orthogonal projection onto $\partial\Omega$. The above expression is well-defined since we only consider the partial wetting case characterized by $0<\theta<\pi$. To construct a local orthonormal basis in three dimensions, one further defines

[TABLE]

Problem formulation:

The motion of a material interface is governed by the kinematic condition

[TABLE]

where $v$ is the transporting velocity field. Here it is assumed that $v\in\mathcal{C}^{1}(\overline{\Omega})$ is divergence free and tangential to the domain boundary, i.e.

[TABLE]

For simplicity, we further assume that the solid boundary is planar. This assumption is not essential and may be dropped to study the contact angle evolution over curved surfaces. Note also that incompressibility of the flow is not necessary to study the contact angle evolution. However, it is assumed here since the VOF method is most commonly applied to incompressible flows.

We emphasize that the assumption of a globally continuously differentiable velocity field is not met in a typical multiphase flow, where the velocity gradient admits a jump which is controlled by the interfacial transmission condition for the stress. However, it has been shown in [5] Lemma 8, that in the case of two spatial dimensions222A similar statement holds in three spatial dimensions. For that case, one can show the continuity property \left\langle\left\llbracket\nabla v\right\rrbracket\alpha,\beta\right\rangle=0 at $\Gamma$, where \left\llbracket\cdot\right\rrbracket denotes the jump over the interface and $\alpha$, $\beta$ are arbitrary vectors in the plane spanned by $n_{\partial\Omega}$ and $n_{\Gamma}$, see [5]. the conditions

[TABLE]

imply continuity of $\nabla v$ at the contact line (under certain regularity assumptions). For simplicity, we therefore assume a globally continuously differentiable field.

In this work, we focus on the Volume-of-Fluid method [10, 1] for advecting the interface, whose reconstruction algorithm must be adapted in order to achieve second-order convergence in the near-wall region. However, the contact line advection problem can be used to verify near-wall advection of any other interface advection method. For example, a brief discussion of the contact line advection problem with the Level Set Method [21, 22] in two dimensions can be found in [23]. An open research data record containing the full C++-implementation of the Level Set Method used in [23] together with a number of computational examples is available online, see [24].

2 Kinematic evolution equation

The kinematic evolution equation (see [5],[6]) allows to compute the rate of change of the contact angle in terms of the transporting velocity field. The evolution of the contact angle is considered along the trajectories (or “characteristics”) of the flow, i.e. along solutions of the ordinary differential equation

[TABLE]

The conditions (4), (5) and (6) imply that both $\operatorname{gr}\Gamma$ and $\operatorname{gr}\Sigma=\operatorname{gr}\overline{\Sigma}\setminus\operatorname{gr}\Gamma$ are invariant for the flow generated by (7), i.e. a solution of (7) starting on the contact line (the fluid-fluid interface) will stay on the contact line (the fluid-fluid interface); see [5] for a proof. Therefore, it is possible to study the time evolution of the contact angle along a solution of (7) starting at the contact line. Along such a curve, one can define the Lagrangian time-derivative of a quantity $\psi\in\mathcal{C}^{1}(\operatorname{gr}\overline{\Sigma})$ according to

[TABLE]

Given a $\mathcal{C}^{1}$-velocity field, the Lagrangian time-derivative of the interface normal vector is given by (see [5], [6])

[TABLE]

Clearly, the evolution of the contact angle can be inferred from the normal vector evolution via the relation (2). Moreover, it has been shown [5], [6] that, for a planar boundary333Note that (8) also holds if $\partial\Omega$ is not planar., the contact angle follows the evolution equation

[TABLE]

where the tangential direction $\tau$ is defined as (see Figure 2)

[TABLE]

The goal of the present work is to verify the numerical solution delivered by the VOF method against an analytical (if available) or numerical solution of (9).

2.1 An analytical solution for linear velocity fields in 2D

We first consider the case of general linear divergence free velocity fields in 2D. In this case, the velocity gradient $\nabla v$ is constant in space and time and the ODE system (7) and (8) is explicitly solvable. Note that this also provides a local approximation for general differentiable velocity fields.

We choose a Cartesian coordinate system $(x_{1},x_{2})$ such that the solid wall is represented by $x_{2}=0$. We consider a velocity field of the form

[TABLE]

The coefficients $c_{1}$ and $c_{2}$ in this formulation have the dimension of $s^{-1}$. Therefore, it is more convenient to choose a length scale $L$ and a time scale $T$ and write

[TABLE]

with the non-dimensional quantities $\hat{x}_{i}=x/L$, $\hat{c}_{i}=c_{i}T$ and $\hat{v}_{0}=(Tv_{0})/L$. In the following, we will use the formulation (11) while dropping the hats. For a field of this form, the (constant) gradient is given by

[TABLE]

Motion of the contact line:

The motion of the contact line is determined by the ordinary differential equation

[TABLE]

The unique solution of the above initial value problem is

[TABLE]

and $x_{1}(t)=x_{1}^{0}+v_{0}t$ for $c_{1}=0$.

Contact angle evolution:

Note that the constancy of $\nabla v$ decouples the system (7) and (8). Hence, the evolution of the normal vector can be solved independently of the evolution of the contact point. To find the solution, we make use of the fact that in two dimensions the normal vector $n_{\Sigma}$ is, up to a reflection, uniquely determined by the contact angle $\theta$. Given a contact angle $\theta$, the two possibilities are

[TABLE]

In the case of a droplet (and for $\theta<\pi$), this corresponds to the two distinct contact points (left and right). The corresponding expressions for $\tau$ are

[TABLE]

This allows to infer the evolution of $\theta$ for the left and the right contact point directly from (9) without the need to solve the system (8).

Inserting the expressions for $\nabla v$, $n_{\Sigma}^{l,r}$ and $\tau^{l,r}$ to (9) yields the nonlinear ordinary differential equation444Note that we now use the Lagrangian formulation and write $\theta(t)$ for the contact angle at $x(t)\in\Gamma(t)$ where $x(t)$ is a trajectory of the flow, i.e. a solution of (7).

[TABLE]

with the “$+$” for the evolution of the right contact point and the “$-$” for the evolution of the left contact point. The unique solution of the above initial value problem is given by the formula (see Appendix for a derivation)

[TABLE]

Remark:

Obviously, the solution is independent of the parameter $v_{0}$. This is to be expected since the two differential equations (7) and (8) decouple and the parameter $v_{0}$ can be eliminated by a change of the frame of reference. Moreover, the evolution of the left and the right contact point is identical if $\theta_{0}^{r}=\theta_{0}^{l}$ and $c_{2}=0$. Finally, we note that (14) has a well-defined limit for $c_{1}\rightarrow 0$ since

[TABLE]

3 Numerical Method

The numerical studies in the present work are carried out with the Volume of Fluid in-house code Free Surface 3D (FS3D) originally developed by Martin Rieber [25]. Since then, FS3D has been further developed at the University of Stuttgart (see, e.g., [26] and more references given there) and at the Technical University of Darmstadt (see, e.g., [15] and references therein).

3.1 Geometrical Volume-of-Fluid Method

The Volume of Fluid method makes use of the phase indicator function for one of the phases $\Omega^{\pm}(t)$, which are separated by the interface $\Sigma(t)$, i.e.

[TABLE]

to describe the interface. Formally, the kinematic condition (4) translates to the transport equation

[TABLE]

for the indicator function, which is its conservative forms reads as

[TABLE]

Note that the right-hand side vanishes for an incompressible flow. However, it is kept in the discretization to enhance volume conservation properties [1]. The latter approach is called “divergence correction”. Integrating (16) over a control volume in two spatial dimensions leads to

[TABLE]

where the volume fraction $\alpha_{ij}$ associated with the control volume $V_{ij}$ at time $t$ is introduced as

[TABLE]

A temporal integration of (17) leads to the exact transport equation of the volume fraction field given as

[TABLE]

The conservative form of (19) allows for the exact conservation of the phase volume also in the discrete case. Algebraic Volume-of-Fluid methods approximate the solution of (19) using a Finite Volume discretization that relies on an interpolation of $\chi$ at $\partial V_{ij}$. The sharp jump in $\chi$ between phases ”$+$” and ”$-$” then causes oscillations in the numerical solution, that are counteracted by adding artificial diffusion which, in turn, leads to artificial smearing of the interface [27]. Geometric Volume-of-Fluid methods are more accurate than their algebraic counterparts, as they approximate the numerical flux, which can be decomposed in a sum over faces $A_{f}$ of the control volume according to

[TABLE]

by reconstructing a sharp geometrical approximation of the indicator function $\chi(t,\cdot)$ and subsequently approximating the integral using geometrical methods.

Within the FS3D code, equation (19) is solved on a structured Eulerian grid using an operator splitting method as described below (see also [1, 25, 26]).

Implementation in FS3D:

FS3D relies on structured Cartesian meshes and therefore decomposes the boundary of each cell into $4$ planar faces (i.e. edges) in two dimensions with normal vectors that are collinear with coordinate axes. The grid points used for the velocity components (face/edge centroids) are shifted with respect to the grid points used for the volume fractions (volume/area centroids) as shown in Figure 3.

The idea of the operator splitting approach [1, 25] is to decompose the full transport problem into a series of one-dimensional transport steps along the coordinate axis. Formally, this is achieved by decomposing the velocity as $v=(v_{1},0)+(0,v_{2})=:\tilde{v}_{1}+\tilde{v}_{2}$. With this notation, equation (19) reads as

[TABLE]

Note that the projected velocity fields $\tilde{v}_{i}$ are no longer divergence free and it is, therefore, important to consider the divergence correction terms in the discretization of (20). The discrete set of equation reads as

[TABLE]

where $\delta V_{i\pm 1/2,j}$ denotes the volume flux over the edge $(i\pm 1/2,j)$ in a $x$-sweep, $\delta V_{i,j\pm 1/2}$ denotes the volume flux over the edge $(i,j\pm 1/2)$ in a $y$-sweep and

[TABLE]

The choices $\beta=0$ and $\beta=1$ correspond to an explicit or implicit discretization of the divergence correction, respectively. For the present study, we choose $\beta=0.5$. The order of the direction of the sweeps is exchanged after each time step to avoid numerical asymmetries [28]. After each directional split transport step, a heuristic volume redistribution algorithm similar to [29] is applied to enforce boundedness of the method, i.e. $0\leq\alpha_{i,j}^{n}\leq 1$.

Approximation of the volume flux:

The volume fraction field $\alpha_{ij}$ is used to reconstruct a planar interface in each cell which is cut by the interface $\Sigma(t)$, in effect approximating $\chi$ with the Piecewise Linear Interface Calculation (PLIC). The reconstructed PLIC indicator $\tilde{\chi}$ is then used to approximate the flux integral

[TABLE]

where $A_{f}$ is a face (edge) of the cell $ij$. The geometrical interpretation of the above integral is a volume of the phase $\Omega^{-}$ (indicated with $\chi$) that passes through the face $A_{f}$. The calculation of this volume is performed by sweeping the face (edge) $A_{f}$ of the cell backwards along the trajectory given by $\tilde{v}_{d}$, and then clipping this swept volume with $\tilde{\chi}$. Note that the operator splitting simplifies the flux calculation because it simplifies the swept volume as a rectangular cuboid, whose intersection with $\tilde{\chi}$ is much simpler to calculate than for the swept volume that is traced along the Lagrangian trajectories given by the full velocity $v$. The simplification comes at the price of performing an additional PLIC interface intersection per splitting step.

Flux at boundary faces:

No special care is necessary to compute the flux over cell faces at an impermeable domain boundary since we assume a vanishing normal velocity leading to zero flux there.

Accurate handling of the contact line motion by the geometrical VOF method, whose verification is outlined in section 2, requires accurate interface reconstruction methods near the boundary of the solution domain. Proposed adaptations of well-known reconstruction algorithms (see Section 3.2) are described in detail in Section 3.3.

3.2 Interface Reconstruction

A large number of methods have been developed for the geometrical reconstruction of the interface from the volume fraction field. An overview of reconstruction algorithms can be found in [30, 31]. In the present work, we consider two reconstruction algorithms, namely the classical method by Youngs [32] and the ELVIRA method due to Pilliod and Puckett [30]. We propose extensions of these two methods to reconstruct an interface close to the boundary in Section 3.3. To keep the formulas simple, assume that the mesh is equidistant in each direction with mesh sizes denoted as $\Delta x_{1}$ and $\Delta x_{2}$.

A natural measure for the interface reconstruction error is the $L^{1}$-error defined as

[TABLE]

where $\hat{\chi}$ is the characteristic function of the reconstructed domain.

Youngs Reconstruction Method:

The idea of the Youngs method is to approximate the interface normal vector by the discrete gradient of the volume fraction field, i.e.

[TABLE]

Then a plane with orientation $n_{\Sigma}^{Y}$ is positioned such that the volume fraction in the local cell is matched (see [33] for details of the positioning algorithm). The gradient in (23) is approximated by weighted central finite differences on a $3\times 3$-block of cells. For an equidistant mesh, the gradient at cell $(i,j)$ is discretized with central finite differences as

[TABLE]

The Youngs Method is known to be one of the fastest methods for interface reconstruction from volume fractions fields. But it is only first-order accurate with respect to the $L^{1}$-norm since it fails to reconstruct all planar interface exactly, see [30]. This may be explained as a consequence of the lack of regularity of the volume fraction field $\alpha(i,j)$, which can be understood as the evaluation at the cell centers of the continuous function

[TABLE]

obtained from averaging the phase indicator function over a control volume. It is well-known [33] that the latter function is only of class $\mathcal{C}^{1}$. One can, therefore, not expect convergence of the finite differences scheme (24). This underlines the need for more advanced methods.

(E)LVIRA Method:

The idea of the Least Squares VOF Interface Reconstruction Algorithm (LVIRA) method proposed by Puckett [34] is to find a planar interface reconstruction that minimizes the quadratic deviations of the volume fractions in a $3\times 3$-block under the constraint that this interface exactly reproduces the volume fraction in the central cell. Hence, one minimizes the functional

[TABLE]

where $\tilde{\alpha}_{i+k,j+l}(n)$ is the volume fraction in cell $(i+k,j+l)$ which is induced by a plane with orientation $n$ satisfying $\tilde{\alpha}_{ij}(n)=\alpha_{ij}$. Due to the nonlinear constraint, the minimization problem cannot be reformulated as a linear system of equations.

Since the minimization of (25) is computationally expensive compared to a direct method like the Youngs reconstruction, Pilliod and Puckett introduced the Efficient Least Squares VOF Interface Reconstruction Algorithm (ELVIRA) [30]. The computational costs are reduced by minimizing (25) only over a finite set of candidate orientations obtained in the following way:

Suppose the interface can be described in the slope-intercept form

[TABLE]

Then the interface normal vector is either

[TABLE]

The slope is approximated by central-, forward- and backward-finite-differences of column sums. The candidates for the slope $m_{x}$ in the cell $(i,j)$ are

[TABLE]

We observe that the slope should be approximated with $+m_{1}^{b}$, $+m_{1}^{c}$ or $+m_{1}^{f}$ if the second component of $n_{\Sigma}$ is positive and vice versa with $-m_{1}^{b}$, $-m_{1}^{c}$ or $-m_{1}^{f}$ if the second component of $n_{\Sigma}$ is negative. This results in $6$ candidates for the normal vector obtained from column sums in the $x_{2}$-direction.

Obviously, it is not always possible to represent the interface as a graph over $x_{1}$. Therefore, one also has to consider the case

[TABLE]

This gives rise to analogous approximations for $m_{2}^{b}$,$m_{2}^{c}$ and $m_{2}^{f}$. This results in $12$ candidates for the interface normal in two dimensions. It can be shown to be sufficient to reconstruct any straight line exactly which makes the method formally second-order accurate with respect to the $L^{1}$-error, see [30].

3.3 Interface Reconstruction close to the Boundary

We propose an adaptation of well-known methods for interface reconstruction in Volume-of-Fluid methods, aiming for accurate reconstruction of the interface close to the domain boundary. The FORTRAN implementations of the developed interface reconstruction methods are available online in an open research data repository, see [35].

For simplicity of notation, we consider in the following the domain boundary at $x_{2}=0$, i.e. we consider an interface cell with index $(i,j=1)$.

3.3.1 Boundary Youngs Method

We consider a $3\times 3$-block of cells and aim at reconstructing the interface in the lower middle cell. We propose to discretize the gradient of the volume fraction field in two space dimensions in the following way (see Figure 4):

1. (a)

Tangential to the domain boundary central finite differences are used. 2. (b)

In normal direction to the domain boundary, weighted forward finite differences are employed.

From Taylor’s formula, one can show that for a $\mathcal{C}^{3}$-function $f$, the first derivative $f^{\prime}(x)$ can be approximated with second-order accuracy according to

[TABLE]

This formula is applied to approximate the derivative of the volume fraction normal to the boundary. Note, however, that the volume fraction $\alpha$ is only $\mathcal{C}^{1}$ (see Section 3.2). One can, therefore, not expect a convergence of the orientation with that method. But we still consider it here since it is a straightforward extension of the widely used Youngs method to the boundary case.

For an equidistant grid in two space dimensions, the Boundary Youngs gradient in a cell with index $(i,1)$ is discretized as

[TABLE]

3.3.2 Boundary ELVIRA Method

In order to allow for mesh convergent results for the contact angle evolution, one needs a reconstruction method which is second-order accurate at the boundary. Therefore, we propose the following adaptation of the ELVIRA method due to Pilliod and Puckett [30]: Minimize the functional

[TABLE]

where $\tilde{\alpha}_{i+k,1+l}(n)$ is the volume fraction in cell $(i+k,1+l)$ which is induced by a plane with orientation $n$ satisfying $\tilde{\alpha}_{i,1}(n)=\alpha_{i,1}$. Here the minimization is performed over a larger stencil of $5\times 3$ cells. This turns out to be necessary to reconstruct every straight line at the boundary exactly. Following the idea of the Efficient Least Squares VOF Interface Reconstruction Algorithm [30], the functional (30) is minimized over a finite set of candidate orientations obtained from finite differences of column sums (see Figure 5, where the column sums in horizontal direction are visualized in red).

The following candidate slopes are computed in the cell $(i,1)$ in normal direction to the boundary

[TABLE]

The following candidate slope is computed from sums tangentially to the boundary

[TABLE]

This yields $16$ candidates for the interface normal. We can demonstrate by numerical experiments that this is sufficient to reconstruct any straight line at the boundary up to machine precision, see [35].

3.3.3 Numerical errors

It is well-known from the literature that the standard Youngs method (24) fails to reconstruct arbitrary straight lines, while the error is typically of the order of a few degrees. As a numerical test, a straight line is moved with a fixed inclination angle on an equidistant grid with $\Delta x_{1}=\Delta x_{2}$. This motion produces volume fractions ranging from [math] to $1$ for the considered computational cell away from the boundary. The reconstructed orientation with the standard Youngs and ELVIRA methods are shown in Figure 6(a). While the ELVIRA method always delivers the correct angle, the Youngs method shows an error of about $1^{\circ}-2^{\circ}$ in the considered example. The situation is much different for the same translation test for a boundary cell, see Figure 6(b). While the Boundary ELVIRA method is still able to deliver the correct orientation, the Boundary Youngs method shows a large error of up to $\pm 20^{\circ}$ that is also highly dependent on the position of the interface. Therefore, one can only expect a very rough estimate of the contact angle from the Boundary Youngs method which cannot converge with mesh refinement.

4 Results

To verify the advective transport of the contact angle, a spherical cap sitting at the boundary is initialized and transported using different velocity fields. To study the convergence in space and time, the time step $\Delta t$ is linked to the grid spacing $\Delta x$ by fixing the Courant number

[TABLE]

The influence of the choice of the Courant number is discussed below (see Figures 17 and 20). The following computational examples are carried out with (unless stated otherwise)

[TABLE]

Fixing the Courant number defines a temporal grid $\mathcal{T}$ (equidistant if $v$ does not depend on time). We report the error for both the contact line position and the contact angle in the maximum norm over all time steps, i.e.

[TABLE]

as a function of $\Delta x/R_{0}$. Note that the error in the contact line position is normalized by the initial radius $R_{0}$. The reference values $\operatorname{x_{cl,ref}}(t)$ and $\operatorname{\theta_{ref}}(t)$ come either from an exact or from a numerical solution of the ordinary differential equations (7) and (9). The numerical values for the contact line position $\operatorname{x_{cl,num}}$ and the contact angle $\operatorname{\theta_{num}}$ are evaluated directly from the reconstructed PLIC element intersecting the domain boundary (see Figure 7). To this end, the point of intersection of the local interface with the domain boundary is computed. If this point lies within the cell, the cell is recognized as a contact line cell and the contact angle and the contact line position are computed. Note that, due to the finite reconstruction tolerance of the VOF method (in this case $10^{-6}$) , irregular cases where no contact point is found may occur. An example is sketched in Figure 7(b), where the point of intersection lies slightly outside the current cell but the volume fraction of the neighbor cell is below the reconstruction tolerance so that it is not recognized as an interface cell. These irregular cases are excluded from the following error analysis.

Computational setup:

For the subsequent examples, we choose the following common setup. The computational domain $\Omega=[0,1]\times[0,0.25]$ is covered by an equidistant Cartesian grid of $N\times N/4$ cells, where $N$ varies from $128$ to $2048$. A spherical cap with dimensionless radius $R_{0}=0.2$ is initialized at the “solid boundary”

[TABLE]

where the flow is assumed to be tangential. The center of the sphere is placed at $(0.4,-0.1)$ yielding an initial contact angle of

[TABLE]

Since we are only interested in the local transport of the interface, we can allow for an artificial inflow boundary to the computational domain away from the contact line. Inflow boundaries are characterized by $v\cdot n_{\partial\Omega}<0$ (see Figure 8). At inflow boundaries we formally apply a homogeneous Neumann boundary condition for the phase indicator function, i.e.

[TABLE]

This condition is straightforward to implement using a simple constant continuation of the volume fraction field into a layer of “ghost cells” (see, e.g., [36]). Here we only consider the case where the interface does not meet the artificial boundary such that (33) simply states that no additional volume is transported into the computational domain. In particular, the boundary condition does not affect the dynamics of the interface and (9) still holds.

Three examples for the transporting velocity field are studied with the Youngs and ELVIRA methods, where these methods are combined with their newly developed boundary versions to treat the contact line advection.

4.1 Vortex-in-a-box test

We start with a classical test for interface advection methods given by

[TABLE]

This particular field called “vortex-in-a-box” has been routinely used to test numerical methods for interface advection; see [1], [2]. In the classical test this velocity field is used to strongly deform a sphere into a spiral. Due to the periodicity in time, it follows that the initial shape at $t=0$ and the final shape at $t=\tau$ would coincide if the problem is solved exactly. This allows to study aspects of the convergence behavior of the advection method even though the solution to the advection problem with the velocity field (34) is not known. The discrete $L^{1}$-error

[TABLE]

is usually used to quantify the rate-of-convergence. Note, however, that this kind of test does not say anything about the intermediate dynamics of the numerical solution555Formally, a numerical method which keeps the volume fractions fixed passes this test with zero error.. Here we revisit this classical test in the presence of a moving contact line. The results for $v_{0}=0.1$ and $\tau=0.2$ are reported in Figure 9. The simulations are carried out with a fixed Courant number of $\operatorname{CFL}=0.2$, where the numerical time step is chosen such that $t=\tau$ is reached after an integer number of time steps. As expected from the case without a contact line, the Boundary Youngs method shows a first-order convergence while the Boundary ELVIRA method is nearly second-order convergent.

Thanks to the kinematic evolution equation, it is also possible to study the dynamics of the advection in terms of the contact line position and the contact angle. The ordinary differential equations (7) and (9) are solved numerically to obtain reference solutions $\operatorname{x_{cl,ref}}(t)$ and $\operatorname{\theta_{ref}}(t)$.

Contact Line Motion:

The numerical evolution of the left (in this case the advancing) contact point reconstructed from the PLIC interface is investigated for the Boundary Youngs and Boundary ELVIRA method. It is found that both the Boundary Youngs and Boundary ELVIRA method deliver at least first-order convergent results for the motion of the contact line, see Figures 10 and 11.

Contact Angle Evolution:

The results for the contact angle over time are shown in Figures 12 and 13. While the numerical solution convergences to the reference solution for the Boundary ELVIRA method, the Boundary Youngs method does, as expected, not deliver mesh convergent results. In fact, one observes a strong oscillation of the reconstructed contact angle with a jump discontinuity when the contact line passes from one cell to the other; see Figure 14, where the reconstructed contact angle is plotted along with the discrete cell index. This behavior might be due to the spatial structure of the reconstruction error as reported in Figure 6. Clearly, the frequency of these jumps increases with mesh refinement leading to the strongly oscillatory behavior. The error in the maximum norm may even increase with mesh refinement, see Figure 12. Therefore, the Boundary Youngs method does not allow for a meaningful evaluation of the contact angle based on the local interface orientation even though it is first-order convergent with respect to the contact line motion and the discrete $L^{1}$-error regarding the initial and final shape comparison.

Following Figures 13 and 14, the evolution of the numerical contact angle for the Boundary ELVIRA method is reasonably smooth even on coarse grids. Some small oscillations are visible which, however, disappear with mesh refinement. In fact, the method shows a first-order convergence with respect to $\operatorname{\mathcal{E}_{\theta}^{\infty}}([0,0.5])$. The maximum error on the finest mesh with $\Delta x/R_{0}=5\cdot 10^{-3}$ is about $0.5$ degrees.

4.2 Linear Velocity Field

We now consider linear velocity fields of the form (11). In this case, the explicit solution (14) is available for verification. We choose the example

[TABLE]

The time evolution is investigated up to dimensionless time $T=0.4$. According to (14), the exact solution for the left contact point is given by

[TABLE]

where $x_{1}^{0}=0.4-\sqrt{0.2^{2}-0.1^{2}}\approx 0.227$ is the initial coordinate of the left contact point.

Contact Line Motion:

Like in the previous example, both methods show first-order convergence with respect to the maximum norm regarding the motion of the contact line, see Figure 15.

Contact Angle Evolution:

The numerical contact angle for the Boundary Youngs method is again subject to strong oscillations ($\pm 10^{\circ}$ in this case) and does not convergence with mesh refinement as visible in Figure 16(a). In contrast to that, the evolution of the contact angle for the Boundary ELVIRA method is first-order convergent and smooth with jumps visible only on a coarse grid, see Figure 16(b).

Influence of the Courant number:

The influence of the Courant number for the convergence of the Boundary ELVIRA method (with respect to the contact angle) and the Boundary Youngs method (with respect to the contact line position) is shown in Figure 17. Apparently, there is hardly any influence in the considered example. Linear convergence is achieved for all reported Courant numbers from $0.1$ to $0.9$.

Both numerical methods show excellent volume conservation. The relative volume error is at most of the order $10^{-10}$. Thanks to the volume redistribution algorithm, the volume fraction fields are exactly bounded up to machine precision.

4.3 Time-dependent Linear Velocity Field

As a third example, we consider the spatially linear, time-dependent velocity field of the form

[TABLE]

As mentioned before, the time-dependent coefficient $\cos((\pi t)/\tau)$ is a classical choice to test advection methods by comparing the phase volumes at $t=0$ and $t=\tau$. Here we consider the full dynamics of the advection by solving the kinematic evolution equation for the field (37) explicitly.

Using the ansatz (40), it is easy to show that the exact solution for the latter velocity field is given by (for $c_{1}\neq 0$)

[TABLE]

where $s(t)$ is defined as

[TABLE]

In particular, the evolution is periodic in $t$ with period $2\tau$. Note that the solution (14) is recovered in the limit $\tau\rightarrow\infty$ since

[TABLE]

As a concrete example, we choose again $v_{0}=-0.2$, $c_{1}=0.1$ and $c_{2}=-2$ together with $\tau=0.2$.

Contact Line Motion:

The contact line motion (see Figure 18) is first-order convergent for both methods.

Contact Angle Evolution:

Like in the previous examples, the numerical contact angle shows strong oscillations for the Boundary Youngs method, see Figure 19(a). The evolution of the numerical contact angle for the Boundary ELVIRA method is reported in Figure 19(b). While the numerical contact angle shows some deviations from the smooth reference curve on coarse grids, the period of the exact solution is still captured correctly. Like in the examples discussed before, refinement of the mesh at a fixed Courant number of $\operatorname{CFL}=0.2$ leads to smoothening of the results and first-order convergence in the maximum norm. The results show the ability of the Boundary ELVIRA method to accurately capture the dynamics of the contact angle evolution.

Influence of the Courant number:

The results for the time-dependent linear field (37) turn out to be much more sensitive to the choice of the Courant number, see Figure 20.

The Boundary Youngs reconstruction fails completely for $\operatorname{CFL}$ greater than $0.5$ due to the appearance of interface cells with $|\nabla_{h}\alpha|\approx 0$. The convergence of the contact angle in the maximum norm for the Boundary ELVIRA method breaks down for $\operatorname{CFL}\geq 0.7$ (see Figure 20(a)). Note, however, that a Courant number as large as $0.7$ is rarely achieved in multiphase flow simulations of systems governed by capillary effects (which is typically the case for wetting problems). In these systems, the numerical time step is usually limited by a stability criterion based on the propagation of capillary waves (see, e.g., [2]) and the (advective) CFL number is small.

5 Summary and Conclusion

The contact line advection problem is studied based on the geometrical Volume-of-Fluid method. Adaptations of the Youngs and ELVIRA methods to reconstruct the interface close to the boundary are introduced (see [35] for the implementations in FORTRAN). This allows to solve the transport equation for the interface without enforcing any boundary condition on the volume fraction field (except for inflow boundary conditions on artificial boundaries). Both the contact line position and the contact angle are evaluated based on the piecewise linear approximation of the interface (PLIC).

The Boundary Youngs method allows to track the motion of the contact line with first-order accuracy. However, a meaningful evaluation of the contact angle in terms of the local interface orientation is not possible. Instead, the numerical contact angle shows strong oscillations resulting from a spatial dependence of the reconstruction error which is already present for a planar interface.

The Boundary ELVIRA method delivers first-order convergent results for the dynamics of both the contact line motion and the contact angle evolution. The results are verified using an explicit and a numerical solution of the kinematic evolution equation (9).

Based on the Boundary ELVIRA method to transport the contact angle, one may develop numerical methods of dynamic wetting, where the contact angle is not prescribed as a geometric boundary condition. Instead, a local force term depending on the interface orientation may be introduced. The latter approach, as discussed e.g. in [19] for an algebraic VOF method, avoids the necessity to manually “adjust” the contact angle after a transport step. The present method allows to evaluate the local interface orientation with high accuracy in geometrical Volume of Fluid methods and may, therefore, serve as a basis for future numerical methods in dynamic wetting.

Acknowledgments:

We kindly acknowledge the financial support by the German Research Foundation (DFG) within the Collaborative Research Centre 1194 “Interaction of Transport and Wetting Processes”, Projects B01 and Z-INF. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt.

Appendix A Appendix

The following definition of a $\mathcal{C}^{1,2}$-family of moving hypersurfaces can also be found in [37], [20] and in a similar form in [38].

Definition 1**.**

Let $I=(a,b)$ be an open interval. A family $\{\Sigma(t)\}_{t\in I}$ with $\Sigma(t)\subset\mathds{R}^{3}$ is called a $\mathcal{C}^{1,2}$-family of moving hypersurfaces if the following holds.

1. (i)

Each $\Sigma(t)$ is an orientable $\mathcal{C}^{2}$-hypersurface in $\mathds{R}^{3}$ with unit normal field denoted as $n_{\Sigma}(t,\cdot)$. 2. (ii)

The graph of $\Sigma$, given as

[TABLE]

is a $\mathcal{C}^{1}$-hypersurface in $\mathds{R}\times\mathds{R}^{3}$. 3. (iii)

The unit normal field is continuously differentiable on $\mathcal{M}$, i.e.

[TABLE]

A family $\{\overline{\Sigma}(t)\}_{t\in I}$ is called a $\mathcal{C}^{1,2}$-family of moving hypersurfaces with boundary $\partial\Sigma(t)$ if the following holds.

1. (i)

Each $\overline{\Sigma}(t)$ is an orientable $\mathcal{C}^{2}$-hypersurface in $\mathds{R}^{3}$ with interior $\Sigma(t)$ and non-empty boundary $\partial\Sigma(t)$, where the unit normal field is denoted by $n_{\Sigma}(t,\cdot)$. 2. (ii)

The graph of $\overline{\Sigma}$, i.e.

[TABLE]

is a $\mathcal{C}^{1}$-hypersurface with boundary $\operatorname{gr}(\partial\Sigma)$ in $\mathds{R}\times\mathds{R}^{3}$. 3. (iii)

The unit normal field is continuously differentiable on $\operatorname{gr}\overline{\Sigma}$, i.e.

[TABLE]

Being the boundary of a submanifold with boundary, the set $\operatorname{gr}(\partial\Sigma)$ is itself a submanifold (without boundary).

Remark 1** (Construction of an analytical solution).**

Equation (13) may be solved with the following Ansatz: We look for solutions of the form

[TABLE]

It is an easy exercise to show that this yields the following ordinary differential equation for $f$:

[TABLE]

The initial condition for $\theta$ translates to

[TABLE]

For $c_{1}\neq 0$, the initial-value problem (41)-(42) has the unique solution

[TABLE]

Hence, we obtain the desired solution (14).