*3.1. Axisymmetric BEM Formulation*

In the previous section the JKR solution of the Guduru contact problem was briefly summarized. Several variants of the Guduru contact problem have been studied by different authors [20–22,24,28], nevertheless all of them assume full contact within the contact patch (a simply connected contact area). To overcome this limitation an axisymmetric Boundary Element Method (BEM) was developed assuming that the rigid sphere and the wavy halfspace interacts with a Lennard–Jones interaction law (LJ in the following, see Figure 2a)

$$\sigma\left(h\right) = -\frac{8w\_c}{3\varepsilon} \left[ \left(\frac{\varepsilon}{h}\right)^3 - \left(\frac{\varepsilon}{h}\right)^9 \right] \tag{9}$$

where *σ* is the traction (*σ* > 0, when compressive), *h* is the gap and *ε* the equilibrium distance (the maximum tensile stress *σ*<sup>0</sup> = − 16*w<sup>c</sup>* 9 √ 3*ε* takes place at separation *h* = 3 1/6*ε*). BEM contact codes that use the LJ interaction law have been derived previously by several authors to solve contact problems similar to the one tackled here. For example Wu [31] solved the adhesive contact problem between a sphere and a longitudinal wavy surface, while Medina and Dini [32] studied the problem of an adhesive sphere squeezed against a rough substrate. Notice that other BEM solution strategies exist that are based on energy minimization [33,34].

The contact problem considered is equivalent to the case of a "rough" axisymmetric rigid Hertzian indenter squeezed against an elastic halfspace (Figure 2b). The Guduru gap function is written as

$$h = -\Delta + \varepsilon + \frac{r^2}{2R} + A\left(1 - \cos\left(\frac{2\pi r}{\lambda}\right)\right) + u\_z\left(r\right) \tag{10}$$

where ∆ > 0 when the Hertzian profile approaches the halfspace and *u<sup>z</sup>* (*r*) is the deflection of the elastic halfspace (see Figure 2b).

**Figure 2.** (**a**) Lennard–Jones interaction law. (**b**) Equivalent representation of the considered contact problem in the undeformed configuration (grey lines) and after applying a certain remote displacement ∆ (black lines).

For axisymmetric frictionless contact problems [35,36]

$$
\mu\_z\left(r\right) = \frac{1}{E^\*} \int \sigma\left(s\right) \mathcal{G}\left(r, s\right) s ds \tag{11}
$$

where *σ* (*s*) is the pressure distribution, *G* (*r*,*s*) the Kernel function

$$\mathbf{G}\left(r,s\right) = \begin{cases} \frac{4}{\pi r} \mathbf{K}\left(\frac{s}{r}\right), & s < r\\ \frac{4}{\pi s} \mathbf{K}\left(\frac{r}{s}\right), & s > r \end{cases} \tag{12}$$

and *K* (*k*) the complete elliptic integral of the first kind of modulus *k*.

After surveying the Literature on the topic [31–37], we developed an axisymmetric BEM inspired by the works of Greenwood [35] and Feng [36]. Assume the radial domain is discretized with *N* elements, so that we have *M* = *N* + 1 discretization points. To solve Equation (10) on a discrete domain, one needs to determine the elastic deflection *u<sup>z</sup>* (*r*). A problem arises in evaluating the integral (11) as the kernel function *G* (*r*,*s*) is singular in *s* = *r*. The common approach is to discretize Equation (11) assuming that the pressure *σ* (*s*) has a simple form over a discrete element. To this end, the simplest approach is to assume that the pressure is constant over each element. Nevertheless, Greenwood [35] reported that this method may lead to suspicious results, particularly in the regions with strong pressure gradients and suggested using the method of the overlapping triangles [29], for which the pressure *σ* (*s*) has a triangular form. Hence the deflection at point *r<sup>i</sup>* due to a triangular pressure distribution being *p<sup>j</sup>* at *r* = *r<sup>j</sup>* and falling linearly to 0 at *r* = *rj*−<sup>1</sup> and *r* = *rj*+<sup>1</sup> is

$$
\mu\_z \left( r\_i \right) = \mu\_{z,i} = \frac{1}{E^\*} G\_{ij} p\_j \tag{13}
$$

