On the singular nature of the elastocapillary ridge

The statics of wetting amounts to finding the state of minimal elastocapillary energy. The substrate deformation is described by a mapping from the reference state prior to deformation, to a current state after deformation. Following standard notation, the mapping is written as

where $\mathbf{X}$ is the position of a material point on the reference domain, mapped onto its current position $\mathbf{x}$. We consider the geometry to be invariant along the contact line, so that the problem is two-dimensional (plane strain elasticity). Hyperelastic solids are described by an elastic energy density $W(\mathbf{F})$, which depends on the deformation gradient tensor $\mathbf{F}=\partial\mathbf{x}/\partial\mathbf{X}$. We now turn to the interface, which in the (plane strain) two-dimensional description is one-dimensional. We define the arc-length material coordinate at the interface as $S$, and the current surface position $\mathbf{x}_{s}(S)=\chi(\mathbf{X}(S))$. The surface stretch, accounting for the change of length of surface elements follows as

which is a scalar in this plane strain description. Now we can express capillarity, usually defined by the excess energy $\gamma$ per unit area of the deformed state, as a free energy $\lambda\gamma$ per unit area in the reference state.

Crucially, elastic media can exhibit a nontrivial capillarity where the surface energy $\gamma(\lambda)$ will itself be a function of the stretch $\lambda$ – this is the Shuttleworth effect Shuttleworth (1950); Andreotti and Snoeijer (2016); Style et al. (2017); Andreotti and Snoeijer (2020). With this, the elastocapillary energy (per unit length along the contact line) takes the form

respectively giving the total (bulk) elastic energy and the (surface) capillary energy. $\mathbf{F}$ and $\lambda$ are the corresponding bulk and surface stretches, and are both defined by the map $\chi(\mathbf{X})$. We anticipate that we will consider incompressible substrates, in which case the constraint of incompressibility will be included in $W(\mathbf{F})$.

Equilibrium configurations of the elastocapillary substrate are found by minimising the functional $\mathcal{E}$ with respect to $\chi(\mathbf{X})$. Considering variations $\delta\mathbf{x}=\delta\chi(\mathbf{X})$, we find

Here we introduced the nominal (or first Piola-Kirchhoff) stress tensor $\mathbf{s}$, and the surface tension $\Upsilon$,

where for the latter we indeed recognise the Shuttleworth relation (1). In addition, we used that $\delta\lambda=\mathbf{t}\cdot\partial\delta\mathbf{x}/\partial S$ along the boundary, where $\mathbf{t}$ is the surface-tangent unit vector in the current configuration.

To study the elastocapillary ridge, we still need to include the pull of the contact line, induced by the liquid drop that is wetting the solid. This can be achieved by making explicit the capillary energy of the drop, via its liquid-vapour surface energy $\gamma_{LV}$. The subtlety here is that one needs to impose a constraint at the contact line  Lubbers et al. (2014); Snoeijer et al. (2018): the position $\mathbf{x}$ of the liquid-vapour interface must (by definition) coincide with that of the solid interface. The effect of this constraint, imposed by a Lagrange multiplier, provides a localised traction on the substrate, pulling with a strength $\gamma_{LV}$ along the direction of the liquid-vapour interface $\mathbf{t}_{LV}$ ¹¹1Note that in the present work the solid interface is treated as part of the substrate, so that the external traction of the liquid-vapour interface is $\gamma_{LV}\mathbf{t}_{LV}$. The force transmitted onto the elastic bulk of the substrate, i.e. after passing the interface, is more intricate as discussed e.g. in Marchand et al. (2012a); Bostwick et al. (2014); Andreotti and Snoeijer (2016). The representation by a local force is indeed commonly used in modelling approaches Style and Dufresne (2012); Limat (2012); Style et al. (2013); Wu et al. (2018); Masurel et al. (2019). Here we therefore treat the contact line as a perfectly localized external traction, with the associated work functional $\mathcal{R}=\gamma_{LV}\mathbf{t}_{LV}\cdot\mathbf{x}(\mathbf{X}_{\rm cl})$, where $\mathbf{X}_{\rm cl}$ is the solid’s material point at which the contact line is acting. During the variation this corresponds to a work

The virtual work principle, $\delta\mathcal{E}=\delta\mathcal{R}$, then gives the equilibrium condition

which should be satisfied for arbitrary $\delta\mathbf{x}$.

Equation (8) defines the elastocapillary equilibrium in the weak formulation. This equilibrium is indeed highly singular. Namely, the forcing on the right hand side appears as a point force, pulling at $\mathbf{X}_{\rm cl}$, while the elastocapillary energies on the left contains only surface and bulk contributions. The debate in the literature precisely revolves around the following question: Do singularities appear in surface (capillarity) or in bulk (elasticity), in order to balance the point force at the contact line?

Our interest pertains to finding equilibrium configurations of the elastocapillary problem, i.e. to minimisers of the energy functional in (4) extended with the work functional $\mathcal{R}$ representing the contact line, subject to appropriate boundary conditions. Specifically, we consider substrates that are flat in the reference configuration, with complete fixation at the bottom boundary and guided fixation (slip) at the lateral boundaries. We allow for the possibility to impose a prestretch $\lambda_{\infty}$, refering to the uniaxial stretch far away from the contact line. Besides the work associated to the point-forcing at the contact line, the top surface is free of traction, as is made explicit in the weak formulation (8) of the minimisation problem. The constitutive relations for the strain-energy density and the surface energy are specified in Section III below. In all simulations, the shear modulus $G$ and the relevant surface energies are chosen such that the wetting ridge is much smaller than the width of the domain, with a typical example given in Fig. 1(a). In that example the domain width and height respectively are $8\gamma_{LV}/G$ and height $\frac{8}{3}\gamma_{LV}/G$, which are representative for all presented results.

