Investigation of Wetting Hydrodynamics Using Numerical Simulations

Meniscus shapes from a simulation of a plate immersing into an inﬁnitely deep liquid bath, for a range of outer length scales, have been obtained numerically. These have been compared with the leading-order prediction from a three-region asymptotic analysis done in the double limit, Capillary number, Ca ! 0, L S / L C ! 0, with Ca ln( L C / L S ) of O (1), where L S and L C represent the slip length and an outer macroscopic length, respectively. For Ca (cid:44) 0.01, the numerically computed and the perturbation solutions show excellent agreement. Within this range of Ca, the meniscus slope at a distance 10 L S from the dynamic contact line is geometry independent, that is, does not vary with changes in the outer length L C . The interface slope at this point can serve as an appropriate material boundary condition for the outer problem. For 0.01 (cid:44) Ca (cid:44) 0.1, the intermediate region solution continues to closely ﬁt the numerically generated solution, while the match in the outer region begins to degrade. By monitoring the pressure difference between the surrounding inviscid gas phase and arbitrarily chosen point in the liquid, we attribute this breakdown to inﬁltration of viscous effects into the outer region, so that static capillarity does not adequately describe meniscus shapes in this regime. For Ca (cid:46) 0.1, there is no match between the numerical and perturbation solutions in both the intermediate and outer regions, indicating that higher-order contributions must be accounted for in the perturbation solutions.¬ © 1996 (cid:64)


