Detailed modelling of contact line motion in oscillatory wetting

The experimental results of Xia and Steen for the contact line dynamics of a drop placed on a vertically oscillating surface are analyzed by numerical phase field simulations. The concept of contact line mobility or friction is discussed, and an angle-dependent model is formulated. The results of numerical simulations based on this model are compared to the detailed experimental results of Xia and Steen with good general agreement. The total energy input in terms of work done by the oscillating support, and the dissipation at the contact line, are calculated from the simulated results. It is found that the contact line dissipation is almost entirely responsible for the dissipation that sets the amplitude of the response. It is argued that angle-dependent line friction may be a fruitful interpretation of the relations between contact line speed and dynamic contact angle that are often used in practical computational fluid dynamics.


INTRODUCTION
A liquid spreading over a dry surface is a phenomenon that is crucial to many natural processes and important in technology. However, the detailed description and understanding of dynamic wetting is still a complex and challenging problem 1,2 . The fact that the continuum equations of fluid mechanics exhibit a nonintegrable singularity 3 of the viscous stress at the contact line (CL) shows that the detailed microscopic and nanoscopic features of the liquid and the surface will be important for the macroscopic flow. This introduces a host of different processes and phenomena that need to be understood to predict and control wetting processes.
In technology, in addition to such examples as spray painting, coating, etc., one particularly important field is microfluidics 4 . A common challenge is to handle small volumes of liquid, often in the form of small droplets, and one means for achieving this is to use wetting phenomena. Another area where surface tension and wetting become dominant is microgravity 5 . In the absence of gravity, surface tension becomes dominant, and wetting will be important in any fluid handling, from liquid fuel to many daily activities and needs of the astronauts.
Dynamic wetting driven by vibration is both of practical importance and a convenient way to study the phenomenon. The dynamic wetting on a glass plate dipping into a tank and oscillated vertically was studied in ref. 6 , and damping of surface waves in a rectangular tank was investigated in ref. 7 , revealing complex dependencies of damping rates on oscillation amplitude. The dynamics of a sessile droplet on a vertically oscillating surface will be sensitive to the detailed conditions at the CL, such as the presence of hysteresis or CL dissipation [8][9][10] . Xia and Steen 9 made careful experiments using droplets on a polydimethylsiloxane (PDMS)-covered substrate, which was oscillated vertically at frequencies near drop resonance. The resulting dynamics was examined through phase plots of the dynamic contact angle, the CL position, and the CL speed. In particular, Xia and Steen used this information to measure the CL mobility.
On earth, the size of droplets that can be used is limited to a radius in the order of millimeters, but in microgravity, a larger parameter space can be investigated. In microgravity a droplet size in the order of centimeters can be used instead, which will be advantageous in several respects; one is that spatial dimensions are larger, and the resonance frequency much lower, allowing for higher both spatial and temporal resolution. The larger droplet size also implies that the droplet and CL dynamics are even more dominated by inertia than in a millimeter-sized droplet on earth. It was the intention of Paul Steen to perform such experiments 11 , and these have now been carried out.
From a strictly thermodynamic point of view, a moving CL should be associated with energetic losses of some kind 12 . The idea of a localized dissipation at the CL has been invoked in different ways over the years. Following Hocking 13 , Xia and Steen introduce the concept of a CL mobility M as a phenomenological parameter that relates the deviation of the dynamic contact angle from equilibrium θ − θ e to the CL speed U CL , see also refs. 8,13 . In computational fluid dynamics, more elaborate phenomenological relations between contact angle and CL speed have been devised 14 , which take the form θ ¼ f ðU CL ; θ e ; ::Þ. In Molecular Kinetic Theory (MKT) 15,16 , dynamic wetting is described as an activated process on the molecular scale, and the line friction ζ is given a phenomenological interpretation on the molecular scale. In its simplest linearized form, this can be written as where γ is the surface tension and ζ is the coefficient of wettingline friction, which in MKT is estimated in terms of molecular quantities and thermal fluctuations.
In the phase field method, the fluid is viewed as a mixture of two immiscible species. The governing equations are derived from the thermodynamic potentials of such a system to yield typically 1 Flow Centre, Department of Engineering Mechanics, The Royal Institute of Technology, 100 44 Stockholm, Sweden. 2 Södertörn University, Alfred Nobels allé 7, 141 89 Huddinge, Sweden. ✉ email: gustav.amberg@sh.se the Cahn-Hilliard equations 17,18 . The interface now becomes a diffuse region separating the two species, which has a definite width ε. The line friction appears as an energy dissipation associated with the CL displacement. Yue and Feng 18 derive the resulting equivalent condition relating CL speed and dynamic contact angle as here Γ is introduced as a rate parameter in the relaxation of the dynamic contact angle boundary condition. It is noted that the parameters in Eqs. (1)-(3) in the approximation of θ À θ e ( 1 can be identified by setting In the last equality, we have introduced the line friction μ f , which is essentially the same as the coefficient of wetting-line friction used in ref. 16 , except that it also absorbs the factor sinθ e . We note that ζ and μ f have the dimensions of viscosity and that γ/ μ f is a velocity.
In many classical treatments, notably the Cox-Voinov law 2 , it is presumed that the static contact angle applies right at the solid boundary and that the angle variations with the speed that are often observed are an "apparent" contact angle, which is attained a short distance away from the wall. It is also often assumed that there is a fluid slip on the wall at the CL, which helps regularize the singularity in stresses. In MKT, and inherent in the introduction of a CL mobility, the contact angle is assumed to be different from the static value right at the wall, when viewed on molecular length scales. In the phase field model, mass diffusion will help regularize the CL and there is no need for a fluid slip. The introduction of dissipation related to CL movement will cause the contact angle at the wall to deviate from the static value.
It is far from clear what the actual conditions are for a given liquid spreading on a particular surface. For a system of decane spreading on a surface covered with a thin layer of PDMS, it was demonstrated experimentally 19 that the microscopic dynamic contact angle is velocity dependent and deviates substantially from the equilibrium value also at very small but nonzero capillary numbers. In ref. 20 , a theoretical model is developed that links the distribution of assumed nanoscopic geometrical surface defects to the line friction dissipation. Recent molecular dynamics (MD) simulations have also shown that for water molecules and a wall with hydrogen bonds with the water molecules, the first layer of water molecules are effectively bound to the surface, a no-slip condition is appropriate, and the contact angle deviates from the equilibrium value 21,22 . It is known that the local molecular arrangements will be different in electrowetting and that this will alter the line friction 23,24 . Also, a microscopic geometrical structure on the surface will change the dynamic wetting and can be described through an effective line friction 25 . In an oil-water system, Rondepierre et al. 26 demonstrated that CL friction was responsible for a decrease in CL speed of three orders of magnitude, as a certain surfactant was added.
The actual nanoscopic cause of local dissipation at the CL can thus be very different depending on the properties of the surface, surface roughness and structure, the liquid properties, the surface chemistry of the wet surface, etc. We conclude that we should not expect any universal answer to the question of what is causing CL friction. Many different nanoscopic or microscopic processes can no doubt have this effect. However, whatever the origin, the effect can be described as a single parameter, the CL friction.
So far, the dependence of line friction on the contact angle has received limited attention, but there is clearly every reason for line friction to vary with the dynamic contact angle. The nonlinear MKT theory gives one possibility and based on MD simulations Johansson and Hess 22 formulated an expression for the angle dependence of line friction for water molecules on a surface with hydrogen bonds to the water. Other surfaces and liquids may certainly show different behavior.
On this note, we evaluate the detailed experimental results of Xia and Steen 9 , using phase field simulations, with the aim of determining precisely what is required in the mathematical model, in order to faithfully reproduce the experiments. We will study how the angle dependence of the CL friction should be chosen. The model and the simulation are then used to draw conclusions on the source of the dissipation that is evident in the experiment.

METHODS
The simulations are made using the Navier-Stokes-Cahn-Hilliard equations 18,27,28 . These describe the two-phase system as two immiscible species and are motivated by the thermodynamics of such a mixture. A phase field variable is introduced that has different values in the two species, and the fluid interface is identified as the steep but continuous transition between those values.
Here u, C, and ϕ are the fluid velocity, the phase field variable, and the associated chemical potential, respectively. C is +1 in the liquid and −1 outside the droplet. S is a reduced pressure such that p ¼ S þ Cϕ is the actual pressure. μ and ρ are the local viscosity and density, respectively, and κ is the phase field mobility. σ and ε are the phase field parameters, where ε denotes the interface thickness, and σ gives the surface tension γ through represents the standard choice in phase field methods and gives a qualitatively reasonable thermodynamic behavior representative of an immiscible mixture.
The simulations are performed in a cylindrical coordinate system that follows the oscillating substrate, giving rise to the acceleration term on the right-hand side of the momentum equation, with a t ð Þ ¼ A d sin 2πft denoting the vertical position of the substrate.
The boundary conditions on the solid wall express that the fluid cannot penetrate through the wall and does not slip on the wall.
One additional boundary condition is needed for the phase field variable, which expresses the wetting conditions.
Here σ and ε are the phase field parameters as above, giving surface tension as γ ¼ σð2 p 2Þ=3. γ 1 and γ 0 are the surface energies of the dry and wet solid surface, respectively, so that the equilibrium contact angle θ e is given by Yue 29 developed a phase field treatment of contact angle hysteresis, where advancing and receding contact angels are introduced in a piecewise continuous function on the right-hand side of the equation corresponding to Eq. (11). The CL is then allowed to move according to whether the value of the dynamic contact angle exceeds (subceeds) the advancing (receding) angle. A treatment inspired by this is also developed for level-set methods 30 .
The left-hand side represents the dissipation associated with CL motion, quantified by the parameter μ f , which we will call CL friction. This has dimensions of viscosity.
By considering solutions to the Cahn-Hilliard equations near equilibrium at the CL, an explicit relation between CL speed U CL and dynamic contact angle θ can be derived from Eq. (11), see ref. 18 and Eq. (3): The dynamic contact angle is equal to the equilibrium angle if μ f = 0, in which case there is also no dissipation at the CL. The last approximate equality holds if θ is near θ e .
The above equations are made nondimensional using R as length reference, the radius of the half-sphere (approximately the initial wet footprint radius), and an inertial capillary velocity The nondimensional parameters that appear are Ohnesorge number Nondimensional oscillation angular frequency. f is the oscillation frequency in Hertz In addition to these, the Cahn-Hilliard equations use the following parameters: Peclet number related to phase field mobility. Set to 100 in the simulations. The simulations are carried out using an adaptive finite element solver, as in ref. 31 . The adaptive grid is automatically refined as needed down to an element size of 0.001. The air viscosity and density are taken larger than those of air but much smaller than the values for the liquid. They are deemed small enough for the air to have a negligible influence on the flow inside the droplet.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Flow fields
Xia and Steen report the most detailed results for their experiment labeled M00, using a 20 µl water droplet oscillated at f = 61 Hz, with an amplitude of A d = 0.1 mm. The substrate is a slightly hydrophobic PDMS-covered silicon surface, with a static contact angle of θ e = 101°. This will be the reference case for the simulations shown here.
The nondimensional numbers for this case are Oh = 0.00256, amplitude A = 0.047, and angular velocity of the driving ω = 4.41. The low value of the Ohnesorge number indicates that the droplet dynamics are inertial and that viscosity plays a minor role. Near the solid surface, we should expect a viscous Stokes layer of thickness $ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi μ=ðρπf Þ p , which is approximately 0.07 mm for the experimental parameters, i.e., much less than the initial droplet radius of 2.1 mm.
The order unity value of the angular velocity shows that the oscillations are reasonably matched with the inertial timescale, as is expected since the experiment aims at being near resonance. The value of the line friction Ohnesorge number Oh f is set to unity, which makes the reference value of the line friction μ f ;ref ¼ ffiffiffiffiffiffiffi ffi ργR p . We will return to the more detailed modeling of the line friction below. Figure 1 shows the simulation results over one cycle for the M00 case. In Fig. 1a, the first frame shows the position near the lower turning point when the droplet is compressed towards the substrate. In the second frame, it is extending upwards and is near its most elongated state in the third frame. In the fourth frame, it is being again compressed as the substrate is moving upwards.
A Weber number based on the oscillation amplitude A d and frequency f can be expressed as A 2 ω 2 , giving the value 0:047 2 4:41 2 ¼ 0:043. We thus expect the flow to be dominated by the surface tension. Figure 1b shows the corresponding pressure fields, always reflecting the curvature of the interface. The pressure is elevated in the droplet due to the capillary pressure and fluctuates as the interface curvature varies. In the most compressed state in the first frame of Fig. 1b, there is a concave shape at the top of the droplet, which causes a local low pressure there. The pressure in the air is almost constant, due to the low air density. Figure 1c shows the velocity fields for the points where the droplet is extending upwards and being compressed (corresponding to frames two and four of Fig. 1a, b).
In Fig. 2 are shown the time signals for the vertical substrate position, together with the CL position and the droplet height. It is noticed that the amplitude of the droplet height is about three times that of the substrate amplitude, as expected for a near resonance situation. All signals quickly go into a periodic motion, with no sign of period-doubling or other more complex dynamics. The responses in both the droplet height and the CL position are nonlinear, however, with shapes departing from sinusoidal.
The CL position signal deviates from a sinusoidal shape, and tends towards a square wave, with flat peaks. This is a signature of a stickslip motion of the CL; it is rather stationary at its extreme values but shifts quickly between them twice per cycle. We will see that this is built into the simulation through angle-dependent line friction.

Model for angle-dependent line friction
In addition to the input data already discussed, angle-dependent line friction has been implemented in this simulation. The line friction in Eqs. (11) and (12) is computed from a regularized well function: and Equation (14) thus describes nondimensional line friction that has a high value μf 0 in the sticking region in the angle range θ e ± G. Amberg dθ, and a low value μf 1 outside of this range. δ smooths the corners of the well function and μ f,ref is a dimensional reference value that is used in the line friction Ohnesorge number Oh f . In the present case, Oh f = 1, and thus μ f ;ref ¼ ffiffiffiffiffiffiffi ffi ργR p . In the simulation in Figs. 1-6, these values are: θ e = 101°, dθ = 7°, δ = 1°, μf 0 = 10, μf 1 = 0.5 and Oh f = 1. We could interpret θ e + dθ = θ a as approximating the advancing contact angle and θ e − dθ = θ r as the receding angle.
Inserting the line friction according to Eq. (14) into Eq. (12), we readily obtain the dynamic contact angle as a function of CL speed (see Fig. 3). As expected, the high line friction near the equilibrium angle creates a region of near stick, with low velocities for angles in a region ≈±7°around the equilibrium. In this region, the well function H(θ) ≈ 0, and μ Ã f θ ð Þ will be constant and equal to the high value μ f0 . For CL speeds up to about 0.3, the angle is approximately constant, between 7 and 10 degrees, establishing a "slip" region. For angles departing >10 degrees from equilibrium, or speeds larger than approximately 0.4, the well function H(θ) is 1 and the line friction μ Ã f θ ð Þ becomes equal to the low value μ f1 . To   G. Amberg summarize: in the window θ À θ e j j< dθ, the line friction is constant, equal to the high value μ f0 , and for angles distinctly outside this window, θ À θ e j j> dθ, the line friction tends to the constant value μ f1 . Connecting these two regions is a velocity range where the dynamic angle is approximately constant, i.e. a region of "free slip".
As seen in Fig. 3a, there is still some variation of angle with CL speed also in the region of "free slip". The slope of the curve in Fig. 3a in this region is related to δ, the width of the transition region in line friction in Fig. 3b, and the line friction far from equilibrium μ f1 . Approximating the contact angles and the CL speeds at the ends of the transition region as θ 0 ¼ θ e þ dθ þ δ, θ 1 ¼ θ e þ dθ À δ, and U 0 ¼ ðθ 0 À θ e Þ=μ f 0 , U 1 ¼ ðθ 1 À θ e Þ=μ f 1 , and assuming that μ f 1 ( μ f 0 , an estimate of the slope in Fig. 3a in the transition or "free slip" region can be obtained as Δθ=ΔU % 2μ f 1 = 1 þ dθ δ À Á (angles in radians). For a very steep transition region (δ ( dθ) the slope will thus approach zero, while it becomes comparable to μ f1 for δ~dθ.
The model for line friction in Eqs. (13)-(15) was constructed to be as conceptually simple as possible, characterized by a high friction region around equilibrium, and low friction further from equilibrium. In order to test this against the experiments of Xia and Steen the simulated results are overlaid on their experimental phase plots (their Figs. 4 and 5), see Fig. 4. Figure 4a shows the graph of dynamic contact angle vs CL speed, the "traditional diagram" (TD) as it is referred to by Xia and Steen. The colored circles are the experiments of Xia and Steen, and the black curve is obtained from the numerical simulation. If the CL friction would be a constant this curve would always be a straight line, thus the angle variation in the CL friction becomes important. The parameters dθ, δ, μ f0 , and μ f1 for the angledependent line friction are chosen to give a fair approximation of the TD, in the following manner. dθ is chosen to capture the magnitude of the sticking region in Fig. 4a, and the sticking region line friction, μ f0 , sets the slope of the curve at the origin in Fig. 4a. The transition region width δ allows for the smooth curvature of the TD at the end of the sticking region. The line friction far from equilibrium μ f1 is then chosen to approximate the slope of the TD in the "free slip" region. It should be noted that the numerical simulation accurately reproduces the theoretical curve in Fig. 3a and that the CL velocities encountered in the experiment are always below about 0.3, so that only the "sticking" and the "free slip" regions are visited. Figure 4b shows a phase plot in terms of the CL speed vs the CL position. In addition to the experimental data shown as colored circles, the black circles denote the position of the substrate. The black line is the trajectory from the simulation. It captures both the elongation and the slight inclination of the experimental loop. The width of the simulated trajectory is slightly less than the experimental one, but overall, the agreement is good. Figure 4c shows the dynamic contact angle vs CL position. Here the stick-slip character of the motions is evident; the shape of the trajectory is a quadrilateral, where the horizontal upper and lower parts show the rapid phase when the CL moves from one almost steady position to another, at a fairly constant contact angle. The vertical sides represent the "stick" phase where the CL is nearly stationary and the angle changes. Figure 4d was introduced by Xia and Steen to highlight the "stick" and the "slip" parts of the motions and quantify those separately. The graph shows dynamic contact angle departure from the equilibrium angle multiplied by CL position r CL À r ref ð Þθ À θ e ð Þ vs CL position multiplied by CL speed r CL À r ref ð Þ U CL . The nonlinearity that is introduced creates two loops, one for the receding and one for the advancing motion of the CL. The simulated trajectory is again in fair agreement with the experiment, even though the experiment extends somewhat G. Amberg further from the origin. In view of the scatter in the experimental data the agreement in Fig. 4 was deemed sufficient for the present discussion.
In terms of the model in Eqs. (13)- (15), we see that the sticking phase is characterized by a high constant line friction μ f0 = 10. This shows up in the central part of Fig. 4d, where the simulated trajectory is steep and near the red experimental circles. The slope of the curve here is proportional to μ f0 . The main source of CL dissipation is at the slip phase, as the uncompensated Young's stress at the receding or advancing contact angle multiplies the CL speed. Here the CL slip speed becomes rather independent of the dynamic angle. In a macroscopic experiment where the overall CL speed is measured along with the contact angle, an overall CL friction for an advancing CL would be identified from Eq. (12) as μ f ¼ γ=U CL ð Þ cosθ e À cosθ a ð Þ =sinθ a . This would not be possible to directly relate to material parameters, since U CL here would be determined from the overall droplet dynamics, and the angle would stay nearly constant, close to the advancing angle. To predict the wetting behavior, the CL friction would need to be modeled and the overall dynamics simulated. For angles further away from the sticking region, μ Ã f θ ð Þ becomes μ f1 and Eq. (12) reduces to μ f 1 U CL ð Þ=γ ¼ cosθ e À cosθ ð Þ =sinθ with a constant line friction. We do not however have any experimental data to verify if this is indeed the case.
The CL mobility parameter M (Eq. (1)), as identified by Xia and Steen is the inverse of the slope of the line formed by the blue circles in Fig. 4d. This is also approximately the same as the inverse of the slope of the straight line obtained by connecting the two regions of blue circles in the "wings" of Fig. 4a.
The agreement with the experiment in Fig. 4 is fair, but there are some differences, so more precise modeling of the experimental points in the TD in Fig. 4a was made. The model in Eqs. (13)-(15) is symmetric around the equilibrium angle, while the experimental points show some differences between the advancing and receding curves. A more complex mathematical expression for the angle dependence of the line friction was designed so that the TD is approximated more closely, see Fig. 5a. Overall the agreement in all four panels is somewhat improved but not dramatically so. The conclusion is that the complexity of Eqs. (13)- (15), with the four parameters dθ, δ, μ f0 , and μ f1 is sufficient to capture the essential features of the flow.

Energy dissipation
It may be asked what the nature is of the dissipation that limits the response amplitude. In the simulation for the standard M00 case, the input energy per cycle was determined from the pressure and speed of the substrate integrated over a cycle. Figure 6b shows a graph of the total vertical force F v from the substrate, vs the substrate position over a cycle. In the same manner, the dissipation at the CL was determined from a graph of 2πr CL cos θ e À cos θ ð Þ vs r CL graph, see Fig. 6a. The total amount of work supplied to the drop by the oscillation over a cycle is now the area inside the loop in Fig. 6b, which is calculated as 0.132914. The area inside the loop in Fig. 6a giving the CL dissipation is 0.131779. As expected, the work input is slightly larger than the CL dissipation, with a relative difference of 0.8%. We should not overinterpret this small difference, in view of other possible sources of inaccuracy, but this still shows that the CL dissipation is the completely dominating cause of damping in this case. There is an almost perfect match between the input energy and the energy dissipated at the CL over a cycle. The other source of dissipation by bulk viscous dissipation is indeed expected to be quite small, given the small value of the Ohnesorge number (Oh = 0.00256).

DISCUSSION
The good qualitative and quantitative agreement between simulation and experiment, in all the different aspects that can be accessed through the phase plots presented by Xia and Steen, indicates that the essential physical processes are well represented in the mathematical model. The experiment is quite sensitive since the response is amplified through the resonance of the droplet motion, and the damping that is present in the system is controlling the resulting periodic CL motion completely. As shown in Fig. 6, the damping is almost entirely due to the CL dissipation.
Simulations were performed where the CL dissipation is removed, either by enforcing the static angle (by setting μ f = 0) or pinning the CL position. The amplitude of the droplet motion then grows quickly, and the droplet will break within a few periods. The most important part of the mathematical model is the angle-dependent line friction in Eqs. (13)- (15). This simple function allows for a narrow high friction region near equilibrium, and much lower friction outside this region, and this is sufficient for capturing the essentials of the flow dynamics.
The simulations presented here only show results for an inertial flow with a low Ohnesorge number. It may be asked if the behavior of a viscous drop, with Oh ≥ 1 could be described in a similar way, with hysteresis modeled as angle-dependent line friction. While this has not been carried out here, I expect this to be the case since the line friction is local to the CL and not directly dependent on the internal bulk flow. It should be noted though that the experimental arrangement of Xia and Steen would probably not be suitable for very viscous drops, since the resonance is needed to have a large enough amplitude of the CL motion.
This also exemplifies how empirical models that postulate a relation between CL speed and dynamic contact angle can be interpreted through Eq. (12) as revealing the angle dependence of the line friction. As commented above, the line friction can have many different micro-or nanoscopic causes, and thus may depend on the dynamic angle in different ways. But in attributing the CL variation to the line friction and hypothesizing it to be reflecting processes that are local to the CL, we could hope to decouple it from the macroscopic flow problem. It would then be possible to restrict the problem to analyzing the micro-or nanospecifics of a particular surface and a particular fluid to find the line friction and its angle dependence 22 and then hope that this line friction could be used with accuracy for all flow situations, for that particular combination of fluid and surface.

DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. Fig. 6 Calculation of CL dissipation. a Uncompensated Young's stress 2πr CL cos θ e À cos θ ð Þ vs r CL over one cycle. The area inside the loop is 0.131779. b Vertical force from the substrate on the drop vs substrate displacement a t ð Þ ¼ A sin ωt. The area inside the loop is 0.132914.