While viscoelastic, adhesive contact rupture of simple indenters is well studied, contact formation has received much less attention. Here, we present simulations of the formation of contact between various power law indenters and an adhesive, viscoelastic foundation. For all investigated indenters, we find that the macroscopic relaxation time *τ* scales approximately with 1/*ρ*^{1.8}, where *ρ* is the range of adhesion. The prolongation of contact formation with Tabor parameter is rationalized by the increased dissipation that short-range adhesion causes on a moving crack.

## I. INTRODUCTION

The modeling of mechanical contacts between elastomers and rigid counterfaces has advanced substantially in the recent past.^{1–3} In particular, the comparison between simulated contact configurations and real-laboratory images has reached new levels, mainly thanks to improved capabilities of computing and measuring detailed features in partial contacts of nominally flat surfaces.^{1,4–8} However, the progress made pertains mostly to situations where either the experiment is conducted quasi-statically or multi-scale roughness is insignificant. Viscoelastic simulations with roughness on separate length scales remain scarce.^{8–10} As a consequence, it remains difficult to ascertain or to predict when contact hysteresis in a given system is mostly viscoelastic in nature^{11} or caused by elastic multistability^{12–14} and whether their effects are deemed additive^{10} or inseparable.^{8} When assessing contact hysteresis, it is as important to describe contact formation as contact rupture.^{15} However, despite continuous progress in the simulation of viscoelastic, adhesive contacts,^{16–20} contact formation is explored rather little although it is similarly important as contact rupture. It could be one reason for the so-called Monday morning problem, which refers to the sticking of valves in production engines after resting over the weekend. Further applications, where (slow) contact formation matters, are hydraulic and pneumatic seals or adhesive gripping devices.

The difficulty of simulating adhesion-driven contact formation involving soft matter is that the short-range nature of adhesion must be accounted for, even for macroscopic objects. First, noncontact puts much greater demands on the range of adhesion than contact:^{15,21} While a Tabor parameter *μ*_{T} ≈ 5 (*μ*_{T} is a dimensionless measure inversely proportional to the range of adhesion and introduced in detail further below) suffices to reproduce the *μ*_{T} → *∞* load–displacement curves under retraction, the energy loss would be less than 50% of the real value, due to a premature jump into contact. Unfortunately, short-range adhesion requires a modeler to use small mesh sizes to avoid discretization artifacts, most notably lattice trapping.^{15} Making matters worse, the errors in energy hysteresis disappear only with the inverse cube of the linear mesh-element size.^{15} Second, short-range adhesion enhances the dissipation of a moving crack or contact line,^{22–26} whereby not only crack opening but also crack closure is impeded, which in turn puts large demands on the computing time.

The discussion above implies that simulating effectively viscoelastic processes using relatively coarse scales might be achievable by reinterpreting the time scales used in the viscoelastic model.This would be possible if multiplying the range of adhesion by a factor *s* accelerated the dynamics by a power of *s*. Testing for the possibility of such a mapping was the original motivation for the research reported in this work. However, the simulations can also serve as a test for theoretical predictions on contact closure and extend the number of geometries, for which gap closure is investigated. In this context, our focus is on indenter shapes, where the height of the indenter is a power law of the distance from the indenter’s symmetry axis. This problem class offers the greatest potential for us to exploit the similarity of solutions.

## II. METHODS

We use Green’s function molecular dynamics (GFMD),^{27} which is a boundary-element method for the simulation of linearly (visco)elastic contact problems assuming elastic bodies to be isotropic in planes normal to their (originally flat) surface and to be periodically repeated in that plane. We focus on normal, frictionless contacts within linear elasticity. Our systems consist of various perfectly rigid indenters and a homogeneous and isotropic, linearly viscoelastic body described in the continuum limit.

*R*has the unit of length and corresponds to a radius of curvature for a Hertzian (

*n*= 2) indenter, while

*r*is the distance from the indenter’s symmetry axis. For a conical indenter (

*n*= 1),

*h*(

*r*) could be written as

*h*(

*r*) =

*r*tan

*φ*, in order to implement different opening angles

*φ*. Investigated exponents are

*n*= 1 for a conical indenter,

*n*= 2 for a Hertzian indenter, and

*n*= 3 as well as

*n*= 4 as crude approximations for a flat punch. A true flat punch has no interesting contact formation dynamics, as its instantaneous contact radius is identical to the relaxed one.

^{15}It is given by an interaction potential density of

*g*refers to the gap or interfacial displacement between elastomer and indenter,

*γ*is the energy gained per unit area when indenter and elastomer touch, while

*ρ*is the range of adhesion. In all investigated systems, a numerical value of

*γ*= 0.01

*E**

*R*is assigned. The numerical value of

*ρ*was adjusted according to the individual Tabor parameters. In principle,

*ρ*could be set to a smaller value for negative than for positive gaps to better mimic hard-wall repulsion. However, since all reported simulations addressed rather short-ranged adhesion and thus small

*ρ*, we used the same value for

*ρ*for attraction and repulsion.

**q**is an in-plane wave vector, or wave number for one-dimensional interfaces,

*q*is its magnitude, while

*E** is the static contact modulus and $u\u0303(q)$ the Fourier transform of the displacement field. However, as we study a viscoelastic system,

*E** is made frequency dependent. Our default model consists of a single Maxwell element in parallel with a spring, yielding a frequency-dependent contact modulus of

*τ*

_{Mxw}, unless stated otherwise. Any deviation from our default model is explicitly pointed out. For example, in one case, the surface modes are coupled to three rather than to one Maxwell element. The real and imaginary parts of

*E**(

*ω*) are shown in Fig. 1 for the three-element model and the default, one-element model. The “weights” $E0(n)$ and relaxation times $\tau Mxw(n)$ in the three-element models are chosen as $E0n+1=6E0n$ and $\tau Mxw(n+1)=\tau Mxw(n)/6$.

While a “brute-force” Fourier method is certainly not the most effective approach to address the questions in this paper, it is sufficiently efficient to obtain satisfying answers. Moreover, by addressing primarily line contacts having a formal interfacial dimension of *D* = 1, much of the computational overhead spent on not exploiting the symmetry of the problem in *D* = 2 is alleviated.

The most critical dimensionless parameter for a (force-free) adhesive contact of a Hertzian indenter is the Tabor parameter *μ*_{T}. It is inversely proportional to the range of adhesion and is designed such that *μ*_{T} ≫ 1 makes a (quasi-static) contact behave as in the limit for zero-range adhesion, which was first solved analytically by Johnson, Kendall, and Roberts (JKR)^{28} for *n* = 2 and *D* = 2. For arbitrary *n* at *D* = 2, the contact problem can be solved using Sneddon’s method.^{29,30} Less relevant for most real purposes, but easier to tackle numerically and theoretically, is long-range adhesion. For *D* = *n*, adhesion effectively acts like an external load that does not depend on *ρ* for *ρ* → *∞* or *μ*_{T} → 0 when *γ* and *R* are fixed, as is readily seen using the Bradley model,^{31} which becomes exact in the long-range limit. Its analysis also reveals that a physically meaningful long-range-adhesion limit does not exist for *n* ≠ *D*, because the offset load vanishes (*n* > *D*) or diverges (*n* < *D*) for a diverging *ρ* at fixed *γ* and *R*.

^{32}generalized the Tabor parameter for arbitrary power law indenters to

*D*= 1. In the following, we will usually report

*μ*

_{T}rather than the range of adhesion

*ρ*. Note that the precise value for

*μ*

_{T}would differ if the ratio

*γ*/