I. INTRODUCTION
In both nature and in industrial processes, moving contact lines are ubiquitous, and they occur over a wide range of length scales-from microns for the movement of oil ganglia through porous rocks to meters in a variety of coating processes. A proper geometry-free characterization of the dynamic wetting process is therefore crucial if predictive models of these phenomena are to be developed. However, several complications arise when such a goal is attempted, largely because of the disparity between the submicroscopic length scales at which wetting occurs and the macroscopic length scales at which typical experimental observations are made. A further complication is that different local models yield the same dynamics at macroscopically observable levels. 1 Therefore, in order to verify refined models, measurements must necessarily be made on a submicroscopic length scale over which this physics remains important. This is a formidable challenge, partly because rapid changes in interface curvature occur near the three-phase juncture, and sophisticated experimental techniques become necessary to produce adequate resolution of data in that region.
Challenges also exist from a continuum modeling perspective. For a pure liquid, the customary hydrodynamic analyses 1 result in a multivalued velocity field at the dynamic contact line. If the liquid is Newtonian, a stress field is produced that is nonintegrable as the dynamic contact line is approached. One consequence of this stress divergence is an unbounded interface curvature at the dynamic contact line, resulting in the inability to specify a true dynamic contact angle. Because this angle serves as a boundary condition for the differential equation governing the shape of the liquidfluid interface, an ill-posed problem results. This difficulty is often resolved by breaking the flow domain into two regions. An inner one, in the vicinity of the dynamic contact line, where slip is permitted ͑thereby removing the source of the double-valued velocity at the dynamic contact line͒, has dimensions characterized by a slip length L S . An outer region, characterized by a much larger geometry-dependent macroscopic length scale L C , utilizes the no-slip dynamic boundary at the solid surfaces. Singular perturbation solutions to the field equations in each region ͑ignoring inertial effects͒ are matched in the overlap region to produce a uniformly valid solution over the entire domain. In the double limit Ca→0, ⑀→0, with Ca ln͑⑀ Ϫ1 ͒ of O(1), where Ca is the Capillary number U/␥, is the liquid viscosity, and ␥ is its surface tension, U is a characteristic speed, and ⑀ϭL S /L C , there is no overlap of the inner and outer regions, 2,3 and a third, intermediate region must exist between the inner and outer regions of expansion. In cylindrical polar coordinates, the inner, intermediate, and outer regions are characterized by dimensionless lengths rЈ/L S , Ca ln(rЈ/L C ) and rЈ/L C , respectively. ͓Variables with a prime ͑Ј͒ superscript are dimensional.͔ Analyses 2-5 of the resulting boundary value problem to O(1) yield a solution for the interface shape in the intermediate region of the form ͑assuming that the displaced fluid has zero viscosity͒ Here (rЈ) is the angle between the moving solid and the tangent to the liquid free surface at a distance rЈ from the a͒ Current address: Clean Sciences, Inc., 163 Whitney Place, Fremont, California 94539. b͒ Author to whom correspondence should be addressed. dynamic contact line, the subscript indicating the value of this angle in the indexed region. In addition, R represents this angle at some distance R located within the intermediate region. As seen in Eq. ͑1͒, the slope of the interface within the intermediate region is independent of macroscopic geometry, that is, it does not depend upon the outer length scale L C . Although the exact magnitude of R depends upon the specific inner model, the form of Eq. ͑1͒ is independent of it.
If the geometry independence of Eq. ͑1͒ could be confirmed, that is, if for a certain Ca, a single value of R can be used to correctly predict complete meniscus shapes for a range of outer length scales, then R could be used as a material boundary condition for the outer problem, replacing the true dynamic contact angle boundary condition. Cutting off the domain at rЈϭR and providing this slope characterization would permit a solution for all the field variables and the free-surface shape in the outer region without resorting to specifying any of the details of the fluid physics in the inner region ͑which are currently not well known͒. Clearly, this would represent a powerful advance in our ability to predict wetting behavior.
Two approaches can be taken to determine if Eq. ͑1͒ is valid. The first is to compare experimentally obtained meniscus shapes in the vicinity of the dynamic contact line, that is, in the intermediate region, with predictions from the perturbation analysis. For a given geometry and Capillary number, R can be chosen by minimizing the difference, in a least squares sense, between the experimental and analytically predicted meniscus shapes. This approach has been utilized by a number of different investigators. 4 -7 The geometry-free nature of this expression for the shape of the interface in the intermediate region for Caр10 Ϫ2 has been established. 6 However, some systematic deviations were observed for CaϾ10 Ϫ2 , which could be attributable either to the difficulty in obtaining accurate meniscus shapes in the proximity of the dynamic contact line or to a possible failure of the O(1) asymptotic analysis in this Capillary number regime. Clearly, to resolve this issue, more accurate experimental meniscus shape measurements are required in the vicinity of the dynamic contact line. Such experiments, using rods entering liquid baths at various immersion angles ͑as a means of varying geometry͒ have now been completed. 7,8 Viscous deformation has been observed at distances of the order of a capillary length scale from the dynamic contact line at Caу10 Ϫ2 , an observation that has direct implications for the measurement of ''apparent'' dynamic contact angles in capillary tubes. 7 Experiments conducted at Caϳ0.45 indicate that the breakdown of the model is caused by the inadequacy of the lowest-order perturbation solution in the geometryfree region, or because of contributions to the shape of the interface from the inner region. 8 Alternatively, numerical simulations of the boundary value problem can be used to obtain complete meniscus shapes, which can be compared with predictions from the perturbation analyses-this approach is used here. Because calculated meniscus shapes are not tainted by experimental errors, interface slopes accurate to a specified tolerance can be obtained over the whole domain, including the immediate vicinity of the dynamic contact line ͑at the expense, of course, of computation time͒. Eliminating gravity in the model, the outer length scale becomes the container length, which can be changed easily. Thus, the geometry independence of Eq. ͑1͒ can be conveniently and rigorously tested. Furthermore, the simulation can be readily extended to high Ca, where experimental measurements might become difficult, clearly identifying the range of validity of the perturbation analyses. These simulations can also be used to determine the limiting Ca value beyond which viscous effects start infiltrating the outer region, and static capillarity becomes inadequate for describing the outer meniscus shape.
The finite element 9-11 and finite difference 12,13 methods have been used for simulation of the steady motion of an interface between viscous liquids in a capillary tube. Comparisons of the simulations with experiments have been restricted to predictions of ''apparent'' contact angles. The goal of our work is to show that numerical simulations are a valid tool for understanding fluid physics in the vicinity of moving contact lines. Therefore, its scope is very different from previous numerical studies of wetting hydrodynamics.

