A Green’s Function Molecular Dynamics Approach to the Mechanical Contact between Thin Elastic Sheets and Randomly Rough Surfaces

Adhesion of biological systems is often made possible through thin elastic layers, such as human skin. To address the question of when a layer is sufficiently thin to become adhesive, we extended Green’s function molecular dynamics (GFMD) to account for the finite thickness of an elastic body that is supported by a fluid foundation. We observed that thin layers can much better accommodate rough counterfaces than thick structures. As a result, the contact area is enlarged, in particular, when the width of the layer w approaches or even falls below the short-wavelength cutoff λs of the surface spectra. In the latter case, the proportionality coefficient between area and load scales is (w/λs)3, which is consistent with Persson’s contact mechanics theory.


Introduction
A central aspect in many engineering applications is the mechanical contact between nominally flat but microscopically rough surfaces. It affects important issues such as friction and wear between solids, energy efficiency, and durability as well as environmental compatibility of many devices. Much progress has been made on the topic in the past few decades. Many analytical, numerical, and experimental approaches have been pursued to determine what occurs, from the macro-to the microscales, when two surfaces come into contact. Most of the attention so far has been focused on surfaces resting on bulk solids, and generally described with a semi-infinite. This is in contrast to many biological surfaces, which often have a hierarchical structure like gecko's feet or consist of thin, elastic layers resting on a fluid-like foundation, such as human skin.
An analysis of the contact mechanics of thin layers may thus constitute a first step towards a better understanding of adhesion in biological systems. Indeed, adhesive devices reproducing, on an industrial scale, the amazing adhesive properties of natural systems, are much sought-after. Examples for the latter span a wide range of scales [1][2][3][4][5][6][7][8][9] from insects such as beetles and flies, to spiders and lizards such as geckos. Independent of its systemic position and even the scale of the animal, we find surprising similarities related to the shape of the terminal elements (setal tips) in all these cases always have a similar thin shape as shown in Figure 1. Basically, all the spatulae are flattened and stuck on a rough substrate, while the free ends are oriented in the distal direction [8]. In order to understand the origin of the super-adhesive properties of these systems and, later, to reproduce them, a variety of studies, starting from the pioneering contributions of Kendall in the sixties [9], have focused on the peeling mechanism, by which a tape detaches from the flat substrate where it is attached [3][4][5][6][7][8]10,11]. Despite the interesting results obtained in these studies, the majority rely on the assumption of perfectly flat substrates. However, in order to reproduce good adhesive properties at macroscales, surface roughness cannot be ignored. This study takes a first step towards including the effect of microscale roughness into the contact mechanics of thin layers. As depicted in Figure 1, our model is a simple, yet explicative model for the contact between a rigid rough surface and an elastic layer of finite thickness supported by a constant pressure. It is important to observe that solving the contact mechanics problem between rough solids is challenging even when using the half-space assumption, which treats any characteristic length in the interface as small compared to the thickness of the contacting bodies. Historically, the so-called multiasperity models [12][13][14][15][16] have been the first attempt to resolve the contact between elastically deformable and microscopically rough surfaces: fundamentally, these methodologies consider the rough surfaces to be composed of asperities-with a certain distribution of heights and radii of curvature-which behave like independent Hertzian punches without any mutual interaction. To include the effect of elastic deformation on all scales, Persson developed a different approach [17,18] based on the assumption that the contact pressure probability distribution is governed by a diffusive process as the magnification at which we observe the interface is increased. The theory is formulated in such a way that it is exact in full-contact conditions, while for partial contact, it provides an approximate solution predicting correctly, for example, a linearity between contact area and load (see, e.g., in [19]) and the logarithmic dependence of the interfacial separation on load (see, e.g., in [20]). Furthermore, a large variety of approaches have been developed to obtain quantitatively and qualitatively accurate results: these techniques include finite element methods (FEM) [21], boundary element methodologies (BEM) [22][23][24][25][26][27][28], molecular dynamics simulations [29][30][31] and hybrid approaches [32,33]. An important issue that these approaches have to deal with is their numerical convergence, which can be quite demanding due to the multi-scale nature of the problem.
Here, we extend the well-known and rather efficient Green's function molecular dynamics (GFMD) technique [31] to the treatment of thin layers. The work is outlined as follows: Section 2 describes the mathematical formulation of the methodology, while results are presented and discussed in Section 3. Final remarks conclude the paper in Section 4.

