On the application of fracture mechanics mixed-mode models of sliding with friction and adhesion

In fracture mechanics, it is well known that mixed mode enhances the toughness observed in pure mode I (Evans and Hutchinson 1989) due to crack faces interlocking and friction resulting for roughness. Unfortunately, these mixed-mode models are not physical laws or general energy principles, but are intrinsically empirical. They are mostly of the form including a mode-mixity function $f\left( \psi\right)$ (Hutchinson and Suo 1992) giving the 'effective toughness' G_c,eff as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{c,ef\!\,f}=G_{Ic}f\left( \psi\right) \label{Gpsi} \nonumber \end{align} \tag{ 1 }$$

where G_Ic is mode I toughness (or surface energy, if we assume Griffith's concept). Also, $\psi$ is phase angle

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi=\arctan\left( \frac{K_{II}}{K_{I}}\right) \nonumber \end{align} \tag{ 2 }$$

where K_II and K_I are, respectively, the mode II and mode I stress intensity factors. This is strictly valid only in the case of Dundurs' second constant equal to zero (see (Barber 2018) or Hutchinson and Suo, p 77 versus p 82). Cao and Evans (1989) experimentally looked at epoxy-glass bimaterial interface, and in general various models for microscopic phenomena affect the interface toughness, such as friction, plasticity and dislocation emission (Hutchinson 1990).

When models like these are applied to a contact mechanics problem in the presence of adhesion and friction, we may expect either the contact area to be largely unaffected by the presence of a mode II loading, in one limit or that $G_{c}\simeq G_{Ic}$, and in this case we effectively expect mode II weakens the mode I condition, so the contact area should decrease in sliding with respect to the pure adhesion case. The case of area enhancement is rather unexpected, as experimentally it is confirmed (Ciavarella 2018, Papangelo et al 2019, Papangelo and Ciavarella 2019, Sahli et al 2019).

Fracture mechanics concepts were firstly applied by Johnson et al (1971) (JKR-theory, 1971) to adhesion between elastic bodies. They are applicable to contact even in the presence of friction, as mixed-mode fracture mechanics problem, but with some special peculiarities. Specifically, the energetic 'JKR-assumptions' correspond to the Griffith criterion, which consists in practice in assuming extremely short range adhesive forces (virtually a delta-function). JKR is the correct limit for soft and large bodies, and the equivalent of the so called 'small-yield' criterion is expressed by the Tabor parameter (Tabor 1977) for the sphere,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_{sphere}=\left( \frac{RG_{Ic}^{2}}{E^{\ast2}\Delta r^{3}}\right) ^{1/3}=\frac{\left( Rl_{a}^{2}\right) ^{1/3}}{\Delta r}=\frac{\sigma_{0} }{E^{\ast}}\left( \frac{R}{l_{a}}\right) ^{1/3}\label{Tabor} \nonumber \end{align} \tag{ 3 }$$

being very large for JKR to apply, where R is the sphere radius, G_Ic is work of adhesion, $\Delta r$ is the range of attraction of adhesive forces, and $E^{\ast}$ the plane strain elastic modulus. $E^{\ast}=\left( \frac {1-\nu_{1}^{2}}{E_{1}}+\frac{1-\nu_{2}^{2}}{E_{2}}\right) ^{-1}$ and E_i, $\nu_{i}$ are the Young modulus and Poisson ratio of the material couple. Also, $\sigma_{0}$ is the theoretical strength of the material, and we have introduced the length $l_{a}=G_{Ic}/E^{\ast}$ as an alternative measure of adhesion—for Lennard-Jones potential of elastic crystals, $l_{a} \simeq0.05a_{0}$, where a₀ is the equilibrium distance, which means that l_a is of the order of angstroms.

The use of energetic criteria extending JKR to the presence of friction was attempted by Savkoor and Briggs (1977) who also conducted experiments between glass and rubber. We write in terms of Irwin stress intensity factors

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{I}=\frac{P_{a}}{2a\sqrt{\pi a}}=\frac{p_{a}}{2}\sqrt{\pi a};\qquad K_{II}=\frac{T}{2a\sqrt{\pi a}}=\frac{\tau_{m}}{2}\sqrt{\pi a}\label{Ki} \nonumber \end{align} \tag{ 4 }$$

where the external normal load $P=P_{H}-P_{a}$ while $P_{H}=\frac{4E^{\ast }a^{3}}{3R}$ is a compressive Hertzian load and P_a is an adhesive load responsible of the contact edge singularity. Writing the energy balance condition in terms of a constant tangential load for which $K_{II}=Y\tau _{m}\sqrt{\pi a}$, where $\tau_{m}$ is the average shear while Y is a geometric factor which contains also some averaging over the perimeter, also of the K_III term, while a is the radius of the circular contact area, Savkoor and Briggs's results in a reduced effective energy for the 'ideally brittle fracture' at equilibrium

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{c}=G_{Ic}=G_{I}+G_{II}=\frac{1}{8}\frac{p_{a}^{2}\pi a}{E^{\ast}} -Y^{2}\frac{\tau_{m}^{2}\pi a}{2E^{\ast}}=\left( \frac{p_{a}^{2}}{8E^{\ast} }-Y^{2}\frac{\tau_{m}^{2}}{2E^{\ast}}\right) \pi a \nonumber \end{align} \tag{ 5 }$$

and therefore the contact area will follow the JKR equation, but following the equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{Ic,ef\!\,f}=G_{Ic}-Y^{2}\frac{\tau_{m}^{2}\pi a}{2E^{\ast}}\label{Savkoor} \nonumber \end{align} \tag{ 6 }$$

where $\tau_{m}$ is the friction average stress. Experiments clearly evidenced a reduction of the contact area when tangential load was applied, but less than expected from assuming $G_{c}\simeq G_{Ic}$, so less than the prediction (6). Experimental findings showed some interference with the development of Schallamach wave which tend to permit slip without affecting the contact.

More recent experiments continue to confirm contact area reduction at both macroscopic and even smaller scales (Sahli et al 2019). Johnson (1996, 1997) and Waters and Guduru (2009) have proposed different models to take into account the interplay between two fracture modes, namely I and II. In particular, Johnson (1997) attempts also to model slip explicitly with cohesive models (as well as the mode I corresponding part), and even in this case, the conclusion remains of the contact area reduction. Various recent other papers (Ciavarella 2018, Papangelo et al 2019, Papangelo and Ciavarella 2019, Sahli et al 2019) have shown that the size and even the elliptical shape of the contact under shear are reasonably found by these LEFM (linear elastic fracture mechanics) models over a wide range of loads and geometrical features, despite the mixed-mode function strictly requires a complex functional form to replicate faithfully the results and hence empirical fitting at least over one set of results. Also, they suggested there is no obvious advantage in trying to model the slip displacements (which correspond to recur to a cohesive model, in the context of fracture mechanics), since this effect is essentially included in the mixed-mode function. Obviously, in the empirical functions models, one could use different empirical functions, and therefore be able to model also enhancements of contact area size.

A recent paper by Menga et al (2018) (MCD, in the following), stemming from some experimental evidence that shear stress should be constant at the interface during sliding at least in rubber versus glass (see (Chateauminois and Fretigny 2008), but also MCD reference list), introduce an interesting and stimulating variant of the friction and adhesion problem suggesting an increase of contact area in sliding—even without the need to postulate an increase of the mixed mode fracture energy. This result would seem in contrast with present experimental evidence on the contact area, but has the advantage to emerge naturally from an apparently thermodynamic rigorous theory, although the range of validity should be discussed. To reconcile the predictions with the experiments, MCD argue that fluctuation in the stresses reduce the effective surface energy, or G_Ic. Notice immediately that assuming the shear stresses are constant at the interface, is equivalent to a fully cohesive model, i.e. a cohesive model like suggested by Dugdale-Barenblatt and Maugis (see (Maugis 2013)), when the size of the cohesive zone fully extends to the entire 'crack' (i.e. contact area). While energy criteria can be still applied with cohesive models, which have been devised to extend LEFM, this is the true limit case which is to be treated with great care.

Indeed, we could define a Tabor parameter for the shear problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_{sphere}^{shear}=\frac{\tau_{0}}{E^{\ast}}\left( \frac{R}{l_{a}^{shear} }\right) ^{1/3} \nonumber \end{align} \tag{ 7 }$$

where $\tau_{0}$ is now a theoretical strength under shear, and in principle, we introduced also a different adhesive length scale $l_{a}^{shear}$. It is hard to estimate $\mu_{sphere}^{shear}$ exactly, but if we assume that $l_{a}^{shear}\simeq l_{a}$, it is possible that $\tau_{0}\ll \sigma_{0}$ and hence as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_{sphere}^{shear}\ll \mu_{sphere} \nonumber \end{align} \tag{ 8 }$$

there is a region for which we can be in the intermediate range for which we can apply essentially LEFM criteria in mode I and a fully cohesive model for mode II, as implicitly, the authors of MCD theory assume. However, this leads to non-obvious results, which pose a number of interesting questions worth examining.

Consider a system in which there is an uncoupled mode I problem (pressures, normal displacements), and a tangential one (shear stresses, shear displacements), as for example in the case of a rigid body against an elastomer with $\nu=0.5$ (or, more generally, Dundurs' second constant equal to zero). We can write the strain energy stored in the body as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U=\frac{1}{2}\int_{A}\sigma vdA+\frac{1}{2}\tau_{0}W \nonumber \end{align} \tag{ 9 }$$

where A the contact area, $\sigma$ is the normal pressure and $v$ the normal component of displacement, W is the displaced volume in tangential direction W  =  Aw where w is the mean tangential displacement which we can write as $w=k\tau_{0}A^{1/2}/E^{\ast}$, where k is a constant factor of the order of 1, which is not important here. Hence, we can write $\tau_{0}=\frac{wE^{\ast} }{kA^{1/2}}$, $W=k\tau_{0}A^{3/2}/E^{\ast}$, W  =  wA, and $\frac{1}{2}\tau _{0}W=\frac{1}{2}\frac{E^{\ast}A^{1/2}}{k}w^{2}$. In the classical formulation of this problem in mixed mode fracture, the state variables would be ($v,w,A$), and one would need to apply a Griffith energy minimum principle, the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial U}{\partial A}\right) _{v,w}=G_{Ic} \nonumber \end{align} \tag{ 10 }$$

