Non-axisymmetric elastohydrodynamic solid-liquid-solid dewetting: Experiments and numerical modelling

We have studied the dewetting dynamics of partially wetting liquid films confined between a soft elastic hemisphere and an elastomer layer by means of systematic experiments. We focused on the experimentally most relevant case of non-axisymmetric dewetting, which initiated at the locations of minimum film thickness near the perimeter of the contact area. We found the contact line speed to be highly anisotropic in this case. It is significantly faster in the azimuthal direction along the perimeter of the contact spot than in the radially inwards direction. We developed a three-dimensional, fully coupled numerical model that reproduces many features observed in the experiments.


Introduction
The presence and stability of thin liquid films between two solid surfaces is of great scientific and technological relevance for many systems and applications. For minimizing wear in bearings, a stable liquid film is desirable [1,2]. To prevent aquaplaning on wet roads, the liquid film should be removed as efficiently as possible to maximize traction [3][4][5][6][7][8][9][10][11][12][13]. Certain carnivorous plants developed a prey capture strategy based on insect aquaplaning [11,14]. Similarly, for residue-free nanoimprint lithography, a swift and complete removal of intervening liquid is required [15][16][17][18][19][20]. The strength of adhesion between organic or polymeric surfaces is considerably reduced by the presence of a water film, as the Hamaker constant for the interaction of e.g. two polystryrene or poly-tetra-fluoro-ethylene halfspaces across vacuum is 7 to 10 times higher than across an ultrathin water layer [21].
Roberts and Tabor reported the spontaneous collapse of thin liquid films confined between a rubber surface and an SiO 2 lens for film thicknesses below 40 nm [22]. Serayssol and Davis studied the influence of surface interactions on the elastohydrodynamic collision of two suspended spheres [23]. Becker and Mugele and de Beer et al. studied the dewetting of molecularly thin liquid films between mica substrates using a surface forces apparatus [24,25]. They observed dewetting transitions and Supplementary material in the form of a .pdf file available from the Journal web page at https://doi.org/10.1140/epje/i2020-11926-3 a e-mail: a.a.darhuber@tue.nl trapped liquid for high approach rates of the two substrates. Persson and Mugele provided an exhaustive review of squeeze-out of molecularly thin wetting and nonwetting liquid films [26]. Bandyopadhyay et al. showed that a micrometer-scale corrugation of the elastomeric punch can effectively suppress the dewetting transition of confined ultrathin polystyrene films [27]. Li et al. found that the addition of the surfactant sodium-dodecyl-sulfate (SDS) can prevent dewetting and induce film stability [28]. Brochard-Wyart and de Gennes were the first to present an analytical model for the dewetting of water films between a rigid solid and a rubber [29]. By balancing the surface energy change with the viscous dissipation in the liquid they derived the scaling law for the growth rate of an axisymmetric dry spot. Martin and Brochard-Wyart presented experiments that confirmed the predicted scaling laws [30]. Persson et al. additionally considered the non-uniformity of the pressure distribution inside the contact area [31]. Brochard-Wyart and Martin extended the model towards materials that are characterized by a static shear modulus, a high frequency modulus, and a relaxation time [32]. Carbone and Persson extended the model of Brochard and de Gennes towards dewetting of liquid films between viscoelastic materials [33].
Brochard et al. focused on axisymmetric dewetting induced by purposely introduced defects or protrusions on one of the confining solids. In contrast, in this manuscript we present systematic experiments of nonaxisymmetric dewetting. The experimental setup is illustrated in fig. 1(a). An elastomeric half-sphere is pressed onto a flat layer of the same material with liquid confined between them. Due to the applied pressure the liquid film thickness h decreases in time. After h has become sufficiently low, dry-spot nucleation may occur usually near defects and surface irregularities, as shown in fig. 1(d), (e). We developed a 3D fully coupled numerical model of the dewetting process, which combines thin film flow with elastic surface and bulk deformation and which reproduces many features observed in the experiments.

