Boundary conditions for dynamic wetting -- A mathematical analysis

The moving contact line paradox discussed in the famous paper by Huh and Scriven has lead to an extensive scientific discussion about singularities in continuum mechanical models of dynamic wetting in the framework of the two-phase Navier Stokes equations. Since the no-slip condition introduces a non-integrable and therefore unphysical singularity into the model, various models to relax the singularity have been proposed. Many of the relaxation mechanisms still retain a weak (integrable) singularity, while other approaches look for completely regular solutions with finite curvature and pressure at the moving contact line. In particular, the model introduced recently in (Lukyanov, Pryer, Langmuir 2017) aims for regular solutions through modified boundary conditions. The present work applies the mathematical tool of compatibility analysis to continuum models of dynamic wetting. The basic idea is that the boundary conditions have to be compatible at the contact line in order to allow for regular solutions. Remarkably, the method allows to compute explicit expressions for the pressure and the curvature locally at the moving contact line for regular solutions to the model by Lukyanov and Pryer. It is found that the solution may still be singular for the latter model.


Introduction
Starting with the work by Huh and Scriven [1], the scientific discussion about singularities became central for the continuum mechanical modeling of dynamic wetting (see, e.g., [2][3][4][5][6][7][8]). While it is generally accepted that the non-integrable singularity in the viscous dissipation introduced by the no-slip condition (see [1]) is unphysical for a viscous fluid, there are different approaches to relax the singularity. For the sake of brevity, we only consider two particular approaches in this paper. See [4,5] and references therein for an overview of existing models. Note, however, that the method of compatibility analysis is general in nature and is applicable to other modeling approaches as well.
One of the most prominent choices is the Navier slip law which allows for tangential slip at the solid boundary according to where λ > 0 is a friction coefficient, w is the velocity of the solid boundary, n is the unit outer normal and D = 1 2 (∇v + ∇v T ) is the rate-of-deformation tensor. The parameter L = η/λ is called slip-length and controls the amount of tangential slip at a given shear rate. For quasi-stationary solutions, it has been shown by asymptotic methods [9] that a finite and positive slip length leads to an integrable singularity with the pressure behaving like where r is the distance to the contact line in a polar coordinate system. In this case, the force balance at the interface implies that the curvature is also infinite at the contact line to oppose the singular force due to pressure. Note that this irregular behavior of the solution might lead to problems when the model is solved numerically [10]. While some authors argue that this "weak" type of singularity has little influence on the macroscopic flow [11], there is also a large body of research looking for continuum mechanical models which do not even show weak singularities (see, e.g., [3,6,12,13]).
One approach chosen to regularize the moving contact line problem is to separately balance the mass contained in an interfacial layer. A framework for this kind of modeling is provided by the non-equilibrium thermodynamics of surfaces [14]. Mathematically, the mass within the interfacial layer is expressed as a density per unit area associated with a sharp interface in the continuum limit. The resulting model, known as the Interface Formation Model (IFM) [4,15], adds another level of complexity to the description since it requires to solve additional balance equations on moving interfaces. Physically, the idea is that the process of formation or disappearance of a piece of interface has a relaxation timescale which leads to dynamic surface tensions. Therefore, the model predicts a dynamic contact angle which is governed by a dynamic version of the Youngs equation, i.e.
where the surface tension coefficients σ i for the respective surface layers depend on the local state of the interface. Notably, the pressure and curvature at the moving contact line is claimed to be regular (see [4] for a detailed discussion).
In the present article, we study a recently proposed model by Lukyanov and Pryer [13], which can be understood as a quasi-stationary adaptation/simplification of the full IFM. The basic idea is to allow for non-zero interfacial mass densities but to require that these are constant in space and time. This assumption allows to substantially simplify the governing equations of the IFM. In particular, there is no need to solve the mass transport equation on the surfaces. Instead, the velocity associated with transport along the surfaces can be eliminated resulting in a modified set of boundary conditions for the stationary Stokes equations. Note that the fundamental difference to a model without interfacial mass is that the impermeability condition (1) at the solid boundary is relaxed because mass can be transported from the bulk to the surface phase. The model is introduced and applied to nanodroplets in [13]. Remarkably, it is stated that there "'a small, albeit natural, change in the boundary conditions is all that is necessary to completely regularize the problem" [13]. However, the authors neither prove the latter claim nor provide a numerical convergence study for the pressure and the curvature at the contact line.
Before we discuss the latter model in more detail, the method of compatibility analysis is introduced and applied to the "standard Navier slip model".