Here we numerically approximate the minimiser of $\mathcal{E}-\mathcal{R}$ by means of a goal-adaptive finite-element method Becker and Rannacher (2001); Oden and Prudhomme (2001). In goal-adaptive methods, the finite-element approximation is locally refined on the basis of an a-posteriori error estimate, in such a manner that an optimal approximation to a predefined quantity of interest (the goal) is obtained. Goal-adaptive finite-element methods generally proceed according to the SEMR (Solve $\rightarrow$ Estimate $\rightarrow$ Mark $\rightarrow$ Refine) process Nochetto and Veeser (2012); van Brummelen et al. (2017b). The SEMR process starts by solving a finite-element approximation on a coarse mesh. Next, the contribution of each element to the error in the goal quantity is estimated, based on a so-called dual problem Becker and Rannacher (2001); Oden and Prudhomme (2001); van Brummelen et al. (2017b). The elements that yield the largest contribution to the error are marked according to a refinement strategy. These marked elements are subsequently refined by subdivision. This process is repeated until a certain threshold for the error estimate is satisfied or a preset number of refinement iterations has been executed. In accordance with our interest in minimisers of $\mathcal{E}-\mathcal{R}$, we take the energy itself as the goal functional. The optimality conditions are resolved by means of the Newton–Raphson method. The goal-adaptive finite-element method for the present problem has been implemented in the open-source software framework Nutils van Zwieten et al. . The optimality conditions (8) are in fact directly derived from an implementation of the energy functional $\mathcal{E}-\mathcal{R}$ via the automatic-differentiation functionality in Nutils.

An illustration of a goal-adaptive finite-element approximation is provided in Fig. 1(a). The approximation is based on 16 refinement iterations. Accordingly, the smallest elements in the adaptive approximation are $2^{16}$ times smaller than the initial element size. The initial mesh comprises $24\times 8$ uniform quadrilateral elements and, correspondingly, the smallest elements are 5–6 orders of magnitude smaller than the elastocapillary length. Importantly, the adaptive procedure automatically introduces the local refinements in the vicinity of the contact line. This refinement pattern is in agreement with the singularity of the pressure towards the contact line, and we extensively verified the numerical convergence of the result. For the result shown in Fig. 1(a), the relative numerical error in the computed value of the solid opening angle $\theta_{S}$ is less than $10^{-6}$.

We first consider the solid angle $\theta_{S}$, as measured at the tip of the wetting ridge in FEM. Figure 1(b) shows $\theta_{S}$ plotted against $\gamma_{LV}/\Upsilon_{S}$, with the value of $\Upsilon_{S}$ taken at the tip. Clearly, $\theta_{S}$ follows a universal curve for all cases considered. The parameters that were varied in our simulations are the contact line force $\gamma_{LV}$ (compared to the value of $\gamma_{0}$ in (10)), while solid surface tensions are with or without Shuttleworth effect ($c_{0,1}=0$, respectively $c_{0,1}=1$). We also considered different amounts of prestretch of the substrate, ranging from $\lambda_{\infty}=1$ (no prestretch) to $\lambda_{\infty}=2$ (extending the length by 100%). The universal curve for $\theta_{S}$ indeed follows Neumann’s law, which for the specific case of identical solid-liquid and solid-vapour energies reads

Here, we emphasise that owing to the Shuttleworth effect, the surface tension $\Upsilon_{S}(\lambda)$ depends on the strain. Since the Neumann balance is to be interpreted as a boundary condition at the contact line, we consider (13) with values of the stretch $\lambda_{\rm cl}$ taken at the contact line. The result of (13) is superimposed as the solid line in Fig. 1(b), providing a perfect description of the FEM results.

We thus reach a first major conclusion: Neumann’s law (based on the local values of the surface tension) universally applies to elastocapillary wetting ridges, irrespective of the large elastic deformations at the contact line. This rejects the recent hypothesis that strong elastic nonlinearity, as encountered for narrow $\theta_{S}$ and large prestretch, would lead to a failure of Neumann’s law Masurel et al. (2019). The universal validity of Neumann’s law has an immediate consequence for measurements of the surface-constitutive relation based on $\theta_{S}$, since we safely conclude that $\theta_{S}$ gives a direct access to the values of $\Upsilon_{S}$. Phrased differently, the experimental observation for PDMS that $\theta_{S}$ increases with $\lambda_{\infty}$ Xu et al. (2017) can, in a macroscopic theory based on hyperelasticity, only be explained via a strong Shuttleworth effect.

To further illustrate this, we closely follow the experimental protocol of Xu et al. (2017) in our simulations, and consider how the geometry of the ridge evolves when stretching the substrate by an increasing amount $\lambda_{\infty}$. Figure 2(a) shows $\theta_{S}$ versus the stretch at the contact line $\lambda_{\rm cl}$. The open circles are FEM results with a Shuttleworth effect ($c_{0,1}=1$, and $\gamma_{0}=\gamma_{LV}$), showing an increase of the solid opening angle $\theta_{S}$. Indeed, the dependence of $\theta_{S}$ is perfectly predicted by Neumann’s law (13), as is indicated by the solid line. In experiments, one of course does not control the stretch at the contact line $\lambda_{\rm cl}$, but rather the global stretch of the substrate $\lambda_{\infty}$. In Fig. 2(b) we therefore plot these two stretches against one another. While $\lambda_{\rm cl}$ is not exactly identical to the imposed stretched $\lambda_{\infty}$, the differences turn out to be minor – consistently with experiments Xu et al. (2017)). As a guide to the eye, the dashed line in Fig. 2(b) indicates $\lambda_{\rm cl}=\lambda_{\infty}$. We expect this near-homogeneity of $\lambda$’s to arise only for nearly symmetric $\gamma_{SL}$ and $\gamma_{SV}$, as asymmetry in general leads to stronger gradients of stretch (cf. Sec V).

