Surface patterns created by spinodal decomposition are the product of the decomposition of two phases interrupted at an arbitrary time and the selective discarding of one of the phases. The Cahn-Hilliard equation governs the spatiotemporal diffusion and evolution of phases, which can be quantified by the local concentration ϕ r , t at position r on the surface and time t (Cahn, 1961), where the magnitude of ϕ ranges from 0 to 1 and 0 and 1 denote distinct material phases, Γ is the penalty parameter for interphases, D is the effective diffusivity, and f ϕ is the free energy functional quantified by the total energy of a particular phase distribution in space. ∂ ϕ r , t ∂ t = D ∇ 2 d f ϕ r , t d ϕ r , t − Γ ∇ 2 ϕ r , t (1) While creating the surfaces via spinodal decomposition, the partial differential equation given in Eq. 1 should be numerically solved for ϕ r , t f where t _f is an arbitrarily long time, and then one of the phases, say close to 0 values, are spatially eliminated. In practice, the elimination process involves a selective material removal process such as chemical etching, and thus perfect separation of phases is not possible, especially around the interphases where 0 < ϕ r , t f < 1 . Therefore, this imperfect etching will result in surface patterns. Here, we assume that the height of those surface patterns will scale with the magnitude of ϕ r , t f , i.e., the spots with ϕ r , t f = 1 will constitute the surface peaks thanks to no etching and the spots with ϕ r , t f = 0 will experience maximum etching and thus form the valleys. The bottleneck in the generation of surfaces in this manner is the long simulation times needed for sufficient spatiotemporal evolution of phase distribution, and the sensitivity of the final patterns to the initial values. Cahn had shown that solution to Eq. 1 can be approximated by a superposition of Fourier components growing or diminishing according to their wavelengths (Cahn, 1961). Assuming slowly-varying amplitudes for those components, spatial variation of phase can be approximated as a Gaussian random field (spinodoid) as done elsewhere (Kumar et al., 2020), where n _i and Γ _i are the direction and phase angles of the ith component of the series, k is the wave number that we set to 2π, and m is the number of waves taken as 1,000 in this study. ϕ r = 2 m ∑ i = 1 m cos k n i ⋅ r + Γ i (2) Obviously, the approximate series representation given in Eq. 2 provides the concentration ϕ(r) and thus surface heights much quicker than the full Cahn-Hilliard equation does. Besides, Kumar et al. (2020) showed that limiting the directions of the components n _i would lead to anisotropic patterns, and thus the creation of various topographies is possible with a single expression. To generate the concentrations as a function of position, we first select the direction vector n _i residing on the surface spanned by e ₁ and e ₂ basis vectors in x and y directions. Direction vector has a unit-norm and can be described by a single angle Θ _i = cos⁻¹ (e ₁ ⋅n _i ) as n _i = (cos Θ _i , sin Θ _i ) in e ₁ and e ₂ basis. Thus, the selection of the direction vector is reduced to the choice of the angle Θ. If that angle is chosen randomly within the range [0, π/2], then the direction of wave components will be randomly distributed in the first quadrant defined by x > 0 ∧ y > 0 surface plane and the resultant surface patterns will be isotropic. Limiting the upper bound of that range to a particular angle θ, in contrast, will lead to strong directionality and thus anisotropy in the concentrations and surface patterns. After selecting the direction vector, phase angle Γ _i is randomly assigned from the range [0, 2π] to each wave. Summing the expressions for all waves as in Eq. 2 then yields the concentration ϕ(x, y) over the discretized surface plane. We then scale the concentrations uniformly to obtain a target of 0.1 [a.u.] for the root-mean-square (rms) of surface heights; i.e., h (x, y) = Cϕ(x, y) such that h _rms = 0.1. The surface plane xy is 10 × 10 [a.u.] in dimension so that surface slopes ≈ 0.1 / 10 are practically sound. Since k = 2π yields wavelengths of unity in each wave component, 10 × 10 surface plane will guarantee about 10 surface peaks per direction and thus enable statistical analyses of the contact parameters, especially in isotropic surface case ( ≈ 100 surface peaks). The xy plane is discretized into 512 × 512 points, which is a compromise between accuracy and computational efficiency in the contact analyses compared to the finer discretizations (e.g., mesh convergence study reveals less than 5% difference from the finest discretization of 2048 × 2048). Figure 1 shows the surface heights h (x, y) for two representative surfaces created by this procedure. All physical lengths in x and y directions and height h (x, y) in Figure 1 are normalized by window size L, 10 [a.u.]. As shown in Figure 1A, limiting the angle Θ ≤ θ = π/180 delivers an anisotropic surface geometry where height variations in y direction are minimal compared to the variations in x. That surface as an extreme example resembles a 2D sine wave profile. When the direction vector is selected randomly in all possible directions; i.e., Θ ≤ θ = π/2, isotropic distribution of the surface heights is achieved (Figure 1B).

