A Hydrodynamic Model of Movement of a Contact Line Over a Curved Wall

TL;DR

This paper develops a two-dimensional hydrodynamic model for the movement of a contact line over curved surfaces, addressing stress singularities and improving boundary condition prescriptions in two-phase flow simulations.

Contribution

It extends the classical hydrodynamic solution to curved walls and proposes a method for applying non-singular boundary conditions near moving contact points.

Findings

The model accurately predicts contact line movement over curved surfaces.

Simulations show improved advection of contact points with the new boundary conditions.

The approach mitigates stress singularities in numerical simulations.

Abstract

The conventional no-slip boundary condition leads to a non-integrable stress singularity at a contact line. This is a main challenge in numerical simulations of two-phase flows with moving contact lines. We derive a two-dimensional hydrodynamic model for the velocity field at a contact point moving with constant velocity over a curved wall. The model is a perturbation of the classical Huh and Scriven hydrodynamic solution [11], which is only valid for flow over a flat wall. The purpose of the hydrodynamic model is to investigate the macroscopic behavior of the fluids close to a contact point. We also present an idea for how the hydrodynamic solution could be used to prescribe macroscopic Dirichlet boundary conditions for the velocity in the vicinity of a moving contact point. Simulations demonstrate that the velocity field based on the non-singular boundary conditions is capable of accurately advecting the contact point.

Figures

Equations

$$\frac{\partial {\text{v}}^{\prime }}{\partial t^{\prime }}+({\text{v}}^{\prime }\cdot \nabla ){\text{v}}^{\prime }+\frac{1}{\rho }\nabla p^{\prime }=\frac{\mu }{\rho }{\nabla }^2{\text{v}}^{\prime }-2{\mathbf{\Omega }}^{\prime }\times {\text{v}}^{\prime }-{\mathbf{\Omega }}^{\prime }\times ({\mathbf{\Omega }}^{\prime }\times r^{\prime }{\text{e}}_r).$$

$$\text{v}=U{\text{v}}^{\prime },r=L{r}^{\prime },t=\frac{L}{U}{t}^{\prime },p=P{p}^{\prime },\mathbf{\Omega }=\frac{U}{R}{\mathbf{\Omega }}^{\prime }.$$

$$\frac{U^2}{L}\frac{\partial \text{v}}{\partial t}+\frac{U^2}{L}(\text{u}\cdot \nabla )\text{v}+\frac{P}{L\rho }\nabla p=\frac{\mu U}{\rho L^2}{\nabla }^2\text{v}-2\frac{U^2}{R}\mathbf{\Omega }\times \text{v}-\frac{L{U}^2}{R^2}\mathbf{\Omega }\times (\mathbf{\Omega }\times r{\text{e}}_r).$$

$$\mathrm{Re}\frac{\partial \text{v}}{\partial t}+\mathrm{Re}(\text{v}\cdot \nabla )\text{v}+\frac{PL}{\mu U}\nabla p={\nabla }^2\text{v}-2\mathrm{Re}\frac{L}{R}\mathbf{\Omega }\times \text{v}-\mathrm{Re}\frac{L^2}{R^2}\mathbf{\Omega }\times (\mathbf{\Omega }\times r{\text{e}}_r).$$

$${\nabla }^4\psi =0.$$

$$v_r=-\frac{1}{r}\frac{\partial \psi }{\partial \theta }$$

$$v_{\theta }=\frac{\partial \psi }{\partial r}.$$

$$\begin{array}{rl}\theta & =\arcsin (\frac{r}{2R})+\pi \text{in A}\\ \theta & =-\arcsin (\frac{r}{2R})\text{in B}\end{array}\}.$$

$$\sin \theta =-\frac{r}{2R}\text{, }\cos 2\theta =1-\frac{r^2}{4R^2}\text{in both A and B}$$

$$\begin{array}{rl}\cos \theta & =-\sqrt{1-\frac{r^2}{4R^2}}\text{in A}\\ \cos \theta & =\sqrt{1-\frac{r^2}{4R^2}}\text{in B}\end{array}\}\begin{array}{rl}\sin 2\theta & =\frac{r}{R}\sqrt{1-\frac{r^2}{4R^2}}\text{in A}\\ \sin 2\theta & =-\frac{r}{R}\sqrt{1-\frac{r^2}{4R^2}}\text{in B}\end{array}\}.$$

$$\text{n}=(-\sin \theta ,\cos \theta )\text{in both A and B}\begin{array}{rl}& {\text{t}}_A=(-\cos \theta ,-\sin \theta )\text{in A}\\ & {\text{t}}_B=(\cos \theta ,\sin \theta )\text{in B}\end{array}\}.$$

$$\text{n}\cdot {\text{v}}_A=\sin \theta \frac{1}{r}\frac{\partial {\psi }_A}{\partial \theta }+\cos \theta \frac{\partial {\psi }_A}{\partial r}=0$$

$$\text{n}\cdot {\text{v}}_B=\sin \theta \frac{1}{r}\frac{\partial {\psi }_B}{\partial \theta }+\cos \theta \frac{\partial {\psi }_B}{\partial r}=0$$