where we have solved numerically the integral in Equation (11) to obtain *Gij* once for all. Notice that the kernel function singularity at *r<sup>i</sup>* = *rj*−<sup>1</sup> and *r<sup>i</sup>* = *rj*+<sup>1</sup> is canceled by the pressure being 0 in *rj*−<sup>1</sup> and *rj*+<sup>1</sup> , instead, for the singular point *r<sup>i</sup>* = *r<sup>j</sup>* , we considered a pressure equal to −*p<sup>j</sup>* at *rj*−<sup>1</sup> and *rj*+<sup>1</sup> rising linearly at 0 at *r<sup>j</sup>* superimposed to a constant pressure ring, equal to *p<sup>j</sup>* , in between the radii *rj*−<sup>1</sup> and *rj*+<sup>1</sup> for which the displacement field is known analytically (see Appendix A). By defining the following quantities

$$H = h/\varepsilon - 1; \qquad u = \frac{r}{\Gamma}; \qquad L = \frac{\lambda}{\Gamma}; \qquad P = \frac{p\mu\varepsilon}{w\_{\varepsilon}}; \qquad \Gamma = \left(\frac{R^2 w\_{\varepsilon}}{E^\*}\right)^{1/3}; \tag{14}$$

Equation (10) is written for the normalized gap *H<sup>i</sup>* (*H* vanishes for *h* = *ε*) at each node *i* as:

$$H\_i = -\Delta^\dagger + \frac{1}{2}\mu u\_i^2 + A^\dagger \left(1 - \cos\left(\frac{2\pi u}{L}\right)\right) + \mu \sum\_{j=1}^N G\_{ij}' P\_j \tag{15}$$

where *µ* = *Rw*<sup>2</sup> *c E*∗2*ε* 3 1/3 is the Tabor parameter, ∆ † = ∆/*ε*, *A* † = *A*/*ε*, *G* ′ *ij* = *Gij*/Γ and

$$P\_{\dot{\jmath}} = -\frac{8}{3}\mu \left[ \frac{1}{\left(H\_{\dot{\jmath}} + 1\right)^{3}} - \frac{1}{\left(H\_{\dot{\jmath}} + 1\right)^{9}} \right] \tag{16}$$

All the results that will be presented below were obtained using *N* = 500 discretization elements with a constant element size, which proved to be sufficient for obtaining converged solutions up to a Tabor parameter *µ* = 5. Lower values of *µ* did not require such a mesh refinement, nevertheless, to avoid confusion, the same discretization was used along all the paper. Unless differently stated, in all the simulations the overall domain length was set equal to the sphere radius *r*max = *R* = 100*ε*. The numerical code was implemented in MATLAB. Similarly to Feng [36] an efficient pseudo-archlength continuation scheme was implemented [38], which is needed to follow the

system solution branches. Furthermore, to make the numerical solution of the nonlinear system of Equation (15) faster to solve, the system Jacobian was provided analytically to the numerical solver ("*fsolve*" in MATLAB) that implements a Newton–Raphson scheme.

#### *3.2. Validation of the Numerical Results*

First to assess the correctness of the numerical implementation, the BEM results are validated against those provided by Feng [36]. In Figure 3 the dimensionless normal load

$$\overline{\mathcal{W}} = \frac{\mathcal{W}}{\pi w\_{\varepsilon} \mathcal{R}} = \frac{2\Gamma^2}{\mu \mathcal{R} \varepsilon} \int\_0^{+\infty} P(u)u du}{}$$

is plotted as a function of the dimensionless approach −<sup>∆</sup> † for *µ* = [1, 2, 3]. Markers refer to the data reported by Feng [36] (its Figure 3) while the solid black lines were obtained numerically using the BEM presented here. Red dots stand for the Bradley [23] rigid solution, which is compared to the numerical solution (solid black line) obtained with *µ* = 0.01. All the curves are in perfect agreement.

**Figure 3.** Dimensionless normal load *<sup>W</sup>* as a function of the approach −<sup>∆</sup> † for *µ* = [1, 2, 3] as reported by Feng [36] (markers) and as obtained here numerically (solid black line). Red dots show the Bradley [23] rigid solution, which is compared to the numerical solution (solid black line) obtained with *µ* = 0.01.

