Effect of Surface Roughness on Adhesive Instabilities for the Elastic Layer

If an incompressible elastic layer bonded to a rigid plane is placed close to a second rigid plane, adhesive interactions between the surfaces can cause elastic instabilities. These lead to spatially non-uniform traction and gap distributions which exhibit a regular pattern with a characteristic wavenumber. However, real surfaces are never completely plane. In this paper, we consider the influence of surface roughness on the instability, with particular reference to the force-displacement relation. With random roughness profiles, the traction distributions are always spatially irregular, so the onset of instability is more difficult to define. One approach is to monitor the amplitude of the power spectrum of the distribution near the characteristic wavenumber. Since surface roughness generally reduces the mean adhesive traction, we might expect it to exert a stabilizing effect. Numerical results confirm this for moderate to large RMS amplitudes, but show that low RMS roughness can actually trigger the instability in ranges where the uniform layer would be stable. The resulting traction-displacement relation is then found to be approximately linear with a slope close to that at the point where the uniform solution loses stability.


INTRODUCTION
If two bodies with plane surfaces are placed close together, they may experience attractive [e.g., van der Waals'] forces, or forces involving both attractive and repulsive ranges (Jones, 1924;Maugis, 2013). Since the attractive forces must eventually decay with increasing separation, they have the character of a "negative spring, " which can trigger an elastic instability. If the bodies are incompressible [Poisson's ratio ν = 0.5], or if a body comprising a thin elastic layer bonded to a rigid plane surface is attracted to another rigid plane surface, the instability may result in a nonuniform [typically periodic] pattern of alternating regions of contact and separation. Patterns of this kind have been observed experimentally (Mönch and Herminghaus, 2001;Gonuguntla et al., 2006a), and predicted theoretically, based on energetic arguments (Shenoy and Sharma, 2001;Sarkar et al., 2004). In particular, the characteristic length scale of the pattern correlates with the unstable wavelength in a linear perturbation of the uniform state. The patterning instability also modifies the mean traction-separation characteristic for the layered system, generally leading to different behavior during approach and separation and consequent hysteresis losses (Ciavarella et al., 2019).
The instability permits self-assembly processes such as elastic contact lithography [ECL], where the pattern in a polymer film is fixed by UV curing or by lowering the temperature (Sarkar and Sharma, 2010;Ghosh et al., 2016Ghosh et al., , 2017. In ECL, the periodicity and size of the pattern are critical parameters and various methods have been proposed to control them, including the use of a curved FIGURE 1 | A rigid body with a plane surface that may contain some surface roughness is placed near to an elastic layer bonded to a rigid foundation. The mean gap isḡ and u(ξ , η) is the local elastic displacement of the layer surface.
Real surfaces are of course never perfectly smooth, and surface roughness generally reduces both the maximum pulloff traction, and the maximum negative slope of the effective traction-separation law (Joe et al., 2018). This should reduce the tendency for patterning in the contact of layered bodies. However, sufficiently small amplitude roughness might also serve as an initial perturbation to trigger an instability. In this paper, we shall therefore use a numerical solution to examine the effect of roughness on both the generation of patterns and the mean traction-separation relation. In particular, we shall examine the extent to which the effect of roughness can be captured by using a modified adhesive traction law developed for the contact of rough elastic half spaces.

DEFORMATION OF A THIN LAYER
We consider an elastic layer of thickness h bonded to a rigid plane, and define Cartesian coordinates (x, y) in the plane of the layer surface and corresponding dimensionless coordinates ξ = x/h, η = y/h. If a second rigid plane is placed a distancē g away from the undeformed surface of the layer as shown in Figure 1, the local gap between the surfaces will then be where u(ξ , η) is the local outward normal elastic displacement of the layer surface.

Interface Energy
We assume that the adhesive tractions between the surfaces can be described by a traction law σ (g), where g is the local value of the gap. We can then also define the mean interface energy per unit area as and a stable final configuration will be one that minimizes the total potential energy = U + Ŵ, where U is the mean elastic strain energy per unit area.

Elastic Strain Energy
The normal traction σ (ξ , η) at the free surface of the layer needed to produce a sinusoidal normal elastic displacement u(ξ , η) = u ζ cos(ζ ξ ) is where (Hannah, 1951), and we recall that ξ = x/h, so ζ is a dimensionless wavenumber.
In this paper, we shall restrict attention to incompressible layers [ν = 0.5], for which the corresponding dimensionless compliance 1/f (ζ ) is shown as a function of dimensionless wavenumber ζ in Figure 2. The curve exhibits a maximum of ∼0.482 at a wavenumber ζ ≈ 2.1, and zero compliance for uniform loading f (ζ ) → ∞ as ζ → 0]. One consequence of this is that with general loading σ (ξ , η), the mean value of u is zero and henceḡ is determined by the rigid-body approach which is a controlled parameter.
The mean elastic strain energy per unit area associated with the deformation (3) is The elastic strain energy for more general displacement distributions can be obtained by writing u(ξ , η) as a double Fourier series or Fourier transform and convoluting the resulting transform with (5).