Analytical models
The stability of thin liquid films confined between solid elastic media is governed by the spreading parameter S ≡ γ AB −γ AL −γ BL [34], which depends on the interfacial energies of the contacting media, γ AB being the surface energy of the dry A-B elastomer interface, γ AL and γ BL the solid-liquid interfacial energies, respectively. For S > 0 the liquid is always stable, but for S < 0 dewetting may occur.
In the case of circular dry spots initiating in the center of the contact area, the growth dynamics may be predicted by comparing the rate of surface energy gain with the rate of viscous dissipation in the rim forming in front of the contact line (see fig. 1(e)). The dissipation rate depends on the shape of the rim, which depends on the surface and elastic properties of the materials in contact. Brochard-Wyart and De Gennes derived the following relation for the time evolution of the dry spot radius [29]: where η is liquid viscosity, Y is Young's modulus of the elastomer and h is the liquid film thickness (assumed uniform). Persson et al. additionally considered the nonuniformity of the pressure distribution inside the contact area [31]. They assumed a Hertzian distribution [35], leading to the following implicit solution for the time t needed for the contact line position to reach a certain position r cl : Here, r cs is the radius of the contact spot, τ the squeezeout time derived from eq. (1) by solving for the time when r cl (τ ) = r cs , and κ is a parameter dependent on the geometrical and material properties of the system where R L is the radius of curvature of the rubber lens and h is also assumed to be uniform. For values of κ typical for our experimental geometry and material systems, the solution of eq. (2) is very well approximated by eq. (1), up to approximately a third of the contact spot radius.

Numerical model
We developed a theoretical model to simulate soft elastohydrodynamic lubrication phenomena and the dewetting of thin films between soft elastic surfaces numerically. The model combines the Reynolds equation for thin liquid film flow and linear elasticity accounting for the stresses and deformations of the soft elastomeric components. The three-dimensional implementation allows us to study nonaxisymmetric geometries. All simulations were performed using the finite-element software COMSOL Multiphysics vs. 5.3.

Linear elasticity
The elastic deformation in the limit of slow deformation and small strain is governed by the Cauchy momentum equation for a stationary system where σ ij is the Cauchy stress tensor and f i are the Cartesian components of an external body force density vector.
Due to the small dimensions, we neglect gravity. The body force density is thus assumed to vanish f i = 0. We assume the elastic material properties to be linear, non-dissipative, isotropic and homogeneous. The corresponding constitutive equation is where are the shear and bulk moduli, Y and ν are Young's modulus and Poisson's ratio. We assumed ν = 0.49 to represent the effective incompressibility of the elastomer at the pressures relevant to our experiments [36]. Moreover, δ ij is the Kronecker delta and e ij is the infinitesimal strain tensor, which depends on the displacement vector u i via For computational simplicity we assume the substrate to be rigid and oriented parallel to the base plane of the hemisphere. In the appendix, we consider the case of elastomeric substrates using an Arbitrary Lagrangian-Eurlerian (ALE) approach. As for boundary conditions (BCs), the base of the hemisphere is assumed to be attached to a rigid surface and thus u i (r, z = 0) = 0. The surface of the solid not in contact with the liquid (i.e. the elastomer-air interface) is assumed to be traction-free and therefore i σ ij n i = 0. In the azimuthal direction the computational domain is bounded by two symmetry planes, ϕ = 0 and ϕ = ϕ 0 , where the normal displacement vanishes, u ϕ = 0. As a consequence, the radial displacement at r = 0 vanishes, u r (r = 0, z) = 0. The considered geometry and boundary conditions (BCs) are illustrated in fig. 2.

