Solvability of a moving contact-line problem with interface formation for an incompressible viscous fluid

We investigate the free-boundary problem of a steadily advancing meniscus in a circular capillary tube. The problem is described using the “interface formation model,” which was originally introduced with the aim of avoiding the singularities that arise when classical hydrodynamics is applied to problems with a moving contact line. We prove the existence of an axially symmetric solution in weighted Hölder spaces for low meniscus speeds.


Introduction
When a liquid moves across a substrate surface, the line along which the free interface intersects with the surface is called the contact line. Free-surface flows with moving contact lines are commonly encountered in our daily lives. However, mathematical descriptions of such flows are beyond the capabilities of classical hydrodynamics. More precisely, if the classical no-slip condition is applied to a moving contact-line problem, a non-integrable stress singularity arises near the contact line (see, e.g., [5]).
The most common approach for avoiding this situation is to replace the no-slip condition with conditions that allow for slipping near the contact line (see, e.g., the review of slip models in Chap. 3 of [25]). This approach removes the above-mentioned stress singularity, but other issues remain unresolved. One such issue concerns the modeling of the dynamic contact angle. In slip models, the contact angle has to be prescribed, and in many cases it is assumed to be a function of the contact-line speed. However, a previous experimental study [2] suggests that the variation in the contact angle depends not only on the contact-line speed, but also on the flow field in the vicinity of the contact line. Another issue concerns the description of the flow kinematics. The experiments reported in [5] imply that a rolling motion arises in an actual flow. However, as reported in [24], two typical slip models (one that prescribes the fluid velocity on a solid surface and another that adopts the Navier slip condition) fail to reproduce the rolling motion.
posed. Following Shikhmurzaev's treatment in [24], a condition whereby the pressure is finite at the contact line is added. This treatment enables us to determine the contact angle as part of the solution without generating any singularity from the motion of the contact line. However, the case in which the contact angle is π/2 is an exception (see the next section in this regard). The pressure diverges at the contact line when the contact angle is in the range (π/2, π), because the contact line is a ridge of the domain. It is shown that the singularity arising in this problem is caused only by the presence of the corner of the domain, and the regularity of the solution near the contact line is determined by the size of the contact angle.
The remainder of this paper is organized as follows. In Sect. 2, we formulate the problem. In Sect. 3, we introduce function spaces and state our main result. In Sect. 4, we give some results for the Stokes equations with the slip boundary condition defined in a sectorial domain. The regularity of the solution to our problem near the contact line is determined from the distribution of eigenvalues of the operator pencil associated with this model problem. In Sect. 5, we derive an estimate of the solution to the linearized problem. The estimate is obtained by combining Schauder's method with an exponential decay estimate for the Dirichlet integral obtained by means of differential inequalities techniques. In Sect. 6, we solve the nonlinear problem using the contraction mapping principle with the estimates established for the linearized problem. We construct a solution in a neighborhood of the rest state of the fluid. We assume that the static contact angle differs from 0 and π , and we do not consider the situation where the dynamic contact angle reaches 0 and π . In our problem, the location of the interface is not known a priori; hence, for technical convenience, we transform the problem onto a prescribed domain with fixed boundaries. However, as the contact angle is unknown in advance, the regularity of the solution near the contact line cannot be determined directly from the transformed problem. The regularity of the solution is determined after the contact angle has been found by solving the transformed problem. Finally, in Sect. 7, the overall conclusions from this study and future tasks are described.