Surfaces generated by amplified instability employ dynamics of a rotating drum to amplify instabilities and pattern formation in thin films of polymer blends such as polydimethylsiloxane (PDMS) and vinylpolysiloxane (VPS-8) (Jambon-Puillet et al., 2021). Here, we simulate this pattern formation by following the physical processes and corresponding equations reported in (Marthelot et al., 2018; Jambon-Puillet et al., 2021). In this process, one first coats a drum of radius R with a thin layer of polymer melt (PDMS + VPS-8) with surface tension γ and density ρ. Then, the drum is rotated at a slow speed Ω₁ to induce Rayleigh-Taylor instabilities in the polymer melt and form an initial hexagonal array of “pendant droplets”. The wavelengths and heights of this initial pattern depend on the characteristic length l c = γ / ( ρ R Ω 1 ) , a measure of capillary length for the melt (Marthelot et al., 2018). The profile of individual droplets can be found by the Young–Laplace equation Eq. 3 that describes the balance of hydrostatic pressure due to the centrifugal acceleration and the effect of the surface tension of the melt (Ku et al., 1968), where s is the arc length, φ(s) is the angle between the tangent vector to s and the vertical axis, and r(s) and z(s) are the radial and height coordinates defining the droplet as shown in Figure 2. d 2 φ s d s 2 = − cos ⁡ φ s ℓ c 2 + d d s cos ⁡ φ s r s , d z s d s = cos ⁡ φ s , d r s d s = sin ⁡ φ s (3) Solving Eq. 3 with properly selected boundary conditions delivers the axisymmetric initial droplet profile, z(s) and r(s). As in (Jambon-Puillet et al., 2021), we use the following boundary conditions: φ(0) = φ(s _f ) = −π/2, r (s _f ) = 0, z (0) = 0, and z (s _f ) = h ₀, where s _f is the total arc length and h ₀ is the maximum height of the initial droplet. The first two of those boundary conditions relate to vanishing slopes of the profile at the edge (s = 0) and symmetry axis (s = s _f ) of the profile, and the last three follow from the coordinate definitions and instability-driven initial height, h ₀ at the axisymmetry axis. Setting h ₀ to a given value, one can use the shooting method to solve Eq. 3 and initially unknown total arc length s _f , and obtain individual droplet profile as shown in Figure 2. Note that r _f denotes a radial extent of the initial droplet. Marthelot et al. (2018) and Jambon-Puillet et al. (2021) show that theoretical solution to Eq. 3 captures well the profile of experimentally manufactured droplet by Rayleigh-Taylor instability. Because we assume the initial droplets are perfectly axisymmetric, we can simply render 3D shape of droplets h (x, y) in discretized xy plane for a single unit: there exists ( x i ) 0 ≤ i < 2 n − 1 and ( y j ) 0 ≤ j < 2 n − 1 such that, for any a = 0, 1, 2, …, n−1 satisfies { h ( x i , y j ) | r a ≤ ( x i − x n ) 2 + ( y j − y n ) 2 0.5 < r a + 1 } = z ( r a + 1 ) − z ( r a ) r a + 1 − r a ( x i − x n ) 2 + ( y j − y n ) 2 0.5 − r a + z ( r a ) , where x _i and y _j denote the ith point of x and jth point of y, respectively, on the discretized xy plane (i.e., linear interpolation). The half-profile of droplet shown in Figure 2 is discretized with n number of points, and r _a denotes the ath point in discretized r, e.g., z (r ₀) = h ₀ and z (r _n−1) = 0. There is no droplet growth between the droplets, and the film thickness of the elastomers after droplet generation is very thin, so we assume the surface patch between the droplets and hairs is perfectly flat; i.e., for any i, j which does not satisfy r a ≤ ( x i − x n ) 2 + ( y j − y n ) 2 0.5 < r a + 1 , h (x _i , y _j ) = 0. Repeating this procedure for many such droplets and distributing them in a hexagonal array over the surface would form the initial pattern. In line with the experimental findings of (Jambon-Puillet et al., 2021), we select h ₀ from a log-normal distribution with 0.034l _c of mean and 0.36l _c of standard deviation and obtain the profile of 115 droplets. We then repeat those droplets with a wavelength of Λ = 2 π 2 l c over the surface with 10Λ × 10Λ area 0.0373 m × 0.0373 m in physical dimensions). The surface is discretized into 1,024 × 1,024 points. This discretization is determined by the comparison with the case of 2024 × 2024 discretization in the same manner for the spinodal surface.A subsequent increase in the rotational speed of the drum leads to a further increase in the heights of the initially seeded droplets. Increasing acceleration promotes this growth while simultaneous curing and thus the formation of elastic resistance balances inertial forces and prevents stripping of partially-cured polymer blend from the drum surface. Resulting droplet profiles look like elongated versions of the initial profiles, “stretched hairs” (Jambon-Puillet et al., 2021) (see Figure 2 for an example). Jambon-Puillet et al. (2021) utilized minimization of energy and provided an estimation for the ratio of final to initial maximum heights of the droplets, ε = h 0 * / h 0 , as a function of effective acceleration defined by g eff = R ( Ω 2 2 − Ω 1 2 ) where Ω₂ is the increased rotation speed of the drum, where β ≈ 0.29 and α = ρg _eff h ₀/G is a nondimensional parameter estimated for each droplet of height h ₀ and shear modulus G. ε α = 1 3 β α + g α + β 2 α 2 g α , g α = β 3 α 3 + 3 / 2 9 + 12 β 3 α 3 + 81 3 (4) Note that α is a measure of acceleration-induced pressure in terms of the elasticity of the partially-cured blend. Assuming lubrication theory, droplet shapes forming under the influence of gravity and surface tension were previously discussed in (Lister et al., 2010). In particular, the profile of an individual droplet at equilibrium can be expressed as Eq. 5, where J₀ is the Bessel function of the first kind of order 0 and r _f is the total radial length. z r = ε h 0 J 0 r f / l c − J 0 r ε / l c J 0 r f / l c − 1 (5) Jambon-Puillet et al. (2021) found that Eq. 5 approximated the measured profiles of the stretched hairs reasonably well. We employ Eq. 5 to obtain the profiles of individual stretched hairs as shown in Figure 2. To obtain the final surface patterns (e.g., Figure 3B), we repeat this procedure for all 115 droplets. The histogram of the peak heights of the stretched hairs shown in Figure 3B agrees reasonably well with the measured histograms reported in (Jambon-Puillet et al., 2021). Thus, our assumptions on the acceleration-induced shape changes of initial droplets deliver surface patterns similar to the ones obtained in practice. Note that the stretch ɛ, the growth factor of the peak height for each droplet Eq. 5 is a strong function of the parameter α and thus effective acceleration. As detailed by Jambon-Puillet et al. (2021), axisymmetric profiles given by Eqs 3–5 are obtained when α ≲ 8. Therefore, in this study, we investigate surface patterns for effective acceleration values residing within that limit. Lastly, the initial and final profiles of the droplets simulated by this process (Figure 3) exhibit clear differences. All physical lengths in x and y directions and height h (x, y) in Figure 1 are normalized by window size L, 0.0373 m. The stretch in heights due to increased acceleration is accompanied by shrinkage in radial extent as illustrated for an individual droplet in Figure 2. Both of those changes in geometry lead to an increase in curvature. Hence, effective acceleration as a surface processing parameter is expected to have a prominent influence on contact mechanics. We will study this more in Section 4.