Thin liquid film flow
Since the typical film thickness (order 100 nm) is much smaller than the radius of the contact spot r cs (order 300 μm), a thin film approximation can be used to describe the liquid flow. Assuming impenetrable liquid-solid interfaces, incompressible laminar flow subject to no-slip BCs, the Reynolds equation can be solved to obtain the pressure distribution p f . Here, t is time, η is the dynamic viscosity of the liquid, and v e and v s are the tangential velocities of the liquid-elastomer interface and the substrate surface (see fig. 2), respectively, which are computed from their elastic displacements. The liquid film thickness h(x, y, t) is calculated from the positions of the boundaries h = z e − z s . Partial wettability in solid-liquid-solid systems can be implemented in an elastohydrodynamic framework by means of the concept of disjoining pressure Π [37]. We assume the following empirical form: with where S Π is an amplitude factor with units of surface tension (N/m), h * is a nanoscopic lengthscale, and n and m are integers (n > m) defining the shape of the disjoining pressure isotherm and the limiting behavior for h → 0 and h → ∞. We typically chose h * = 10 nm. The value of Π(h) is different from zero essentially only for ultrasmall film thicknesses (h 10h * ), whereas it is immeasurably small at microscopic distances, Π(h h * ) ≈ 0. The divergence of Π for h → 0 enforces a non-zero minimum film thickness, which helps alleviate the hydrodynamic stress singularity associated with moving contact lines [38].
The two-dimensional computational domain of the liquid phase has three boundaries. At the azimuthal bounding planes ϕ = 0 and ϕ = ϕ 0 , we prescribe symmetry conditions ∂h ∂ϕ = 0 and At the external boundary, i.e. outside of the contact spot, the pressure is set to zero.

Interfacial coupling conditions
At the elastomer-liquid boundary a traction T i = j σ ij n j is exerted by the liquid on the solid where δ ij is the unit tensor and n i are components of the unit normal vector of the liquid-elastomer interface. The second term in eq. (12) corresponds to a pressure load applied to the solid-liquid interface, which is a superposition of the total fluid pressure calculated from eq. (8) and the disjoining pressure Π(h).

Initial and kinematic conditions
In order to induce dry spot nucleation at a predefined location, we introduce a defect with shape as part of the surface topography of the rigid substrate.
Here, Δh d and w d define the height and width of the defect, r d is the defect position and f hs is a smooth Heaviside function with two continuous derivatives The dynamics of the system is generated by a vertical motion of the rigid substrate, while the base plane of the elastomeric hemisphere is held fixed. The vertical position of the liquid-elastomer interface is where u (e) z is the z-component of the displacement vector of the elastomer at the liquid-elastomer interface and r 2 ≡ We prescribe the time dependence of the vertical position of the liquid-substrate interface as where b 0 and b 1 are the initial and final positions of the substrate and t c is the approach time constant. The value of the time constant t c is restricted by the requirements of slow elastic deformation and by the timescale of the dewetting process. We chose t c much shorter than the time required for dry spot nucleation and found that the precise value of t c then effectively has no influence on the liquid film thickness evolution. At t = 0 the displacement field of the solid is set to zero. The initial pressure distribution is assumed to be zero in the entire liquid. The initial condition for h is determined by eqs. (15) and (16).

Axisymmetric variant of the model
If a dry spot nucleates in the center of a circular contact spot, dewetting can evolve in an axially symmetric manner. In this case we simplified the model and considered a one-dimensional liquid and a two-dimensional solid domain using cylindrical coordinates. Equation (7) Equation (8) becomes The boundary conditions remain in effect. At the symmetry axis r = 0 we have