Formulation of the problem
Consider the situation in which a meniscus is moving at a constant speed W > 0 in an infinite circular tube = {|x | ≡ x 2 1 + x 2 2 < 1}. We assume that the total flux across an arbitrary cross-section of the part filled with the liquid is prescribed and given by ρ|B|W , where ρ is the density of the liquid and |B| is the area of the cross-section. In the formulation below, we assume that the flow is steady in the coordinate system moving with the meniscus.
First, we establish the governing equations. Let the region h = {x 3 < h(x ), x ∈ B} (B : |x | < 1) be filled with a liquid and h = {x 3 = h(x ), x ∈ B} be the free interface of the liquid.
We assume that the flow over h is governed by the steady Navier-Stokes equations where v is the velocity, p is the pressure, and ν is the kinematic viscosity.
For h , where n is the unit outward normal to h , we consider the following equations: where v 1 and ρ 1 are the velocity and surface density of the liquid-gas interface, respectively; T(v, p) = 2νD(v) -pI is the stress tensor, where is the strain-rate tensor; H is twice the mean curvature of h , which is assumed to be negative where h is convex; is the projection operator onto the tangent plane of h ; ∇ is the gradient operator defined on h ; and ρ 1,e , τ 1 , χ , γ , and ρ 0 are positive constants.
In particular, ρ 1,e is the equilibrium surface density. For , we assume that v · n = - where v 2 and ρ 2 are the velocity and surface density of the liquid-solid interface, respectively; W = (0, 0, -W ); n, , and ∇ are defined for in the same way as for h ; and ρ 2,e , τ 2 , α, and β are positive constants. In particular, ρ 2,e is the equilibrium surface density. In regard to the interfaces, (2) and (6) represent the momentum balance, (3) and (7) represent the mass balance, (5) and (9) are equations of state, and (4) and (8) are phenomenological laws between fluxes and thermodynamic forces arising in the entropy production rate of the interface (for details of the derivation, see [22,23,25]).
Remark 1 In the present paper, we assume that the right-hand sides of (4) 1 and (8) 1 vanish. In the original version of the interface formation model, Shikhmurzaev derived these conditions by assuming that the effect of mass exchange on the bulk is negligible (see, e.g., [23]). However, as revealed in [24], these terms play an essential role in describing the flow in the immediate vicinity of the contact lines more accurately. We have not yet succeeded in proving the solvability of the same problem when the right-hand sides of (4) 1 and (8) 1 do not vanish. The main difficulty lies in the derivation of an inequality corresponding to (54), which is essential in proving the solvability of the linear problem.
On M, we assume the conservation of mass, i.e., where e 1 and e 2 are unit normal vectors to M that are tangential to h and , respectively. We also assume the following force balance: where θ d is the contact angle determined by cos θ d = e 1 · e 2 , and σ sg is the solid-gas surface tension, assumed to be a nonnegative constant. Now, let us formulate our problem. After eliminating the surface velocity v 1 with the aid of (4) 2 and v 2 with the aid of (8) 2 , the above equations can be written as follows: Further, we assume that the liquid-solid interface reaches its equilibrium state and that the flow is fully developed in the far field, and we impose the following condition: Here, V is the Poiseuille flow given by the pair (V, P), where This Poiseuille flow was derived in [34]. As seen in (13) 3 and (14) 2 , the surface densities ρ 1 and ρ 2 satisfy second-order elliptic partial differential equations. Therefore, from a mathematical point of view, we need two boundary conditions on M to determine the surface densities using these equations. However, between conditions (15) and (16) on M, we can use only (15) for that purpose, because condition (16) is used to determine the contact angle. To use (16) as a boundary condition to determine the surface densities, we need to prescribe the contact angle, but doing so would eliminate the main advantage of this model. This "lack of a condition" was first pointed out in [1].
When the contact line is moving, the term βW arises in the generalized Navier condition (6). Shikhmurzaev [24] adds condition (50) (see [24]) to remove a pressure singularity caused by the above term and constructs an asymptotic solution without a pressure singularity. We follow the above treatment and assume his condition (50) in the following equivalent form: We can also find this condition from the right-hand side of conditions (13) 1,4 and (14) 1,3 by setting a 1 = γ ∇ ρ 1 · e 1 , b 1 = -( γ 2 ∇ρ 2 + β (v -W)) · e 2 , a 2 = b 2 = 0 and applying condition (26) to these terms (in this case the problem domain is not a sector, so the term ∂a 2 ∂τ is replaced by ∂a 2 ∂τ -v · ∂n ∂τ ).
Remark 2 Another approach for adding a condition at the contact line can be found in [1], where a boundary condition is added by accounting for the entropy produced at the contact line. Shikhmurzaev disagrees with this modeling approach from a physical point of view (see Sect. 4.3.5.3 in [25]).
Condition (20) ensures that the compatibility condition corresponding to (34) is satisfied for problem (91), which makes it possible to construct a solution without a singularity at the contact line when θ d < π/2. However, when θ d > π/2, the pressure becomes singular because the model problem (23) has singular eigensolutions. When θ d = π/2, in addition to (20), the condition that which comes from condition (30), has to be imposed for the pressure to be regular at the contact line. However, if we further assume condition (21), our problem becomes overdetermined. In addition, it is impossible to impose a condition for a specific contact angle in advance for a problem where the contact angle is not known. For these reasons, we cannot construct a regular solution when θ d = π/2. We hereafter denote the problem consisting of equations (12)- (20) for unknown functions (v, p, ρ 1 , ρ 2 , h) as problem (P).