The scenario changes dramatically when the substrate does not exhibit a Shuttleworth effect (i.e. $c_{0,1}=0$). In that case, both $\theta_{S}$ and $\lambda_{\rm cl}$ take on a constant value, that is totally independent of the imposed $\lambda_{\infty}$. This is indicated in Fig. 2(a) by the closed circle – which in fact corresponds to various simulations with $\lambda_{\infty}$ ranging from 1 to 2. This invariance of $\theta_{S}$ with respect to $\lambda_{\infty}$ is easily understood from the Neumann balance. Namely, surface tensions are constant when $c_{0,1}=0$, and since we consider $\gamma_{0}=\gamma_{LV}$ we find that $\theta_{S}=120^{\circ}$. By contrast, the invariance of the stretch at the tip comes as a surprise and its explanation calls for a better understanding the nature of the elastic singularity. Below, we will derive analytically that without the Shuttleworth effect, $\lambda_{\rm cl}=\sqrt{\pi/\theta_{S}}$, irrespective of the externally imposed prestretch $\lambda_{\infty}$ of the substrate.

Measurements of the stretch at the contact line thus provide important additional information on the Shuttleworth effect, that till date has not yet been explored. Namely, experiments by Xu et al. (2017) reveal an increase of stretch at the contact line upon a global stretching of the substrate. From the above it is clear that such a dependence can, in a macroscopic theory based on hyperelasticity, not occur when there is no Shuttleworth effect.

To further analyse the vicinity of the tip, we now turn to the elastic stress measured along the free surface. In Fig. 3(a,b) we plot the pressure $p$ as a function of the distance to the contact line $x$, on a semilogarithmic scale. In all cases the FEM simulations exhibit a weak singularity of the pressure, diverging logarithmically with the distance to the tip.

Panel (a) corresponds to a case without Shuttleworth effect $(c_{0,1}=0)$, for different ratios $\gamma_{LV}/\gamma_{0}$. With this, we cover a broad range of $\theta_{S}$ down to very narrow angles with $20^{\circ}$. The prefactor of the logarithmic pressure singularity is larger for narrow $\theta_{S}$. The pressure plotted in Fig. 3(a) is scaled by $G\times\left(\frac{\pi}{\theta_{S}}-\frac{\theta_{S}}{\pi}\right)$, which indeed captures the $\theta_{S}$ dependence of the prefactor of the singularity. We remark that for very narrow angles the logarithmic asymptotic only emerges at distances much below the elastocapillary length $\gamma_{0}/G$; this illustrates the challenge of accurate numerical resolutions for small $\theta_{S}$. Panel (b) corresponds to the case with a strong Shuttleworth effect $(c_{0,1}=1)$, for different amounts of substrate prestretch $\lambda_{\infty}$ (the corresponding $\theta_{S}$ are in Fig. 2). Figure 3(b) again reveals a logarithmic singularity of pressure, with a weak variation of the prefactor with $\lambda_{\infty}$.

Interestingly, the Shuttleworth effect allows for a new phenomenon induced by gradients of surface tension. For liquid interfaces, gradients in surface tension arise due to gradients in composition or in temperature – this is known as the Marangoni effect, and leads to tangential interfacial stress. For the elastic interfaces considered here, the gradients in surface tension are due to gradients of $\lambda$ along the interface. Given this analogy, we refer to this as the elastic Marangoni effect.

Figure 3(c) indeed reveals the emergence of elastic (Cauchy) shear stress $\sigma_{nt}$ along the interface, which we will refer to as elastic Marangoni stress. Somewhat surprisingly, the Marangoni stress is not singular, but converges to a constant value upon approaching the contact line. This elastic Marangoni stress can be positive or negative, depending on the prestretch that is imposed. Without prestretch ($\lambda_{\infty}=1$), the contact line region will have the largest surface tension, giving a Marangoni stress that is oriented towards the contact line ($\sigma_{nt}<0$). Conversely, when the imposed $\lambda_{\infty}$ is large, the contact line region has the smallest surface tension and the Marangoni stress is directed away from the contact line ($\sigma_{nt}<0$). This is further quantified in Fig. 4, where the change of direction of the Marangoni effect is observed to be close to $\lambda_{\rm cl}\approx 1.2$. Indeed, this nearly coincides with the point where $\lambda_{\rm tip}\approx\lambda_{\infty}$ [cf. Fig. 2(b)]. So the orientation of the Marangoni stress depends on whether the stretch at the tip is larger or smaller than the stretch imposed at large distance.

We will now pursue a fully analytical theory for the numerical observations above. We have seen that the elastic singularity is weak, only logarithmic in the stress, so we first try to split off the singularity. For this, we perform an integration by parts on (8) by writing

The integral over the third term indeed gives point-like contributions

where the sum runs over all possible discontinuities along the contour. The term $\mathrm{Div}(\mathbf{s}\cdot\delta\mathbf{x})$ can be brought to the surface using the divergence theorem. For a smooth domain of integration, the divergence theorem holds for any vector field which is in ${\mathcal{L}}^{1}$ and whose spatial derivatives are in ${\mathcal{L}}^{1}$ Willem (2013.). This is allowed as long as the corresponding singularity is weaker than $1/|\mathbf{X}|$.

