Role of the rate of surface dilatation in determining microscopic dynamic contact angle

The phenomena of wetting play an important role in many industrial and natural processes, such as oil recovery,^1,2 lubrication,^3,4 electrowetting,⁵ microfluidics,^6–9 biosurfaces,¹⁰ and emulsifications.^11,12 The extent to which a fluid wets a solid surface is commonly characterized by the contact angle (angle formed by the fluid-fluid interface with the solid). For the case of a static droplet, Young¹³ determined the contact angle by balancing the forces acting at the contact line as

where θ_s is the static contact angle and the subscripts AB, AW, and BW denote the interface between fluid A-fluid B, fluid A-Wall, and fluid B-Wall, respectively. This relation has also been derived by minimizing the total free energy of the system,¹⁴ validated using statistical mechanics,^15,16 and recently shown to be valid in the microscopic scale¹⁷ as well.

In comparison, for the dynamic case, there are several competing dynamic contact angle models in the literature that are based on a diverse set of physical mechanisms with applicability to length scales ranging from macroscale to molecular scale. The three common theories that dynamic contact angle models are based on are the hydrodynamic theory, the molecular kinetic theory (MKT), and the interface formation theory (IFT). Most dynamic contact angle models fall under the hydrodynamic theory,^18–22 which is based on classical continuum fluid mechanics, where the singularity of the moving contact line is alleviated by relaxing the no-slip boundary condition^22–25 or by artificially truncating the solution at a molecular length scale where the solution breaks down.¹⁹ Usually in these models, the contact angle is observed at the macroscale, with the viscous bending of the interface occurring at the mesoscale, while at the microscale it is assumed that the contact angle is the same as its static equilibrium value. However, in contrast, recent findings^26,27 have shown the dependence of the microscopic contact angle on the contact line velocity.

In the case of models based on the molecular kinetic theory,²⁸ the dynamic contact angle is attributed to the local frictional force acting at the contact line. Here, the motion of the contact line is based on the adsorption-desorption of the fluid molecules on to the surface within the three-phase zone. In this case, there is no viscous bending and the observed angle is the microscopic dynamic contact angle, which is velocity dependent.

Finally, in the interface formation theory proposed by Shikhmurzaev,²⁹ the contact line advances by a rolling type motion. This results in the transport of fluid particles from the fluid-fluid interface to the fluid-wall interface, through the contact line and vice versa. Since the properties of the two interfaces are not the same, there is a transition region in the vicinity of the contact line, where a surface tension gradient exists. The gradient results in a local surface tension value different from that of the equilibrium value. This is described by a constitutive equation proposed by Shikhmurzaev, in which surface tension is related to the surface (interface) density. Assuming that Young’s equation is still valid, IFT suggests that the gradient in surface tension results in a change in dynamic contact angle.

Here, we listed only some of the dynamic contact angle models and the theories on which they are based; for a more detailed list, the reader is referred to the review articles by Snoeijer and Andreotti,³⁰ Chau,³¹ and Blake.^32,33 All these dynamic contact angle models, some of which are based on competing mechanisms, show reasonable agreement with experiments.^34–36 This makes it difficult to ascertain the underlying mechanism governing the dynamic contact angle, especially the microscopic dynamic contact angle. In this manuscript, we aim to resolve this issue by using interfacial hydrodynamics to analytically derive a model for the microscopic dynamic contact angle and validate it using molecular dynamics (MD) simulations.

To this end, we find the functionality of the microscopic dynamic contact angle by writing the force balance at the contact line by considering Gibbs’ representation of the interface as a dividing surface and the contact line as the intersection of these dividing surfaces. Gibbs’ representation has the benefit of maintaining the modeling simplicity of continuum mechanics seen in hydrodynamic contact angle models, while obtaining the microscopic resolution of the MKT and IFT models, discussed above. The force balance is given by the balance of surface stress of the respective interfaces, which is usually assumed to be equal to the thermodynamic surface tension. Here, we shall show that this is true only in the static case or far from the moving contact line, where slip velocity is a constant. In the dynamic case, the surface stress at the moving contact line is a function of both the rate of surface dilatation and the thermodynamic surface tension. Hence, from the force balance, it follows that the microscopic dynamic contact angle is dependent on the one-dimensional rate of surface dilatation [∂u^(s)/∂s] in addition to the thermodynamic surface tension, namely,