Main result
To state our main result, we first introduce the function spaces used throughout this paper. Let D be a domain in R n , l be a nonnegative integer, and α ∈ (0, 1). We use C l+α (D) to denote the space of functions f (x), x ∈ D with the norms Below, we introduce weighted Hölder spaces. Let s ∈ R and M ⊂ ∂D be a closed set. Using For s < 0, we assume that C l+α s (D, M) = We also introduce the space C l+α,μ s ( , M) of functions defined on the lateral surface , which is defined in the same manner as C l+α,μ s ( h , M), except that h is replaced by . Finally, let W l 2 (D) be the function space with the norm The following is our main result.
, and that their supports are compact. Additionally, assume that the condition is satisfied, and, if λ 0 > 1, further assume that the following compatibility condition is satisfied: for a constant C > 0 that is independent of the data.
Rewrite problem (23) with ν = 1 in the polar coordinates (r, φ) defined by x 1 = r cos φ, x 2 = r sin φ, then apply the Mellin transform with respect to r: Thus, we have the following system of ordinary differential equations for u r = u 1 cos φ + u 2 sin φ, u φ = -u 1 sin φ + u 2 cos φ: We denote the operator of this boundary value problem by U(λ). If λ is not a root of (24), this problem has a unique solution. The roots of (24) are called eigenvalues of problem (28). A nontrivial solution corresponding to each eigenvalue is called an eigenvector.
and let λ 1 , . . . , λ n be the eigenvalues of problem (28) that exist in the range (s 1 , s 2 ). Then the solutions (u 1 , where c j are constants depending on the data, and (u(λ j ), q(λ j )) is the eigenvector of problem (28) corresponding to eigenvalue λ j .
, and that their supports are compact. Further, assume that compatibility conditions (25), (26), and for a constant C > 0 that is independent of the data.

Linear problem
In this section, we consider the following problem: where ν, γ , κ 1 , κ 2 , β, χ 1 , and χ 2 are positive constants, and is defined for a given function h ∈ C 3+α s+1 (B) in the same manner as h was defined in Sect. 2. Similarly, , , M, n, e 1 , e 2 , ∇ , and ∇ are defined in the same way as in Sect. 2. We assume that the domain and the data f, g, a, a 3 , b, b 3 , h 1 , h 2 , m 1 , and m 2 do not depend on the rotation angle φ around the x 3 -axis, and the φ-components of f, a, and b vanish in the cylindrical coordinates (r, φ, z) defined by (x 1 , x 2 , x 3 ) = (r cos φ, r sin φ, z). We denote the contact angle between and as θ and assume that θ = 0, π .
The following is the main result in this section.
Theorem 2 Let 0 < α < 1 and s be a constant satisfying the conditions s ≤ 2 + α, s / ∈ Z, and 0 < s < min(λ 0 , π/θ ), where λ 0 is the constant given in Lemma , and m 1 , m 2 ∈ C s (M). In addition, the following condition holds: For the case θ < π 2 , we also assume the compatibility condition at the contact line M: where τ is given by 2 ) in the cylindrical coordinates introduced above. Then there exists a constant μ > 0 such that problem (32) has a unique axisymmetric solution with the following properties: , and ρ 2 ∈ C 3+α,μ s+1 ( , M), and these functions satisfy the estimate for some constant C > 0 that is independent of the data.
We begin with the following lemma.