with imposed displacements $v,w$. This would not satisfy the requirement that shear stress distribution is constant in the contact area: the solution has the well-known LEFM square root singularities, in both the normal pressure distribution, and also in the shear stresses. Indeed, it is the full stick solution which incidentally is used with success in previous papers (Ciavarella 2018, Papangelo et al 2019, Papangelo and Ciavarella 2019, Sahli et al 2019), even for this problem under transition to sliding, but using 'empirical' mixed-mode functions to take into account of the complex effect of the influence of mode II into mode I.

In MCD, the authors explore what a pure energy condition implies without any mixed-mode empirical function. By applying the Legendre transform to change the state variables from ($v,w,A$) to ($v,\tau_{0},A$), the trivially obtain a new thermodynamic potential H,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H =U-\left( \frac{\partial U}{\partial W}\right) _{v,A}W=U-\left( \frac{\partial U}{\partial w}\frac{\partial w}{\partial W}\right) _{v,A}W=U-\tau_{0}W \nonumber \end{align} \tag{ 11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle =\frac{1}{2}\int_{A}\sigma vdA-\frac{1}{2}\tau_{0}W \nonumber \end{align} \tag{ 12 }$$

as imposing the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial H}{\partial A}\right) _{v,\tau_{0}}=G_{Ic} \nonumber \end{align} \tag{ 13 }$$