$$\text{At the interface, i.e. θ=ϕ:}{v}_{A\theta }=\frac{\partial {\psi }_A}{\partial r}=0$$

$$\text{At the interface, i.e. θ=ϕ:}{v}_{B\theta }=\frac{\partial {\psi }_B}{\partial r}=0$$

$$\text{At the interface, i.e. θ=ϕ (r>0):}{v}_{Ar}=v_{Br}\iff$$

$$\frac{\partial {\psi }_A}{\partial \theta }=\frac{\partial {\psi }_B}{\partial \theta }$$

$$\text{At the interface, i.e. θ=ϕ (r>0):}{\tau }_{Ar\theta }={\tau }_{Br\theta }\iff$$

$${\mu }_A\frac{{\partial }^2{\psi }_A}{\partial {\theta }^2}={\mu }_B\frac{{\partial }^2{\psi }_B}{\partial {\theta }^2}$$

$${\text{t}}_A\cdot {\text{v}}_A=\cos \theta \frac{1}{r}\frac{\partial {\psi }_A}{\partial \theta }-\sin \theta \frac{\partial {\psi }_A}{\partial r}=U.$$

$${\text{t}}_B\cdot {\text{v}}_B=-\cos \theta \frac{1}{r}\frac{\partial {\psi }_B}{\partial \theta }-\cos \theta \frac{\partial {\psi }_B}{\partial r}=-U.$$

$$\text{v}=v^0+\frac{r}{R}{v}^1+\dots .$$

$$\psi (r,\theta )$$

$$+\frac{r^2}{R}(e+f\cos 2\theta +g\sin 2\theta +h\theta )+\dots ,$$

$$v_r^0=(b+d\theta -c)\sin \theta -(a+c\theta +d)\cos \theta$$

$$v_{\theta }^0=(a+c\theta )\sin \theta +(b+d\theta )\cos \theta$$

$$v_r^1=2f\sin 2\theta -2g\cos 2\theta -h$$

$$v_{\theta }^1=2(e+f\cos 2\theta +g\sin 2\theta +h\theta ).$$

$$\begin{array}{rl}\theta & =\pi +\frac{r}{2R}+\mathcal{O}(\frac{r^3}{R^3})\text{in A}\\ \theta & =0-\frac{r}{2R}+\mathcal{O}(\frac{r^3}{R^3})\text{in B}\end{array}\}$$

$$\sin \theta =-\frac{r}{2R}\text{, }\cos 2\theta =1+0+\mathcal{O}(\frac{r^2}{R^2})\text{in both A and B}$$

$$\begin{array}{rl}\cos \theta & =-1+0+\mathcal{O}(\frac{r^2}{R^2})\text{in A}\\ \cos \theta & =1+0+\mathcal{O}(\frac{r^2}{R^2})\text{in B}\end{array}\}\begin{array}{rl}\sin 2\theta & =0+\frac{r}{R}+\mathcal{O}(\frac{r^3}{R^3})\text{in A}\\ \sin 2\theta & =0-\frac{r}{R}+\mathcal{O}(\frac{r^3}{R^3})\text{in B}\end{array}\}.$$

$$M^1{Z}^1=F^1.$$

$$M^1=\left(\begin{array}{cccccccc}1& 1& 0& \pi & 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 0& 0\\ 1& \tilde{C}& \tilde{S}& \varphi & 0& 0& 0& 0\\ 0& 0& 0& 0& 1& \tilde{C}& \tilde{S}& \varphi \\ 0& 2\tilde{S}& -2\tilde{C}& -1& 0& -2\tilde{S}& 2\tilde{C}& 1\\ 0& {\mu }_A\tilde{C}& {\mu }_A\tilde{S}& 0& 0& -{\mu }_B\tilde{C}& -{\mu }_B\tilde{S}& 0\\ 0& 0& -1& -1/2& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 1/2\end{array}\right),$$

$$Z^1=\left(\begin{array}{c}{e}_A\\ f_A\\ g_A\\ h_A\\ e_B\\ f_B\\ g_B\\ h_B\end{array}\right)\text{ and }{F}^1=\left(\begin{array}{c}\frac{1}{2}(a_A+\pi c_A+d_A)\\ \frac{1}{2}(a_B+d_B)\\ 0\\ 0\\ 0\\ 0\\ \frac{1}{2}(b_A-c_A+\pi d_A)\\ \frac{1}{2}(c_B-b_B)\end{array}\right).$$

$$\tilde{S}=\sin 2\varphi$$

$$\tilde{C}=\cos 2\varphi .$$

Full text

††journal:

1 Introduction

Consider two immiscible and incompressible fluids. If the fluids are in contact with a solid there is a line where the interface separating the two fluids meet the solid. This line is called a contact line, and when it evolves with time we have a moving contact line problem. Such problems form an important class of two-phase flows and appear both in nature and in many industrial applications REW . Examples of phenomena in nature are raindrops falling on a window or water bugs resting on water surfaces LIU . Industrial applications include lubrication, inkjet printing, gas and oil recovery in porous media ZAHEDI ; ROCCA ; 1ROCCA ; 2ROCCA ; 3ROCCA ; 4ROCCA and the development of microfluidic devices such as micropumps and lab-on-a-chip devices MARTIN .