where in the last step we transformed the result to the current domain, using the definition of the true stress (or Cauchy stress) according to $\boldsymbol{\sigma}=\mathbf{s}\cdot\mathbf{F}^{T}/\det(\mathbf{F})$.

The condition of equilibrium, $\delta\mathcal{E}=\delta\mathcal{R}$ obtained from (7),(16), then splits into bulk, surface and point conditions:

where $\mathcal{D}$ denotes the current domain of the deformed state. Besides the classical elastic stress equilibrium in bulk (17), the interface condition (18) gives the Marangoni effect where $\sigma_{nt}\equiv\mathbf{t}\cdot\boldsymbol{\sigma}\cdot\mathbf{n}$ balances gradients in surface tension $\partial\Upsilon/\partial s$, while the normal component of elastic stress $\sigma_{nn}\equiv\mathbf{n}\cdot\boldsymbol{\sigma}\cdot\mathbf{n}$ balances the Laplace pressure. Finally, the Neumann condition appears at the contact line, equation (19), expressed as a discontinuity of the surface tangents. The only assumption made in the derivation above is that the stress singularity is sufficiently weak for the divergence theorem to be applicable, as is the case for a logarithmic singularity.

We now analytically establish the nature of the elastic singularity, through an asymptotic analysis near the contact line. For this we express the mapping $\chi(\mathbf{X})$ in polar coordinates, $(r,\varphi)$ and $(R,\Phi)$, respectively for the current and reference state. The contact line is located at $r=0$ and $R=0$, and without loss of generality the initially flat free surface is chosen to be along the lines $\Phi=0$ and $\Phi=\pi$. We make use of the fact that the boundary condition (19) forces the solid into an angle $\theta_{S}$, which is defined by the property

As is common with singularities Eggers and Fontelos (2015), we expect the asymptotics to be scale-invariant, so we propose a similarity Ansatz

Imposing (20) one finds that $\beta=0$. A critical feature of soft elastic solids is that these are basically incompressible, i.e. $\mathrm{det}(\mathbf{F})=1$. Combined with $\beta=0$, this then dictates $\alpha=1$, which implies that the radial stretch $\lambda_{r}=dr/dR$ remains finite and is independent of $R$. In the azimuthal direction, incompressibility implies a relation between the functions $g_{1,2}$, which can be accounted for by writing

so that the solid angle follows as $\theta_{S}=\int_{0}^{\pi}d\Phi\,f(\Phi)$. The deformation gradient tensor of this mapping reads

which indeed satisfies $\mathrm{det}(\mathbf{F})=1$ for arbitrary $f(\Phi)$. The corresponding Finger tensor reads

This defines the most general scale-invariant incompressible map that generates a corner.

For the special case where $f^{\prime}=0$, one recovers the classical solution by Singh & Pipkin Singh and Pipkin (1965). However, that solution is shear-free (i.e. $\mathbf{F}$ and $\mathbf{B}$ are diagonal) and therefore cannot be universally valid. Here we derive the most general corner solution that satisfies mechanical equilibrium, $\mathrm{div}(\boldsymbol{\sigma})=0$. We focus on a neo-Hookean material defined by (9), which has a Cauchy stress $\boldsymbol{\sigma}=G\mathbf{B}-p\mathbf{I}$, so that (17) becomes

This implies that $\textrm{div}(\mathbf{B})$ must be irrotational, i.e. $\mathrm{curl}\left(\mathrm{div}(\mathbf{B})\right)=0$, which here takes the form

where $K$ is an integration constant. Inserting (24) and bearing in mind that $\partial/\partial\varphi(\cdots)=(\cdots)^{\prime}/f$, we find

This is a nonlinear second order ODE for $f(\Phi)$. As boundary conditions we impose the stretch at each of the boundaries, which will subsequently give the shear stress via the connections

We note that $\lambda_{r}$ in the similarity solution is independent of $R$, and can therefore be identified to the stretch at the contact line $\lambda_{r}=\lambda_{\rm cl}$. The constant $K$ can be adjusted to accommodate the desired $\theta_{S}$. Explicit solutions will be presented below, and compared directly to FEM simulations.

Once a solution is found, one can explicitly integrate (25) to obtain the pressure

up to an integration constant. This completes the analytical description of incompressible corner solutions in the fully nonlinear regime.

The similarity solutions derived above capture all FEM results of Sec. III, in the vicinity of the contact line. First, we consider the stress, which for a neo-Hookean solid is given by $\boldsymbol{\sigma}=G\mathbf{B}-p\mathbf{I}$. Our theory explains the FEM result that the normal stress diverges logarithmically, following the singularity of pressure (29), and offers a way to compute the prefactor $K$. Furthermore, the corner solution shows that $\mathbf{B}$ as given in (24) remains finite at the contact line. This explains why the Marangoni stress $\sigma_{nt}=\sigma_{r\varphi}$ remains finite at the contact line.