leads to the solution under mode I displacement control, and mode II 'strength' control.

This leads to their equation (26) (both under force, or under displacement control in mode I)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{c,ef\!\,f}=G_{Ic}+\Delta G_{Ic}\left( a\right) =G_{Ic}+\frac{4\tau_{0}^{2} a}{\pi E^{\ast}}\quad\label{menga} \nonumber \end{align} \tag{ 14 }$$

i.e. the effective surface energy (or toughness of the interface) is increased, rather than decreased in Savkoor's theory (6), and curiously of a very similar quantity, perhaps even exactly the same since $\tau_{0}=\tau_{m}$, if $Y=\frac{2}{\pi}$. This should not be confused from the result of the theories using the mode-mixity functions of the type (1), like (Ciavarella 2018, Papangelo et al 2019, Papangelo and Ciavarella 2019, Sahli et al 2019), since in the latter case, there is no size-effect associated to the contact area a, as there is in (14), and this has profoundly different implications, as we shall explore.

Indeed, the first reactions to this result are that

• (i)  
  even in the limit of no surface energy $G_{Ic}\rightarrow0$, there would be an 'effective adhesion', as $G_{c,ef\,\!f}\rightarrow\frac{4\tau_{0} ^{2}a}{\pi E^{\ast}}$, which, under sufficiently large compressive normal forces, would imply an unbounded increase of the effective energy;