The standard mathematical model is the incompressible Navier-Stokes equations influenced by the interface through surface tension. The interface in its turn is advected by the fluid velocity. A main challenge in numerical simulations of moving contact line problems is that the adherence, or no-slip boundary condition leads to a non-integrable stress singularity at the contact line REW ; DUSSAN1 . With a no-slip condition an infinite force would be required to move the contact line. However, the surface tension force, which is believed to cause the contact line movement, is finite RENARDY . To accurately model contact line dynamics atomistic phenomena must be taken into account in a small region near the contact line.

In their paper Huh and Scriven HUH argue that even though the ultimate resolution of a moving contact line must rest on molecular considerations on a microscopic length scale, the problem can be approached through kinematics and dynamics of fluids. They mean that it is still instructive to use a hydrodynamical model to investigate the macroscopic behavior of the fluids close to, but not exactly at, a contact line. In HUH Huh and Scriven construct such a hydrodynamical model for the case of steady, two-dimensional flow with a perfectly flat fluid interface. The model is aimed at describing the flow at an intermediate length scale $L$, which is much smaller than the large scale features of the flow but still large enough for the atomistic phenomena to be negligible. The model is based on the assumption that at this intermediate scale the viscous effects dominates the convection and the creeping flow approximation of Navier–Stokes equations is valid. The analysis results in a family of solutions, which depends on three parameters: the velocity of the contact point (in 2D the contact line is reduced to a contact point), the contact angle (the angle between the interface and the wall), and the viscosity ratio of the two fluids.

The Huh and Scriven solution is referred to extensively in the fluid dynamics community, and it is recognized to be a useful tool for describing flow in an intermediate region near contact points, see for example the recent review paper by Snoijer and Andreotti Review2013 , and also Yulii ; DUSSAN1 ; DUSSAN . It is also well recognized that this solution has shortcomings. Firstly, the Huh and Scriven solution is singular at the contact point. It is by now well known that at the contact point the Navier–Stokes equations do not model the physics. In this small region atomistic phenomena come into play, and these are not included in the standard Navier–Stokes system. Even if the Huh and Scriven solution were regular there, it would not be physically relevant. The Huh and Scriven solution is relevant at the intermediate length scale, much smaller than that of main features in the flow, but larger than the atomistic scale. Secondly, the assumption of a planar interface is unrealistic. There is a jump in the pressure over the interface, which can only be balanced by surface tension of a curved surface. However, if the surface tension effect is strong, which is the case for capillary dominated flow, a small curvature suffices. In chen1997 a modified Huh and Scriven solution is presented, which reveals that at an intermediate length scale the curvature decreases away from the contact point, and the flow approaches the flow of the planar case. There is also experimental agreement. Further away from the contact point, the assumption of a flat fluid interface, as well as the assumption of a low Reynolds number, are in general no longer valid. The conclusion is that the planar interface and the Huh and Scriven solution are appropriate for matching an outer solution with an inner solution, at some intermediate distance from the contact point. Thus, important to note is that the Huh and Scriven solution on its own does not suffice for determining the movement of the contact point, it must be combined with a model taking the atomistic processes into account. Examples of such models are the phenomenologically based phase field model, or a molecular dynamics based model.

The main focus of this work is the modification of the hydrodynamical model derived by Huh and Scriven to be valid for a different geometry. In HUH the solution is derived for a contact point moving over a flat wall. A flat wall is often considered in numerical simulations of moving contact line problems because of its simplicity LIU . However, in many applications more complex geometries need to be taken into account LIU .

In the first part of this paper, we extend the hydrodynamical model presented in HUH to a circular solid boundary with radius of curvature $R$. If the curvature of a more complex boundary varies sufficiently slow, the boundary can be well approximated to have a circular shape locally at the contact point, and our model can be used. Further, by assuming $R>L$ we can neglect centrifugal and Coriolis forces. Since we are interested in the dynamics close to the contact point we use perturbation analysis over the small quantity $r/R$, where $r<L$ is the distance to the contact point. In the first part of this paper we determine the first two terms in the expansion. However, it is possible to extend the theory and calculations presented here to account for more terms in the expansions if required.

In the second part of this paper we show one possible example of how our hydrodynamical model could be used in numerical simulations of moving contact point problems. The focus is on how to incorporate small-scale effects on the macro-dynamics. The small-scale effects are assumed to be given by a relation between the macroscale wall contact angle and the contact point speed. We present numerical experiments designed to test how well a velocity field with boundary conditions based on the extended hydrodynamic solution advects the contact point.

2 PART @slowromancapi@: Derivation of the Hydrodynamic Model

In this section we will construct a hydrodynamical model for the velocity field close to a contact point moving over a curved wall. In Sect. 2.1 we start by describing our two-dimensional model problem, and discuss the conditions under which the creeping flow approximation (also called the Stokes model) can be used in the vicinity of the contact point. In Sect. 2.2 we introduce a perturbation expansion, and solve for the first two terms. Finally, Sect. 2.3 discusses resulting velocity fields.