To simulate the elastic response of patterned surfaces under normal contact, we employ TAMAAS, an open-source C++-based Python Library developed by the Computational Solid Mechanics Laboratory at EPFL (Frérot et al., 2020). We simulate non-adhesive elastic contact and so utilize the conjugate gradient-based contact scheme (Polonsky and Keer, 1999) and boundary integral method. For that, as explained in the previous section, we first generate surfaces simulating the patterns created by the spinodal decomposition and amplified instability. Patterned surfaces are treated as rigid punches repeating periodically outside the xy domain to be pressed on elastic half-space. We assign 1 [a.u.] for Young’s modulus E and 0.5 for Poisson’s ratio ν of the elastic half-space. We then impose the nominal mean contact pressures p ₀ ranging from 0.001E* to E* and solve for the contact pressures at each grid point, p (x, y). All pressure units are normalized to the contact modulus E* = E/(1 − ν ²).

Figure 4 shows representative contact pressure distributions obtained from the isotropic spinodal patterned surface and the initial surface generated by the amplified instability. Both surfaces are imposed to nominal mean contact pressures p ₀ of 3.16 × 10⁻² E*. One can see that only a minute portion of the nominal contact area is load-bearing and local contact pressures can attain very high values up to 2.5E* at the peaks of both surface patterns hinting at inevitable stress concentrations and potential failure. Both features are signatures of randomly rough surfaces. We will show in Section 3 that surface patterning via the two methods mentioned in the previous section can lead to a reduction in contact areas and stress concentrations comparable to randomly rough surfaces. To quantify the real area of contact A _r , we estimate the sum of grid areas with non-zero pressures, and obtain the real mean contact pressure p _r as p 0 A 0 A r . Since peak heights and curvatures vary within nominal contact areas of both patterns, there will be local areas with much larger pressures than p _r . To study stress concentrations at those locations, we define p conc = p r + 3 p rms ( x , y ) p 0 , where p _rms denote the rms of contact pressures. This definition captures stress concentrations well without potential amplification and error due to discretization (especially under light loading).