Stability Criterion
If both surfaces are plane [i.e., smooth], the state u(ξ , η) = 0, g(ξ , η) =ḡ, σ (ξ , η) = σ (ḡ) is clearly an equilibrium state, but it will be unstable to small sinusoidal perturbations of dimensionless wavenumber ζ if there exists any ζ such that The critical wavenumber is defined by the maximum of the curve in Figure 2, from which we deduce that the uniform solution will be unstable if and only if −σ ′ (ḡ) > E/0.482h.  (6), wavenumbers in the range ζ A < ζ < ζ B are unstable.}.

Solution Method
More general displacement distributions for a square domain 0 < x < L, 0 < y < L can be written in the form and equation (5) can then be used to obtain the mean strain energy per unit area U as Adding the corresponding interface energy from equation (2), we obtain the total potential energy , which must be a minimum at a stable equilibrium state. We used the gradient descent method to identify the values of the Fourier coefficients A mn for a local energy minimum. Notice that the upper limit N must be chosen so as to provide an adequate number of integration points in the evaluation of Ŵ.

RESULTS FOR SMOOTH SURFACES
If the two surfaces are smooth, the uniform state is always an equilibrium solution at which energy gradients are zero, and since the material is incompressible [ν = 0.5], this corresponds to u(ξ , η) = 0. Even when the criterion (6) is satisfied, the numerical solution may remain at the uniform state unless some small perturbation is introduced.
We used the traction law (Maugis, 2013) derived from Lennard-Jones molecular force law (Jones, 1924), where ε is an interatomic length scale and γ = γ (ε) is the interface energy per unit volume at the equilibrium spacing g = ε. The maximum tensile traction occurs at g = 3 1/6 ε and is of magnitude σ 0 ≈ 1.06 γ ε.
If the maximum negative slope [−σ ′ (g)] max satisfies the condition there will exist a bounded range g 1 <ḡ < g 2 in which the uniform solution is unstable. Also, for a value ofḡ strictly within this range, there will exist a range of unstable wavenumbers ζ A < ζ < ζ B , including but not limited to ζ ≈ 2.1. This range can be identified by drawing a horizontal line at the height in Figure 2 as shown, and finding its intersection with the curve. We start the solution procedure with a value ofḡ outside the unstable range and then changeḡ by small increments, using the solution at the previous step as an initial guess for the gradient descent solution. This is expected to mimic the behavior of the physical system under controlled displacement conditions. Numerical noise might also be expected to emulate the effect of noise [e.g., vibration] in an experimental system.
We characterize the inverse thickness of the layer by the dimensionless parameter size was chosen so as to define a fundamental wavenumber [e.g., ζ 01 ] equal to 0.25, so that L = 8πh. Results are shown for both separation dḡ/dt > 0 and approach dḡ/dt < 0. In each case, the uniform state is preserved up to the appropriate stability boundary [denoted by A and B, respectively], but the non-uniform solution then persists significantly beyond the point at which the uniform solution reverts to stability. We deduce that even in the stable range there exist local energy minima corresponding to non-uniform states, and that these states are separated from the lower energy uniform state by energy barriers. Notice that both the approach and separation curves in Figure 3 are approximately straight lines with slope close to the critical slope defined by (6).
The non-uniform deformation states are characterized by the development of regular patterns. Figures 4A-C shows contours of local gap g(x, y) corresponding to the points on the approach curve labeled (a), (b) and (c) in Figure 3. A "labyrinth" [i.e., a connected system of passageways (high g) separated by walls of "contact" (low g)] develops at the onset of instability (a) and Frontiers in Mechanical Engineering | www.frontiersin.org then morphs into an inverse labyrinth at (b) and into an array of isolated regions of separation ["holes"] at (c).
Figures 5A-C shows the same results in Fourier transform space. In all cases we see high values clustered near the most unstable wavenumber ζ = 2.1 and the distribution is axisymmetric within the limits of statistical variance, indicating that the pattern is statistically isotropic.