Materials and methods
The elastomer used was a silicone-based, heat-curable, two-component polymer resin (Smooth-On, Encapso K). Samples were prepared by mixing the two components in a 1:1 weight ratio and depositing a droplet on a borosilicate glass slide. The droplet spreads into a layer approximately 200 to 300 μm thick and 2 cm wide. In order to obtain a spherical silicone cap, resin was deposited on the flat side of a fused silica half-ball lens (Edmund Optics, product number 67-395, radius of curvature 1.5 mm), which was glued to a borosilicate glass substrate using a UV-curable adhesive (Norland Optical Adhesive 65). The procedure was similar to the one described in [39], except that we used Gibbs confinement [40] to define the shape rather than a flat hydrophobic substrate. The radius of curvature of the spherical cap can be adjusted by the deposited volume of resin. After curing for at least 24 hours at room temperature, an optically transparent elastomer was obtained, with Young's modulus given by the manufacturer as 1.365 MPa. Compared to using micro-machined molds this contact-free fabrication procedure yields molecularly smooth surfaces. The liquid used in the experiments was a perfluoropolyether (Solvay, Fomblin Y LVAC 14/6, average molecular weight 2500). It is essentially non-volatile at room temperature, has a viscosity η of 140 cSt at 20 • C, a density of 1890 kg m −3 at 25 • C, and a refractive index n D of 1.298. With the chosen materials, swelling of the polymer was never observed. The spreading parameter |S| can vary with contamination and non-uniformity of the polymer surfaces. An estimate derived from the analysis of the shape of trapped liquid droplets, as discussed in ref. [41], is close to (10 ± 4) mN/m.
A detailed sketch of the experimental setup is presented in fig. 3. An incoherent blue LED light source (Roithner) is used with a center wavelength of 450 nm. A non-polarizing beam splitter (Thorlabs, product number CCM1-BS013/M) is used to direct the light through the objective (Olympus, LUCPlanFL N 20×/0.45NA) towards the sample and the reflected light from the sample to the CCD camera (Basler, model piA1000-60gm). A borosilicate glass cover slip with dimensions 1 mm × 25 mm × 75 mm is used as a rigid support for the silicone layer. The latter is covered with liquid and brought in wet contact with the silicone hemispherical lens. A grey filter with a drop of an index matching liquid is used to lessen unwanted reflections from the glass-air interface. The top glass slide is pivoting on a rubber platform on one side. On the other side it is supported by a force gauge (Strain Measurement Devices, model S100 0.2N) with measurement range 0.2 N and spring constant 89 N/m. By using a vertical translation stage the silicone lens can be lowered onto the silicone substrate. The load on the lens can be regulated from 0 to roughly half of the weight of the top glass slide.
Light reflected from the silicone-liquid interfaces gives rise to interference fringes ( fig. 1(c)-(e)), which allow to determine the film thickness distribution of the liquid. Details of the film thickness evaluation are provided in the electronic supplementary material (ESM).

Nucleation inside the contact spot
In fig. 4 the positions of the contact line are plotted (symbols) as a function of time for the experiment shown in fig. 1(d), (e), where dewetting started 91 μm from the center of the contact spot. The contact line positions were extracted in the directions radially inward and outward relative to the center of the contact spot. The dashed lines represent fitted powerlaw relations r cl ∼ t 3/4 according to eq. (1), which are excellent approximations during the early stage of the dewetting process. Later, the measured data deviate from this behavior, because the contact line is accelerated by the radially increasing pressure gradient in the contact spot. The dash-dotted lines correspond to eq. (2), where κ was treated as a fitting parameter, which shows excellent agreement over the entire data range.

Nucleation at the edge of the contact spot
Most frequently dry spots nucleate at the edge of the contact spot because the liquid film thickness is the smallest there. An example is shown in fig. 5(a). In such a case, the contact line velocity exhibits a high degree of anisotropy as shown in fig. 5(b)-(f). The dewetting speed is approximately twice as fast along the perimeter of the  contact spot compared to the direction radially inwards. In fig. 5(a) two dry spots nucleated roughly simultaneously at the edge of the contact spot. Both dry spots grow at the same rate and the complex contact line morphology evolves in an almost perfectly mirror-symmetric fashion. This is a clear indication that the dewetting dynamics and the contact line morphology are governed by elastocapillary and elastohydrodynamic effects rather than by randomly distributed surface irregularities. Figure 6 illustrates the time dependence of the contact line position measured both along the edge of the contact spot and radially inwards, as indicated with arrows in fig. 5(e). In the very early stage of dewetting the dry spot expands roughly isotropically. Given enough time, the speed towards the center (crosses in fig. 6) conforms to a power law s ∼ t 0.62 (solid line), whereas the speed along the perimeter becomes a constant, s ∼ t (dash-dotted line), and is significantly higher. We note that near the perimeter the pressure gradient has its strongest position dependence. Therefore, it is not to be expected that the time dependence of the contact line position should follow eq. (1), which assumes uniform pressure. Furthermore, we note that in our experiments the scaling exponents for the initial azimuthal and radial dewetting typically vary between 0.6 and 1.3, and 0.6 and 1, respectively, depending on the exact distance of the dewetting spot from the edge of the contact spot.
To investigate this phenomenon further, fig. 7 presents high magnification images of the region near the apex of the contact line moving along the perimeter. After the initial growth of a small dry spot ( fig. 7(a)) the liquid film thickness distribution is constant in the immediate vicinity of the apex of the contact line. Effectively the local morphology is time-invariant and propagates with a constant shape and azimuthal velocity. Thickness nonuniformities ("bumps") form on the neighboring liquid rim ( fig. 5(c)-(f) and fig. 7(c)-(f)). However, these do not influence the edge dewetting speed, as they occur sufficiently far away. The azimuthal dewetting speed along the perimeter is higher for several reasons. There is relatively less liquid to displace compared to the center of the con-  tact spot, as the film is the thinnest at the edge. Moreover, approximately half of the liquid is displaced into the liquid bulk volume outside of the contact spot, which has a much smaller viscous flow resistance. Furthermore, the outwards escape of liquid does not increase the elastic nor the surface energy of the system. Therefore, the liquid rim width and height at the apex are relatively small and do not increase in time.
We have systematically measured the dewetting speed both along the perimeter of the contact spot and radially inwards as a function of the liquid film thickness. We determined the film thickness close to the center of the contact spot. Because the liquid film thickness is essentially uniform everywhere within the contact spot (except for the region close to the edge), the exact location is immaterial. The results are presented in fig. 8. Since we saw in fig. 6 that r cl ∼ t β , the dewetting speed in the radi-ally inwards direction is not a constant but decreases as t β−1 . Consequently, a characteristic lengthscale or position needs to be chosen, at which to evaluate the inwards velocity. We chose the radius of the contact spot. We fitted the experimental data for r cl (t) with Accordingly, the contact line velocity v cl at a radial distance equal to r cs becomes v cl (r = r 0 ) = α 1 β βr which is represented by the orange rectangles in fig. 8. The dewetting speed increases for thinner films. The solid lines represent the power relation v cl,edge ∼ h −3/2 . Interestingly, the ratio between the radial and the edge dewetting speed is approximately constant for h ≤ 120 nm.