II. THE MODEL
A prototype problem that captures the essential features of a dynamic wetting process is chosen for the numerical simulation, and is illustrated in Fig. 1. It is comprised of an incompressible liquid of infinite depth in a container of length L C , with the left sidewall moving into the liquid with constant velocity U. Gravity is ignored in the simulations, so that L C is the appropriate outer length scale. The upper boundary of the liquid, hЈ(xЈ), is a free surface, so that its location is not known a priori, but is obtained as part of the solution to the transport equations and applicable boundary conditions. A hypothetical lower boundary ͑shown by the dashed line, yЈϭ0͒ restricts the computational domain. At this boundary, the flow is assumed unidirectional, with zero net volumetric flow, and is sufficiently well removed from the free surface so that any further displacement has no impact on the shape of the free boundary. The volume/width of liquid within this computational domain is A. The static contact angle is S .
For a single component Newtonian liquid, the field variables ͑velocities, pressures, free surface heights, and gas pressure͒ depend upon the Reynolds number, Re, and the Capillary number, Ca. Inertial effects will be ignored in the numerical simulations presented ͑Reϭ0͒, conforming to the assumptions used in the perturbation analyses.
The dimensionless steady-state governing equations and boundary conditions, including the slip model, are provided in Appendix A.

III. NUMERICAL SOLUTION TECHNIQUE
Following a well-developed procedure, the Galerkin Finite Element technique is used to discretize the continuity and conservation of linear momentum equations, as well as the kinematic condition. 14 -16 The latter is distinguished in this scheme as the equation to be used for computing the free surface shape. The details leading to the appropriate weak forms of each of these equations are presented elsewhere. 16 The contact angle conditions replace the kinematic condition for calculating free surface heights at the contact lines. The unknown pressure in the inviscid gas phase, relative to a datum pressure at the lower right corner in the liquid is obtained using the volume constraint.
The liquid domain is subdivided into N x ϫN y elements. Vertical spines originating at the lower hypothetical boundary form element borders, while the ends of these spines are used to represent the free surface. Element corner nodes are located proportionally along these spines. The velocity fields within each element are approximated by nine-node Lagrangian biquadratic basis functions, while the pressure field is approximated by four-node bilinear basis functions. The free surface location is expanded using one-dimensional quadratic basis functions. These choices have proven convergent for Newtonian flows, and provide C 0 continuity for all variables across interelement boundaries. As is customary in finite element practice, the basis functions are developed on a square parent element in an ͑,͒ Cartesian coordinate system. This parent element is transformed onto the deformed quadrilateral element in the real domain through the use of isoparametric mapping. With this mapping, the free boundary coincides with ϭϩ1 for each element bordering the surface. Isoparametric mapping also facilitates evaluation of the unit normal and tangent vectors. 14 The residuals are calculated using the four-point tensor product Gaussian quadrature.
The discretization results in as many nonlinear algebraic equations as the number of unknowns ͑the unknowns consist of values of the two components of the velocity, the pressure, and free surface heights at appropriate finite-element nodes, and the gas pressure͒. This algebraic equation set is solved by Newton's method. The linear equation set to be solved at the nth Newton iteration is corr n J M n ϭϪR n ,¬ ͑3͒ where corr n ϭs nϩ1 Ϫs n .¬ ͑4͒ The elements of vector R n are the residuals of the weak forms of the equations, the elements of vector s n are the values of the unknowns at the nth Newton iteration, and the Jacobian matrix J M n ϭ‫ץ‬R n /‫ץ‬s n . The elements of J M n are obtained numerically using a one-sided forward difference scheme. 15,16 At each Newton iteration, Eq. ͑3͒ is solved by frontal elimination. 17 The Newton iterations are stopped when the L 2 norm and the L ϱ norm of R n are below 10 Ϫ6 . All the unknowns are updated at each Newton iteration, producing nearly quadratic convergence. Starting from initial guesses where only the known essential boundary conditions are specified, along with a static meniscus shape and zero values for all variables at all other nodes, convergence is achieved within nine iterations.
The convergence of numerically generated solutions with grid size has been carefully monitored, and details are provided below. This is especially critical in moving contact line problems, because large changes in velocity occur over dimensions comparable to the slip length, and an adequate tessellation must be provided to capture these features. Furthermore, roundoff errors associated with discretization of the transport equations and boundary conditions automatically produce slip, so that the force singularity referred to earlier is evident only when a plot of the force on the moving plate versus mesh size diverges as the mesh size is reduced. 16 Therefore the mesh must be highly refined in the region around the dynamic contact line; we use a minimum of five nodes within a slip length.
The lower hypothetical boundary is displaced successively larger distances from the free boundary ͑through adjustment of the volume/width, A͒ until the L 2 norm and L ϱ norms of the slopes of the free surface from two successive values of A are Ͻ10 Ϫ8 ͑note: because the piecewise quadratic expansion used for the free surface shape only guarantees C 0 continuity at interelement boundaries, slopes are calculated at Gauss points on the element surfaces͒. Our results indicate that this distance should be a minimum of 4L C . We have used Aϭ4L C 2 in all the simulations. The algebraic grid generation scheme, using proportionally spaced spines, proves adequate for generating converged solutions for ⑀у10 Ϫ4 . The requirement that at least five nodes be placed within a distance of one slip length from the dynamic contact line implies a large disparity between element sizes in different parts of the computational domain. For ⑀Ͻ10 Ϫ4 , when the free surface nodal locations are perturbed in preparation for evaluation of the appropriate entry in the Jacobian matrix, the smallest elements bordering the free boundary get sufficiently distorted, so that the mapping to the isoparametric domain becomes singular. Furthermore,