Here, u^(s) is the tangential component of surface fluid velocity, ∂s is in the direction tangent to the interface, θ_d is the microscopic dynamic contact angle, μ_s is the surface shear viscosity, and κ_s is the surface dilatation viscosity. In addition, we shall show that one of the key factors responsible for the difference in contact angle at the advancing and receding edges is the rate of surface dilatation. We argue that the rate of surface dilatation results from varying slip velocity in the vicinity of the moving contact line²³ and is inseparable from it. Stated differently, we cannot have varying slip velocity without having a rate of surface dilatation and vice versa.

We also demonstrate that mechanical surface tension is related to thermodynamic surface tension by the rate of surface dilatation as $γ¯=γ+κs∇(s)⋅v(s)$⁠, which is obtained using the constitutive equation, commonly referred to as the Boussinesq surface fluid model. This is similar to the relation between thermodynamic and mechanical pressure in fluid dynamics.^37,38 The above relation shows that in the vicinity of the moving contact line, where the rate of surface dilations is not negligible,^23,39–41 the mechanical and thermodynamic surface tension cannot be assumed to be equal, as is almost always the case.⁴² Based on a phenomenological explanation, IFT suggested the existence of a thermodynamic surface tension gradient, but here we can show analytically that it is only the mechanical surface tension that has a gradient, while the thermodynamic surface tension still remains constant, despite having assumed no surface mass flux at the moving contact line.

Hence, in our work, we discover that the rate of surface dilatation is the fundamental mechanism governing the microscopic dynamic contact angle besides the thermodynamic surface tension. We reiterate that only the mechanical surface tension has a gradient in the vicinity of the contact line as a result of this surface dilatation, while the thermodynamic surface tension remains constant. Also, it is found that the difference in the sign of surface dilatation (expansion and compression) in the vicinity of the advancing and receding contact lines is one of the primary causes for the difference in advancing and receding contact angles in the case of forced wetting.

This manuscript is organized as follows: In Sec. II, an expression for the microscopic dynamic contact angle model is derived from the force balance at the moving contact line. Using Boussinesq’s surface fluid model and the mass conservation at the interface, the model is presented in three different forms, namely, surface dilatation form, mechanical form, and surface mass flux form. The numerical setup with details of molecular dynamics simulations and methods of computation of respective variables are presented in Sec. III. Finally, in the results Sec. IV, various findings are validated and discussed in greater detail.

Our objective here is to obtain the microscopic dynamic contact angle by finding the correct matching condition for the force at the moving contact line. First, we briefly outline the steps undertaken to achieve this task. Following Gibbs’ interpretation, the conservation of momentum is written for a control volume encompassing the bulk media, the interface, and the contact line. This is presented for a general multiphase system, which we then simplify for the specific case of a contact line formed by two immiscible fluids, steadily moving over a planar surface. The force balance thus obtained at the moving contact line is expressed in terms of surface stress, which in itself is not in a closed form and hence of little use. Therefore, using a constitutive relation presented by Boussinesq for a surface fluid and considering the surface mass balance at the interface, we find the microscopic dynamic contact angle model. This model is presented in three different forms, each having its own advantages. Throughout this manuscript, we refer to the mass and force balance in the bulk media as bulk mass balance and bulk force balance, respectively. Similarly, we refer to the balances in the surface as surface mass and force balance, while the balances in the contact line are referred to as line mass balance and line force balance.

With the brief outline discussed, we now present the force balance for a general multiphase system,⁴³ by considering a control volume encapsulating the bulk media, the interfaces, and the contact line as illustrated in Fig. 1. In the force balance, the terms are grouped into volume, area, and line integrals, corresponding to the bulk media, the interface, and the contact line, respectively; see the following equation:⁴⁴ 