*σ*

_{max}were used instead of

*ρ*in Eq. (5).

*r*

_{c}(

*t*), which we define to be the distance from the indenter’s symmetry axis to the point where the tensile stress takes its maximum, specifically, the distance in a direction parallel to a main axis of the simulation cell rather than its diagonal, though the two measures are very close to each other and thus show identical scaling. A small relaxation run is done first, using “regular” rather than viscoelastic GFMD, during which the elastomer is allowed to relax to the shape that it would have under the assumption that $E\u221e*$ rather than $E0*$ was its static contact modulus. The estimates of quasi-static contact properties are obtained from similar calculations with the proper static contact modulus $E0*$. This way, initial contact radii

*r*

_{in}≡

*r*

_{c}(

*t*= 0) and quasi-static contact radii

*r*

_{qs}≡

*r*

_{c}(

*t*→

*∞*) can be determined with high accuracy. Results for

*r*

_{c}(

*t*) are reported in terms of a relaxation function defined as

While we do not yet have experimental reference data, it might be important to stress that it might be difficult to rigorously define *r*_{in} in real experiments, since contact with a rigid counter body cannot be simply switched on as in a computer simulation, while inertial effects would prevent an immediate contact at *t* = 0^{+} with the high-frequency modulus from occurring. Thus, we expect deviations between our idealized model and real experiments in the very early stages of contact formation to be unavoidable.