304¬
Phys. Fluids, Vol. 8, No. 2, February 1996¬ Finlow, Kota, and Bose the calculations presented here are restricted to surface shapes that are single valued in this simple representation using spines. Therefore this simple discretization scheme will certainly be inadequate at ''low'' or ''high'' contact angles, or where viscous deformation causes some part of the free surface to become tangent to a vertical spine. Alternate mesh generation schemes will then be required, but that remains outside the scope of this paper.

IV. RESULTS
As is customary in numerical simulations of transport problems, we demonstrate first that the numerical solution presented is not dependent upon the mesh size. Figure 2 illustrates the location and values of v x maximum , v x minimum , and v y maximum for successively refined meshes ͑23ϫ23, 25 ϫ25, and 26ϫ26-the last mesh is shown in Fig. 3͒. These variables approach constant values as the mesh is refined. We have used the 26ϫ26 nonuniform mesh shown in Fig. 3 for the rest of our computations. This results in 6401 unknowns.
Since the location of the intermediate region is not known a priori, the numerically generated solution is compared with a composite solution from the perturbation analysis. The outer solution to O(1), characterized by a constant curvature, is given by where 0 ϭg Ϫ1 ͓g͑ R ͒ϩCa ln͑L C /R͔͒. The uniformly valid composite solution is obtained by summing Eqs. ͑1͒ and ͑5͒ and subtracting the common part 0 , and becomes For given inner model parameters ⑀ and D , and Capillary number, the simulations have been run for a range of outer length scales 100L S ϽL C Ͻ2000L S . Choosing R to be 10L S ͑R must be a point within the intermediate region͒, the value of R for each L C is obtained by a single parameter least square minimization of the difference between the numerically generated solution and Eq. ͑6͒ for r 1 ϽrЈϽr 2 . The upper and lower bounds on rЈ are intended to limit the fitting to the intermediate region because R is an intermediate region parameter. In our simulations, r 1 ϳ5L S and r 2 ϳ0.25L C ͑we have varied r 2 by Ϯ20% and found no significant change in the value of R ͒. The quality of the fit, using the L 2 norm of the numerical and composite solutions, is used to establish an upper bound for the Capillary number, above which this fitting procedure is no longer meaningful.
Computed meniscus shapes for a wide range of Ca and a fixed value of container length, L C ϭ2000L S , are shown in Fig. 4. At low Ca, the interface shape is dominated by surface tension, and is marked by near constant curvature. As Ca increases, viscous effects on the meniscus shape become more apparent-a rapid change in shape in the immediate vicinity of the moving contact line is followed by a region of low curvature. In order to obtain a more detailed understanding of the nature of these shapes, meniscus slopes are calculated at Gauss points along the free surface, and compared with the predictions from the perturbation analyses. Figures 5͑a͒-5͑d͒ show the numerical solution, the best fit composite solution, as well as the outer and intermediate solutions for L C ϭ2000L S and a range of Ca. As illustrated in Figs. 5͑a͒ and 5͑b͒, the fit over the full container length is excellent for Caр10 Ϫ2 . Figure 5͑c͒ reveals that at Caϭ10 Ϫ1 , the fit within the intermediate region is good, while the numerically generated meniscus shape shows a lower curvature in the outer region than that predicted from the analysis. The leading-order perturbation analysis appears to adequately describe the strong viscous effects in the proximity of the contact line, but viscous effects, not accounted for in the leading-order solution, are starting to infiltrate the outer region. At Caϭ3.5ϫ10 Ϫ1 , there is a complete breakdown of the fit between the numerically generated and the leadingorder perturbation solution over the entire domain, as shown in Fig. 5͑d͒. At this Ca, higher-order corrections to both the intermediate and outer solutions must be included in order to adequately describe meniscus shapes. The quality of the fit for the numerical and best fit composite solutions ͑fitted in the intermediate region͒ for both the intermediate and outer regions at each Ca, is assessed by calculating the L 2 norms, using where m is the number of data points used in each region. The results are shown in Fig. 6. For Caр10 Ϫ2 , the difference  The infiltration of viscous effects into the outer region is further illustrated in Fig. 7, where the pressure in the inviscid gas phase P g , relative to a datum pressure in the liquid at the lower right corner of the domain, is plotted as a function of Ca. Our calculations have revealed that this pressure is invariant with the outer length scale L C for all Ca. At low Ca, surface tension dominates the shape of the meniscus. Our scaling then requires that P g ϳCa Ϫ1 ͑cos D ϩcos S ͒ for CaӶ1, where D and S represent the contact angles at the moving and stationary plates, respectively. In our simulations, D ϭ S ϭ45°. At large Ca, viscous effects dictate the shape of the free boundary and P g becomes independent of the Capillary number. The limiting value is P g ϳϪ6y for Caӷ1, where yϭA/L C 2 . Both the low and high Ca asymptotes are shown, along with the results from the numerical simulation ͑the constant difference between the numerically calculated values of P g , and the asymptote at high Ca is of the order of Ϫ2, and represents the pressure drop caused by the presence of a bounding liquid surface͒. Viscous stresses start becoming important beyond Caϳ10 Ϫ2 , producing the observed lack of fit with the leading-order solution in the outer region. Note that because Aϭ4L C 2 ͑A will always be some multiple of L C 2 ͒, these results are independent of L C , implying that this limiting value of Ca is independent of macroscopic geometry.
To assess the geometry independence of Eq. ͑1͒, we have plotted, in Fig. 8, the best fit values of R vs Ca for a range of outer length scales between 100L S and 2000L S . Here R appears to be independent of the outer length scale L C for Caр10 Ϫ1 . Within this range of Ca, the first term in Eq. ͑6͒ is at least one order of magnitude larger than the second at rϭ10L S , confirming that the meniscus shape at this point is dominated by the geometry-free viscous contribution. For CaϾ10 Ϫ1 , this plot shows an apparent geometry dependence. However, as discussed above, the fit has also degraded substantially, so that those R values may no longer be meaningful.
The geometry-independent meniscus shape in the intermediate region should be independent of the specifics of the  7. Variation of the pressure in the inviscid surrounding phase ͑-͒, relative to an arbitrary datum pressure of 100 at the lower right corner of the domain, with Capillary number. The low Ca asymptote, shown as the inclined dashed line, is given by P g ϳCaϪ1͑cos S ϩcos D ͒, while the high Ca asymptote, shown as the horizontal dashed line, is given by P g ϳϪ6A/L C 2 . The small constant difference between the simulated and asymptotic results at high Ca is the end effect caused by the presence of the bounding liquid surface. Both asymptotes are independent of L C . The calculated P g values are also independent of L C , indicating that these results are independent of macroscopic geometry. Viscous effects infiltrate into the outer region for Caу10 Ϫ1 . inner model. In order to test this directly, we have completed a set of simulations for L C ϭ200L S and 2000L S using an entirely different slip behavior, 18  The geometry independence of R for Caр0.1 implies that it can be used to replace the true dynamic contact angle as a material boundary condition when modeling fluid flow problems containing dynamic contact lines within this range of Ca. Essentially, the troublesome region containing the dynamic contact line can be eliminated, and the modeling can proceed on a truncated domain, to which geometry-free boundary conditions along with the usual hydrodynamic conditions can be applied.