Figure 4 shows the pull-off force *W pull*−*o f f* (panel (a)) and the approach <sup>−</sup><sup>∆</sup> † *pull*−*o f f* (panel (b)) at pull-off as a function of *<sup>µ</sup>* (the pull-off force *<sup>W</sup>pull*−*o f f* is defined as the minimum of the *<sup>W</sup>*(<sup>∆</sup> † ) curve). Blue stars have been obtained from Feng [36], while red squares have been obtained using the BEM developed here. The results we obtained for both load and approach at pull-off compare very favorably with those in [36].

**Figure 4.** (**a**) Pull-off force *W pull*−*o f f* and (**b**) approach <sup>−</sup> <sup>∆</sup> † *pull*−*o f f* at pull-off as a function of the Tabor parameter *µ* = [10−<sup>2</sup> , 5]. Red squares were obtained with the presented Boundary Element Method (BEM) scheme, while blue asterisks were obtained from Feng [36].

#### **4. Rigid Limit**

The majority of the authors, who have tackled the Guduru contact problem, have focused their attention on the JKR limit, where it was clear since the early papers by Guduru [17] and Guduru and Bull [19] that substantial enhancement could be obtained, with few exceptions, as the works of Waters [22] and that of Ciavarella [24]. Waters et al. [22] developed a Maugis–Dugdale cohesive model for the Guduru problem that showed adhesion enhancement is mostly limited to the JKR regime. The cohesive model clearly depended on an additional parameter with respect to the JKR model, the Tabor parameter *µ*. Strictly speaking, Waters et al. [22] used the parameter introduced by Maugis [39], which anyway differs only by a small multiplicative factor from *µ*. Waters et al. [22] analysis showed that for small *µ* the pull-off detachment force converged to the Bradley rigid solution for the sphere, i.e., *<sup>W</sup><sup>B</sup>* <sup>=</sup> 2. Nevertheless, this holds only for a smooth sphere in contact with a flat halfspace. For example, let us assume *λ* ≈ *R* ≈ *A*, then the pull-off force in the rigid regime could be estimated by considering the contact between the first crest of the halfspace waviness and the sphere. The radius of curvature of the crest is

$$R\_2 = \frac{\lambda^2}{4\pi^2 A} \approx \frac{R}{4\pi^2} \tag{18}$$

and the composite radius

$$R^\* = \left(\frac{1}{R} + \frac{4\pi^2}{R}\right)^{-1} \approx \frac{R}{4\pi^2} \tag{19}$$

hence the pull-off force of the sphere would be reduced by about factor 1/4*π* <sup>2</sup> <sup>∼</sup> 0.025. Indeed, using the Rumpf–Rabinowicz model [25–27] Ciavarella [24] recognized this. Although analogous to the Guduru problem, the Rumpf–Rabinowicz model refers to a different geometry. It considers the contact between a large sphere of radius *R* and a rigid small hemisphere of radius *r*<sup>2</sup> placed on a rigid plane (see Figure 5).

**Figure 5.** Left panel: simplified sketch of the Rumpf–Rabinowicz model. Right panel: simplified sketch of the Guduru model.

Two competitive mechanisms for adhesion take place: while the radius of the hemisphere increases macroscopic adhesion increases due to the interaction with the hemisphere but decreases as the rigid plane gets further away from the countersurface, which, using *r*<sup>2</sup> = *<sup>λ</sup>* 2 4*π*2*A* , can be written as [24]

$$\left|\overline{W}\right|\_{pull-off} = \frac{1}{1 + 4\pi^2 a} + \frac{1}{\left(1 + \frac{R^\dagger}{4\pi^2 a}\right)^2} \tag{20}$$