## III. RESULTS

*t*≲ 10

*τ*

_{Mxw}, a regime where the contact radius depends approximately logarithmically on time, i.e., for times 10 ≲

*t*/

*τ*

_{Mxw}≲ 10

^{4}, and a final regime at

*t*≳ 10

^{4}

*τ*

_{Mxw}, which has an exponential time dependence according to

*τ*. While the relaxation in the intermediate time curve is somewhat flattened for the three-element model compared to the one-element model, the contrast in different

*r*

_{c}(

*t*) curves in Fig. 2 appears somewhat less noticeable to the eye than those of the frequency-dependence of the elastic modulus shown in Fig. 1. Interestingly, though perhaps not surprisingly, the early-time (later-time) dynamics of the two models superimpose reasonably well when the relaxation times of the three-element model are scaled such that the high-frequency (low-frequency)

*E*″(

*ω*) or |

*E*(

*ω*)| exhibits similar behavior. In the following, we will merely consider the one-element model, widely known also as the standard linear solid (SLS).

Similar crack closure dynamics as those just discussed are maintained when extending the range of Tabor parameters. To this end, *C*(*t*) rather than *r*_{c}(*t*)/*R* is shown for different *μ*_{T}, one time for two-dimensional interfaces in Fig. 3 and one time, as before, for one-dimensional interfaces in Fig. 4. Similar behavior is revealed in both figures. A double logarithmic representation was chosen for them as to highlight differences in the functional form of *C*(*t*) during the early-time dynamics. Note that the time is expressed in units of *τ*, which allows us to reveal that *C*(*t*/*τ*) is quite insensitive to *μ*_{T} at sufficiently large *t*. The collapse of curves when represented as *C*(*t*/*τ*) improves with increasingly large *μ*_{T}. Factors and related information helpful to reconstruct the full *r*_{c}(*t*/*τ*) are collected in Table I, however, each time only for the largest Tabor parameter.

D
. | n
. | r_{in}
. | r_{qs}
. | τ_{1/2}
. | τ
. | μ_{T}
. | s . |
---|---|---|---|---|---|---|---|

1 | 1 | 0.0021 | 0.062 | 130 | 1060 | 4.1 | 19 |

2 | 0.11 | 0.46 | 23 | 910 | 41 | 5.9 | |

3 | 0.28 | 0.67 | 9.8 | 670 | 65 | 4.1 | |

4 | 0.41 | 0.78 | 6.6 | 750 | 79 | 3.5 | |

2 | 1 | 0.013 | 0.10 | 23 | 110 | 0.30 | 24 |

2 | 0.17 | 0.52 | 3.7 | 83 | 3.0 | 7.2 | |