We now turn to a fully quantitative analysis, by solving (27) for various boundary conditions. Typical (symmetric) similarity solutions are represented graphically in Fig. 5 denoting the Lagrangian grid both in undeformed (grey) and deformed (black) configurations. The three panels each correspond to $\theta_{S}=120^{\circ}$, with different amount of stretch imposed on the free surfaces. In Panel (a) we report the solution without shear stress, for which $f^{\prime}=0$ for all $\varphi$. In this case, (22) reduces to the classical solution by Singh & Pipkin Singh and Pipkin (1965), with the constant $f=\theta_{S}/\pi$. In the context of elastocapillary ridges, the absence of shear corresponds to a substrate without a Shuttleworth effect. This explains why in the absence of a Shuttleworth effect, the stretch at the contact line $\lambda_{\rm cl}$ was found to be independent of $\lambda_{\infty}$ in our FEM simulations: in a shear-free corner, the stretch takes on a specific value that depends only on the solid angle, as $\lambda_{r}=\lambda_{\rm cl}=\sqrt{\pi/\theta_{S}}$. The stretch at the contact line is therefore locally determined by $\theta_{S}$, irrespective of the conditions imposed at large distance. Furthermore, in this specific case without shear stress, we find an analytical expression for the strength of the pressure singularity, the constant $K$ in (29). Inserting $f=\theta_{S}/\pi$ in (27) gives $K=\frac{\pi}{\theta_{S}}-\frac{\theta_{S}}{\pi}$. Indeed, this was exactly the scaling used in Fig. 3(a), necessary to account for the $\theta_{S}$ dependence. This demonstrates that the corner solutions are fully quantitative and provide the correct asymptotics observed in FEM, valid in the strongly nonlinear regime.

The Shuttleworth effect dramatically changes the physical picture. Now, a variety of surface stretches $\lambda_{r}$ is possible, as shown in Fig. 5(b,c). Each of these solutions comes with its own value of the elastic Marangoni stress. Figure 4 illustrates this point, where the prediction of the similarity solutions is shown as a solid line and compared directly to the Marangoni stress in FEM. For the symmetric surface tensions considered in our simulations, the corresponding similarity solution is naturally symmetric and can be found without any adjustable parameters: it follows directly from the surface constitutive relation (10), which in combination with Neumann’s law determines the appropriate combination of $\theta_{S}$ and $\lambda$. The perfect prediction of the elastic Marangoni stress in Fig. 4 confirms that the corner solutions indeed offer the correct asymptotic description of the singularity – also in the presence of the Shuttleworth effect.

As a conclusive remark, we emphasise again that the observation in PDMS that $\lambda_{\rm cl}$ increases upon varying $\lambda_{\infty}$ Xu et al. (2017) cannot be explained in a hyperelastic theory without a Shuttleworth effect.

So far we have considered an isolated contact line, at some prescribed position $\mathbf{X}_{\rm cl}$, pulling vertically with perfectly symmetric wetting conditions. In a real wetting problem, however, a droplet will spread dynamically until it reaches its equilibrium liquid angle – simultaneously the contact line reaches an equilibrium material position $\mathbf{X}_{\rm cl}$, which is not known a priori but which needs to be found self-consistently. Hence, the full equilibration involves a free exploration of the contact line over the substrate. Technically, such an equilibrium without pinning implies that the change of material coordinate is energetically neutral. Naturally, this is the case when the substrate is perfectly homogeneous in its reference state. Indeed, in contrast to the rigid case, there are various examples where well-prepared soft polymeric substrates are basically free from pinning and contact angle hysteresis Xu et al. (2017); Schulman et al. (2018); Snoeijer et al. (2018); Lhermerout et al. (2016).

Here we take the opposite perspective and consider the possibility that the presence of the contact line itself induces heterogeneity in the material – in its reference state. Even when the originally prepared soft polymeric substrate does not exhibit permanent defects that provide a frozen surface energy landscape, the substrate can develop heterogeneities dynamically, due to the presence of the contact line. Indeed, a large-stress region builds up at small scale, which can lead to irreversible plastic flow, like in the “Fracture Process Zone” that forms at a crack tip. Although wetting-induced damaging processes have been evidenced in experiments where a soft gel exhibits fracture by wetting Bostwick and Daniels (2013), we focus here on non-damaging plastic deformations in the near-surface region – so that the bulk reference is not affected. Plasticity typically occurs in situations where there is multistability, where multiple stable configurations coexist, which can lead to a hysteretic response upon contact displacement Caroli and Nozières (1998). The large strain may indeed provide a configurational plasticity, without damaging the material. When chains between cross-linkers are long enough to produce entanglement, strain may trigger changes of glassy chain conformation. As an alternative mechanism, the contact line may lead to a local strengthening associated with the elongation of polymeric chains, producing a highly dissipative zone when the contact line explores its environment. Below we derive the consequences of a non-damaging, plastic process zone induced by the presence of the contact line. By analogy with fracture mechanics, or with defects in crystalline solids, such a plastic process zone can be described by a defect singularity in the theory of elasticity. The singularity then represents the effect of the plastic process zone on the elastic “outer” region. We reveal how the strength of such a defect directly relates to contact angle hysteresis.

The consequence of a defect, representing the effect of a process zone on the outer region, can be computed from the change in energy associated to a global displacement of the solution. This is illustrated in Fig. 6, showing such a displacement $\delta\mathbf{X}_{\rm cl}=\delta U\mathbf{T}$ on the reference domain (panel a), and on the current domain (panel b). The change in elastic energy associated to the displacement of a defect is known as the Eshelby force Eshelby (1975),

where the integral encloses the defect and we define the Eshelby’s energy-momentum tensor

The Eshelby force reduces to the J-integral in small deformation (linear) elasticity, where it finds an interpretation as the energy release rate in fracture mechanics Rice (1968). To derive the capillary energy released by moving a defect, it is instructive to follow the derivation of the classical result (30), which is based on the application of Noether’s theorem in the space of material coordinates Eshelby (1975).

On the reference domain, the displacement simply amounts to a translation $\delta\mathbf{X}_{\rm cl}=\delta U\mathbf{T}$ of the contact line force, as in Fig. 6(a). The corresponding translation on the current domain is sketched in Fig. 6(b). The idea of deriving the elastic energy released by displacing a defect is to interpret the translation $\delta U$ as a variation $\delta\mathbf{x}$, which can be expressed as