where *R* † = *R*/*ε*. Similar mechanisms are expected to be at play in the Guduru problem.

Assume that the rigid sphere and the wavy halfspace interact with a Lennard–Jones interaction law (9). From Equation (15), neglecting the elastic deformations, the dimensionless interfacial gap is

$$H(\theta) = -\Delta^{\dagger} + A^{\dagger} \left[ \frac{\theta^2}{8\pi^2 a} + (1 - \cos\left(\theta\right)) \right] \tag{21}$$

where *θ* = <sup>2</sup>*π<sup>r</sup> λ* has been introduced. Using Equations (16) and (17) the total load is

$$\overline{\mathcal{W}}\_{rigid} = \frac{2\pi}{\pi w\_{\varepsilon}R} \int\_{0}^{+\infty} \sigma(H) r dr = -\frac{4A^{\dagger}}{3\pi^{2}a} \int\_{0}^{+\infty} \left[ \left( \frac{1}{H(\theta) + 1} \right)^{3} - \left( \frac{1}{H(\theta) + 1} \right)^{9} \right] \theta d\theta \tag{22}$$

which clarifies that at a given approach ∆ † the rigid solution depends only on two parameters *α*, *A* † .

Notice that for a smooth sphere-flat contact one can use in the first integral of Equation (22) *dH* = *<sup>r</sup> R dr* and obtain that the rigid solution depends only on the adhesion energy and not on the shape of the potential. This cannot be done for the Guduru gap function, which implies that the Guduru rigid solution will be slightly affected by the shape of the interaction law used.

In Figure 6 the loading curves are shown for log<sup>10</sup> *α* = [−4, −2, −1, 0, 1] and *A* † = 1. The solid black lines are the theoretical predictions based on Equation (22), while the red markers have been obtained numerically using the BEM with *µ* = 10−<sup>4</sup> . Numerical and theoretical predictions are in perfect agreement, with the curve log<sup>10</sup> *α* = −4 corresponding to the Bradley [23] rigid solution for the smooth sphere (Equation (22))

$$\overline{\mathcal{W}}\_B\left(\Delta^\dagger\right) = -2\left(\frac{4\left(\Delta^\dagger - 1\right)^6 - 1}{\Im\left(\Delta^\dagger - 1\right)^8}\right) \tag{23}$$

The curves dimensionless normal load *<sup>W</sup>* versus dimensionless approach −<sup>∆</sup> † plotted in Figure 6 show that the pull-off force is not monotonically decreasing with *α* and that the critical approach at detachment is close to ∆ † <sup>≃</sup> 0 only for the smallest values of *<sup>α</sup>*. Figure 7 shows the pull-off force as a function of *α* for log<sup>10</sup> *A* † = [−3, <sup>−</sup>0.5, 0, 1, 2, 4] (solid lines obtained as the minimum of (22), while red markers obtained numerically *µ* = 10−<sup>4</sup> ). All the curves start at small *α* from the Bradley rigid solution *WB*,*pull*−*o f f* <sup>=</sup> 2 (smooth sphere), while increasing *<sup>α</sup>* the pull-off force decays and, after a transition zone around *α* ≈ 1, reaches a limit value at large *α*. It is shown that at large *α*

and *A* † the pull-off force can decrease by more than three orders of magnitude with respect to the smooth case. To allow a comparison the predictions of the Rumpf–Rabinowicz model are reported for log<sup>10</sup> *R* † = [0, 1, 2, 3, 4] (blue dashed lines in Figure 7). First notice that, while the rigid limit of the Guduru problem depends on *α*, *A* † , the Rumpf–Rabinowicz model depends on *α*, *R* † . Both the models show a similar decay with *α*, but they give two different limits for large *α*. In the Rumpf–Rabinowicz model, a large *α* implies a very small hemisphere (*R*/*r*<sup>2</sup> = 4*π* <sup>2</sup>*α*), hence the case of a large sphere interacting with a smooth plane is retrieved. In the rigid Guduru model, *A* † and *λ* † = *λ*/*ε* are not coupled, hence increasing *α* leads to a vanishing wavelength of the sinusoid but does not affect *A* † , which gives the observed adhesion reduction.

**Figure 6.** Loading curves for the rigid model. Solid lines have been obtained from Equation (22), while red markers are numerical solutions for *µ* = 10−<sup>4</sup> . The curves are obtained for log<sup>10</sup> *α* = [−4, −2, −1, 0, 1] and *A* † = 1.

**Figure 7.** Pull-off force as a function of the parameter *α* (log scale) obtained using the rigid solution of the Guduru model (Equation (22), solid black line) and the Rumpf–Rabinowicz model (Equation (20), blue dashed line), while red squares are numerical solutions for *µ* = 10−<sup>4</sup> . For the Guduru rigid model the curves are obtained for log<sup>10</sup> *A* † = [−3, <sup>−</sup>0.5, 0, 1, 2, 4], while for the Rumpf–Rabinowicz model log10(*R* † ) = [0, 1, 2, 3, 4].