3 | 0.35 | 0.71 | 1.6 | 45 | 4.8 | 4.9 | |

4 | 0.48 | 0.81 | 1.0 | 36 | 5.8 | 4.2 |

D
. | n
. | r_{in}
. | r_{qs}
. | τ_{1/2}
. | τ
. | μ_{T}
. | s . |
---|---|---|---|---|---|---|---|

1 | 1 | 0.0021 | 0.062 | 130 | 1060 | 4.1 | 19 |

2 | 0.11 | 0.46 | 23 | 910 | 41 | 5.9 | |

3 | 0.28 | 0.67 | 9.8 | 670 | 65 | 4.1 | |

4 | 0.41 | 0.78 | 6.6 | 750 | 79 | 3.5 | |

2 | 1 | 0.013 | 0.10 | 23 | 110 | 0.30 | 24 |

2 | 0.17 | 0.52 | 3.7 | 83 | 3.0 | 7.2 | |

3 | 0.35 | 0.71 | 1.6 | 45 | 4.8 | 4.9 | |

4 | 0.48 | 0.81 | 1.0 | 36 | 5.8 | 4.2 |

The increase in relaxation time compared to the times associated with the Maxwell element can certainly be related to the increased dissipation that short-range adhesion causes in a propagating crack.^{22,24,25} Initial dynamics are relatively fast, when the crack-tip radius is large and thus dissipation small^{10,26} but then slows down as the crack-tip radius becomes smaller, while the slopes of the displacement field increase. For the investigated tip shapes, we find the relaxation times to scale approximately with $\mu T1.8$. This is true for the ultimate relaxation time *τ* as well as intermediate relaxation times *τ*_{x} defined by *C*(*τ*_{x}) = *x* with *x* ≳ 0.5. To make the different data sets appear in a rather narrow range, the Tabor parameter was scaled (multiplied) with a scaling variable *s*, while the relaxation time was divided by *s*^{1.8} (Fig. 5).

It remains to be discussed to what degree our results depend on the precise form of the cohesive-zone model (CZM). To address this issue, we repeated some simulations using a Dugdale model^{33} and found the same $\tau \u221d\mu T1.8$ scaling law as when using Eq. (2). Prefactors of the scaling laws cannot be compared directly because the length scale used in the definition of the Tabor parameter can be either the literal range of adhesion *ρ* or the ratio of surface energy and maximum tensile stress. Taking the geometric mean of those two options makes the Dugdale potential relax roughly 15% more slowly than the default CZM with the same Tabor parameter. We attribute the slight increase in relaxation time and thus dissipation by the Dugdale model to the discontinuity of the CZM at the cutoff. It leads to a diverging second derivative, which in turn induces enhanced local surface velocity for a moving crack and thus enhanced velocity gradients in the material and increased total dissipation compared to more smoothly evolving CZMs.

## IV. CONCLUSIONS