To verify, compare and contrast our numerical results with existing contact theories, we will adopt the following models for different cases.

• For comparison to the anisotropic spinodal surface, we will use 2D-sinusoidal contact formulations (Johnson, 1985). This choice is apparent from Figure 1A, which is created by limiting the direction vector to within θ ≈ 0° and thus surface heights vary along the x-axis only. The ratio of real to nominal areas of contact A _r /A ₀ and the contact pressures p(x) for 2D-sinusoidal surface h (x, y) ≈ z(x) = Δ cos (2πx/λ) are given as Eqs 6, 7, where a is the half contact width at each peak, and p* = πE*Δ/λ;

Figures 5A, B show the real area of contact normalized by the nominal area A _r /A ₀ and real mean contact pressure p _r as a function of mean pressure p ₀ for two extreme cases: θ ≈ 0° and θ = 90°. Note that those cases correspond respectively to the most anisotropic (2D sinusoidal) and isotropic patterns shown in Figure 1. For each case, 20 randomly sampled surfaces are realized with randomly selected direction vectors and phase angle per wave components Eq. 2. The mean and ± 1 standard deviation of the responses of those 20 surfaces are shown in Figure 5 and the subsequent figures. As expected for both cases, the real area of contact is less than a few percent of the nominal area at light loading (p ₀ ≈ 10⁻³ E*). Around that light-loading asymptote, discrete GW simulations overlap nearly perfectly with the numerical solutions, which further verifies the numerical modeling. 2D sine contact theory correlates with the numerical solutions of the anisotropic case even at heavier loads reaching up to p ₀ ≈ E*. At that level of loading, complete contact condition is attained; i.e., A _r ∼ A ₀. 2D sine contact theory also predicts that complete contact is reached when p ₀ = p* = πE*Δ/λ, which is ≈ 0.55 E * for the 2D surfaces rendered. At p ₀ = 0.55E*, the mean of numerical solutions is A _r /A ₀ ≈ 0.85 suggesting 15% deviation from the theory, which is pretty close to the complete contact. This comparison with theory confirms our expectation on anisotropic cases resembling 2D-sinusoidal wavy patterns whereas isotropic cases delivering contact performance similar to randomly rough surfaces.