2.1 The Contact Point Model Problem

Consider a model problem in two space dimensions with constant curvature of the solid wall and constant contact point speed $U$, see Fig. 1. A polar coordinate system $(r,\theta)$ is used, with the origin fixed to the contact point position, and with $\theta=0,\,180$ along the tangent of the wall at the contact point. This coordinate system is rotating with angular velocity ${\bf\Omega}$ of magnitude $U/R$, where $R$ is the radius of curvature of the wall. In this rotating frame of reference the contact point appears to be at rest and the wall to move with constant speed $U$ (but with opposite direction to that of the contact point in a fixed frame). The fluid interface between phase A and B is assumed to be perfectly flat and inclined at an angle of $\phi$ from the line $\theta=0$.

2.1.1 Creeping flow approximation

In a rotating frame of reference centrifugal and Coriolis forces appear in the Navier-Stokes equations:

[TABLE]

Here $\textit{{v}}^{\prime}$ is the fluid velocity, $p^{\prime}$ is the pressure, $\mu$ is the viscosity, and $\rho$ is the density. Introduce dimensionless quantities by

[TABLE]

Now the equations are

[TABLE]

With the Reynolds number $\mathrm{Re}=LU\rho/\mu$ we have (after multiplying with $\frac{\rho L^{2}}{\mu U}$)

[TABLE]

Close to the contact point viscous effects dominate over convection, and $\mathrm{Re}\ll 1$. Thus we neglect inertia terms (first two terms above). In addition, since $R>L$ we similarly neglect the centrifugal and Coriolis forces (last two terms). The result is that we can use the same creeping flow equations as in the flat wall case.

The creeping flow equations can be formulated to take the form of a two-dimensional biharmonic equation for the stream function $\psi(r,\theta)$ BIHARMONIC :

[TABLE]

In terms of the stream function the fluid velocity components in plane polar coordinates are:

[TABLE]

In our setting we will have one stream function in each fluid, $\psi_{A}$ and $\psi_{B}$, and corresponding velocities.

2.1.2 Boundary and Interface Conditions

To close the system (2) we need boundary and interface conditions. Before applying the boundary conditions, we parametrize the circular wall as:

[TABLE]

When applying the boundary conditions, the following reformulations of the parametrization in (5) will also be needed:

[TABLE]

[TABLE]

Further, the normals and tangents to the circle are described in plane polar coordinates by

[TABLE]

Note that the components of these vectors are in the radial and angular directions, respectively.

Now, the following boundary and interface conditions can be applied:

1. (i)

Kinematic conditions: vanishing normal component of velocity at the solid wall and fluid interface.

[TABLE] 2. (ii)

Kinematic condition: no slip, i.e. continuity of velocity at the interface.

[TABLE] 3. (iii)

Dynamic condition: continuity of tangential stress at the interface.

[TABLE]

Equations (10) and (11) have been used to simplify the second expression. 4. (iv)

Kinematic condition: no slip, i.e., no tangential relative motion of the fluid at the solid wall.

[TABLE]

2.2 Approximate Solution by Perturbation Analysis

In this section we use perturbation analysis to find an approximate expression for the velocity field close to a contact point at rest at a circular wall, which rotates with constant speed. We are looking for an expansion of the velocity field in terms of different powers of $r/R$,

[TABLE]

The starting point is the biharmonic equation (2) in plane polar coordinates for each fluid, together with the boundary and interface conditions (8) - (15). The velocity is given by derivatives of the stream function according to (3) and (4). The general solution of the biharmonic equation in plane polar coordinates was derived in a paper by Michell MICHELL by a separation of variables. The solution is in the form of an expansion containing terms with different powers of $r$. Here, just as in HUH we require that the velocity is bounded as $r\rightarrow 0$, and therefore all terms with negative powers of $r$ in the Michell solution are excluded. In HUH the expansion was to be used also far from the contact point, and therefore the velocity was required to be bounded also for large values of $r$. The consequence was that all terms with powers of $r$ higher than one in the expression for the stream function were also omitted. In our setting the expansion is only used close to the contact point. Further away it needs to be matched to an outer solution, and therefore there is no reason to omit terms with higher powers of $r$.

Based on MICHELL we use the following expansion of the stream function

[TABLE]

for each of the fluids A and B. For the stream function expansion to correspond to (16), we have scaled the constants $e,f,g$ and $h$ in the second term above by $1/R$, compared to the solution given in MICHELL (i.e. the constants $e,f,g$ and $h$ here equals the constants from the Michell solution multiplied by $R$). Further, in the paper by Michell the term $r^{2}\theta$ was not included due to an assumption of theta-periodicity, which is not assumed here.

Using the expansion for the stream function (17), the zeroth and first order components (in $r/R$) of the expansion for the velocity (16) are given by

[TABLE]