Here, ρ is the density, v is the velocity, T is the stress tensor, and b is the body force. Terms without any subscript or superscript are associated with bulk media, (·)_(s) and (·)^(s) are terms associated with an interface, while (·)_(l) and (·)^(l) are associated with the contact line. R is the region encompassed by the control volume, Σ is the area formed by the interface, and c is the contact line (normal to the plane of Fig. 1). The volume integral is the same as the commonly used momentum equation for a bulk medium. The surface and line integrals are analogous to the bulk momentum equation, except for the additional jump terms denoted by $〚⋅〛$⁠. The first of the jump terms accounts for the momentum associated with mass transfer across the interface, $〚ρ(v−v(s))(v−v(s))⋅n^s〛$⁠, and the contact line, $〚ρ(s)(v(s)−v(l))(v(s)−v(l))⋅n^l〛$⁠.⁴⁵ The second jump term accounts for the jump in bulk stress across the interface, $〚T⋅n^s〛$⁠, and the jump in surface stress across the contact line, $〚T(s)⋅n^l〛$⁠.⁴⁶ The importance of the jump terms lies in the fact that it is through these jump terms that a higher dimension is connected and communicates with the lower dimension and vice versa. Hence, from Eq. (2), we obtain the force balance at the moving contact line in its most general form

It must be noted that although Eq. (2) is discussed in the context of a multiphase flow, the force balance itself is applicable to any problem with a surface and line discontinuity. Xia and Mohseni⁴⁷ had used a similar force balance to study the formation of vortex sheets at the trailing edge of an airfoil. They applied conservation laws to a “Y”-shaped control volume encompassing the two boundary layers at the trailing edge that merge to create the free shear layer. DeVoria and Mohseni⁴⁸ also used similar equations to define an inviscid model for boundary and shear layers, where the mass and momentum in the layers were confined to a two dimensional surface.

We now simplify this general force balance [Eq. (3)] for a multiphase system for the specific case of a steadily moving contact line by employing the following assumptions:

1. The contact line is in dynamic equilibrium, d_(l)v^(l)/dt = 0.

2. The contact line is subject to no body force, b^(l) = 0.

3. The contact line is a material boundary with no mass transfer across it,

1. There is no gradient in the line stress along the contact line (normal to the plane of Fig. 1), ∇_(l) · T^(l) = 0.

Using these assumptions, the force balance [Eq. (3)] reduces to

The force balance in its current form is not closed, and therefore, a surface constitutive law presented by Boussinesq is used to rewrite Eq. (4) in terms of surface tension (surface pressure) and surface velocity. Thereafter, we present the above force balance in two additional forms, each having its own advantages.

The Boussinesq surface fluid model^49,50 describes the surface stress as T^(s) = [γ + (κ_s − μ_s)∇_(s) · v^(s)]P + 2μ_sD^(s). Here, $P=I−n^n^$ is the surface projection tensor and $D(s)=1/2(P⋅∇(s)v(s)+(∇(s)v(s))T⋅P)$ is the rate of surface deformation tensor. Since we assume a two-dimensional problem with the domain periodic in the binormal direction, there is no surface shear and the surface stress simplifies to T^(s) = γP + (κ_s + μ_s)(du^(s)/ds)P. Hence, substituting the surface stress, the force balance at the moving contact line is written as

By expanding and projecting the force in the direction tangent to the wall, we obtain

Here, γ, κ_s, and μ_s are surface material properties that are known a priori. The rate of surface dilatation [∂u^(s)/∂s] can be computed from the flow, leaving the microscopic dynamic contact angle as the only unknown. In the field of surface rheology, a substantial effort has been devoted to computing the surface shear and dilatation viscosity,^51–54 though this has been primarily focused on a fluid-fluid interface rather than a fluid-wall interface. Hence, we demonstrate that the microscopic dynamic contact angle, θ_d, is dependent on the surface dilatation in addition to the thermodynamic surface tension. In order to eliminate one of the material properties, which is required as an input, namely κ_s, we rewrite the force balance in Eq. (6) by substituting the thermodynamic surface tension at the respective interfaces with the mechanical surface tension as is shown next.

From the constitutive equation of a bulk fluid stress, it is known that the mechanical pressure (⁠$p¯$⁠) is related to the thermodynamic pressure (p) as $p¯=p+κ∇⋅v$⁠, where κ is the bulk dilatational viscosity. Considering that the Boussinesq surface fluid model [T^(s) = [γ + (κ_s − μ_s)∇_(s) · v^(s)]P + 2μ_sD^(s)] is analogous to the constitutive equation for the bulk fluid stress, it follows that the mechanical surface tension is related to the thermodynamic surface tension in the same way as the two forms of pressure,