• (ii)  
  this equation also implies a distortion of the JKR solution load versus area of contact which instead in the AFM experiments by Carpick et al [13], was a nearly perfect fit, even during sliding, leading to the conclusion that the force was simply proportional to the contact area (see figure 1).

Zoom In Zoom Out Reset image size

Therefore, this area enhancement theory is more difficult to reconcile with experiments than models like (Ciavarella 2018, Papangelo et al 2019, Papangelo and Ciavarella 2019, Sahli et al 2019), which of course remains semi-empirical, but also are more robust in the predictions—particularly not having any size effect in the increase of toughness.

The authors of the MCD theory recur to suggesting how large pressure oscillations may compensate for this effect, but these are not an outcome of the theory, and this should be further investigated. However, there seem to be simpler explanations as to why the implied effects are not measured, as we shall explain in the next paragraph with a simpler example, which leads to an even more surprising conclusion.

2.1. A mode I 'nanoscale' paradox?

The MCD theory that the contact area should increase upon application of tangential force may be due to their use of a energy condition for linear elastic fracture mechanics, while being in the limit of fully developed cohesive zones in mode II. If we now speculate of an imaginary adhesive problem in which 'constant tensile pressure $\sigma_{0}$' is observed the way MCD suggest, we have the strain energy is more simply

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U=\frac{1}{2}\sigma_{0}V \nonumber \end{align} \tag{ 15 }$$

where $V$ is the displaced volume $V=Av$, where $v$ is the mean normal displacement which we can write as $v=k\sigma_{0}A^{1/2}/E^{\ast}$, and as usual k is a constant factor of the order of 1, which is not important here. Hence, $\sigma_{0}=\frac{vE^{\ast}}{kA^{1/2}}$, $V=k\sigma_{0}A^{3/2} /E^{\ast}$, $V=vA$, and finally $U=\frac{1}{2}\frac{E^{\ast}}{k}A^{1/2}v^{2}$.

By applying the Legendre transform to change the state variables from ($v,A$) to ($\sigma_{0},A$) in order to have the 'pressure-control', we obtain a new thermodynamic potential H,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H =U-\left( \frac{\partial U}{\partial V}\right) _{v,A}V=\frac{1} {2}\sigma_{0}V-\left( \frac{\partial U}{\partial v}\frac{\partial v}{\partial V}\right) _{v,A}V \nonumber \end{align} \tag{ 16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle =-\frac{1}{2}\sigma_{0}V=-\frac{k}{2}\frac{\sigma_{0}^{2}}{E^{\ast}}A^{3/2} \nonumber \end{align} \tag{ 17 }$$

as indeed $-\sigma_{0}V$ is the potential energy associated with the uniform stress distribution $\sigma_{0}$. Notice that now $v$ is no longer prescribed, but the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial H}{\partial A}\right) _{\sigma_{0}}=G_{Ic} \label{conclusion} \nonumber \end{align} \tag{ 18 }$$

is never satisfied, as the energy decreases without limit and the minimum is clearly at infinite contact area. It is a similar conclusion than MCD theory—to the extreme limits the increase is unbounded.

In fact, this example instead may not be such an ideal situation at all: it is known that at nanoscale, the adhesion problem becomes controlled by the theoretical strength $\sigma_{0}$ (Gao and Yao 2004), and if we have a flat indenter, the correct adhesion pull-off is simply

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{po}=\sigma_{0}A \nonumber \end{align} \tag{ 19 }$$

where A is the size of (flat ended) indenter, a principle which is found in biology of nanoscale fibrillar structures to maximize their adhesion on their feet, below sizes on the order of 100 nm. In the case of an indenter having a shape, including the spherical classical one, one would need very small sizes to reach this limit (Greenwood 2009) unless the shape is the 'optimal one' involving elliptic integrals. In any case, the energetic formulation simplified from MCD does not seem to provide a meaningful result, as (18) seems to imply always infinite contact area, regardless of shape of indenter. It is a limit case whose range of validity we are not able to identify, but suggests a warning also on the more general mixed-mode corresponding result.