V. CONCLUSIONS
Numerical simulations of a model problem containing a dynamic contact line have been utilized to establish the geometry-independent nature of the slope of the free surface in the immediate vicinity of the moving contactline for Caр10 Ϫ1 . Within this region, local viscous effects have a strong impact on the shape of the free surface, so that the macroscopic geometry does not play a role. The slope of the interface at any point within this region can be used as a material boundary condition for simulation of dynamic wetting problems. For 10 Ϫ2 ϽCaр10 Ϫ1 , viscous effects infiltrate into the outer region, significantly diminishing the quality of fit between the numerically computed and lowest-order perturbation solutions in this region. For CaϾ10 Ϫ1 , higher-order corrections to the perturbation solutions appear necessary for both the intermediate and outer regions to adequately described meniscus shapes.

ACKNOWLEDGMENTS
This work was supported by a grant from the Pittsburgh Supercomputer Center ͑CBT 910043P͒. We are grateful to E. Ramé, S. Garoff, and E. B. Dussan, V, for several insightful discussions.

APPENDIX: DIMENSIONLESS GOVERNING EQUATIONS
All lengths are scaled with L C , velocities with the plate speed U, and the pressure and viscous stresses with U/L C . All variables presented here are dimensionless. Continuity: "-vϭ0.¬ ͑A1͒ Here vϭv x e x ϩv y e y , where e x and e y are unit vectors in the x and y directions. Conservation of momentum (creeping flow, no gravity): Ϫ"Pϩ"-ϭ0, ͑A2͒ ϭ"vϩ͑"v͒ T is the viscous stress tensor.