Here, the mechanical surface tension is defined as the average of the diagonal terms of the surface stress tensor, in the same manner, the mechanical pressure is defined as the average of the diagonal terms of the bulk stress tensor. This definition for mechanical surface tension can be readily shown to be equivalent to the classical definition, $γ¯=∫−∞∞(PT−PN)dn$⁠, where P_T and P_N are components of the bulk pressure, tangential and normal to the interface. It must be noted that surface tension could be interpreted as the surface pressure. By rewriting the force balance at the contact line [Eq. (6)] in terms of the mechanical surface tension, we show that the microscopic dynamic contact angle model is

In the limit of a static case, the surface flow is surface divergence free and $T(s)=γP=γ¯P$⁠; hence, the above equation reduces to Young’s equation for the static contact angle γ_AB cos θ_s = γ_AW − γ_BW. Next, we show that, using the surface mass balance at the interface, the surface dilatation term can be rewritten in terms of the mass flux of bulk fluid into the interface, thereby eliminating the need of computing the rate of surface dilatation.

One of the challenges in the surface dilatation form of the force balance is computing the rate of surface dilatation (∂u^(s)/∂s) in the macroscale. Hence, it is replaced using the surface mass balance at the interface. The surface mass balance for a steady incompressible surface flow is given as⁴³ 

where $v$ is the component of bulk velocity normal to the interface. Substituting in the dynamic contact angle model, Eq. (6) is rewritten as

where $λs=(κs+μs)1ρ(s)$ is called the total surface kinematic viscosity. Therefore, in this form, the only surface quantities required to compute the microscopic dynamic contact angle model are surface material properties, which are known a priori.

This model is further simplified by using Young’s relation for static contact angle, θ_s, to eliminate thermodynamic surface tension at the wall-fluid interfaces to obtain

Here, we assume that the mass flux into the fluid-wall interface is only from the bulk fluid side, as the wall side has no relative velocity normal to the interface. Hence, (ρ$v$)_A and (ρ$v$)_B correspond to the mass flux into the interface from bulk fluid A and B, respectively. In addition, we have also assumed the fluids to be immiscible; therefore, there is no mass flux across the fluid-fluid interface. Hence, we infer that the change in contact angle between the static and the dynamic case is a result of the difference in mass flux at the two fluid-wall interfaces in the vicinity of the moving contact line.

The findings are validated and discussed in Sec. IV. First, however, the details of numerical setup and methods used for the evaluation of surface quantities are presented in Sec. III.

Here, we describe the problem geometry, provide details of the numerical simulation, and describe the methods used to evaluate various surface quantities from MD simulations.

The moving contact line is simulated by modeling a two-phase Couette flow for immiscible fluids, where the walls move in the opposite directions with a constant speed U (see Fig. 2). Here, periodic boundary conditions are imposed along the x and z directions. The contact angle (θ) is defined as the angle formed by the wall and fluid-fluid interface, as measured in fluid A. Due to the symmetric nature of the problem, θ₁ = θ₄ and θ₂ = θ₃ are the leading and trailing edge, respectively, in the frame of reference of the wall.

As previously eluded to, one of the major difficulties faced in studying the dynamic contact angle is that direct measurement is limited by the ability of an experimental technique to resolve the interface. With the aid of molecular dynamics simulations,⁵⁵ we are able to overcome this hurdle. The pairwise interaction of molecules, separated by a distance r, is modeled by the Lennard-Jones (LJ) potential

Here, ϵ and σ are the characteristic energy and length scales, respectively. The potential is zero for r > r_c, where r_c is the cutoff radius, which we set to r_c = 2.5σ, unless otherwise specified.