Reduction in the real load-bearing areas leads the inevitable increase in real mean contact pressures. Figure 5B shows the real mean contact pressures p _r obtained from the simulations and contact theories. The real mean contact pressures at light loading (p ₀ < 0.01E*) are constant around 0.3E* for the isotropic patterns. Discrete GW model predicted similar values. The continuous form of GW model with randomly rough surfaces also predicts a constant real mean contact pressure independent of the applied load. This is thanks to the proportional increase of load-bearing areas with increasing load. For the anisotropic case, real mean contact pressures increase nonlinearly with the applied load and 2D sinusoidal contact solution captures that change reasonably well. At heavier loads (p ₀ > 0.1E*), both anisotropic and isotropic patterns deliver similar real areas of contact and mean contact pressures. At that asymptote, sensitivity to the initial patterns is lost through long-range interactions, shape alterations, and mergers among contact patches. Similar observations were reported elsewhere with detailed numerical analysis of sinusoidal contacts (Yastrebov et al., 2014).

Figures 6A, B show the real area of contact and the real mean contact pressure for the surfaces generated by amplified instability. Numerical simulations and corresponding discrete GW model predictions are plotted for three different instability cases. The first case is the initial surface rendered at Ω₁ = 41.9 rad/s, and the second and third cases are created with incrementally increased rotating speeds corresponding to effective accelerations of g _eff = 137 m/s² and g _eff = 271 m/s², respectively. The second and third surfaces are originated from the same initial surfaces. Note that initial heights are selected randomly from an experimentally-measured distribution, which then leads to a slightly different set of peaks for each surface patterns. Twenty surfaces are rendered at each case to determine statistically-relevant contact parameters numerically. The discrete GW model is obtained for one of those surfaces delivering mean values of the contact parameters shown in Figure 6. The shaded area represents the mean ± 1 standard deviation of numerical results for 20 surfaces per case. Similar to spinodal surfaces, all of the cases deliver contact areas much smaller than the nominal area under light loading. Besides, the real area of contact seems proportional to the applied load, which leads to constant real mean contact pressures at that extreme (p ₀ < 0.01E*). The real area of contact is consistently the smallest for high effective acceleration cases and highest for the initial surfaces, and vice versa is true for the real mean contact pressure. As the effective acceleration increases, so do the stretch of initial droplets and the height of individual hairs. Besides, the curvature for stretched hairs increases with the effective acceleration as shown in Figure 2. Both of those changes in geometry lead to less contact area and more stress concentration.

The real mean contact pressure p _r for the initial surfaces and the surfaces created at low g _eff is almost constant from 0.35E* to 0.45E*, while the surfaces created at high g _eff show decreasing real mean contact pressure for p ₀ < 0.1E*. This will be discussed in detail in Section 4 with the asperity curvatures and asperity heights. In Eq. 4, the stretch ɛ increases as α increases, and α is linearly proportional to the initial droplet height h ₀. Consequently, the variance in the initial surface’s asperity heights may result in a wider range of variance in amplified surfaces. We can see the shaded area is getting larger as the increased rotating speed Ω₂ is getting higher.

Figure 7 shows the evolution of stress concentration p _conc on anisotropic and isotropic spinodal surfaces as a function of applied nominal pressure p ₀. As defined in Section 2, p _conc is essentially the ratio of maximum and real mean contact pressures. The stress concentration for both the anisotropic and isotropic surface patterns starts at values ranging from 200 to 800, decreases monotonically with increasing nominal pressures, and converges to about 2 under heavy loading. The factor of 2 is common to any wavy pattern at the inception of complete contact. Except for that heavy-loading asymptote, isotropic patterns experience higher stress concentration than anisotropic patterns.

Figure 8 shows the evolution of stress concentration p _conc on the surfaces created by amplified instability as a function of applied nominal pressure p ₀. Three cases of amplified instability are investigated as above: after initial rotation, after low, and after high effective acceleration. Similar to the spinodal surfaces, stress concentration starts at high values ranging around from 10³ to 10⁴ under low nominal pressures, decreases monotonically, and approaches 2 under high nominal pressures. The patterns created at high effective acceleration (g _eff = 271 m/s²) exhibit higher stress concentration than the other patterns. This is because the initial droplets get taller and sharper as effective acceleration increases (see, for example, Figure 2).