Boundary conditions
No slip and no penetration of liquid on all stationary solid walls,

vϭ0.¬ ͑A3͒
At the free surface yϭh(x), n-TϭCa Ϫ1 2HnϪnP g ,¬ ͑A4͒ n-vϭ0.¬ ͑A5͒ Equations ͑A4͒ and ͑A5͒ represent the stress balance and kinematic conditions, respectively. Here TϭϪPIϩ, I is the identity tensor, n is the unit outward pointing normal from the free liquid surface, H is the local mean curvature of the interface 2HϭϪ" II -n, where " II is the surface divergence operator ͑IϪnn͒-", and P g is the unknown pressure in the surrounding inviscid phase, relative to an arbitrarily chosen datum pressure at the lower right-hand corner of the domain.
At the contact lines: At the dynamic contact line, the microscopic dynamic contact angle is ϭ D and at the static contact line, ϭ S .
Slip model: The liquid speed along the moving plate is v y ϭϪ1ϩexp͕Ϫ͓h͑0͒Ϫy͔/⑀͖, ͑A7͒ where h(x) is the dimensionless height of the free surface and ⑀ϭL S /L C . At the hypothetical lower surface, yϭ0, v x ϭ0,¬ ͑A8͒ v y ϭϪv y0 ͓͑1Ϫx͒͑1Ϫ3x͔͒.¬ ͑A9͒ Here v y0 is the dimensionless liquid speed at (x,y)ϭ(0,0), and is defined using the slip model ͑A7͒. Note that the imposition of Eqs. ͑A8͒ and ͑A9͒ is equivalent to using the boundary condition ‫ץ‬v y /‫ץ‬yϭ0 at that location. The mass of liquid within the computational domain must be conserved. For a liquid of constant density, this reduces to the volume contraint, where A is the liquid volume/width within the computational domain.