for each of the fluids A and B. We will in the following determine the coefficients in the two leading order terms of the expansions for each of the fluids by applying the boundary conditions (8) - (15). Note that it would be straight forward to continue the expansion to higher powers of $r/R$.

To consider the boundary conditions along the circular wall to different powers of $r/R$, we need to taylor expand the parametrization of the circular wall given in (5) and (6) around $r/R=0$:

[TABLE]

[TABLE]

[TABLE]

We are now ready to apply the boundary and interface conditions (8) - (15) to the expression for the stream function (17), and consider separatly different powers of $r/R$. The zeroth order equations involve $[a_{i},b_{i},c_{i},d_{i}]$, $i=A;B$, and are exactly the same as in HUH . It follows that the velocity field derived in HUH for a contact point at rest at a flat wall, which translates with constant velocity, is the unperturbed solution $v_{0}$. For clarity we have stated the equations for the coefficients $[a_{i},b_{i},c_{i},d_{i}]$, $i=A;B$, and their solutions, in the appendix. We note that the system is singular for $\phi=0,180$, and non-singular for other angles.

Next we collect terms that are linear in $r/R$ and get the following system for the coefficients $[e_{i},f_{i},g_{i},h_{i}]$, $i=A;B$ from the second term of (17):

[TABLE]

Here

[TABLE]

[TABLE]

We have used the notations

[TABLE]

The rows in system (23) originate from the boundary and interface conditions (8)-(15), respectively. This system is also singular for precisely $\phi=0,180$.

When all coefficients $[a_{i},b_{i},c_{i},d_{i},e_{i},f_{i},g_{i},h_{i}]$, $i=A;B$ have been calculated the resulting approximated velocity field is obtained using the expressions (18) - (21).

2.3 Resulting Velocity Fields and Verification

We look at the resulting approximated velocity field for three different cases. In all cases the contact point is stationary while the circular solid of radius $R=1$ is rotating with angular velocity $U/R$ where $U=1$. In the first case the angle between the fluid interface and the line $\theta=0$ is $\phi=45$ and the viscosity ratio is $\mu_{A}/\mu_{B}=1$. In the second case the viscosity ratio is increased to $\mu_{A}/\mu_{B}=100$. In the third case the higher viscosity ratio is used again, but the contact angle is now increased to $\phi=170$. The resulting velocity fields for the three cases are shown in Fig. 2, Fig. 3 and Fig. 4 respectively.

Increasing the viscosity ratio from $\mu_{A}/\mu_{B}=1$ to $\mu_{A}/\mu_{B}=100$ leads to a significant change in the flow pattern, compare case one in Fig. 2 and case two in Fig. 3. In the first case, when the two phases have the same viscosity, the flow along the circular boundary in phase A and along the interface both are directed inwards towards the contact point. This leads to the formation of a jet out from the contact point in phase A (directed upwards in Fig. 2). In the second case, when the viscosity ratio is increased to $\mu_{A}/\mu_{B}=100$, the direction of the flow along the interface has changed sign compared to the equal viscosity case. Consequently the flow along the circular boundary in phase B and along the interface are now both directed outwards from the contact point, which leads to a jet inwards towards the contact point in phase B (the region with lower viscosity) instead. In case three, when the contact angle is increased to $\phi=170$ but the viscosity ratio is kept at $\mu_{A}/\mu_{B}=100$, the same flow pattern with an ingoing jet towards the contact point in the phase with lower viscosity, phase B, is observed. This behavior agrees with the observations in HUH .

To verify the perturbation expansion approximation we consider how well the first two terms satisfy the boundary conditions. The velocity component tangential to the circular wall should be unity, $U=1$, due to the no slip boundary conditions (14) and (15). Further, the velocity component normal to the wall should be zero according to the boundary conditions (8) and (9). In Fig. 5 the error in the normal and tangential velocity components along the circular wall for case one are shown for both fluid phases A and B. The error in the velocity components decrease quadratically with decreasing $r$ ($R$ constant), which agrees with the perturbation expansion approximation (since we have neglected all terms containing powers of $\frac{r}{R}$ of two and higher). Similar second order behavior of the errors in case two and three are shown in Fig. 6 and Fig. 7 respectively. However, for case three, i.e. Fig. 7, a small deviation of the second order behavior of the error in the velocity component normal to the circular wall in phase B is observed for the larger values of $r$. The reason for the deviation in the second order behavior is that the error changes sign between $r=0.2$ and $r=0.3$. This can bee seen in Fig. 8 where the behavior of the velocity component normal to the circular wall in phase B is plotted. For the smaller values of $r$, approximately $r<0.25$, the normal velocity component is negative, i.e. too small, while for the larger values of $r$ it is too big. For $r>0.25$ the error again follows a second order behavior.

The systems for determining the coefficients in the stream function expansion is singular for $\phi=0$ and $\phi=180$. In case three when $\phi$ is close to $180$ an increase in the error is observed. Compare upper right figures in Fig. 6 and Fig. 7 for example. The error in the tangential velocity component along the circular wall for large values of $r$ in phase A is also clearly visible in Fig. 4.