MCD assume there is a steady state solution, and probably make the error to consider the potential energy of a constant tangential load and at the same time, impose that shear stress at interface is constant. But in this case, the problem is: there is no equilibrium. We have an imposed tangential load T₀, an imposed shear strength $\tau_{0}$, and the problem is solved in terms of equilibrium in just one line

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{0}=\frac{T_{0}}{\tau_{0}} \nonumber \end{align} \tag{ 20 }$$

independently on anything else. The rest of the energy calculation is in conflict with the former equation in general.

Notice that this is not limited to the frictional problem. Even in the pure mode I problem, the same problem persists. Indeed, there is no equilibrium. We have an imposed normal load P₀, an imposed pressure strength p ₀, and the problem is solved in terms of equilibrium in just one line

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{0}=\frac{P_{0}}{p_{0}} \nonumber \end{align} \tag{ 21 }$$

independently on anything else. The rest of the energy calculation is in conflict with the former equation in general.

We consider the system of figure 2. A sphere is loaded under constant forces, and s_x is the halfspace mean displacement, u_x0 is the relative sphere-halfspace displacement, so that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_{x}=u_{x0}+s_{x}\label{0}. \nonumber \end{align} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

We assume there is a steady state⁵.

The stiffnesses and the entire problem are uncoupled in the $x,y$ directions as the sphere is rigid and the substrate is incompressible, and hence Dundurs' second constant is zero (Barber 2018).

4.1. Energy balance at imposed external forces

The external forces $F_{x},F_{y}$ are imposed and dead loads. The variation of work of active external forces is equal to the increase of internal energy of the system, plus the increase of the dissipated energy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta L=dU+\delta Q\label{1} \nonumber \end{align} \tag{ 23 }$$

where dU is localized in the halfspace as the sphere is rigid. Hence, the various terms are equal to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta L =F_{x}du_{x}+F_{y}du_{y}\label{2} \nonumber \end{align} \tag{ 24 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle dU =d\left[ U_{x}\left( F_{x},A\right) +U_{y}\left( F_{y},A\right) \right] \label{3} \nonumber \end{align} \tag{ 25 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta Q =\delta Q_{fr}+\delta Q_{coes}\label{4} \nonumber \end{align} \tag{ 26 }$$

where the term $\delta Q$ is here the sum of a frictional work $\delta Q_{fr}$ and the work to increase of cohesive area $\delta Q_{coes}$. In general, notice that $\delta Q_{coes}$ depends on the path to internal equilibrium of the adhesive stratum, which in this treatment is external to the thermodynamical system

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta Q_{coes}=-R_{C}\left( \frac{s_{shear}}{s_{opening}}\right) \delta A \nonumber \end{align} \tag{ 27 }$$

where $s_{shear},s_{opening}$ are, respectively, shear and direct stresses acting on the cohesive zones, and as yet not precisely defined. In the special limit case of path independent (reversible) energy, we return to the classical Griffith JKR Signorini type of contact,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta Q_{coes}=-G_{c}\delta A \nonumber \end{align} \tag{ 28 }$$

where G_c is surface energy.

The equilibrium of the contact area size is obtained when the virtual change with respect to contact area is zero in (23), or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{dU}{dA}-\frac{\partial L}{\partial A}+\frac{\partial Q_{fr}}{\partial A}+\frac{\partial Q_{coes}}{\partial A}=0 \label{5} \nonumber \end{align} \tag{ 29 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{dU}{dA} =\left( \frac{\partial U_{x}}{\partial A}\right) _{F}+\left( \frac{\partial U_{y}}{\partial A}\right) _{F}\label{6} \nonumber \end{align} \tag{ 30 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial L}{\partial A} =F_{x}\left( \frac{\partial u_{x}}{\partial A}\right) _{F}+F_{x}\left( \frac{\partial u_{y}}{\partial A}\right) _{F}\label{7} \nonumber \end{align} \tag{ 31 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial Q_{fr}}{\partial A}\right) =F_{x}\left( \frac{du_{x0}}{dA}\right) _{F}=F_{x}\left( \frac{d\left( u_{x} -s_{x}\right) }{dA}\right) _{F}\label{8} \nonumber \end{align} \tag{ 32 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial Q_{coes}}{\partial A}\right) _{F} =-R_{c} \label{9} \nonumber \end{align} \tag{ 33 }$$

where we used (22).

Hence, substituting (30)–(33) in (29), we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \left( \frac{\partial U_{x}}{\partial A}\right) _{F}+\left( \frac{\partial U_{y}}{\partial A}\right) _{F}-F_{x}\left( \frac{\partial u_{x}}{\partial A}\right) _{F}-F_{y}\left( \frac{\partial u_{y}}{\partial A}\right) _{F}\nonumber \\&+F_{x}\left( \frac{d\left( u_{x}-s_{x}\right) }{dA}\right) _{F} -R_{c}=0 \nonumber \end{align} \tag{ 34 }$$

or, cancelling terms due to the work of frictional forces,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial U_{x}}{\partial A}\right) _{F}+\left( \frac{\partial U_{y}}{\partial A}\right) _{F}-F_{y}\left( \frac{\partial u_{y}}{\partial A}\right) _{F}-F_{x}\left( \frac{ds_{x}}{dA}\right) _{F}-R_{c}=0. \nonumber \end{align*}$$