Materials and Methods
We simulate a flat, elastic manifold of finite width that touches on one side a rigid, rough substrate fixed in space and, on the other side, a fluid reservoir, or, alternatively, a solid with negligible shear modulus in such a way that a uniform pressure is applied to the elastic layer (as shown in Figure 1). This is done using the GFMD method [31,34,35] as described in [36]; with Green's functions [37] appropriate for our problem. The central aspects of the approach are described as follows. We denote the position of surface points facing the substrate by their lateral and normal coordinates, r i , which are equally spaced on a two-dimensional surface, representing N = n × n grid points, n being typically 4096. In-plane periodic boundary conditions are employed and the system is treated as being homogeneous within the plane. The coordinate system is chosen with the highest point of the substrate set to zero. The (undeformed) elastic manifold touches the substrate in just this one point, when we apply an infinitesimally small force. The mean displacement at this load is said to be zero. Contact is complete when all points of the manifold touch the substrate. The displacement at that point is normalized to one.
The grid points are propagated in time according to Newton's equation of motion, which is achieved with the Verlet algorithm. The total normal force, F i , acting on a grid point is the sum of an external force, the elastic force, a damping force, and the (unknown) constraint force preventing the rigid surface from penetrating the thin layer. The external force is simply the external/fluid pressure σ 0 times the area element ∆A = a 2 associated with the grid point. The elastic restoring on each atom is computed from the Fourier transform of the stress in Fourier space. The latter is obtained as where the coefficientsz(q) are the Fourier transform of the grid point normal to the surface positions z(r i ) and q are in-plane wavevectors. For the given system, marked by a width w , the Fourier transform of the Green's function is where E * is the contact modulus of the elastic manifold and a finite-width correction factor [37] of where t = qw is a dimensionless width. We note in passing that f (t) tends to one for t 1 and scales as t −3 for t 1. This is, indeed, consistent with the equation governing thin plates (see [38,39] for more details).
The damping force is linear in particle velocity, or, when dynamics are solved in Fourier space, linear in the velocity of a given mode. The proportionality prefactor is best chosen such that the slowest mode of the system is slightly underdamped. For medium loads, it often turns out to be proportional to the load. The time step should be chosen as large as possible before dynamics become unstable, since we only attempt to identify mechanical equilibrium rather than true dynamics. The constraint is established by setting each material point exactly onto the substrate if z i < h(r i ) after the Verlet time step, where h(r) denotes the height of the substrate. It is generated from its Fourier transform, which are complex random numbers with the following expectation values Real and imaginary components ofh(q) are drawn independently from each other. To avoid rare but potential freak fluctuations, we do not draw from Gaussian but instead a uniform random number with zero mean and a second moment according to the target spectrum. The surface spectra satisfy where q r = 2π/λ r and q s = 2π/λ s are the roll-off and the cut-off wavenumber, respectively. As proposed in [40], under the small-slope approximation method, we add an adhesive stress, which acts normal to the surfaces and, in real space, is equal to: where γ 0 is the energy gained per unit area when the two solids touch mechanically, z(r) − h(r) is the local gap between substrate and top solid, and ρ is the characteristic range of the attraction. We chose ρ sufficiently small so that adhesion can be characterized as short ranged. The implementation of the approach is summarized schematically in Figure 2. The GFMD method can handle a very large amount of elements on a single node. This exceeds the size of most experimental surface topography measurements. Equilibration can be achieved within a few thousand time steps for most loads in most cases. The method only becomes inefficient for relative contact areas of less than 0.1% , in which case several hundred thousand iterations might be needed. Here, we only solve for static equilibrium. However, if the method was extended to dynamics, the accessible time span T acc would depend on the level of discretization. As a rough rule of thumb, T acc would be a few µs times the distance between two grid points expressed in units of nanometers.