For larger values of $r$ than included in the error plots presented here the assumption of a small $\frac{r}{R}$ is no longer valid, and the perturbation expansion approximation fails. Furthermore, far away from the contact point the viscous effects no longer dominate the convection, and the creeping flow approximation also fails. For small values of $r$ on the other hand, the velocity field contains large gradients and is therefore very sensitive to small changes in $r$ or $\theta$.

2.4 Change of Frame of Reference

The approximated velocity field can be transformed to a frame of reference where the contact point moves over a stationary wall. This is done by subtracting a rigid body rotation from the velocity field derived in the previous section. The rigid body rotation is centered at the center of curvature of the solid wall, and given by $r_{d}\,\Omega\,\textit{{e}}_{\Theta}$. Here $\Omega=U/R$ is the magnitude of the angular velocity, $r_{d}$ is the distance to the center of curvature of the wall, and $\textit{{e}}_{\Theta}$ is the unit vector in the rotational direction. The two latter quantities are

[TABLE]

As demonstrated in Fig. 9, the tangent to the circular wall at the contact point, i.e the line $\theta=0$, will not remain horizontal when the contact point moves over the wall (with the frame of reference of a stationary wall). Therefore we also need to take into account the angle of which the tangent is shifted from a horizontal line when evaluating the polar coordinate $\theta$. The requested $\theta$ is the angle to the horizontal line plus the angle of shift from a horizontal line of the tangent, here denoted by $\varphi$.

3 PART @slowromancapii@: Numerical Simulation of Moving Contact Points

In this section we demonstrate how the hydrodynamic model derived in previous section could be used in a multiscale model for simulation of moving contact point problems. In the standard model a non-vanishing contact point velocity is in conflict with the standard no-slip boundary condition for the fluid velocity. We will instead impose boundary conditions for the fluid at a small distance from the solid wall. The precise conditions are given by the hydrodynamic solution.

A basic assumption is that there is both a spatial and temporal scale separation between the local microscale contact point behavior and macroscopic fluid flow MARTIN2 . When the dynamics of the moving contact point is driving the flow, the assumption that the flow on the microscale is reacting much faster than the flow on the macroscale is justified MARTIN2 . This temporal scale separation implies that the microscopic dynamics is in equilibrium on the time scale of the dynamics of the macroscale. In particular, as the macroscale contact angle evolves, there is an equilibrium microscale contact point speed for each macroscale wall contact angle. Therefore the macroscale contact angle is the only information from the macroscopic model required to determine a contact point velocity. Consequently, there is a relation between the slip velocity of the contact point and the macroscopic wall contact angle and it is this relation that incorporates the microscopic effects at the macroscale.

Examples of micro-models that can be used for determining such relations between the contact point velocity and the wall contact angle are molecular dynamics REN1 ; REN2 or phase field models MARTIN2 . In order for the multiscale model to be efficient, the microscopic model can for example be used to pre-compute the contact point velocity $U$ for a set of contact angles. The focus here is on how to incorporate the microscopic effects, given by this relation, on the macroscopic dynamics. We will therefore not further discuss the microscopic origin of this relation.

In the proceeding sections we start by discussing the multiscale model in Sect. 3.1 and then present numerical experiments and results in Sect. 3.2.

3.1 Multiscale Model for the Moving Contact Point

A schematic illustration of the different scales of a moving contact point problem is given in Fig. 10. We introduce an intermediate region of length scale $L$, the red region in Fig. 10. If a typical length scale of features in the flow is $M$, then $L\ll M$. In this region standard continuum flow is still assumed, but with viscous effects dominating over convective effects, i.e. the Reynolds number $\mathrm{Re}\approx 1$. Therefore, the creeping flow approximation of Navier-Stokes equations is valid. Further, at this length scale, the fluid interface is to first approximation flat, which is illustrated in Fig. 10. Consequently, the hydrodynamic model is a good approximation in this region, and we can use the analytic expressions for the velocity derived in the first part of this paper. The assumptions on which the hydrodynamic model are based are not valid closer to the contact point, see the microscopic view in Fig. 10. There viscous bending becomes important EGGERS , and molecular phenomena such as diffusive mass transport come into play. These microscopic phenomena should be included in the coupled micro-model discussed above.

We will impose boundary conditions for the fluid velocity at a small distance from the solid wall, in the intermediate region. The precise conditions are given by the hydrodynamic solution.

3.1.1 Contact Point Boundary Conditions

The analytical velocity field given by the hydrodynamic model, is used to formulate boundary conditions for the velocity components. Since the hydrodynamic model is not valid exactly at the contact point, the innermost part of the intermediate region will not be included in the simulation. The solution will not be computed in the excluded region. The simplest approach is to use a computational domain that is $\delta$ smaller in the direction perpendicular to the solid wall, see Fig. 11. Along the new, artificial boundary, which is $\delta$ inside the physical boundary, we impose the analytical velocity from the hydrodynamic model as a Dirichlet boundary condition for the velocity components. The analytical velocity depends on the contact point velocity $U$, which is given by the relation between the macroscopic contact angle and the contact point velocity mentioned above. Thus, the information concerning the movement of one single point (in 2D), is transformed into a velocity boundary condition along the whole boundary. For the domain in Fig. 11 with $\delta=0.05$, the magnitude of the velocity at the artificial boundary is plotted in Fig. 12. When implementing this velocity boundary condition, care has to be taken when choosing the value of $\delta$. The creeping flow approximation is valid in a region with length scale $L$ (the intermediate region), and $\delta<L$ is required. However, the smaller $\delta$ the sharper is the peak in the boundary velocity function (see Fig. 12). Therefore, when discretizing in space the grid must sufficiently resolve the features of the velocity at the boundary, that is $h<\delta$ is required.