The associated change in elastic energy can be computed from this variation, as

Importantly, in the last step one uses that the (reference) substrate is homogeneous everywhere except at the defect Eshelby (1975). When in the vicinity of the singularity $\mathbf{\Pi}\sim 1/|\mathbf{X}-\mathbf{X}_{\rm cl}|$, the integral is finite and can be expressed as (30). When the material is homogeneous everywhere, i.e. no defects, the Eshelby force uniformly vanishes as a consequence of translational invariance.

We now follow the same scheme for the capillary energy, upon replacing $W$ by $\lambda\gamma$, and the deformation gradient tensor $\mathbf{F}$ by its vectorial surface analogue $\mathbf{F}_{s}=\partial\mathbf{x}_{s}/\partial S$. Subsequently, we define the surface-equivalent of the Eshelby tensor (31), which now is a scalar, and which takes the form:

which is the chemical potential anticipated in (12). Indeed, the associated change in capillary energy reads

where the integral runs over an infinitesimal domain across the singularity. It is clear that a finite capillary defect-energy appears only when $\mu$ exhibits a discontinuity at the contact line, i.e. $[\mu]_{-}^{+}\neq 0$.

We thus conclude that the total energy release rate $\Gamma$, liberated upon displacing the elastocapillary defect at the contact line, takes the form

Given that the defect represents a process zone, this indicates a loss of energy $-\Gamma\delta U$, dissipated inside the process zone during the translation. For the special case where there is no pinning defect and the contact line is free to move, the variation of the contact line position should be energetically neutral, so that $\Gamma=0$.

The notion of the (elastic) Eshelby force in wetting was recently proposed in Masurel et al. (2019), where it was argued that the formation of a ridge would already be sufficient to induce an elastic Eshelby force. However, from the above it is clear that this is not the case when the substrate is perfectly homogeneous in its reference state, so that there is a translational invariance of the space of reference coordinates: applying Noether’s theorem to this translational invariance Eshelby (1975), one finds $\partial\mathcal{E}_{\rm el}/\partial U=0$. This vanishing of the Eshelby force is indeed confirmed by our FEM results and analytical solutions: the stress is only logarithmically singular, so that for an infinitesimal integration volume around the contact line (30) gives $f_{\rm Esh}=0$. Therefore, for homogenous substrates, the condition $\Gamma=0$ reduces to $[\mu]_{-}^{+}=0$. The continuity of $\mu$ across the contact line can be interpreted as an “equality of chemical potential”, necessary for a free exchange of material points across the contact line. This condition of no-pinning was previously derived within the strong restrictions of linear elasticity Snoeijer et al. (2018) – but it turns out to be valid also when deformations are large.

For a non-damaging process zone, i.e. the reference state remains intact, we expect the Eshelby force to vanish owing to translational invariance. Nonetheless, a capillary defect $\Gamma$ could still emerge, associated with the interfacial microstate of the polymer.

Up to here we have considered properties of the solid, and did not discuss explicitly the liquid. Yet, the liquid contact angle $\theta_{L}$ is the prime feature that characterises the wetting of a liquid drop. To complete the theory, we now show how the equilibration determines $\theta_{L}$ on homogeneous substrates – and how the maximum strength of a contact line defect can be related to contact angle hysteresis on elastic substrates.

We restrict ourselves to the case of a sufficiently large drop, so that far away from the contact line one encounters a flat substrate (Fig. 6). At a large distance from the contact line, the substrate respectively has a solid-liquid energy $\gamma_{SL}(\lambda_{\infty})$ and a solid-vapour energy $\gamma_{SV}(\lambda_{\infty})$. The usual argument leading to Young’s law for the contact angle amounts to the global horizontal displacement de Gennes et al. (2002). In the present case the (Eulerian) displacement reads $\lambda_{\infty}\delta U$, so that the solid capillary energy increases by $(\gamma_{SL}-\gamma_{SV})\lambda_{\infty}\delta U$, the value of which has to be taken far away from the contact line. This balances the work $-\gamma_{LV}\cos\theta_{L}\lambda_{\infty}\delta U$ performed by the liquid-vapour interface, which together gives Young’s law. The situation is modified by the presence of a defect: as described above, such a displacement also involves a dissipation inside the process zone, indicating a loss of energy $-\Gamma\delta U$. By consequence, we find a modification of Young’s law

where in the second line we anticipate that $f_{\rm Esh}=0$ (owing to the weak logarithmic elastic singularity). For homogeneous substrates $\Gamma=0$, and we recover Young’s law for the liquid contact angle. We remark that, $\theta_{Y}$ is based on the surface energies corresponding to $\lambda_{\infty}$.

The analysis above, in particular (37), can be verified by the FEM simulations. In the numerics, we fix a priori the material position $\mathbf{X}_{\rm cl}$ of the pulling force, so that we effectively work with a pinned contact line. For symmetric surface tensions and pulling vertically, this is equivalent to the unpinned case, but we can consider any liquid angle $\theta_{L}$, by changing the pulling direction $\mathbf{t}_{LV}=(-\cos\theta_{L},\sin\theta_{L})$ in (8). We then measure the jump $[\mu]_{-}^{+}=[\mu]_{SL}^{SV}$ across the contact line obtained for the corresponding solution, as a function of $\theta_{L}$. We consider two cases: (i) symmetric surface energies $\gamma_{SL}=\gamma_{SV}$ (so that $\theta_{Y}=90^{\circ}$), and (ii) asymmetric surface energies $\gamma_{SL}\neq\gamma_{SV}$ (here with $\theta_{Y}=113.6^{\circ}$).