Results and Discussion
Results reported in this study relate to the contact between an elastic layer and a rigid fractal surface numerically generated by means of the spectral method as described in [31]. In particular, here we have generated fractal surfaces with λ s = 4a for the short wavelength cut-off, λ r = 64λ s for the roll-off wavelength, and a system length of 8λ r . The Hurst roughness exponent H is equal to 0.8. Furthermore, we scale the surfaces in such a way that the root mean square gradient √ < ∇h 2 > is equal to 1.
Firstly, we focuse our attention on the contact mechanics between the rigid surface and an elastic layer with a negligible value of the surface adhesive energy γ: at this stage, our aim is to isolate the effect of the layer thickness and, specifically, how the contact properties change in comparison with thick substrates where the half-space assumption is valid. In Figure 3, we plot the relative contact area as a function of the dimensionless load σ 0 /(E * √ < ∇h 2 >) with E * being the composite Young modulus E * = 1 Pa and for different values of the dimensionless width W = w/λ s , where w and λ s are the layer width and the short wavelength cut-off, respectively. We observe that, by decreasing the ratio W, the contact compliance is increased: a normal load is fixed and a larger contact area is obtained. Indeed, as established by Persson's theory [17,18] and generally accepted by the scientific community ( [23,41]), for an half-space, the area/load relation is described by the equation < ∇h 2 > with k being a proportionality coefficient that, when adhesion can be neglected, is approximately equal to k = 2 [19,31]. When dealing with thin layers, we observe an increase of such a coefficient of several orders of magnitude.  In more detail, in order to better understand the role played by the thickness, in Figure 4 in a log-log plot, we report, again for different values of W, the quantity A(A 0 − A)/A 2 0 as a function of the load, thus showing that the influence of a finite value of the layer width appears also with very low loads. Indeed, this entails a marked difference of the system under investigation (i.e., an elastic layer sustained by a constant pressure), not only in comparison with the classic half-space regime but also with other types of boundary conditions and, in particular, with the case reported in [42], where the deformable layer is bonded to a rigid half-space. In this last case, at low values of the contact area, the system behaves like an half-space and, only when the contact load is increased, we obtain a transition to a different regime where the thickness acquires an important role. This is consistent with the cut-off effect produced by a corrective coefficient Θ(qw), the main effect of which is to remove the low frequency contribution to the contact solution [42]. In contrast, under the conditions we investigate here, also at very low loads, the system shows a marked different trend and, in particular, the contact compliance results increased. This is due to the different mathematical form of the corrective coefficient f (qw) and may have important implications in many systems, including biological membranes and human skin, where very large contact areas may be produced with relatively low loads.  To appreciate the quantitative importance of a correct estimation of the finite thickness effects, in Figure 5a, we collapse the curves representing A(A 0 − A)/A 2 0 into a master curve coinciding with the half-space solution by dividing the load for a corrective coefficient c that will be, as expected, a function of the dimensionless width W. Consequently, for this kind of system, we can obtain the following generalized area/load relation: Indeed, beyond some differences for small contact area, the curves show a good overlap in Figure 5a. At the same time, we observe that the corrective coefficient c has a variation range spanning several orders of magnitude (see Figure 5b): this dramatically illuminates how important it is to account for the actual width of the layer. Furthermore, if we observe the asymptotic trend for very small values of W, we notice that c (W) ∼ W −3 as expected given the form of the Green's function coefficient f (qw) .
The increased compliance of the system can be noticed also when looking at the displacement averaged in the contact area as a function of the load. In detail, as shown in Figure 6, at relatively high loads, in agreement with many theoretical and numerical predictions [17,19,43], a logarithmic dependence between the quantitys = (1− < u > /u max ) and σ 0 /E * is found; as expected when running numerical simulation, at smaller loads, such a logarithmic trend is lost due the finiteness of the rigid surface employed in the computations. Indeed, we notice that, when the thickness W is reduced, fixed the load σ 0 /E * , we obtain a much smaller separations , consistently with a system that is more compliant.  Results shown so far refer to a case with negligible adhesion; however, the contact solution also shows a similar trend when adhesive effects are accounted for. Indeed, such a behavior is directly correlated with the elastic energy stored in the deformable layer. In our analysis, the dominant contribution to the system deformation is assumed to be due to the bending and, indeed, for W << 1, the mechanics of our layers coincides with that of thin plates [39]. The convention implies that any contribution related to the stretching energy is neglected: such an approximation is fully justified in a large variety of conditions, including, for example, the case of adhesive contact of insects and geckos [38].
Elastic energy stored in the layer is, indeed, a crucial quantity that determines the adhesive performance of the structure. More specifically, we define the parameter θ as θ = U el /U ad with U el and U ad equal to the elastic and to the adhesion energy, respectively. This is a qualitative measure of the competition between the elastic energy, which has to be stored to deform the contacting bodies, and the adhesive term, which corresponds to the bodies inclination to come into contact and create contact areas: the smaller θ is, the easier it is to obtain large contacts results. Consequently, θ can be considered a good way to estimate the adhesion capability. In Figure 7, for γ = 5 10 −4 and µ = 2, and for three different values of the contact area, θ dramatically decreases with W . This demonstrates that, indeed, thin layers are preferred to thick structures when good adhesion properties are desired due to the reduced compliance. Interestingly, as expected, for small values of W, θ scales again as ∼ W 3 . It is readily argued that the results we have just obtained are consistent with Persson's contact mechanics theory [17]. Persson predicts that the contact area follows er f (σ 0 /∆σ √ 2), where σ 0 is the apparent contact pressure perceived macroscopically, while ∆σ can be interpreted as the width of the pressure distribution within the microscopically fully resolved contact points. This relation does not depend on the thickness of the elastic manifold. However, the layer width would affect the degree of broadening of the pressure distribution within the contact. Specifically, according to Persson theory, adding ∆C(q) to the height spectrum at wavevector q leads to an increase of ∆σ 2 that is proportional to G(q)∆C(q). Adding up all roughness would then lead to a value of ∆σ 2 scaling proportionally to (t 3 E * ) 2 ∇h 2 (see also, for example, Equation (31) in [36]). For self-affine rough surfaces with a Hurst exponent 0 < H < 1, the mean-square surface-height gradient is dominated by the short-wavelength terms. For truly thin substrates resting on a fluid foundation, the integrated broadening ∆σ therefore scales as (t/λ s ) 3 times the Green's function of the semi-infinite manifold, as we descrived in Section 2, and as is revealed in our numerical results.

Conclusions
In the present study, we develop a numerical methodology based on the GFMD to study the contact between a rigid fractal surface and an elastic layer sustained by a uniformly applied pressure.
Indeed, such a scheme provides a model for a large variety of biological systems, spanning from the human skin to biological membranes. All these cases can be studied by the GFMD method that has been modified multiplying the Green's function G(q) for a corrective parameter f that is related to the dimensionless width t = qw . Our results show that, in comparison with solids with large values of the thickness, where the half space approximation may be adopted, thin layers show a much larger compliance that marks both the contact area and the separation when analyzed as functions of the load. In particular, when looking at the contact area, we observe that the classic area/load relation can be corrected by introducing a coefficient c . Interestingly, c scales as W −3 with W being the dimensionless width of the layer: this is related to the mathematical form of the the Green's function corrective coefficient f (qw) and is consistent with Persson's contact theory. Ultimately, we show why thin layers should be preferred when good adhesive properties and large contact areas are desired.