3.2 Numerical Experiments

The purpose of the numerical experiments is to investigate how well the velocity field based on the boundary conditions from the new hydrodynamic model advects the level set function and hence the contact point. We will consider low Reynolds number 2D flow in a symmetric channel with variable cross section, with an interface separating two identical fluids. We will assume the interface at each moment is in a quasi-steady state shape in the form of a circular arc. Such shapes are formed for instance when liquid rises in a narrow tube due to capillarity.

3.2.1 Two-Phase Flow Model

The motion of the two immiscible fluids is given by the incompressible Navier–Stokes equations for velocity u and pressure $p$ in non-dimensional form,

[TABLE]

Here, $\mathrm{\textbf{F}}_{st}$ is the surface tension force at the fluid interface. Further, $\nabla^{\mathrm{s}}\mathbb{u}=\frac{1}{2}(\nabla\mathbb{u}+\nabla\mathbb{u}^{T})$ denotes the rate of deformation tensor and the parameters $\rho$ and $\mu$ denote the density and viscosity measured relative to the parameters of fluid 1,

[TABLE]

The standard level set method presented in LEVELSET is used to keep track of the fluid interface and the moving contact point. The level set function $\phi(\textbf{x},t)$ is a signed distance function and the fluid interface $\Gamma$ is given by the zero level set of $\phi$, i.e. $\phi=0$ on the interface. The subdomain $\Omega_{1}$ occupied by fluid 1 is given by $\phi>0$ and the subdomain $\Omega_{2}$ occupied by fluid 2 is given by $\phi<0$. The level set function is advected in time by the fluid velocity according to the following Hamilton–Jacobi equation

[TABLE]

After advecting the fluid interface, the surface tension force $\mathrm{\textbf{F}}_{st}$ is calculated,

[TABLE]

where We is the Weber number and $\delta_{{\Gamma}}$ is a Dirac delta function with support on $\Gamma$. The normal and curvature of the interface are computed using the level set function,

[TABLE]

Over time the level set function will loose its signed distance property due to discretization errors and non-uniform velocity fields. To smooth the level set function and prevent the formation of large gradients, $\phi$ has to be reinitialized with a regular interval. In section Sect. 3.2.3 we present our model problem and a reinitialization method specially adapted for this problem. The standard technique for reinitialization (Sussman et.al REINIT ) consists of solving the following partial differential equation to steady state

[TABLE]

where $\phi_{0}$ is the level set function before reinitialization and $\tau$ is a pseudo time step. However, solving (27) may result in unphysical volume changes in the fluid phases and does not guarantee preservation of the contact point position. Accurate reinitialization for moving contact lines is a research area itself and is not addressed here.

3.2.2 Discretization and Implementation

We use the two-phase flow solver developed in MARTINBABA with suitable modifications to account for moving contact points. The solver is implemented in the C++ based Finite Element open source library deal.ii DEAL1 ; DEAL2 . Piecwise continuous linear shape functions on quadrilaterals, $Q_{1}$ elements, are used for the level set function and Taylor–Hood elements, $Q_{2}Q_{1}$ elements, are used for the Navier–Stokes equations (see MARTINBABA ). For time stepping, each of the level set equation and Navier–Stokes equations are discretized separately using the second order accurate implicit BDF–2 scheme. In order to avoid an expensive coupling between the incompressible Navier–Stokes part and the level set part (via the variables u and $\phi$) a temporal splitting scheme is introduced. After time discretization and linearization of the Navier–Stokes equations linear systems need to be solved. For the level set equation, a BiCGStab solver is used due to non-symmetry. The resulting system after discretization of the Navier–Stokes equations is of saddle point structure and solved by an iterative GMRES solver SAAD . For preconditioning, a block-triangular operator constructed using the so called Schur complement of the block system is applied from the right ELMAN . For more details about the two-phase flow solver we refer to MARTINBABA .

3.2.3 The Test Problem

The test problem consists of a circular interface moving in a channel with circular walls, see Fig. 13. The non-dimensional width of the domain is 6, the smallest distance between the walls is 2 and the radius of curvature of the circular walls is 20. Both fluids are assumed to have viscosity $\mu=0.7$ and density $\rho=1$, i.e. $Q=1$ in the hydrodynamic model.

The multiscale model does not depend on what microscopic model is used and we do not use a specific relation between the contact angle $\phi$ and velocity $U$, from a specific micro model. Instead we construct hypothetical examples for the relation between $U$ and $\phi$.