Lemma 4
Under the assumption stated in Theorem 2 for g, a 3 , and b 3 , there exists a vector field w ∈ C 2+α,μ and the estimate for a constant C > 0 that is independent of the data.
Proof We construct the function in the form w = ∇ , where is a solution of We first construct a function satisfying conditions (38) 2,3 and the estimate Then we introduce the new unknown -to reduce problem (38) to one in which a 3 = b 3 = 0. In the following argument, we again denote -and g -by and g, respectively.
Now, let us set and letḊ 1 2 ( ) be the space of all equivalence classes [u]. As shown in [6], the spaceḊ 1 2 ( ) is a Hilbert space with the scalar product We prove the weak solvability of problem (38) with a 3 = b 3 = 0 in classḊ 1 2 ( ) for arbitrary g ∈ L 2 ( ). Let ζ k (k ∈ N) be the cut-off function such that Set g k = ζ k g and let k be the corresponding solution of the problem with g = g k . We multiply the equation k = g k by φ ∈Ḋ 1 2 ( ) and integrate both sides to obtain The Riesz representation theorem then implies that there exists a unique solution k ∈ D 1 2 ( ) of (39), and we obtain the estimate Thus, by taking the limit k → ∞ in (39), we have a solution ∈Ḋ 1 2 ( ) satisfying the estimate For problem (38) with a 3 = b 3 = 0, we can obtain a decay estimate that is similar to (56). Based on this estimate, we can obtain estimate (37) using a method of localization similar to that used to obtain estimates (71) and (73). Because of the axial symmetry of the data and the domain, the solution constructed above is axially symmetric. Therefore, we can perform the above localization procedure for a two-dimensional problem using estimate (41) given below. In Lemma 5, d θ , γ θ , γ 0 , and M are defined in the manner specified in Sect. 4. As the result is well known, the proof is omitted.
, and that their supports are compact. Additionally, assume that the condition is satisfied. Then problem (38) with = d θ , = γ θ , and = γ 0 has a unique solution ∇ ∈ C 2+α s (d θ , M) that satisfies the estimate for a constant C > 0 that is independent of the data.
Thus, we have proved Lemma 4.
for a constant C > 0 that is independent of the data.
Proof The proof is similar to that of Lemma 4. The existence of a weak solution of problem (42) is easily shown. In addition, for the problem where χ ≥ 0 is a constant, we can obtain the following estimate in the same manner as in the proof of Theorem 7.1 in [30] under the assumption that the supports of the data are compact: Using this estimate, in addition to a decay estimate for s 1 , s 2 similar to that in (56), we achieve the desired result in a similar manner to that used to obtain estimates (71) and (73).
Let us prove the weak solvability of the problem. After multiplying (32) 1 by v , (32) 5 by ρ 1 , and (32) 8 by ρ 2 , we use integration by parts to obtain Let J( ) be the function space defined by We now prove the following inequality. The main difficulty is in the estimation of v 2, . A similar inequality in an infinite strip domain is derived in [18]. Because the domain is three-dimensional in the present case, more complicated arguments are necessary.

Lemma 7
For arbitrary v ∈ J( ), the following inequality holds: Proof For any v ∈ J( ), with the aid of Korn's inequality, we have As is compact, using a similar argument as in the proof of Lemma 4 in [33] implies that there exists a constant C( ) satisfying the inequality for arbitrary > 0. Combining these inequalities, we have Now, from the conditions ∇ · v = 0 and v · n| = 0, B v 3 (x , x 3 ) dx = 0 holds for arbitrary By integrating the identity with n = 2, where v = (v 1 , v 2 ). We integrate (48) and (50) with respect to x 3 to obtain Next, we integrate (49) with n = 3 in the domain 0,1 and find that and we note the estimate v(·, -1) 2,B ≤ C v (1) 2, 1 from (51) and (52) to obtain Estimate (46) is obtained by combining inequalities (47) and (53).
We now introduce the Lax-Milgram theorem.
Then, for every linear functional f : H → R, there exists a unique element u ∈ H such that Lemma 7 indicates that the bilinear form defined by the left-hand side of (45), which we denote as B[(v, ρ 1 , ρ 2 ), (v , ρ 1 , ρ 2 )], satisfies the condition Thus, from Theorem 3, we have the following.