Notice therefore that the work of frictional forces does not affect the energy balance, since we assume the corresponding term (32) is uncoupled from the term due to cohesion energy (33).

An important consideration is that, being the contact problem uncoupled, i.e. Dundurs' $\beta=0$, normal displacements do not affect shear tractions and viceversa, then the two modes of cohesion are naturally separated. Moreover, local cohesive laws do not distinguish between mode II and mode III shear modes, we shall indicate with G_s the strain energy release rate under shear, whereas G_I is the more standard strain energy release rate under mode I

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{I} =\left[ \left( \frac{\partial U_{y}}{\partial A}\right) _{F}-F_{y}\left( \frac{\partial u_{y}}{\partial A}\right) _{F}\right] \nonumber \end{align} \tag{ 35 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{II} =\left[ \left( \frac{\partial U_{x}}{\partial A}\right) _{F}-F_{x}\left( \frac{\partial u_{x}}{\partial A}\right) _{F}\right]. \nonumber \end{align} \tag{ 36 }$$

Moreover, we shall assume that

• the shear force F_x is, for a given contact area A, a non-linear function of u_x
• we can write the elastic strain energy in the tangential direction in the form
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U_{x}=\frac{1}{2}F_{x}^{2}C_{x}\left( A\right) \nonumber \end{align} \tag{ 37 }$$
  where C_x is the tangential compliance of the system, andwe can write
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s_{x}=C_{x}\left( A\right) F_{x}. \nonumber \end{align} \tag{ 38 }$$

Hence, we can rewrite (33) in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \left( \frac{\partial U_{y}}{\partial A}\right) _{F}-F_{y}\left( \frac{\partial u_{y}}{\partial A}\right) _{F}\right] -\frac{F_{x}^{2}} {2}\left( \frac{dC_{x}}{dA}\right) =R_{c} \nonumber \end{align} \tag{ 39 }$$

which is a quite general result where we do not have to specify the distribution of shear stresses at the interface

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{I}-\frac{F_{x}^{2}}{2}\left( \frac{dC_{x}}{dA}\right) =R_{c}. \nonumber \end{align} \tag{ 40 }$$

In our case, $F_{x}=\tau_{0}A$, and hence, similarly to (5), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{I}-\frac{\tau_{0}^{2}}{2}\left( A^{2}\frac{dC_{x}}{dA}\right) =R_{c}. \nonumber \end{align} \tag{ 41 }$$

In our case, it is obvious that $\left( \frac{dC_{x}}{dA}\right) _{F}<0$, and probably this holds even in much more general cohesive law cases, so that the presence of a tangential force will always decrease the contact area. We do not search here for how general this result is, as this is outside the scope of the present note.