If no reinitialization is used the gradient of the level set function, $\nabla\phi$, gets distorted at the interface close to the contact point. We use the following simple reinitialization procedure:

The contact point position is evaluated by cubic interpolation using the level set function values at the two degrees of freedoms closest to the contact point in each direction along the wall. 2. 2.

The level set function is redefined to represent a distance function to a circular arc interface with the contact point position calculated in previous step.

We emphasize that this reinitialization procedure is adapted to our idealized test problem and is not useful in most other, more realistic settings.

3.2.4 Results

We perform two sets of simulations. In the first set the wall contact angle is $\phi=140$ and the initial contact point position is $x_{initial}=3$. In the second set $\phi=170$ and $x_{initial}=0$. We hypothetically relate both contact angles to a contact point velocity of $U=0.02$. We investigate spatial grid convergence by performing simulations using a set of structured meshes with different sizes. The smallest element length for each mesh is $h=1/24,1/32,1/40,1/48$ and $1/56$ respectively. The mesh where $h=1/24$ is depicted in Fig. 14. For $\phi=170$ the coarsest mesh $h=1/24$ is not used.

As explained in Sect. 3.1.1, when choosing $\delta$ (the distance between the artificial and physical boundary) we need to make sure $\delta<L$, where $L$ is the characteristic length scale of the intermediate region. For the creeping flow to be valid we need $\mathrm{Re}=\frac{\rho UL}{\mu}\ll 1$ or $L\ll\frac{\mu}{\rho U}$. For this model problem this implies the following condition on the distance to the physical boundary: $\delta<L\ll 35$, and we use $\delta=0.05$. The velocity functions that are used as boundary conditions at the artificial boundaries are shown in Fig. 15 and Fig. 16 for $\phi=140$ and $\phi=170$ respectively. We use homogenous Neumann boundary conditions for the velocity at the left and right boundaries.

The simulations are run for a non-dimensional time of $T=10$ with a time step size of $\Delta t=0.0001$ for the case when $\phi=140$ and $\Delta t=0.005$ when $\phi=170$. We are interested in the spatial errors when investigating the velocity boundary conditions and with this time step sizes the spatial discretization errors are assumed to dominate the temporal. The resulting velocity fields at $T$ for the finest mesh ($h=1/56$) are plotted in Fig. 17 and Fig. 18. It can be seen that away from the interface the flow profile is a regular poiseuille profile with velocity very close to zero at the artificial boundary. Close to the interface and the contact point however, the velocity is far from zero at the artificial boundary due to the velocity boundary condition.

The resulting contact point velocity in each time step is plotted as a function of time in Fig. 19. The period of the oscillations corresponds to the time it takes for the contact point to pass one grid cell, except for the coarsest mesh when $\phi=170$ (i.e. $h=1/32$). For this case it seems that the mesh is not capable of accurately resolving the velocity boundary function. The error in the velocity is higher for the case with the larger contact angle $\phi=170$. This agrees with the results presented in Sect. 2.3 where it can be seen that the error in the velocity field from the hydrodynamic model is higher for the case of a larger contact angle.

We measure the relative error in the contact point position at time $T$ compared to the correct position $UT$ for the different grid sizes. In Fig. 20 it can be seen that grid convergence is obtained with a rate of convergence of at least $2$ for both contact angles.

4 Summary and Conclusions

We have derived a two dimensional hydrodynamic model for the fluid flow near a contact point moving along a curved boundary with local radius of curvature $R$ at the contact point. The model is valid at length scales $L$ and by assuming $R>L$ we have neglected centrifugal and Coriolis forces. The model is based on the Stokes equations and was derived using perturbation analysis. As for the original flat wall hydrodynamic model HUH , the model is not valid in the immediate vicinity of the contact point, at lengths scales where molecular phenomena come into play.

We have also discussed a new idea for imposing non-singular Dirichlet boundary conditions for the velocity near a moving contact point (in a two-phase flow simulation). Numerical experiments are presented for a two dimensional model problem of two-phase flow with a moving contact point. Here the non-singular velocity boundary conditions are used together with a level set description of the interface. The results from these simulations demonstrate that the velocity field based on the non-singular boundary conditions is capable of accurately advecting the contact point.

The numerical simulations are based on the idea of excluding a small area around the contact point, and impose the non-singular Dirichlet boundary conditions at an artificial boundary. In this paper we have demonstrated and investigated the idea by excluding a thin region along the whole boundary. However, in a more realistic setting the artificial boundary should only exist along a $\mathcal{O}(L)$ part of the physical boundary. Further, we have applied a special reinitialization, which would not be applicable in most realistic cases. Investigating possible reinitilization strategies is a future research topic. Finally, it is important to note that the model for the moving contact point can also be combined with an interface tracking method as well.

Appendix

The equations for the coefficients in the first term in the stream function expansion (17) are written in matrix form as

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The rows in system (28) originates from the boundary conditions (8)-(15) respectively. The solution is

[TABLE]

where

[TABLE]

Note that the system is non-singular for $0<\phi<180$.