Numerical results
We now turn to the predictions of our numerical model presented in sect. 2. For the purpose of validation we start with the axisymmetric case and compare our results with the analytical predictions of eqs. (1), (2). Figure 9(a)-(c) shows a typical time evolution of the liquid film thickness for defect-induced dewetting. For t ≤ 0.15 s the elastic hemisphere remains essentially undeformed and the minimum film thickness occurs on the symmetry axis at r = 0. For later times deformation becomes apparent and the minimum film thickness is assumed at a non-zero radial position defining the contact spot radius a. Dewetting initiates at t d ≈ 15 s, after which a well defined contact line is seen to propagate radially outwards. We identify the contact line position r cl with the location at which h = 1.5 h * . Figure 9(d), (e) shows the time evolution of the pressure distribution and its (negative) radial derivative. The pressure distribution converges to the Hertzian profile [35] with decreasing film thickness. Figure 10 presents the time evolution of r cl along with fits according to eqs. (1), (2). The time t cs is defined as the instant when r cl reaches the edge of the contact spot for each numerical simulation. The numerical result for (n, m) = (3, 2) differs significantly from the analytical solutions. This is caused by the somewhat unrealistically long range of the surface forces present in the model for m = 2.

Axisymmetric dewetting -model validation
We repeated the simulations using higher exponents (n, m) equal (4, 3) and (5,4), which converge towards the curve corresponding to eq. (2). We note that the values of n and m are especially important for dry spot nucleation at the edge of the contact spot. A value of m = 2 induces rapid, spontaneous dewetting and a large number of dry spots all around the perimeter. This is in contrast to the  experimental observation of a low dry spot density of typically one over the entire contact spot area. Higher values of m ≥ 3 correspond to a faster decay of the disjoining pressure isotherm with film thickness, which implies that dewetting only initiates at the introduced defect. In these cases the dry spot density in the simulation is in accordance with our experimental results.
We now investigate the dependence on the various material parameters appearing in eq. (1) and compare with our numerical simulations. By differentiating eq. (1) we obtain the contact line speed v cl The time t 0 required for the contact line to reach a certain radial position r 0 is Substituting t 0 for t in eq. (22) gives the contact line speed at dry spot radius r 0 v cl (r cl = r 0 ) = 3 4η The liquid film thickness h is coupled with the parameters Y and |S|. Changing one of them will modify h at the time of dewetting as shown in fig. 11(a). According to ref. [29], a dry spot will grow and initiate the dewetting process, if its radius is larger than a critical value R crit , In our model, dewetting is initiated by a defect with halfwidth of w d = 5 μm when the liquid film thickness falls below a critical value h crit . Assuming R crit ∼ w d , we obtain from eq. (25), where E is a dimensionless constant. We tested the applicability of this equation in our model by investigating the relation between the average critical liquid film thickness h crit ≡ 2r −2 cs rcs 0 h(r) r dr evaluated at the time of dewetting and eq. (26). Figure 11(b) shows the result. The solid line represents a linear fit, which indeed represents the data very well. The proportionality factor is determined to be E = 18.9. Consequently, we can use eq. (26) to eliminate h from eq. (24) and obtain the following scaling law:  The solid lines correspond to power laws with fitted exponents −0.505 and 1.508, respectively, which are consistent with eq. (27). Figures 9 to 11 provide an extensive validation of our numerical model.