EFFECT OF SURFACE ROUGHNESS
We assume that the surface is rough with a power spectral density [PSD] of the form where B is a constant, H is the Hurst exponent, here taken as H = 0.2, and ζ 1 , ζ 2 define the range of dimensionless wavenumbers in where u 0 (x, y) describes the deviation of the surface from the mean plane in the undeformed state, and The magnitudes of the coefficients |B mn | were chosen so as to ensure that the resulting surface PSD was of the form (13) and the corresponding arguments [phases] were chosen randomly. With this definition, the gap can be written where u(x, y) is given by (7). Notice that we arbitrarily assign the value B 00 = 0, so that the roughness makes no contribution to the mean gap. The interface energy is then determined from (2) and the constants A mn are chosen so as to minimize as in the smooth surface case.

A Two-Scale Approximation
We anticipate that the wavelengths in the roughness spectrum will generally be significantly smaller than the layer thickness, and this suggests the possibility of a scale-separation approach. Compared with the roughness scale, the thickness of the layer is large, so local effects can be approximated by those in a corresponding half space. The effect of the roughness, as compared with a corresponding smooth surface, can therefore be described in terms of a modified traction law. An inductive method for estimating this law is described in (Joe et al., 2018). If the modified traction law σ M (g, ζ 0 ) is known for a surface with spectral content only in ζ 0 < ζ < ζ 2 , where ζ 0 > ζ 1 , the corresponding law for σ M (g, ζ 0 −δζ ) can be obtained by convoluting σ M (g, ζ 0 ) with the additional roughness tranche ζ 0 − δζ < ζ < ζ 0 . Successive applications of this technique allow us to determine the modified law for the entire spectrum. In Joe et al. (2018), this procedure was implemented using discrete tranches of the spectrum, but the same approach can be used to develop a partial differential equation for σ M (g, ζ ), following the methodology of Persson (2001). The modified traction law for the complete roughness spectrum is then defined by σ M (g) = σ M (g, ζ 1 ).
On the scale of the layer thickness, the effect of the surface roughness can then be approximated by using σ M (g) in place of equation (1), and treating the layer surface as smooth.

Results
The two-scale approximation is likely to be most accurate when the most unstable wavenumber ζ ≈ 2.1 is much smaller than the smallest wavenumber ζ 1 in the PSD [recall that wavenumbers are normalized by the layer thickness h]. However, this degree of scale separation is difficult to achieve in a direct numerical simulation. Figure 6 shows the relation between mean traction and mean gap for an incompressible elastic half space with roughness of the form (13) with ζ 1 = 6.5, ζ 2 = 8.0 and height variance Frontiers in Mechanical Engineering | www.frontiersin.org Barber (2018). In this figure, we compare two methods of solution: (i) a direct numerical solution of the problem [solid line] by minimizing the total energy for a gap defined by (16) using the Lennard-Jones traction law and (ii) a solution in which the layer is smooth [u 0 (x, y) = 0], but the traction law is modified to describe the effect of roughness [circles]. The agreement is clearly extremely good. The advantage of using the modified traction law is that the deformation of the layer can then be adequately described on a much coarser grid, or equivalently, by a more severely truncated Fourier series, and this is particularly useful if the spectral range [ζ 1 , ζ 2 ] is broad.