Theorem 4
For arbitrary (f, a, b) ∈ L 2 ( ) × L 2 ( ) × L 2 ( ), there exists a unique solution for a constant C > 0 that is independent of the data.
For the solution (v, ρ 1 , ρ 2 ) obtained above, the pressure p is determined from the following equation defined for arbitrary η ∈ H( ) ≡ {f ∈ W 1 2 ( )|f · n| ∪ = 0} with a compact support: If we take the vector η such that where ⊂ is a bounded domain (for the construction of η, see, e.g., [13]), we have the estimate Here, p is normalized by the condition p dx = 0. The axial symmetry of the solution constructed above readily follows from the symmetry of the data and the domain. In addition, when f φ = a φ = b φ = 0, if we take v φ e φ (e φ = (-x 2 /|x |, x 1 /|x |, 0)) as v in (45) and take 0 as ρ 1 and ρ 2 , we have the equality The following decay estimate is obtained in a similar manner to the proofs of Theorem 5.3 in [27] and Lemma 1 in [19]. that f, a, for constants μ > 0 and C > 0 that are independent of the data.
We integrate (57) with respect to t over the interval (η -1, η)(1 < η) to obtain where The terms on the right-hand side are estimated as follows: where > 0 is an arbitrary constant. Thus, from (58), we have Multiplying this inequality by e -η/c 1 and integrating over the interval (2, -λ/2) with respect to η gives Now, we note the following inequalities: v (1) 2, (λ,1) (ii) the radius of U i (U i ∩ ∂ = φ) is sufficiently small that U i divides the surface ∂ into two connected parts; and (iii) there exists an integer N 0 such that the intersection of N 0 + 1 arbitrary different balls is empty.
We further introduce a partition of unity {ζ i } corresponding to the covering {U i }.
For U i with ξ i / ∈ M, the following estimates are obtained (for the derivation, see [11]): where D i and D i denote the domains D ∩ U i and D ∩ U i , respectively, for D = , , , and U i denotes a subset of U i on which ζ i = 1 holds. To estimate the solution defined on U i with ξ i ∈ M, we use the axial symmetry of the solution and rewrite problem (32) as the following two-dimensional problem in cylindrical coordinates: where on , (1, 0) on , τ = (τ r , τ z ) = (n z , -n r ), and f r and f z denote the r-component and z-component, respectively, of the term f.
We rewrite the above problem in the coordinates (y 1 , y 2 ), which are related to the original coordinates by y 1 y 2 = R r-1 z , where R is the matrix of rotation θ -(3/2)π around the origin. In the coordinates {y}, is given by the line y 2 = y 1 tan θ (or y 1 = 0 when θ = π/2), and we represent as the curve y 2 = g(y 1 ).

Nonlinear problem
In this section, we prove Theorem 1.
Thus, the functions (v, p, ρ 1 , ρ 2 , h) = (0,p, ρ 1,e , ρ 2,e ,h) describe the rest state of problem (P). Here,h is a solution of the equation where ∇ = (∂ x 1 , ∂ x 2 ) and ν is the unit outward normal to ∂ with the starting point located on M.
We seek a solution of the form where ζ is a smooth cut-off function satisfying In terms of the new unknowns in the above formulation, conditions (13) 2 and (16) can be written as and A x ∇ δh · ν| M = σ 3ρ 0 + ρ 2,e (ρ 0ρ 1,er 1 )(ρ 0ρ 1,e ) r 1 + 1 where A(x ) and B (j) (ξ ) are matrices whose (k, l)th components are ∂ ξ l F k | ξ =∇ h and ∂ 2 ξ k ξ l F j (ξ ), respectively, and F j (ξ ) = ξ j / 1 + |ξ | 2 (ξ = (ξ 1 , ξ 2 )). We now prove the following theorem, which will be applied to the problem consisting of (76) and (77). -1 (B, M), and g ∈ C(M), and thath, f , and g are axially symmetric with respect to the x 3 -axis. Further, assume that the compatibility condition is satisfied. Then the problem has an axisymmetric solution satisfying the estimate for a constant C > 0 that is independent of the data. This solution is uniquely determined if u(P) = 0 is satisfied for a point P ∈ M.
Proof When s > 1, the above result is readily obtained with the aid of a standard theory of second-order elliptic partial differential equations. Therefore, we consider only the case 0 < s < 1.
As the higher-order derivatives of the solution are estimated using Schauder's method, we derive only the following estimate: In the polar coordinate system, problem (79) is written as where H = {1 + ( ∂h ∂r ) 2 } 3/2 . The following function u satisfies the above problem: From this expression, we have the estimate The first term on the right-hand side is estimated as follows: The second term is estimated in the same way. Next, we derive the estimate u (0) The relation u (r) = ( H(r) r ) r 0 ξ f (ξ ) dξ + H(r)f (r) can be used to obtain The first term on the right-hand side is estimated as Now, let us estimate [u ] (s) I . For this purpose, we use the following two relations: and u (r + τ )u (r) = u (t)τ for a number t between r and r + τ .
Let X and Y ,δh be function spaces defined by where and K are positive constants. It is easy to show that the following lemmas hold. Note that an additive constant of q is determined from the compatibility condition of the problem consisting of (76) and (77). Thus, the following estimate holds for q: |q|