The result is presented in Fig. 7(a). It is clear that both cases, symmetric and asymmetric, are in perfect agreement with (37) with $\Gamma=-[\mu]_{SL}^{SV}$. This implies that $f_{\rm Esh}=0$, consistent with the weak logarithmic singularity. Hence, $\theta_{L}$ can be different from its equilibrium value $\theta_{Y}$ by the presence of a non-damaging process zone, represented by a capillary defect. In that case, interfacial plasticity could be associated with a contact angle hysteresis. A typical asymmetric similarity solution is shown via the grid representation in Fig. 7(b), for which there is a jump in stretch at across the contact line. We remark that in all cases, Neumann’s law was still observed to be valid, irrespective of the defect.

First boundary condition. In this macroscopic description, it was found that the stress singularity associated with the contact line ridge is weak (i.e. logarithmic) and therefore integrable, under all conditions that were considered. Hence, the singularity does not behave analogously to an elastic disclination defect and no qualitative difference emerges when the substrate is globally stretched. As a consequence, in this limit where $\gamma/Ga\to\infty$, the Neumann tension balance at the contact line is strictly valid. In the scheme of energy minimisation, Neumann’s law emerges as auxiliary condition (19), and as such serves as a first boundary condition at the contact line. We have no explanation why previous continuum simulations suggested a deviation from Neumann’s law Wu et al. (2018); Masurel et al. (2019). We emphasise, however, that the present numerics are based on an adaptive method, which was necessary to fully resolve the elastic singularity, and that we extensively verified that the results are fully converged. Furthermore we derived new analytical solutions of nonlinear elasticity that describe the singularity – these are indeed perfectly recovered by the numerics.

How can one understand the deviation to Neumann’s law observed in molecular dynamics simulations of wetting on cross linked polymer networks Liang et al. (2018)? This deviation finds its origin in the lack of scale separation between $\gamma/G$ and the molecular scale $a$, which is inevitable in molecular simulations – the scale $a$ there enters as a molecular cutoff of the continuum and also gives a finite width of the interface. As argued in Style et al. (2017); Andreotti and Snoeijer (2020), the elastic contributions near the contact line can be computed by integrating the elastic stress over a small but finite region – in molecular simulations, the smallest possible size for this region would be $a$. In the present work we have demonstrated that the stress singularity is always logarithmic, $\sigma\sim G\log(r/\frac{\gamma}{G})$. Hence, the integral over stress gives an elastic contribution

which needs to be compared to the surface tensions. In molecular simulations, where typically $\gamma/Ga$ is of order 1 to 100, a measurable elastic correction to Neumann’s law indeed appears. We refer to Andreotti and Snoeijer (2020) for a quantitative test of (38), as the elastic correction to Neumann’s law. However, in typical experiments performed with polymeric gel, where the $\gamma/G$ is well above the micron scale, the scale-separation correction would be $10^{-4}$, and one approaches the macroscopic continuum limit. In such experiments, one safely concludes that Neumann’s law holds.

Second boundary condition versus pinning. While much theoretical work focussed on the validity (or not) of Neumann’s law, very little attention was given to the implications of contact line pinning Andreotti and Snoeijer (2020). In many experiments on soft polymer networks contact line pinning is virtually absent, as quantified by a very small hysteresis Xu et al. (2017); Schulman et al. (2018); Snoeijer et al. (2018); Lhermerout et al. (2016). This implies that the contact line can freely move, exchanging the substrate’s material point touching the liquid-vapour interface without any energetic cost. Here we demonstrated that such a free motion occurs only under a very specific condition. Namely the chemical potential defined by $\mu=\lambda^{2}d\gamma/d\lambda$ must be continuous across the contact line. This is the second boundary condition that needs to be imposed when there is no contact line pinning. Such a condition was previously derived under the restrictive assumption of linear elasticity Snoeijer et al. (2018) – here we demonstrate this to be valid also at large deformation, and explored its consequences in numerical simulations. In particular, we have confirmed numerically that Young’s law is only recovered when the second boundary condition, $\mu_{SV}=\mu_{SL}$, is satisfied at the contact line. For asymmetric surface energies, the second boundary condition in general implies a jump in stretch across the contact line, so that in the presence of a Shuttleworth effect one generically expects large deformations.

The possibility of pinning is interesting in itself. Depending on the material strength, large deformations might lead to fracture, as observed in Bostwick and Daniels (2013), or local plasticity. We demonstrated how such a local “process zone” can be accounted for by introducing a defect in the elastocapillary continuum theory. With the defect, one can accommodate a range of angles $\theta_{L}$ by adjusting the strength of the defect at the contact line. In our simulations we only encountered a weak logarithmic singularity of elastic stress, which implies that the strength of the defect received no contribution from elasticity (i.e. the elastic Eshelby force vanishes). The defect strength, in fact, was found to be equal to the discontinuity in chemical potential at the contact line, i.e. $\Gamma=-[\mu]_{-}^{+}$, giving rise to a modified Young’s law (37). In practice, one would expect the defect to exhibit a “toughness”, just like in fracture, which is the maximum value that can be sustained before depinning occurs. Given (37), we immediately infer that this implies a contact angle hysteresis, with advancing and receding angles $\cos\theta_{r}-\cos\theta_{a}=2|[\mu]_{-}^{+}|_{\rm max}$. Future theoretical work should be dedicated to a more detailed description of the interior of the process zone.

Experiments that probe the strain-dependence of surface tension have so far been based on wetting experiments, with a key role to the contact angles of the solid and of the liquid. Having established the elastocapillary continuum framework for soft wetting, for the first time consistently accounting for large deformations and the Shuttleworth effect, we can now critically assess the experimental situation.