Compatibility analysis for the standard slip model
For the following calculations, we consider the geometrical configuration depicted in Figure 1, where for simplicity we only consider the case of two spatial dimensions. We choose Cartesian coordinates in a reference co-moving with the contact line, i.e. the solid boundary is moving with velocity (V 0 , 0) to the right, where V 0 is the speed of the contact line relative to the wall. Moreover, the interface normal and tangential vectors at the contact line have the form n 1 (0, 0) T = (− sin θ, cos θ), τ 1 (0, 0) T = −(cos θ, sin θ), where θ is the contact angle. To simplify the analysis, we assume that the solid boundary is flat, i.e. κ 2 = −∇ Γ2 · n 2 = 0.
Idea of the compatibility analysis: The basic idea of the compatibility analysis is to consider a local linear expansion of the (divergence free) velocity field at the contact line, i.e.
where V 0 = 0 is the contact line speed and L > 0 is the slip length. Therefore, the unknown coefficients c i are dimensionless. The crucial property is: If the solution is regular it has to satisfy all boundary conditions at the contact line, where the free surface and the solid boundary meet (see [16]). This requirement leads to a system of equations for the unknown coefficients c i . Note that the velocity itself and the velocity gradient at the contact line can be expressed as Clearly, higher order terms in (4) do not contribute to ∇v(0, 0). As we will see below, most of the boundary conditions only involve the symmetric part of (∇v), the latter tensor being the rateof-strain tensor D. At the contact line, it can be written in terms of the coefficients according to Due to the structure of D(0, 0), it is convenient to formally introduce a new unknown Application to the standard slip model in the free surface formulation: We consider all boundary conditions evaluated at the contact line.
(i) The kinematic boundary condition 1 at the contact line implies (in a co-moving frame) (ii) Due to the impermeability condition v · n 2 = 0 on Γ 2 (t) (evaluated at the contact line) it follows that c 4 = 0.
applies at the free boundary in the absence of surface tension gradients. Evaluating the latter condition at the contact line yields 1 Here V Γ 1 denotes the normal speed of the free surface. (iv) The Navier slip condition evaluated at the contact line reads as (v) The normal stress condition (for p ext = 0), i.e.
provides no additional information since both p and κ 1 are unknown.
Hence we obtain a system of 4 linear equations for 4 unknowns (c 1 , c 2 , c 35 , c 4 ). The linear system reads as  The determinant of the system matrix is given by For θ ∈ {0, π/2, π} the solution reads as Note that the coefficient c 2 becomes singular for θ → π/2. It is easy to show that no solution of the linear system exists for θ = π/2. The latter case is distinguished mathematically since the tangential stress and Navier slip conditions produce linearly dependent equations which are incompatible (note that c 1 = 0). For the limiting cases θ ∈ {0, π}, it is straightforward to show that (9) has a family of solutions given by Qualitative behavior of regular solutions: The evaluation of the boundary conditions at the contact line only provides information about v(0, 0) and D(0, 0). Since the impermeability condition v, n 2 = 0 applies along the whole solid boundary Γ 2 (t), it follows by differentiation with respect to τ 2 (since κ 2 = 0) that (∇v)τ 2 , n 2 = 0 on Γ 2 (t).
The latter condition evaluated at the contact line shows So we found the local linear expansion of the velocity field for a regular solution (if existent) to the standard slip model. It has been shown in [16], that the rate-of-change of the contact angle can be computed from ∇v at the contact line according tȯ Application of (11) to the solution derived above yields (for θ / ∈ {0, π/2, π}) The relation (12) shows that a (nontrivial) quasi-stationary solution of the model cannot be regular since V 0 = impliesθ = 0. Moreover, the relation (12) is unphysical since the contact angle cannot relax to equilibrium if the thermodynamic condition for the contact line velocity, given by (see [4,16,17] for a discussion of entropy production at the contact line) is satisfied. This result shows that regular solutions to the standard model are unphysical and a singularity must be present at the contact line even if the contact angle is variable (see [16] for a detailed discussion).