Lemma 9
Assume that (u, ∇q, r 1 , r 2 ) ∈ Y ,δh . Then, for δh ∈ X , the following inequality holds: where C is a positive constant that is independent of and K .
With the aid of the estimates in Lemma 10 and Theorem 2, we have the following lemma. The proof is almost the same as that of Lemma 7 in [11], and is thus omitted.
In the above equations,f andf denote the functions f • δh and f • δh , respectively;ē 1 and ∇h are defined for the surface y 3 =h(y ) in the same way as e 1 and ∇ δh , respectively, are defined;ˆ is the operator defined byˆ f = f -(n · f)n;∇ and∇ δh are the operators defined by∇ = ( t J δh ) -1 ∇ and∇ δh =∇ -(n ·∇)n;D(û) denotes a tensor with elements where A ij is the (i,j)th element of ( t J δh ) -1 ; and the operators˜ ,∇, and∇ δh and the tensorD( u ), which arise as a result of the change of variables x = δh (y), are defined in the same way.
For the terms on the right-hand sides of (102)-(111), we obtain the estimates given in the following lemma.
With the aid of the above estimates, applying Theorem 2 to the problem given in (102)-(111) yields the following estimate for sufficiently small and W : We now turn to the problem given by (76)-(77). For notational simplicity, we represent the problem as L(δh) = R(u, q, r 1 , r 2 , δh).
Let the mapping F map each δh ∈ X to the solution δh of the problem where (u, q, r 1 , r 2 ) ∈ Y ,δh is a solution to problem (91) defined for δh given above. Then, with the aid of Theorem 5 and estimate (93), we have Thus, if we choose C 1 as K and choose such that (K + K 2 ) < 1/2, we obtain | δh| (3+α) s +1,B,M ≤ (3/2)K . This implies that F maps X into itself.
Next, we consider the problem where (u , q , r 1 , r 2 ) ∈ Y ,δh is the solution to problem (91) with δh = δh . Then, with the aid of estimate (120), we have Thus, if we choose and W such that C( , W ) < 1, the estimate indicates that F is a contraction in X . Thus, by the contraction mapping principle, we have proved the existence of a solution to our problem.

Conclusion
A free boundary problem describing a steadily advancing meniscus in an infinite circular tube has been investigated. The problem was formulated using the interface formation model, and it was proved that an axisymmetric solution exists in weighted Hölder spaces when the velocity of the meniscus is low. The problem was closed by adding a condition stating that the pressure must be finite at the contact line. Our analysis shows that this treatment determines the contact angle as part of the solution and removes the singularity caused by the motion of the contact line (with the exception of the case in which the contact angle is π/2). It was also shown that the singularity of the solution is caused by the presence of a corner of the domain, and the regularity of the solution near the contact line is determined by the size of the contact angle. For a technical reason, this paper neglected the term representing the mass exchange between the bulk and the interface arising in the normal component of the velocity on the interface. Shikhmurzaev [24] has suggested that the corner of the domain causes no singularity if this term is not neglected-a rigorous verification of this point will be considered in future work.