Different series of experiments have been performed with stretched PDMS gels, for static Xu et al. (2017, 2018) and dynamical wetting Snoeijer et al. (2018). They consistently show a change of the solid angle $\theta_{S}$ under stretching. Similarly, the solid angle was found to change in dynamical experiments on PVS van Gorcum et al. (2020). Till date, these experimental observations have not received any other explanation than via a surface tension that depends on the strain (or, on the history of strain). Hence, they offer a convincing case for a nontrivial surface constitutive relation in soft polymer networks, at least for two different systems.

Another direct piece of evidence for the Shuttleworth effect is that experiments in Xu et al. (2017) reveal an increase of stretch at the contact line upon a global stretching of the substrate. This information was previously not used to interpret results in the context of a Shuttleworth effect. However, our numerical and analytical results show that such a variation of the stretch at the contact line can only occur in the presence of elastic Marangoni stresses, induced by a Shuttleworth effect – if surface tension were constant, the stretch at the contact line would take on a constant value.

This evokes an important question that remains to be resolved: What is the microscopic origin of the coupling between surface energy and strain? The polymer is expected to be liquid-like at small scale, where surface tension is exerted: What can produce the coupling between the microscopic scale and the deformation of the network of reticulation (or entanglement) points? A possible scenario is that the coupling emerges from a superficial layer where the mechanical and structural properties are different from bulk Andreotti and Snoeijer (2016); Schulman et al. (2018). Related to this open question is the experimental observation that, in contrast to solid angle $\theta_{S}$, the liquid contact angle $\theta_{L}$ turns out not to depend on stretching Schulman et al. (2018) – a property that was confirmed for 6 different liquid-substrate systems in Schulman et al. (2018), and which also holds for PDMS Xu et al. (2017); Snoeijer et al. (2018), for substrates stretched up to 100%. This is surprising, since Young’s law for the liquid angle should be valid for sufficiently large drops, but with surface energies $\gamma_{SV}-\gamma_{SL}$ based on the externally imposed stretch $\lambda_{\infty}$ Schulman et al. (2018). This interpretation of Young’s law is confirmed in Sec. V, in an analysis where the large deformation elasticity and the Shuttleworth effect are explicitly accounted for. The implication of the experimental invariance of $\theta_{L}$ (within an experimental resolution of $\pm 1^{\circ}$) is that, for all imposed $\lambda_{\infty}$, the strain dependence $d\gamma/d\lambda$ must be nearly the same on both sides of the contact line. While there is no understanding of the microscopic/mesoscopic origin of the strain dependent surface energy, there is a fortiori no real understanding of this property, observed to be valid for many different pairs of liquid and reticulated polymers.

Another assessment of the Shuttleworth effect makes use of an elastic Wilhelmy plate, where a polymeric wire is partially immersed in a liquid reservoir – allowing one to measure the stretch discontinuity across the contact line. In Chen et al. (2019), it was found that the strain remains very small and no discontinuity was observed – implying once more that $d\gamma/d\lambda$ is equal on both sides. In the initial experiment in Marchand et al. (2012a), conversely, a strong discontinuity of strain was observed at the contact line, implying a jump in $d\gamma/d\lambda$. Given that strains remain very small in these experiments, we can assume that the measured strain reflects the actual strains close to the contact line. Therefore, one can interpret these experiments using the no-pinning condition of Sec. V, i.e. $[\mu]_{-}^{+}=0$, which at small strains implies the continuity of $d\gamma/d\lambda$ across the contact line – in perfect agreement with the observation in Chen et al. (2019). It was argued in Chen et al. (2019) that discontinuous strains could be an artefact due to swelling. As an alternative interpretation, we note that in Marchand et al. (2012a) a strong contact angle hysteresis was observed, which in the Shuttleworth-interpretation would also be consistent with a breakdown of the no-pinning condition $[\mu]_{-}^{+}=0$.

In conclusion, this research opens the promising perspective of identifying different conditions or different preparation protocols to get, or not, polymer networks with intricate surface properties. The main open question is to understand the microscopic origin of the Shuttleworth effect, which in the present understanding is unambiguously confirmed for at least two different systems. We emphasise that mechanically, none of the experimental observations are in contradiction with the presence of a Shutttleworth effect, in particular since $\theta_{L}$ and the elastic Wilhelmy plate only probe the difference of strain-dependence on either sides of the contact line. By contrast, the independent measurements of both the solid angle and the stretch at the contact line Xu et al. (2017, 2018) cannot be explained by a hyperelastic theory without explicitly accounting for a strong Shuttleworth effect. Future experiments on a broad class of soft materials should therefore simultaneously explore both contact angles $\theta_{L}$ and $\theta_{S}$, as well as the strains near the contact line. Combined with the fully nonlinear numerics as presented here, this will offer a systematic quantification of the capillarity of soft solids. A next step is to extend the numerical method to a ridge travelling at constant velocity, including the substrate’s bulk viscoelasticity, and possibly history-dependent surface rheology van Gorcum et al. (2018).

Acknowledgments.  We thank J. Eggers, M. van Gorcum, T. Salez and R. Style for discussions. We acknowledge financial support from ERC (European Research Council) Consolidator Grant number 616918, (to A.P. and S.K.), from ANR (French National Agency for Research) grant SMART (to B.A.), from NWO through VICI Grant No. 680-47-632 (to J.H.S.) and an Industrial Partnership Program (a joint research program of Canon Production Printing, Eindhoven University of Technology, University of Twente, and NWO (to E.H.B.).