The model
The model introduced in [13] is based on the Interface Formation Model (IFM) due to Y. Shikhmurzaev.
The main simplifying assumptions are that the flow is quasi-stationary (i.e. the interface is fixed in a co-moving frame) and can be described by the stationary Stokes problem The interfacial mass densities are assumed to be constant, i.e.
This allows to eliminate the interface velocity v Σ from the model leading to a modification of the boundary conditions for the Stokes equations (see [13] for a derivation). In particular, the presence of interfacial mass allows for a non-zero normal component of the fluid velocity at the solid wall.
In the following, we abbreviate this model with the term "CSM model" (constant surface mass model). Note that, in contrast to the IFM, the CSM model does not assume the presence of surface tension gradients at the contact line. Instead, the dynamic contact angle is an input parameter for the model.
The CSM model introduced in [13] consists of the incompressible stationary Stokes equations, i.e.
and the boundary conditions given by v, −p + 2η n 1 , where P i = 1 − n i ⊗ n i is an orthogonal projection operator, L is the slip length, V 0 is the velocity of the solid wall in the co-moving reference frame of the contact line and Note that the ratios α i have the dimension of a length. Therefore, it is convenient to introduce the dimensionless quantitiesα Note that the normal stress condition (20) is formulated for a constant ambient pressure p ext = 0.

Compatibility analysis for regular solutions
We consider again the linear expansion of the velocity field at the contact line given by (4). The set of boundary conditions (15)-(20) leads to a system of algebraic equations that is discussed below.
Mass balance at the interfaces: Note that the mass balance equations (15) and (17) involving the surface divergence operator ∇ Γi · require some clarification. A short computation shows the relation Making use of the incompressibility condition 0 = ∇ · v = ∇ Γi · v + n i , (∇v)n i , we conclude that It follows that the mass balance equations (15) and (17) can be expressed as whereκ i = κ i L is the dimensionless curvature. We assume the solid boundary to be flat, i.e.κ 2 = 0.
Note that (21) implies that the system of equations (15) Mass balance at the contact line: Equation (19) expresses the fact that the contact line itself cannot store mass. In terms of the unknown coefficients it translates to (α 2 cos θ +α 1 ) c 1 +α 1 c 4 sin θ = −α 2 .
Zero stress condition and Navier Slip condition: Both the condition on the tangential stress (6) and the Navier slip condition (7) remain unchanged with respect to the standard model. Note, however, that in this case c 1 may be non-zero since the free surface is not a material interface.
Note that equation (23), expressing the mass balance in the free surface phase, is non-linear, while equations (24)-(27) constitute a linear subsystem Note also that, in the limit of vanishing surface mass (α 1 ,α 2 → 0), the nonlinear equation becomes (23) linear and the set of equations (23)-(26) reduces to (9) while (27) becomes obsolete. In this sense, the CSM model is a generalization of the standard slip model discussed in Section 2.

Solution of the non-linear system
The determinant of the system matrix A CSM is given as det A CSM = 2 (α 1 cos θ +α 2 ) sin 2θ −α 1α2 sin θ cos 2θ.
In the special case θ = π 2 this simplifies to det A CSM =α 1α2 . Clearly, the linear part of the problem is uniquely solvable in case det A CSM = 0. If, moreover, the solution satisfieṡ (30) Remarks: (i) It is easy to show that the condition is sufficient for det(A CSM ) = 0 on (0, π) for givenα 1 > 0.
(ii) A short calculation using (30) shows that equation (23) can be simplified according to The latter equation is central for the regularity of solutions to the CSM model as we will discuss below.
(iii) Provided that 0 <α 2 ≤α 1 2 , there is a unique contact angle given by which makes the right-hand side of equation (32) equal to zero. In this case, it follows that Otherwise (i.e. for θ = θ * andṁ 1 (θ * ) = 0), the above equation (32) becomes obsolete and the curvature at the contact line is not determined by the compatibility conditions.

Singularities in the model by Lukyanov and Pryer
We will now show that for any choice of surface mass densitiesα 1 ,α 2 > 0 satisfying the invertibility condition (31), there is a choice of the contact angle θ s < π 2 such that no regular solution exists for the parameters {α 1 ,α 2 , θ s }.