Dry-spot nucleation at the edge of the contact spot
Three-dimensional simulations of non-axisymmetric dewetting were performed, where dry-spot nucleation was initiated by a defect at the perimeter. Figure 12(a)-(c) shows snapshots from the experiment depicted in fig. 5. Figure 12(d)-(f) shows snapshots from a simulation where the dry-spot had the same radial extension as in the experiment. The simulation reproduces the dry-spot morphology and the effect of dewetting speed being higher along the perimeter. In both fig. 12(c) and (f) we note that the liquid rim developed bulges near the perimeter, as indicated by the blue triangles.
In fig. 12(g) we present the azimuthal and radial contact line displacements s a ≡ r cs ϕ cl and s r ≡ r cs − r cl , respectively, for the simulation shown in fig. 12(d)-(f). Here ϕ cl is the azimuthal position of the azimuthal apex of the dry spot as indicated in fig. 12(b). The solid lines correspond to power law relations s ∼ (t − t 0 ) ζ . The exponent found for the azimuthal motion matches very well with the early time behaviour (ζ = 0.65) in fig. 6. However, due to computational cost, the simulation did not

Conclusion
We have studied the dewetting of thin liquid films confined between a soft solid hemisphere and an elastomeric layer. For the most frequent case of dry-spot nucleation occurring at the perimeter of the contact spot, we found the contact line speed to be significantly higher in the azimuthal than in the radial direction. The spontaneous formation of dry spots is highly sensitive to the local film thickness h and exhibits a much higher rate and probability in regions of lower h, which predominantly occur near the margin of the contact spot. We found that the height of the dewetting rim in the azimuthal direction along the perimeter of the contact spot remains essentially constant, whereas it grows both in height and width for motion in or against the radial direction. Moreover the film thickness is the lowest at the edge of the contact spot. Both these effects account for the observed strong anisotropy of the contact line speed for non-axisymmetric dewetting.
We developed a fully coupled three-dimensional numerical model of the dewetting process based on a disjoining pressure formalism. We successfully validated our model using the available analytical models for axisymmetric dewetting. The model reproduces the morphology of the dewetting spots as well as their anisotropic growth dynamics very well.
The results we found for a liquid perfluoropolyether are representative of the behaviour of water and aqueous polymer solutions [12,13], which are relevant for biological systems. where σ (liquid) ij is the stress tensor, defined as The BC at the outlet is where p out is automatically adjusted to prevent backflow.

Appendix A.2. Impact of substrate elasticity
In fig. 13 we compared simulations for a rigid and soft substrate. In both cases, the film thickness in the center scales asymptotically as h(0, t) ∼ t −1/2 . However, fig. 13 shows that the liquid film thickness is more non-uniform for rigid substrates. In the latter case, h is thicker at the center and thinner at the edge of the contact spot. The higher degree of uniformity for soft substrates is consistent with the experimental observation that the grayscale value of the interference pattern does not vary significantly in the contact spot (see, e.g., fig. 1(d)). The displacement of the substrate support was kept constant at 45 μm in the simulations shown in fig. 13, resulting in a 19% larger force for the rigid substrate. Moreover, the contact spot radius is slightly larger for the soft substrate.
While the ALE simulations more faithfully represent the experimental geometry, they imply a higher computational cost. The introduction of a non-zero disjoining pressure caused convergence issues such that simulations of dewetting unfortunately were not possible.

Open Access
This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.