Each wall is comprised of at least two layers of molecules oriented along the (111) plane of a face centered cubic (fcc) lattice, with the molecules fixed to their respective lattice sites. The fluid molecules are initialized on a fcc lattice, whose spacing is chosen to obtain the desired density, with initial velocities randomly assigned so as to obtain the required temperature. The fluid in its equilibrium state has a temperature T ≈ 1.1k_B/ϵ and number density ρ ≈ 0.81σ⁻³. The temperature is maintained using a Langevin thermostat with a damping coefficient of Γ = 0.1τ⁻¹, where $τ=mσ2/ϵ$ is the characteristic time and m is the mass of the fluid molecule. The damping term is only applied to the z direction to avoid biasing the flow. The equation of motion of a fluid atom of mass m along the z component is therefore given as follows:

Here, ∑_j≠i denotes the sum over all interactions and η_i is a Gaussian distributed random force. The value of dynamic viscosity of the fluid is μ ≈ 3.3ϵτσ⁻³ and the Reynolds number is Re ≈ 0.95, which corresponds to the case with maximum wall velocity (U = ±0.075στ⁻¹). For all the cases considered, the numbers of fluid A and fluid B atoms are 225 552 each and the numbers of top and bottom wall atoms are 354 56 each.

The LJ coefficients and relative number density (scaled by the number density of the wall, ρ_$w$ = 0.81σ⁻³) of the various cases simulated are listed in Table I. The immiscibility of the two fluids is modeled by choosing appropriate LJ interaction parameters, such that the interatomic forces between them is predominantly repulsive. This is done by choosing a small value for epsilon (which represents the depth of the potential well). For the results presented here, the LJ parameters for interactions between fluid A and fluid B are ϵ^ff = 0.2ϵ, σ^ff = 3.0σ. Since σ^ff = 3.0σ, a cut-off radius of r_c = 5.0σ is used for the fluid A-fluid B interactions. For simplicity, the two fluids are assigned identical fluid properties and only the properties of the fluid-wall interactions are changed. The fluid channel measures 153.0σ × 27.4σ × 144.0σ. The channel height is determined such that the individual moving contact lines are isolated from each other. This is ensured by choosing a channel height such that the fluid-fluid surface tension reaches its equilibrium value along the interface, before it is affected by the MCL on the opposite wall.

The equations of motion were integrated using the Verlet algorithm^56,57 with a time step Δt = 0.002τ. The equilibration time is determined by the time it takes for the fluid-fluid surface tension to reach a steady value, which is 60 000τ (supplementary material). The simulation is initially run until the flow equilibrates, after which spatial averaging is performed by dividing the fluid domain into rectangular bins of size ∼0.52σ × 0.54σ along the x–y plane, and extending through the entire width of the channel. In addition to spatial averaging, time averaging is done for a duration of 20 000τ for the moving contact line problem.

Studies by Priezjev⁵⁸ and Pahlavan and Freund⁵⁹ have shown that the stiffness of thermal walls affects the slip length and its dependence on the shear rate. It is also well known that the property of the secondary fluid in a two-phase flow governs the contact angle, which in turn could affect the local stresses in the vicinity of the contact line in the primary fluid. In this paper, for simplicity and in order to isolate these effects, the test cases are modeled with wall molecules fixed to the lattice site.

In order to validate the contact angle model, we need to evaluate the respective surface quantities and determine the microscopic contact angle from MD simulations. We start by describing how the microscopic contact angle is evaluated, for which the shape of the interface needs to be defined first.

The distribution of number density of fluids A and B, along the length of the channel (x direction), is presented in Fig. 3(a). It is observed that the number densities of both fluids A and B drop sharply and the profiles intersect within the interface region. Hence, the fluid-fluid interface is defined as the loci of these intersections (ρ_A = ρ_B) along the height of the channel, see Fig. 3(b). By fitting a cubic polynomial function through these discrete points, a continuous interface curve is obtained. Now, that the function describing the shape of the interface is determined, the contact angle is evaluated by computing the slope of this function relative to the wall at the edge of the contact line region (CLR). Next, we define the interface and contact line region.

The width of the interface region is determined by locating its two boundaries. The boundary of an interface region is defined as the location where the local number density deviates by 1% of its bulk value. This is defined based on Gibb’s definition of an interface, where he defines an interface as a hypothetical dividing surface separating two homogeneous media. Hence, the interface region and in turn the interface boundary are defined as where the fluid properties deviate from its bulk value by 1%. In the case of the wall-fluid interface region, the location of the boundary on the fluid side is given by the height at which density layering in the fluid begins, depicted in schematic shown in Fig. 4. On the other hand, the boundary on the wall side is defined at a distance of 0.25σ away from the wall lattice site. Similarly, the two boundaries of the fluid-fluid interface region are determined by the location where the density decreases by 1% of its bulk value. Now, the contact line region (or three phase zone) can be defined as the intersection of these three interface regions. The surface tension forces and the force due to surface dilatation are evaluated at the edge of this contact line region.