Pressure at the moving contact line
We can now evaluate the normal stress condition (20) at the free surface. If the curvature κ 1 (0, 0) is uniquely determined by the compatibility conditions, we also obtain the pressure at the contact point. The non-dimensional form of equation (20) readŝ Here we defined the non-dimensional quantities At the contact point we have Therefore, the dimensionless pressure at the contact point is given bŷ Ca − 2(c 2 cos 2θ + c 35 sin 2θ).
Hence the relation for the pressure (relative to the ambient pressure) reads as In particular, the pressure converges to the Laplace pressure p L = −κ 1 σ as V 0 → 0.
3.6 Special case θ = π/2 The system of equations (23)-(27) is substantially simplified in the case θ = π/2. In this case, the nonlinear equation (23) reads as and the matrix A CSM is given by The solution of the linear system is given by According to (36), it follows that the curvature at the contact line is So the curvature is always negative, i.e. the free surface is always convex locally at the contact line for regular solutions with θ = π/2. Moreover, the curvature becomes singular asα 1 → 0, which can be understood as the transition to the standard model. Interestingly, both the curvature and pressure at the contact line do not dependent on the surface mass densityα 2 in the liquid-solid phase.
For the dimensionless pressure at the contact line we obtain, according to (35), the relation Therefore, the pressure (relative to the gas phase) is zero for θ = π 2 and Ca = − 1 4 independently ofα 1 andα 2 . For the pressure in physical units we find Making use of the expression (37) for the mean curvature, we conclude p = −σκ 1 (1 + 4 Ca).

Comparison with results by Lukyanov and Pryer
In the following we revisit some examples given in [13]. The reported values for the interfacial mass densities areα 1 = 1.6,α 2 = 0.51. The latter values are obtained from Molecular Dynamics (MD) simulations (see [13] for details). The matrix is invertible for all θ ∈ (0, π) since (31) is satisfied. Therefore, the non-dimensional mass fluxesṁ i at the contact line can be computed from the general solution (30) (see Figure 3). It is observed that bothṁ 1 andṁ 2 become singular for θ → 0, π. Moreover, the mass fluxṁ 1 has two roots, corresponding to (see Section 3.4) The values for the dimensionless curvature and pressure at the contact line for the examples given in [13] are summarized in Table 1. For the cases (i)-(ii) and (iv) the sign of the curvature does not agree with the macroscopic form of the interface as reported in [13]. Hence, the present mathematical analysis predicts a bending of the interface close to the contact line even for the nanodroplet considered here. The absolute value of the curvature for the case (iii) is much larger then the macroscopic curvature of the interface reported in [13]. In fact, the contact angle θ = 65 • is close to the singular point θ 1 s ≈ 66.1 • . This is also the reason for the extremely low dimensionless pressure ofp ≈ −711 (measured relative to the ambient pressure). This value is about three orders of magnitude lower than the pressure reported in [13] (being approximately −0.9).

Conclusion
The method of compatibility analysis is applied to two continuum mechanical models of dynamic wetting, namely the standard Navier slip model and the model introduced by Lukyanov and Pryer [13]. It is shown that no quasi-stationary regular solutions with a moving contact line exist for the standard model. Even if the contact angle is allowed to vary, we showed that regular solutions behave unphysically (see also [16]). Therefore, some kind of singularity must be present at the contact line for the standard model.
The model [13] allows to store mass on the liquid-gas and liquid-solid interfaces (in the general framework of [15]). However, the surface mass density is assumed to be constant in space and time leading to modified boundary conditions. It is shown that the standard model is recovered in the limit of vanishing surface mass densities (i.e.α 1 ,α 2 → 0). The compatibility analysis for the model shows that regular solutions with finite pressure and curvature at the contact line are possible. In fact, the curvature and the pressure at the contact line can be computed from the compatibility conditions provided that the mass fluxṁ 1 to the free surface at the contact line is non-zero. On the other hand, if for certain model parameters the mass fluxṁ 1 goes to zero, then the solution becomes singular very much like in the standard model (whereṁ 1 is always zero). It is shown that such a singular point exists for every pair of interfacial mass densities (α 1 ,α 2 ) satisfying (31). Hence the model does not always "cure" the singularity completely. Explicit expressions for the pressure and the curvature at the contact line are derived for the special case of θ = π/2. Interestingly, the latter values for θ = π/2 do not depend on the surface mass densityα 2  [13].
of the liquid-solid interface.
The numerical values for the local curvature and the pressure at the contact line are compared with the continuum mechanical simulations of the model reported in [13] showing a quantitative and even qualitative disagreement. A possible explanation is a strong bending of the interface close to the contact line which is not resolved by the simulations in [13].
We emphasize that the method of compatibility analysis discussed here does only provide information locally at the contact line (or, in two dimensions at the contact point). There is no statement about the change in curvature or pressure in the vicinity of the contact line, which might explain the discrepancy to the results in [13]. However, the mathematical method of compatibility analysis is quite general and, therefore, applicable to a variety of models, possibly allowing for new insights into the complex problem of dynamic wetting.