In this work, we confirmed numerically that the crack closure dynamics of indenters with power law profile in contact with an adhesive, linearly viscoelastic foundation has a universal behavior, at least at the late stages of the contact formation process when the contact radius has grown to a size that adhesion can be labeled short-ranged. Analysis of the final stages of the contact formation reveals that the relaxation times increase approximately with the inverse range of adhesion to the power of 1.8. Since Tabor parameters can be quite large in practical applications, this means that crack closure can last rather long times. To give an admittedly extreme example, assuming *E** = 10 MPa, *γ* = 40 mJ/m^{2}, *R* = 20 *μ*m, and a range of adhesion typical for Lennard-Jones interactions, say *ρ* = 3.5 Å, yields *μ*_{T} = 195. Thus, using the data for *D* = *n* = 2 in Table I, we obtain *τ*/*τ*_{Mxw} ≈ 83 · (195/3)^{1.8} ≈ 0.15 · 10^{6}. Making the link to real viscoelastic materials is challenging, because they have a broad distribution of relaxation times, which are quite sensitive to temperature and materials composition. This in turn makes it difficult if not impossible to state a generally valid numerical value for *τ*_{Mxw}. However, if a dynamical analysis is known, an upper estimate for *τ*_{Mxw} should be given by the equation *E*″(1/*τ*_{Mxw}) = *E*_{0} on the small-frequency branch of the loss modulus.

Unfortunately, it is not clear to what degree the scaling might be affected by small-scale roughness. In the extreme case, there will be contact line pinning due to elastic multistability, which can be caused by structural heterogeneity.^{14} In this case, our estimates would only imply a crude lower bound for the contact formation time. A hint in this direction comes from previous simulations.^{8} They mimicked successfully dynamical experiments on the adhesion between an elastomer and a flat punch to which (single-sinusoidal) small-scale roughness was added. In that work, we had identified a steeper dependence of relaxation times on the range of adhesion than in the present work, although our old estimate was based on an analysis, which was much less systematic than the one presented here.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**C. Müller**: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (equal); Software (lead); Writing – original draft (supporting); Writing – review & editing (equal). **M. H. Müser**: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: SCALING ANALYSIS AND GENERALIZED TABOR PARAMETER

While the contact mechanics of both short-range^{30} and finite-range^{32} adhesive power law indenters have been solved, it might be beneficial to have simple arguments allowing one to estimate pull-off forces, work of adhesion, and Tabor parameters in simple terms up to prefactors of order unity without having to embark on the high level of complexity pioneered by Maugis.^{34} For this purpose, it is helpful to realize that the crack closure dynamics for a given *E*(*ω*) can only depend on the dimensionless numbers describing the problem, i.e., *n* and *μ*_{T}. In the following, we derive an expression, which, up to constants of order unity, gives the generalized Tabor parameter. To this end, we define *μ*_{T} such that *μ*_{T} ≡ 1, when simple estimates for the stress-standard deviation in the limit of short-range adhesion match the maximum tension for finite-range adhesion. Our treatment also allows us to identify scaling relations for power law indenters. While similar relations can be deduced from existing literature, we believe our derivation to be original while requiring close to the least possible amount of prior background on contact mechanics and graphite, chalk, ink, or toner in order to arrive at generally valid scaling relations.

*U*

_{el}in areal contacts (

*D*= 2) or its line density

*u*

_{el}≡ Δ

*U*

_{el}/Δ

*L*

_{y}in line contacts (

*D*= 1) can only be of the form

*a*is the contact radius, while the $gn,Dshp$ are constants that depend on the exponent

*n*, the interfacial dimension

*D*, and on the shape of the displacement field. The corresponding total surface energy

*U*

_{s}or line density

*u*

_{s}gained on making contact are

*h*(

*a*) with the help of Eq. (1) and minimizing the total energy in the load-free case w.r.t.

*a*yield a zero-load radii of

*F*

_{po}and the work of adhesion

*W*

_{po}in

*D*= 2 and their corresponding line densities using lower-case letters in

*D*= 1 to be estimated. To this end, we assume that the mean contact stress at pull-off scales linearly with the stress-standard deviation at zero load. Thus, $Fpo\u221da02E*g\u0304c$ for

*D*= 2 while $fpo\u221da0E*g\u0304c$ for

*D*= 1 so that

*F*

_{po}∝

*γR*for a regular Hertzian geometry (

*n*= 2,

*D*= 2) or $Fpo\u221dE*\gamma R3$ for a regular flat punch (

*n*→

*∞*,

*D*= 2). The work of adhesion of a power law indenter can only scale as the product of surface energy and (zero-load) contact area so that

*D*= 1 and

*D*= 2 alike is the same when conducted very slowly or very quickly, i.e., when probing it with the high- or the low-frequency modulus, but it would be large at intermediate pull-off velocities, as argued, for example, in Ref. 26.