Before we describe how the individual surface quantities are evaluated we refer back to Gibbs’ definition of a dividing surface (interface). Gibbs describes the interface as a hypothetical dividing surface that separates two homogeneous phases. In the homogeneous phase, all variables, such as density, viscosity, and stress, have uniform values and the constitutive equations apply uniformly. Any excess properties or fluxes not accounted by the homogeneous phase are assigned to the dividing surface. This analogy can be directly extended to the contact line, formed by the intersection of two or more dividing surfaces. Here, any excess surface quantities are assigned to the common line.

In order to compute the surface tension, the components of stress tensor need to be evaluated first. The stress tensor is computed by accounting for the forces acting on the fluid element.^60–62 The fluid element corresponds to the rectangular bins formed by dividing the fluid domain. By spatially averaging the stress, we obtain the average stress tensor in each bin. In Fig. 5, we plot the stress components (σ_xx and σ_yy) across the fluid-fluid and at the fluid-wall interfaces, where the anisotropy of stress around the interface is observed. Knowing the components of the stress tensor, we now calculate the surface tension as given by its mechanical definition, $γ=∫L1L2[PN−PT]dl$⁠.⁶³ Here, P_N and P_T are the normal and tangential components of the local stress tensor with respect to the interface and dl is in a direction normal to the interface. The limits of integration are chosen such that at l = L₁ and l = L₂ the stress is isotropic (P_N = P_T). In order to evaluate the surface tension for the fluid-fluid interface, the limits L₁ = −5σ and L₂ = 5σ are used. Since we use inert walls,⁶⁴ the wall-fluid interface surface tension can be computed using the same mechanical definition as given above. The limits of integration for the wall-fluid interface extend from the wall to the middle of the channel. Next, we detail the steps involved in evaluating the surface density and surface velocity.

Similar to surface tension, the surface density [ρ^(s)] is evaluated by integrating the bulk density of the bulk fluid across the width of the interface region. Surface velocity [$v$^(s)] on the other hand is evaluated by first computing the surface momentum by integrating the bulk momentum across the width of the interface region and then dividing it by the surface density. Finally, we present the steps used to compute the surface viscosities.

The surface dilatational viscosity is evaluated by using the relation between thermodynamic and mechanical surface tension presented in Eq. (7), $γ¯=γ+κs∇(s)⋅v(s)$⁠. The thermodynamic surface tension (γ) is found by evaluating the mechanical surface tension (⁠$γ¯$⁠) far away from the contact line, where ∇_(s) · $v$^(s) = 0 and $γ=γ¯$⁠. Then, the surface dilatational viscosity (κ_s) is found by fitting the data of $γ¯−γ$ to κ_s∇_(s) · $v$^(s), along the vicinity of the receding contact line. The surface shear viscosity is computed by using the momentum balance at the interface, Eq. (2). From looking at the MD results, it is found that inertia and the jump in momentum due to mass flux terms are negligible in comparison to the surface stress and the jump in bulk stress terms. Hence, the surface momentum balance reduces to $∇(s)⋅[γ¯−μs∇(s)⋅v(s)]P+2μsD(s)+〚T⋅n^s〛=0$⁠. Using this relation, we can compute the surface shear viscosity (μ_s) as it is the only unknown. Here, the surface dilatational and shear viscosities are assumed to be a constant for a given interface. They are computed for the top wall in the vicinity of the contact line and for a wall velocity of u = 0.05στ⁻¹. The list of surface dilatation and surface shear viscosities corresponding to various wall-fluid interfaces is tabulated in Table II. The surface viscosities are not computed for the fluid-fluid interface as the rate of surface dilatation is found to be negligible as a result of the fluids being immiscible.