Identifying Instability Effects
The argument of the previous section suggests that pattern instabilities should occur if and only if the maximum slope of the modified traction law [σ M (ḡ)] satisfies the condition (11) -i.e., The modified traction law σ M (ḡ) is independent of h, so instabilities are more likely for thicker layers, always assuming the finite linear dimensions of the interface are sufficient to accommodate a wavelength in the unstable range.
In some cases, pattern instabilities can be detected by examining the corresponding contour plots for g(x, y). For example, Figure 7A shows the contour plot for the actual rough surface during separation at a mean gapḡ = 1.4, and Figure 7B shows the corresponding pattern predicted for a smooth surface with the modified traction law. There is some blurring of the observed patterns, but the overall morphology is clearly similar.
However, if the roughness amplitude is larger, patterns become more blurred, a typical example being shown in Figure 8A. In this case, detection of instability from gap contours is more difficult, but the corresponding Fourier plot of Figure 8B clearly shows substantial spectral content in the unstable range near ζ = 2.1. This suggests that we might quantify the extent of pattern formation in the numerical solution for the rough surface using the dimensionless parameter m u 0 /ε 2 , where m u 0 is the variance of that part of the gap PSD that lies in the unstable range in the "smooth" solution, identified in Figure 2. We assume here that the roughness PSD has no content in this wavenumber range, since otherwise it would be difficult to distinguish the separate effects of instability and roughness.
In Figure 9 we plot [m u 0 /ε 2 ] max , obtained from the numerical solution, as a function of the maximum negative slope of the modified traction law [−σ ′ M (ḡ)] max , normalized by the corresponding expression for the Lennard-Jones law. Each set of points corresponds to a different value of the lower wavenumber of the roughness ζ 1 = 4.5 : 0.5 : 6.5 and a range of roughness variances −2.7 < log 10 (m 0 /ε 2 ) < −1.6. Results were obtained under both approach and separation conditions, but no significant difference was observed. The vertical dashed line in Figure 9 corresponds to the criterion (18), below which the twoscale approximation would predict m u 0 = 0, though the direct numerical results exhibit a level of noise as one might expect. However, the results exhibit a remarkable level of consistency, showing that [−σ ′ M (ḡ)] max is a very good indicator of the effect of roughness on the instability, and more generally that the two-scale approach to the layer problem defines a good approximation to important features of the system behavior.

CONCLUSIONS
If a smooth elastic layer is placed close to a plane surface, elastic instabilities due to adhesive tractions lead to the development of patterns, and to a modification of the traction-separation law. However, roughness with RMS amplitude comparable with the range of the adhesive force law can have a significant effect on this process. Here we have described a model for analyzing the contact of both rough and smooth surfaces using a double Fourier series.
We also developed a two-scale approximation to the rough contact behavior, by (i) estimating the effect of roughness on the mean traction between two half spaces, using a previously published method (Joe et al., 2018), and then (ii) using this modified traction law in the analysis of a smooth elastic layer. Results show that this gives a very good approximation to the traction-separation law obtained by direct numerical simulation. In particular, the development of patterns is predicted if the maximum slope of the modified traction law satisfies the inequality (18) and the corresponding results correlate extremely well with a criterion based on the spectral content in the unstable range from the numerical solution.
Local layer deformations decay spatially at a rate linked to the layer thickness, so this method is expected to give good predictions for bodies of finite size sufficient to support wavelengths in the unstable range.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/ supplementary files.

AUTHOR CONTRIBUTIONS
JJ wrote the first draft of the manuscript, performed the statistical analysis and organized the database. JB wrote sections of the manuscript. All authors contributed to conception and design of the study, manuscript revision, read, and approved the submitted version.