Turning back to the energy balance (23), let us consider the most general mixed boundary condition on the boundaries

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_{k} =\overline{F}_{k},\qquad k=1,...,n_{F} \nonumber \end{align} \tag{ 42 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_{k} =\overline{u}_{k},\qquad k=n_{F+1},...,2 \nonumber \end{align} \tag{ 43 }$$

where $k=x,y$ and n_F is the number of DOFs where we impose the forces.

Hence

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta L =\overline{F}_{k}du_{K},\qquad k=1,...,n_{F} \nonumber \end{align} \tag{ 44 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle dU =dU\left( s_{x},u_{y},A\right) \nonumber \end{align} \tag{ 45 }$$

and the condition for contact area equilibrium is obtained by imposing the virtual variation of (23) satisfying b.c.s, with respect to A, independently on the variation path

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{dU}{dA}-\frac{\partial L}{\partial A}+\frac{\partial Q_{fr} }{\partial A}+\frac{\partial Q_{coes}}{\partial A}\right) _{Path} =0\label{path} \nonumber \end{align} \tag{ 46 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \left( \frac{dU}{dA}\right) _{Path} =\left( \frac{dU}{dA}\right) _{s_{x},u_{y}}+\left( \frac{\partial U}{\partial s_{x}}\right) _{u_{y} ,A}\left( \frac{\partial s_{x}}{dA}\right) _{Path}\nonumber \\ &+\left( \frac{\partial U}{\partial u_{y}}\right) _{s_{x},A}\left( \frac{\partial u_{y}}{dA}\right) _{Path}\label{20} \nonumber \end{align} \tag{ 47 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{dL}{dA}\right) _{Path} =\overline{F}_{k}\left( \frac {du_{K}}{dA}\right) _{Path}\label{21} \nonumber \end{align} \tag{ 48 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{dQ_{fr}}{dA}\right) _{Path} =\overline{F}_{k}\left( \frac{du_{ax}}{dA}\right) _{Path}=\overline{F}_{k}\left( \frac{d\left( u_{ax}-s_{x}\right) }{dA}\right) _{Path}\label{22} \nonumber \end{align} \tag{ 49 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{dQ_{coes}}{dA}\right) _{Path} =R_{c}\label{23}. \nonumber \end{align} \tag{ 50 }$$

Additionally, equilibrium is imposed by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial U}{\partial s_{x}}\right) _{u_{y},A} =\overline {F}_{x}\label{24} \nonumber \end{align} \tag{ 51 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( \frac{\partial U}{\partial u_{y}}\right) _{u_{x},A} =\overline {F}_{y} .\label{25} \nonumber \end{align} \tag{ 52 }$$

Substituting (47)–(52) into (46), we get

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle &\left( \frac{dU}{dA}\right) _{s_{x},u_{y}}+\overline{F}_{k}\left( \frac{ds_{x}}{dA}\right) _{Path}+\overline{F}_{y}\left( \frac{du_{y}} {dA}\right) _{Path}\nonumber \\ &-\overline{F}_{k}\left( \frac{du_{K}}{dA}\right) _{Path}+\overline{F}_{k}\left( \frac{d\left( u_{x}-s_{x}\right) } {dA}\right) _{Path}-R_{c}=0 \nonumber \end{align*}$$

i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G=\left( \frac{dU}{dA}\right) _{s_{x},u_{y}}=R_{c} \nonumber \end{align} \tag{ 53 }$$

where in general with coupled stiffness, we cannot split G into mode I and shear mode components. However, when they are indeed uncoupled as we are assuming in the present MS,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G=\left( \frac{dU_{x}}{dA}\right) _{s_{x}}+\left( \frac{dU_{y}}{dA}\right) _{u_{y}}=G_{I}+G_{s}=R_{c}. \nonumber \end{align} \tag{ 54 }$$

Hence, we have demonstrated that even with generalized boundary conditions and non-linear elasticity, the energy release rate G is independent on the boundary conditions. For uncoupled problems, further

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U\left( s_{x},u_{y},A\right) =U_{x}\left( s_{x},A\right) +U_{y}\left( u_{y},A\right) \nonumber \end{align} \tag{ 55 }$$

and obviously $G_{I},G_{s}>0$ for any opening or shear load.

5.1. Role of toughening mechanisms

To clarify the role of cohesive resistance R_c which itself is a generalization of surface energy G_Ic, we can write (for the uncoupled case)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{I}+G_{S}=R_{c}\left( \frac{s_{shear}}{s_{opening}}\right) =R_{c}\left( \frac{G_{s}}{G_{I}}\right) \nonumber \end{align} \tag{ 56 }$$

where cohesive resistance should be modelled with some surface mechanism which may be material-dependent, or directly from experimental data.

In general, we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{I}=G_{Ic}+\left( R_{c}\left( \frac{G_{s}}{G_{I}}\right) -G_{Ic}\right) -G_{s} \nonumber \end{align} \tag{ 57 }$$

and therefore, adding a tangential load to an equilibrium contact area, leads to a reduction of the contact area if the resistance increases with respect to G_Ic faster than the contribution to mode II does

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( R_{c}\left( \frac{G_{s}}{G_{I}}\right) -G_{Ic}\right) -G_{s}<0. \nonumber \end{align} \tag{ 58 }$$

In general, we cannot rule out that this condition is possible, although ample literature shows that this is very unusual if ever observed ([3–5]).

For path-independent cohesive energy, this condition is satisfied as $\left( R_{c}\left( \frac{G_{s}}{G_{I}}\right) -G_{Ic}\right) =0$ and hence contact area must decrease.

We have suggested that some hidden problems in the otherwise rigorous thermodynamic theory of MCD may be due to the uncertain range of applicability. Let us summarize what a well known cohesive model obtains in the context of contact mechanics, for the mode I problem, namely the Dugdale–Maugis solution for the contact problem of a sphere (Maugis 2013). The idea is to postulate a cohesive zone having an outer constant stress $\sigma_{0}$ in an annulus a  <  r  <  c outside of the contact, where d  =  c  −  a is the size of the cohesive zone. With this 'trick', energetic methods can still be applied to the problem, even beyond the LEFM formulation, since we know the energy release rate of the cohesive zones. However, in the limit of very low Tabor parameter, when the cohesive zone is extremely large, $m=c/a\rightarrow \infty$, is very subtle, since the cohesive stresses are constant in the annulus, but they must be zero, i.e. $\sigma_{0}=0$. A limit solution, the so-called DMT-M solution for the sphere (Maugis's version of the DMT solution, see (Ciavarella 2017)), it does not appear possible to obtain it with the MCD procedure, despite there is no reason why the same Legendre transform idea could be applied in this simpler problem.

Hence, although we share with the authors of MCD theory the fascination for the elegance of the Legendre pure thermodynamic formulation, and despite we really enjoyed it for the number of stimulating discussion it generates, we suggest it is problematic for various reasons:

• mixed mode problems have hardly been solved with simple energy formulations without additional 'empirical' constants and criteria, which is why (Hutchinson 1990) and the other references given in the introduction paragraph, devised them for the problem of mixed mode fracture;
• with respects to the empirical formulations which have found some validation in experiments, it seems that MCD theory leads to a size-effect of the surface energy/toughness which is not just giving an area increase, as MCD noticed in the paper, but also contrasts with the excellent fits of JKR theories done by Carpick et al (1996) and many others;
• there is a problem associated with the concept of an experiment neither under force nor under displacement control, but rather in 'shear stress' control. As we have to decide how we are going to cause the body to move, and the only ways we can postulate a proper problem are (i) pushing at a constant force, (ii) pushing at a constant speed or (iii) an intermediate case where a spring is connected to the body, the end of which is pushed at constant speed, it is hard to imagine an experiment (possibly even a numerical one), that corresponds to the assumed conditions of the thermodynamic Legendre transform theory in MCD's theory;
• a fully cohesive developed zone in mode II, in the spirit of fracture mechanics, corresponds to a very low Tabor parameter in shear, while the energetic treatment can treat at most intermediate Tabor conditions. It is hard to imagine a rigorous case in which the Tabor parameter in shear should be extremely low, while the Tabor parameter in pressure should be very high.
• similarly, trying to predict the contact area changes in a JKR (short -range) theory due to a very long-range adhesion effect under shear appears also possibly a problem.

The thermodynamics treatment of Menga et al (MCD theory) leads to a number of paradoxes, not only in mixed mode conditions, but even just in mode I. The reason was identified in that the energy method obtained by the Legendre transform does not satisfy equilibrium. We have obtained more consistent treatments which show the contact area is always decreasing with tangential load, consistent with experiments. Although experimentalists have indeed been measured directly in soft materials that shear stresses appear constant during sliding, specifically in glass versus rubber, more recent deconvolutions considering the high strain gradients reached (Nguyen et al 2011) start to find deviations from the perfect constant shear distribution; this may be another reason not to insist too much on this model showing a very low Tabor number in mode II while we have a very high Tabor number in mode I.

This work was supported by the Italian Ministry of Education, University and Research under the Programme "Department of Excellence" Legge 232/2016 (Grant No. CUP-D94I18000260001).