So far in this manuscript, we have presented the force balance at the moving contact line and in turn presented a model for microscopic dynamic contact angle. This is presented in three different forms. In the first form, it is shown that the microscopic dynamic contact angle is a function of the rate of surface dilatation and the thermodynamic surface tension, using the Boussinesq surface fluid model. Next, in the mechanical form, the thermodynamic surface tension and the rate of surface dilatation are replaced by the mechanical surface tension using the Boussinesq surface fluid model. Finally, in the surface mass flux form, surface mass balance is used to replace the rate of surface dilatation with mass flux of bulk fluid into the interface. In this section, we validate these findings and the microscopic dynamic contact angle by performing an MD simulation of multiphase Couette flow problem with different velocities and different wetting properties. Based on the results from MD, these findings are discussed in greater detail.

We start by validating the surface mass balance, presented in Eq. (9). From the surface mass conservation, it is found that the rate of surface dilatation is balanced by the mass flux of the bulk fluid into (and out of) the interface. For a two-dimensional problem, the rate of surface dilatation term becomes ∂u^(s)/∂s and the equation for mass conservation reduces to

Here, s is in the direction tangent to the interface. Unlike the bulk fluid, where the bulk mass balance under the assumption of incompressibility (constant density) guarantees a divergence free flow or no bulk dilatation, in the interface, the surface mass balance under the assumption of incompressibility (constant surface density) does not ensure zero surface dilatation. In fact, the rate of surface dilatation is rather predicated on whether there is a bulk mass flux into (or out of) the interface or not. In Fig. 6, we present the distribution of bulk mass flux and the rate of surface dilatation along the wall, in the vicinity of the advancing and receding contact lines. The data presented are for interfacial properties corresponding to case 3 and for three different cases of wall velocities (u = 0.025στ⁻¹, 0.050στ⁻¹, and 0.075στ⁻¹). The results show good agreement between the terms corresponding to the convection of surface fluid and the mass flux of bulk fluid into the interface. It is observed that, near the receding contact line, the deceleration of the surface fluid results in a mass flux out of the wall-fluid interface and into the bulk fluid. In contrast, the acceleration of the surface fluid, in the vicinity of the advancing contact line, leads to the bulk fluid flowing into the interface.

The rate of surface dilatation can be interpreted as a spatial acceleration or linear strain rate, which is nonzero only in the vicinity of a moving contact line. This is because the streamline on approaching the contact line is forced to turn and hence decelerate along the direction parallel to the interface. The presence of linear strain rate is in line with previous numerical^22,23,39,40 and experimental⁴¹ findings, where a sharp reduction in the magnitude of tangential component of fluid velocity, in the vicinity of the moving contact line, is observed. Here, the velocity magnitude reduces from a value corresponding to no-slip (or finite slip), far away from the contact line, to that corresponding to perfect slip, at the moving contact line. This variation in slip is inseparable from the rate of surface dilatation, in that we cannot have varying slip velocity without having a rate of surface dilatation and vice versa. The importance of this linear strain rate (or spatial acceleration) in accurately defining the slip boundary condition for a moving contact line problem has been discussed by Thalakkottor and Mohseni²³ and Qian et al.⁴⁰ Although surface dilatation is not often discussed in the context of wetting, this has been well established in the context of foam, emulsions, and thin films.⁵¹ 

As previously discussed, in the vicinity of the contact line, the surface flow is not divergence free; hence, the mechanical surface tension does not equal the thermodynamic surface tension in this region. This is different from the usual assumption of them being equal.⁴² As presented in Eq. (7), the mechanical surface tension equals the sum of thermodynamic surface tension and surface dilatation. In Fig. 7, the distribution of mechanical surface tension along the wall is compared to the sum of thermodynamic surface tension and the rate of surface dilatation. It is compared for interfacial properties corresponding to case 3 and for three different cases of wall velocities (u = 0.025στ⁻¹, 0.050στ⁻¹, and 0.075στ⁻¹). The results show good agreement, thereby validating the relation presented in Eq. (7). Here, the thermodynamic surface tension is evaluated by computing the mechanical surface tension far away from the contact line, where ∇_(s) · v^(s) = 0 and $γ=γ¯$⁠. The surface dilatation viscosity (κ_s) is tabulated in Table II. In addition, by comparing the distribution of surface tension in the vicinity of the contact line, for a static and a dynamic case (Fig. 8), it is clearly seen that as we approach the moving contact line, the magnitude of mechanical surface tension for the dynamic cases begins to deviate more from that of the static case and hence the existence of mechanical surface tension gradient is demonstrated. From Eq. (7), it can be concluded that the cause of this deviation between the mechanical and thermodynamic surface tension is a result of the increased rate of surface dilatation observed in the vicinity of the moving contact line.

