How short-range adhesion slows down crack closure and contact formation Open Access

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¹¹ or caused by elastic multistability^12–14 and whether their effects are deemed additive¹⁰ or inseparable.⁸ When assessing contact hysteresis, it is as important to describe contact formation as contact rupture.¹⁵ 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.¹⁵ Making matters worse, the errors in energy hysteresis disappear only with the inverse cube of the linear mesh-element size.¹⁵ 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.

The remainder of this work is organized as follows: The methods used are laid out in Sec. II, results are presented in Sec. III, and conclusions are drawn in Sec. IV. Some scaling arguments allowing the Tabor parameter to be generalized to adhesive power law indenters are given in the  Appendix.

We use Green’s function molecular dynamics (GFMD),²⁷ 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.

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)²⁸ 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,³¹ 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.

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.

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.

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 $μ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³³ and found the same $τ∝μ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.

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², 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⁶. 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₀ 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.¹⁴ 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.⁸ 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.

The authors have no conflicts to disclose.

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).

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

While the contact mechanics of both short-range³⁰ and finite-range³² 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.³⁴ 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.