The microscopic contact angle model is validated by comparing the contact angles from several different cases of wall-fluid properties and for the varying wall velocities listed in Table III. The local mechanical surface tension values and the rate of surface dilatation are directly computed at the contact line (edge of the contact line region) and substituted in Eq. (8) (the mechanical form) to obtain the contact angle. The surface shear viscosity is tabulated in Table II. The prediction made by the microscopic dynamic contact angle model agrees well with the results of the MD simulations. For most cases, the error is under 5%. For the limiting case of a stationary wall (U = 0), the angle θ₁ ≈ θ₂ as dictated by Young’s equation for static contact angle. This shows the consistency of both the MD simulation setup and the methodology used to evaluate the local surface tension force for the respective interfaces. Here, even though the results are only presented for fluid A, the model is also validated for fluid B as they share the same interface. The contact angle for fluid B is 180 − θ, where θ is the contact angle for fluid A. The local mechanical surface tension values at the contact line for respective interfaces are tabulated for various test cases in Table IV.

Our findings appear to be contrary to the recent claim that Young’s equation can describe the microscopic dynamic contact angle,^26,66 particularly the microscopic dynamic contact angle computed at the edge of the contact line region at the molecular scales. The reason for this apparent discrepancy is because in the work by Fernandez et al. for a two fluid system, they computed the cumulative tangential force exerted by the wall on the fluid within the three phase zone (contact line region), instead of decomposing the forces into respective force contributions from thermodynamic surface tension and the rate of surface dilatation.

The gradient of mechanical surface tension, apart from making an important contribution in determining the microscopic dynamic contact angle, also provides an explanation for the difference in contact angle at the advancing and receding edges of a steadily translating droplet. This difference in the contact angles is verified in the results presented in Table III, which also include the limiting case of zero wall velocity so that θ₁ ≈ θ₂, which is in agreement with Young’s equation. This is also consistent with the literature as the wall modeled here is homogeneous and has no macroscopic roughness.^67–69 Physically, the difference in contact angles is explained by referring back to Fig. 8(a) and Table IV, from which it can be seen that the mechanical surface tension along the wall-fluid interface decreases at the receding edge, while it increases at the advancing edge, resulting in a momentum imbalance. Smaller values of the mechanical surface tension at the receding edge indicate an increased local wettability and hence the fluid-fluid interface tends to form a smaller contact angle. In the same way, at the advancing edge, there is a larger value of the mechanical surface tension, which results in a larger contact angle. As previously discussed, this is directly attributed to the spatial acceleration of the flow along the interface of the wall towards the moving contact line.

In this manuscript, following Gibbs’ interpretation of an interface, we found the force balance at the moving contact line and in turn presented a model for the microscopic dynamic contact angle. The model states that in the dynamic case the microscopic contact angle is dependent on the rate of surface dilatation in addition to thermodynamic surface tension of respective interfaces. We also found that as a result of this rate of surface dilatation there is a disparity between the mechanical and thermodynamic surface tension, and they cannot be assumed to be equal in the vicinity of the moving contact line. In addition to finding and interpreting the correct force balance at the contact line, we demonstrated that in the case of forced wetting one of the key factors responsible for the difference in contact angle at the advancing and receding edges is the rate of surface dilatation. We showed that in the vicinity of the moving contact line the rate of surface dilatation is a result of varying slip velocity and as such inseparable from it, that is, we cannot have varying slip velocity without having a rate of surface dilatation and vice versa. All of these findings were validated using molecular dynamics simulations.

See the supplementary material for the time evolution of surface tension of the fluid-fluid interface.

This research was partially supported by the Office of Naval Research and the National Science Foundation. The authors thank Dr. A. C. DeVoria and Dr. P. Zhang for their help.