Wrinkles, rucks, and folds formed in a heavy sheet on a frictional surface

Abstract

Soft elastic sheets resting on rigid surfaces develop wrinkles, rucks, and folds due to the combined influence of elasticity, gravity, and contact interactions. Despite their ubiquity, the principles governing their morphology and transitions remain unclear. We introduce a minimal experiment in which the center of a gravity-loaded sheet is gradually lifted from the supporting plane. This operation generates a clear sequence of shapes: an axisymmetric uplift, a finite number of wrinkles, system-spanning rucks produced by global buckling, and folded states that can arise from ruck collapse upon unloading at larger lifts. Combining experiments, finite-element simulations, and Föppl–von Kármán theory, we establish a unified physical picture of this morphology sequence. In the frictionless case, elasticity and gravity alone govern the response, leading to a universal wrinkling threshold: the wrinkle number is fixed and the onset displacement scales linearly with the sheet thickness. With interfacial friction, the wrinkled state is described by introducing an additional nondimensional parameter that compares frictional and elastic–gravitational forces. These results suggest a simple route to programmable sheet morphogenesis via friction and gravity.

---

1 Introduction

Soft, thin structures resting on external surfaces frequently develop visually striking wrinkles, rucks, and folds under their own weight. A familiar example is the cartoon-like “ghost” shape that appears when a soft sheet is draped over a small object on a table [Fig. 1(a)]. Examples of such gravity-induced morphologies span a wide range of systems, from everyday drapery^1–5 to solar sails^6,7 and even geophysical plates.^8–10 In thin yet heavy systems, contact interactions can generate unpredictable surface patterns that pose practical concerns.⁷

Fig. 1 (a) Halloween ghost formed by draping a soft, thin fabric over a convex object. (b) Schematic of the experimental setup. A thin elastic sheet of radius a is indented vertically at its center by a distance d from beneath the sheet on a rigid flat substrate. A region with a certain radius R detaches from the substrate. The indentation force F is measured as a function of the indentation displacement d by a load cell attached to the indenter. (c) Representative images of three distinct lifted shapes for increasing d: (i) axisymmetric, (ii) wrinkled, and (iii) globally buckled. (d) and (e) Color maps of the vertical displacements z = w(x, y) for the indicated values of d. (d) The point-cloud data acquired by a 3D scanner in our experiments and (e) those in our FES. The parameters used both in experiment and FES are E = 477 kPa, ν = 0.47, ρ = 1127 kg m⁻³, h = 0.27 mm, a = 104 mm and μ = 0.32.

Wrinkles may be deliberately introduced through in-plane loading,^11–13 but more commonly emerge from the uncontrolled sticking and sliding of thin materials against their surroundings.^14–16 When a thin sheet interacts with external surfaces, the contact forces and friction significantly complicate the buckling process, producing patterns that are far more diverse than those arising from simple in-plane compression alone.^17–22 These effects are particularly important in gravity-loaded sheets, where the self-weight naturally brings the material into contact with the supporting surface.^23–27

Although the mechanics of elastic beams or rods in contact with external surfaces can be analyzed in considerable detail,^25–30 the mechanics of two-dimensional sheets remain far less tractable. In sheets, the curvature and in-plane stress are coupled through geometric compatibility,³¹ making the influence of the contact forces far more complex than in slender one-dimensional systems.^16,32–36 In this context, indentation tests offer an appealing minimal model with a simple axisymmetric geometry. They have been widely used to probe thin films–from biological membranes³⁷ to nanoscale sheets³⁸ – and readily reveal symmetry breaking into wrinkles or folds.^{35,36,39–41}

To investigate the mechanism by which gravity and contact interactions shape the morphology of heavy sheets, we introduce an indentation experiment in which an elastic sheet resting on a rigid substrate is lifted gradually from its center [Fig. 1(b)]. This operation produces a sequence of patterns [Fig. 1(c)–(e)]: initial axisymmetric uplift, a finite number of wrinkles, and then system-spanning rucks produced by buckling at the periphery (global buckling). Folded states can also arise upon unloading after large lifts [Fig. 13 and SI Movies]. Using a combination of tabletop experiments, finite-element simulations (FES), and Föppl–von Kármán theory, we characterize the mechanical response across this sequence. In the frictionless case, the wrinkled state shows remarkable universality: the wrinkle number is fixed, and the critical onset displacement scales linearly with the sheet thickness. When friction is present, it changes the hoop stress distribution, and both the wrinkle number and its onset are described by introducing an additional dimensionless parameter measuring the relative importance of frictional and elastic–gravitational forces.

The remainder of this paper is organized as follows. Section 2 describes experimental setup, including sample fabrication, apparatus, measurement protocols, and the finite-element methodology. Section 3 presents the force response and growth of the lifted region, and identifies a characteristic length scale that organizes these observations. Section 4 develops a theoretical framework based on the Föppl–von Kármán equations and compares its predictions with experiments and simulations. Section 5 presents the numerical results for the in-plane displacements and stress distributions in the contact region. Section 6 examines the wrinkling transition through additional experiments, FES, and theoretical considerations, including the effects of friction. Section 7 derives a scaling law for the onset of global buckling and validates it against experimental and numerical results. Section 8 highlights the hysteresis and fold formation observed during unloading. Finally, Section 9 summarizes the main findings and discusses future research directions.

2 Indentation test

2.1 Experiments

In our experiments, we used circular sheets with a uniform thickness h in the range 0.14 mm ≤ h ≤ 4.96 mm, and radius a cut in the range of 40 mm ≤ a ≤ 118 mm. To obtain sheets of uniform thickness, we performed spin-coating with addition-cure silicone rubbers (Elite Double 8 (Zhermack, Italy), HTV-4000 (Engraving Japan), Mold Star 15 SLOW and Ecoflex 00-20 (Smooth-On, USA)). Their Young's modulus E, Poisson's ratio ν, and mass density ρ are summarized in Table 1. Talc was sprinkled to the sheet surfaces to prevent adherence. We also used a commercially available urethane sheet with h = 2.0 mm, E = 3780 kPa, and ν = 0.29. All the h values were measured using a laser displacement sensor (LK-G3000, KEYENCE, Japan). E and ν are measured using cantilever bending tests or tensile tests.

Table 1 Properties of the elastic sheets fabricated by authors, including: the material, Young's modulus E, Poisson's ratio ν, mass density ρ, and the coefficient of static friction between the sheet and substrate μ

Material	E [kPa]	ν	ρ [kg m^−3]	μ
HTV4000	756	0.34	1142	0.34 ± 0.04
Mold Star 15 SLOW	477	0.47	1127	0.32 ± 0.01
Elite Double 8	226	0.49	1030	0.44 ± 0.03
Ecoflex 00-20	44	0.43	1053	0.56 ± 0.03

---

We used a stainless-steel substrate with a diameter of 300 mm and a thickness of 1.5 mm, which can be treated as a semi-infinite rigid plane [Fig. 1(b)]. The substrate had a circular hole with a diameter of 11 mm at its center. The center of the sheet was pushed through this hole using a cylindrical indenter with a length of 30–70 mm and radius r_ind = 0.69–1.70 mm. Because the hole was small but larger than the indenter, airflow could be generated among the hole, indenter, and deformed sheet, preventing any vacuum effects during indentation. The coefficient of static friction between the talc-coated elastomer and the stainless-steel substrate was measured to be μ = 0.32–0.56 [Table 1] using a slip-angle measurement experiment. Before placing the sheet, we spread a cationic surfactant (an antistatic agent, MonotaRO, Japan) on the substrate to suppress the static electricity effects. After evaporation, we placed a sheet and blew it with a hair dryer to remove the fine pre-wrinkles.

The indenter was positioned beneath the substrate as shown in Fig. 1(b). The vertical motion of the indenter was controlled using a stepping motor (ARM46AC, ORIENTAL MOTOR, Japan). The indenter moved sufficiently slowly upward (0.1 mm s⁻¹), thereby imposing a vertical displacement d on the sheet. The reaction force F exerted at the center of the sheet was measured using a load cell (LTS-2KA, KYOWA, Japan) attached to the indenter, and recorded as a function of d.

To quantify the lifted shape, we obtained point-cloud data w(r, θ) using a 3D scanner (EinScan-SP, SHINING 3D, China). The scanning experiment was conducted separately from the force measurements. We placed the indentation apparatus on a desktop 3D scanner and carefully raised the indenter by Δd ≈ 1 mm using a hand-controlled labjack; a point-cloud dataset was acquired at each step [Fig. 1(d)]. From these data, we extracted the mean lifted radius R for axi-symmetric and wrinkled states. After global buckling, where axisymmetry was completely lost, we defined the effective radius through the lifted area S_lift as .

Our method is a type of indentation test, often referred to as a “blister test,” which is traditionally used to measure the strength of adhesion in thin films.⁴² Several studies have examined pattern formation in adhered films (that can still slide laterally).^39,43,44 Although motivated by similar considerations, our system differs in that gravity, rather than adhesion, prevents the lifting, and our focus is on the influence of dry friction.⁴⁵

As h becomes significantly larger than the indenter radius r_ind, the local contact near the indenter is expected to approach that predicted by punch-indentation theory for a Boussinesq-type problem.⁴⁶ Because the sheets used in our experiments were sufficiently thin, we did not analyze the detailed 3D contact deformation; instead, we focused on the overall shape evolution of the sheet.

2.2 Finite-element simulations (FES)

To complement our experimental results, we performed FES using the commercial package Abaqus (Dassault Systèmes, France). A linear elastic circular sheet was modeled using quadrilateral linear shell elements with reduced integration and finite membrane strain (S4R). The element size Δx was chosen to be sufficiently small to resolve azimuthal variations such as wrinkle wavelengths, ensuring λ ≫ Δx. The stainless-steel substrate was modeled as a rigid shell, and normal and tangential contact interactions were included.

To prevent the undesired initial penetration, the sheet was initially placed slightly above the substrate and then we dropped onto it without friction. After equilibration, Coulomb friction was introduced, and the central circular region of radius r_ind was raised at 0.1 mm s⁻¹, with r_ind/a ≲ 0.01. Both the dropping and indentation steps were computed using an implicit dynamic analysis (ABAQUS step type: Dynamic, Implicit). Throughout the simulation, the kinetic-to-strain energy ratio was typically ∼10⁻⁷–10⁻⁵, independent of μ. Even in frictionless runs exhibiting global buckling, it peaked at only ∼10⁻², confirming effectively quasi-static behavior. The geometric and material parameters were almost the same as those used in the experiments, while μ was varied over a broader-than-typical range, 0 ≤ μ ≤ 2.0.

When we focused on large displacements (d/h > 10), R was obtained from the lifted area S_lift using , as in the experiments. Note that the definition of R is valid when the pattern exhibits the axisymmetry and quasi-axisymmetry such as in wrinkled one; however, in the case of globally buckled morphologies, the assumed symmetry is generally broken and R merely serves as a guide to the size of the lifted region. Still, by using the common definition of R, it is informative for quantifying the magnitude of the lifted area in the different regimes shown in Fig. 2. For smaller deformations (d/h ≲ 10), the lifted area became difficult to identify accurately because the out-of-plane displacement was very small. In this regime, we instead extracted R from the azimuthally averaged profile .

Fig. 2 Lifting force and lifted radius vs. indentation height. The left axis and red points indicate the lifting force F, normalized by the total weight of the sheet Mg. The lifted radius R, normalized by the sheet's full radius a, is plotted on the right axis (blue points). Filled and open symbols represent data obtained from experiments and FES, respectively. The data are taken from the experiment and simulation presented in Fig. 1. Inset figures show lifted shapes of the sheet obtained from our FES. Axisymmetry breaks at d = d_w, where m-fold wrinkles emerge (m = 8 and d_w ≈ 4.9 mm). The critical displacement for global buckling is d_c ≈ 9.2 mm.

To investigate wrinkling instability, we introduced a small random vertical imperfection at the sheet nodes (|δw|/h < 1%) to initiate symmetry breaking. Simulations performed with and without this imposed imperfection showed nearly identical force–displacement curves and wrinkle onsets, suggesting that the intrinsic numerical and geometrical imperfections already present in the model were likely sufficient to trigger symmetry-breaking.

3 Typical morphology and force response

Fig. 2 shows the typical force–displacement and lifted radius–displacement curves from both the experiments and FES. The sheet material was Mold Star 15 Slow, with thickness h = 0.27 mm and radius a = 104 mm [material properties in Table 1]. We integrate some results for the different parameters into Fig. 3, which shows that both F and R exhibit power law behaviors with respect to d when plotted on a log–log scale. It also appears that the exponents change when d/h ≪ 1 and d/h ≫ 1, which will be rationalized by the scaling argument given below.

Fig. 3 Dimensionless (a) lifted radius R/ℓ_g and (b) lifting force F/(ρghℓ_g²) plotted against normalized central displacement d/h. The characteristic length ℓ_g is defined in eqn (2). Filled and open symbols denotes experimental and FES data, respectively. The scaling behaviors predicted by eqn (1) and (3) are confirmed in both the shallow (d ≪ h) and large indentation (d ≫ h) regimes, as shown by the dotted and solid lines with the coefficient in eqn (16), (17), (23) and (24).

Although axisymmetry is lost and wrinkles appear at d = d_w (Fig. 2), the F(d) and R(d) curves remain monotonic and are unaffected by symmetry breaking. By contrast, at the onset of global buckling (d = d_c), the force exhibits a discontinuous drop, and R(d) curve shows a sharp change in the slope. After global buckling, F and R increase more gradually.

We now focus on the power-law behavior of F and R in the axisymmetric and wrinkled phases. For a lifted region of radius R and height d, the typical strain and curvature scales are ε ∼ (d/R)² and K ∼ d/R². The corresponding stretching and bending energies are as follows: and , where B ≡ Eh³/[12(1 − ν²)] is the bending modulus.³¹ Both decrease with R, and their ratio [scr E, script letter E]_s/[scr E, script letter E]_b ∼ (d/h)² indicates that bending dominates for d ≪ h, while stretching dominates for d ≫ h. Gravity adds an energy cost [scr E, script letter E]_g ∼ ρghR²d, which prefers a smaller R. Minimizing the total energy [scr E, script letter E]_tot = [scr E, script letter E]_s + [scr E, script letter E]_b + [scr E, script letter E]_g with respect to R yields

wheredefines an “elasto-gravitational length,” and c_b and c_s are numerical prefactors determined later. Although the term “elasto-gravitational length” is often referred to as the characteristic wrinkling wavelength, λ ∼ (B/ρ_lg)^1/4, of an elastic sheet floating on a liquid of density ρ_l,^47–49 the density ρ considered in this study corresponds to that of the sheet itself. The former length arises from a balance between the bending forces and buoyancy, whereas the latter is determined by the competition among bending, gravitational, and stretching forces.

The difference between the well-known gravito-bending length^1,50 ℓ_gb = (B/(ρgh))^1/3 and the ℓ_g introduced here is significant. The former is the characteristic length that appears when the bending elasticity balances gravity such as drapes or beams bent substantially by their own weights. In our situation, however, stretching elasticity prevents large-amplitude deformations (until global buckling occurs). Therefore, the dominant length scale in our problem is given by ℓ_g, which includes stretching elasticity.

Using R in eqn (1), we find that the total energy of the system scales ρgℓ_g²h^1/2d^3/2 when d ≪ h, and ρgℓ_g²h^−1/2d^5/2 when d ≫ h. They must balance the work done during indentation, Fd. Therefore we obtain

where k_b and k_s are additional dimensionless prefactors.

Note that the scaling discussion here implicitly assumes that the sheet can be treated as a thin structure and that the indenter size is sufficiently small. A detailed examination of the case where d is extremely small compared to h reveals that the effect of the indenter shape becomes more pronounced, and the thin-plate approximation may no longer be valid. In the present study, we therefore focus on the regime in which the sheet can reasonably be regarded as thin (d/h ≳ 10⁻¹). To experimentally and numerically verify eqn (1) and (3), we restrict attention to the region where the indenter size is smaller than the representative lengths ℓ_g(d/h)^1/4 and ℓ_g(d/h)^3/4. For the elastomer used in our experiments (ℓ_g ∼ 1 mm), to observe the behavior in Fig. 3 for 10⁻¹ < d/h ≲ 1, it is desirable to use an indenter smaller than ℓ_g(d/h)^1/4 ∼ 1 mm. Accordingly, we employed an indenter with radius r_ind = 0.69 mm in the experiments. In the FES, we used r_ind = 0.1 mm when d/h < 10, and a larger indenter of radius r_ind = 1.0 mm for d/h > 10. For d ≫ h, the relation R ∼ ℓ_g(d/h)^3/4 ensures that r_ind/R ≪ 1, justifying the approximation of a point indenter.⁵¹ Additional FES results showing the influence of r_ind are provided in the SI Section I.

The dashed and solid lines in Fig. 3 represent the predicted scaling laws (1) and (3). Up to the global buckling threshold d = d_c, a wide range of data for elastic modulus, thickness, and even friction coefficient appears consistent with the predicted lines.

4 Theoretical analysis of the axisymmetric state

Our aim here is to analytically determine the pre-factors in the scaling relations eqn (1) and (3). We define a cylindrical polar coordinate (r, θ, z) and the origin is set to coincide with the center of the bottom surface of the sheet. We introduce the displacement vector of the material points on the middle surface of the sheet as u(r, θ) = u_r(r, θ)e_r + u_θ(r, θ)e_θ + w(r, θ)e_z. w(r, θ) represents the out-of-plane displacement, and the middle surface of the deformed sheet is represented by z(r, θ) = h/2 + w(r, θ).

The mechanical equilibrium of a lifted sheet is described by the Föppl-von Kármán (FvK) equations:³¹

B∇⁴w − hσ_αβK_αβ + ρgh = 0, (4)

h∇·σ = 0, (5)

where σ_αβ is components of the in-plane stress tensor, and the Greek indices run over (r, θ). The stress and strain are related by Hookean linear constitutive equations, whereas the strain includes geometric nonlinearity in w as ε_rr = u_r,r + w_,r²/2, ε_rθ = ε_θr = u_r,_θ/(2r) +u_θ,r/2 − u_θ/(2r) + w_,rw_,θ/(2r), and ε_θθ = u_r/r + u_θ,θ/r + w_,θ²/(2r²), where f_,r ≡ ∂f/∂r, and f_,θ ≡ ∂f/∂θ.³¹K_αβ is the curvature tensor, given by K_rr = w_,rr, K_rθ = K_θr = ∂_r(w_,θ/r), and K_θθ = w_,r/r + w_,θθ/r². Eqn (4) and (5) represent the vertical and lateral force balances, respectively.

Assuming axisymmetry, the FvK equations can be reduced to the following ODEs:

where ψ(r) denotes so-called the derivatives of the Airy stress function defined as hσ_rr(r) = ψ/r, hσ_θθ(r) = dψ/dr.⁴¹ The in-plane equilibrium (5) is automatically satisfied by this representation, whereas ψ(r) must obey the compatibility relation (7). The reaction force exerted at the center of a sheet equals to F, so thatUsing this relationship, we integrate eqn (6), we obtain thatNext, we analytically solve the coupled eqn (7) and (9).

4.1 Small displacement: d ≪ h

We introduce the following dimensionless variables, motivated by the scaling relations eqn (1) and (3) for d ≪ h:With this nondimensionalization and assuming R = c_bℓ_g(d/h)^1/4, eqn (9) and (7) becomewhere f′ ≡ df/dξ and . Focusing on d ≪ h, the second term in (11) can be neglected, and the equations for W and Ψ, i.e., eqn (11) and (12) are decoupled. The vertical displacement can then be determined from a single ODE as follows:which can be solved exactly.^47,52 We solve this by the following boundary conditions:

W(0) = 1, W′(0) = W(1) = W′(1) = W″(1) = 0. (14)

The final condition is the moment-free condition at the detachment points r = R.^25,28,31 Solving this, we obtain

W(ξ) = 1 − ξ⁴ + 4ξ² ln ξ, (15)

c_b = 2^3/2, (16)

k_b = 4π. (17)

In Fig. 3, we compared eqn (1), (3), (16) and (17) with the experimental and FES data. Despite the simplifying assumption of the point-load indentation, we observe that our analytical predictions show reasonable agreement with numerical and experimental data obtained for finite size indenters.

4.2 Large displacement: d ≫ h

Next, we introduce the following dimensionless variables:where R = c_sℓ_g(d/h)^3/4. For these variables, eqn (6) becomesand considering d ≫ h, eqn (9) and (7) are reduced towhere α ≡ c_s⁴/(12(1 − ν²)).

Exact analytical solutions of eqn (20) and (21) are known in α = 0, and take a particularly simple closed form for the special value of the Poisson's ratio ν = 1/3.^39,53 However, for α > 0, they do not seem to admit any concise analytical solutions. As we can estimate α ∼ 0.3 from our experimental and FES results, we treat the second term in the left-hand side of eqn (20) as a perturbative term. Therefore, we use the closed-form solution at α = 0 as the unperturbed state and compute O(α) corrections by expanding W, Ψ, and [scr F, script letter F] in α, subject to the boundary conditions

with ν = 1/3. The last condition corresponds to a clamped radial displacement at r = R, i.e., u_r(R) = 0. Although this condition is not satisfied exactly in our system (see below), it is still reasonable for our present purpose, because the radial displacement at the detachment radius remains small compared to other length scales (for example, Fig. 5 shows |u_r(R)| < h).

The perturbation calculation yields the approximate solutions for W, Ψ, and [scr F, script letter F] as a series in α. We then determine the unknown constant α by minimizing [scr E, script letter E]_s + [scr E, script letter E]_g with respect to R, which yields

The full derivation is provided in SI Section II. Eqn (1) and (3) with the prefactors (23) and (24) are plotted as solid lines in Fig. 3, and show excellent agreement with both the experiment and FES.

From this analysis, we also obtained the in-plane stresses σ_rr and σ_θθ in the lifted region (see SI Section II). However, the discrepancy between these analytical stress profiles and the FES becomes significant near r = R [see Fig. 4]. Although the analysis based on the clamped boundary condition, u_r(R) = 0, provides good estimations of F and R, it fails to predict a detailed stress distribution.^39,54 Because the center of the sheet is pushed upwards, the entire region of 0 < r < R must be stretched both radially and azimuthally under the clamped boundary condition. Consequently, the present theoretical analysis predicts the positive stress field over the entire considered region. In contrast, both in the indentation test and FES, the material is allowed to slide on the substrate at r = R. As d increases, material points undergo a small inward radial displacement, and the azimuthal direction contracts, leading to the appearance of the circumferentially compressive region σ_θθ < 0. Thus a more accurate analytical description of the stress field requires matching the solution of the FvK equations for r < R with a planar solution for r > R.^39,41,44 We leave this matching problem for future work, and we here will focus on the fundamental properties of the stress field within the contact region.

Fig. 4 The profile of the azimuthally averaged stress is plotted normalized by the characteristic stress scale B/ℓ_g². Dots and solid lines represent a FES result based on the simulation described in Fig. 1 (with d = 4.66 mm) and our theory in Section 4, respectively. The shaded region corresponds to the contact region.

5 Displacement and stress in the contact region

In previous studies on the indentation-induced wrinkling of thin sheets,^39,41,49 the emergence of wrinkles has been attributed to the destabilization of axisymmetric configurations once the compressive azimuthal in-plane stress exceeds a buckling threshold. Here we are particularly interested in the role of contact forces in wrinkle formation. Therefore, we focus on the in-plane displacement field and the stress distribution within the contact region, and we reveal how friction affects the hoop stress distribution.

In the following, we often use azimuthally averaged quantities defined as , for example, when characterizing the displacement and stress fields. Numerically, f̄(r′) is computed by discretizing r and averaging f over all nodes within a small radial interval r′ − Δr < r < r′ + Δr. When the system exhibits a high degree of axisymmetry, f̄(r) provides a useful representation of the spatial variation. However, it should be noted that for globally buckled structures, in which axisymmetry is broken, f̄(r) should be regarded only as the quantity for reference.

5.1 Displacement

We first examine the in-plane displacement using the simulation data shown in Fig. 1 (with a finite value of μ = 0.32). The displacement vector u(r, θ) in the contact region is shown in Fig. 5(a). The color of the arrows represents |u(r, θ)|/Δx, where Δx is the representative mesh size. We observe inward slip within the contact region as the indentation proceeded, most prominently near r ∼ R.

Fig. 5 Displacement field in the contact region for (i) small indentation height (d = 0.86 mm) and (ii) large indentation height (d = 4.66 mm). (a) Quarter of the sheet is shown as the shaded area. The color of the inward-pointing vectors represents the magnitude of the displacement u(r, θ) in the contact region, normalized by the representative mesh size Δx (here equals to 1 mm) of the finite elements. (b) Radial profiles of the azimuthally averaged radial displacement, defined as , plotted as a function of r/R. The solid and dashed reference lines indicate and −10⁻³, respectively. (c) Radial displacement at the edge of the sheet, , as a function of d. The values of are extracted from the red circle markers in panel (b). indicates that the entire sheet has slipped. In the parameter ranges explored in this study, almost all of the sheet slip completely before the onset of wrinkling. These results are obtained from the finite-element simulations shown in Fig. 1.

For sufficiently small d, a region of moderately large displacements (|u|/Δx ≳ 10⁻³) exists, whereas the displacement in the surrounding area is much smaller than the mesh size Δx [Fig. 5(a-i)]. The fact that |u| is extremely small compared to Δx may suggest that the shear stress between the substrate and the sheet does not exceed the maximum static frictional force (per unit area). Based on this observation, we assume the threshold for determining whether or not the sheet slides with respect to the substrate as |u|/Δx = 10⁻³. Fig. 5(b-i) shows the azimuthally averaged radial displacement , which reveals that the interface between regions of “large” and “small” displacement lies near r/R ∼ 5. As d increases, this interface moves outward, and eventually the entire contact region undergoes a substantial inward slip [see Fig. 5(a-ii) and (b-ii); see also SI Section III]. We refer to this state as complete slipping.

We also plot the azimuthally averaged radial displacement at the sheet edge, r = a, in Fig. 5(c). This indicates that the contact region has already reached a completely slipping state at the onset of wrinkling. This scenario, in which wrinkle formation occurs after the onset of complete slipping, is observed for a wide range of parameter values considered in this study. The critical displacement at the onset of complete slipping is further discussed in SI Section III. Accordingly, in the following analysis, we restrict our theoretical investigation to the stress field at the contact surface in the case of complete slipping.

5.2 Stress

The indentation speed in our experiments is quite low (ḋ = 0.1 mm s⁻¹); therefore we proceed the theoretical analysis with the quasi-static limit ḋ → 0. In the Abaqus simulations, we use the default Coulomb friction model, in which the interfacial shear traction saturates at the Coulomb threshold μf_n and is treated as independent of sliding velocity once slip has initiated,⁵⁵ where f_n represent the normal force. Although in principle, rate- and state-dependent friction can lead to a small difference between the friction coefficient at the onset of sliding and that during slow steady sliding,^45,56 these effects are expected to be negligible under the present loading protocol. It is therefore reasonable, both in our experiments and in the corresponding Abaqus modeling, to treat the interfacial tangential traction as being everywhere close to the Coulomb threshold, with an effective magnitude f_r ≈ μρgh, which we regard as spatially uniform over the contact surface. Under this assumption, the axisymmetric in-plane force balance⁵⁷ within the contact region isThis equation can be solved once the boundary values σ_rr(a) = 0 and σ_rr(R) are specified. In the frictionless case (μ = 0), we immediately obtain a solution of σ_θθ(r) ∼ −σ_rr(R)r²/R². Such a negative stress field cannot be matched smoothly at r = R with the clamped solution for r < R derived from Section 4. Instead, Fig. 4 shows that the FES results exhibit σ_θθ < 0 near r ≈ R, consistent with the above prediction. The negative hoop stress at r = R arises because the vertical deflection pulls material radially inwards (u_r(R) < 0), and it can destabilize the axisymmetric state.^39,41,43,49

The analytical solutions of eqn (25) are provided in eqn (S45) and (S46) in the SI; here, we present a simpler approximation valid for r ∼ R and R ≪ a:

We fit the FES data using the exact expressions (S45) and (S46), treating σ_rr(R) as a fitting parameter. The fitted value of σ_rr(R) was found to be approximately half of that predicted from the clamped solution for r < R in Section 4. The results are presented in Fig. 6. Our solution reproduced the FES stress profiles well for r > R, both in the frictionless case (μ = 0) and with friction (μ = 0.32), even after wrinkle formation.

Fig. 6 Azimuthally averaged in-plane stress profiles obtained from FES. Open and closed symbols represent the radial and hoop components of the azimuthally averaged stress tensor, respectively. The solid line indicates the planar solution of the force balance eqn (25) with σ_rr(R) determined from fitting to the FES data. The shaded region corresponds to the contact region. (a) The results without friction μ = 0, with all the other parameters identical to those in Fig. 1. (b) The results obtained from the same parameter set as that shown in Fig. 1 with μ = 0.32. The critical displacements d_w in (a) and (b) are 4.4 and 4.9 mm, respectively.

Eqn (26) and (27) show that, near r ∼ R, the radial stress σ_rr is only weakly affected by friction, whereas the magnitude of the compressive hoop stress σ_θθ is reduced by the μρga term. The frictional dependence of σ_rr(R) is shown in SI Section IV.

6 Onset of wrinkling instability

Building on the previous sections, in which we clarified the stress fields in both the lifted and contact regions, we now investigate how these in-plane stresses contribute to wrinkle formation through a linear stability analysis of the FvK equations.^49,58

We consider a perturbative solution of the FvK equations describing a sinusoidal out-of-plane variation with m wrinkles and a small amplitude f(r), added to the axisymmetric base state w⁽⁰⁾(r):

w(r, θ) = w⁽⁰⁾(r) + f(r)cos(mθ). (28)

Substituting eqn (28) into eqn (4), and assuming axisymmetric stress fields, the equation for f is given byFig. 6 indicates that |σ_θθ| reaches its maximum near r = R. Given that this large compressive stress is essential for wrinkle formation, we estimate the magnitude of each term in eqn (29) in the vicinity of r = R.⁵⁸ The first term in the l.h.s of eqn (29) comes from the bending elasticity ∼Bm⁴f(R)/R⁴, whereas the third term ∼h|σ_θθ(R)|m²f(R)/R² represents the azimuthal compression. The second term, which is the product of the radial tension and curvature, ∼hσ_rr(R)f(R)/R² acts as a Laplace-pressure-like restoring force when |∂_r²f| > 0. This is analogous to that of a compressed elastic beam resting on an elastic foundation^47,48 with an effective stiffness K_eff ∼ hσ_rr(R)/R². By balancing the restoring force from the foundation (∼K_efff(R)) with the bending term, we obtain the optimal wavelength as λ ∼ (B/K_eff)^1/4.¹¹

6.1 Frictionless case (μ = 0)

In the frictionless case (μ = 0), according to eqn (1), σ_rr ∼ |σ_θθ| ∼ E(d/R)² ∼ B/(hℓ_g²)(d/h)^1/2. The balance between the restoring and azimuthal compressive force in eqn (29) yields m ∼ (σ_rr(R)/|σ_θθ(R)|)^1/2 ∼ O(1). Using m = 2πR/λ ∼ (K_eff/B)^1/4R, as derived in the previous paragraph, together with the scaling relations for σ_rr(R) and R, we obtain m ∼ (d/h)^1/2. With m ∼ O(1), this gives d/h ∼ m² ∼ O(1), i.e., the wrinkling threshold depends only on the sheet thickness, d_w ∝ h.

To verify the above prediction, we performed FES for the case of μ = 0 using the parameters listed in Table 2. We investigated the critical displacement at the onset of the wrinkling, d_w and the number of wrinkles, m using a fast Fourier transformation of the node profile, w(R, θ) (using numpy.fft in Python). The results are presented in Fig. 7. When the sheet radius a is sufficiently larger than the elasto-gravitational length ℓ_g, the number of wrinkles typically ranges from seven to eight, with m = 7 being the most frequent, and d_w is proportional to h with a slope of 16.33. By combining our theoretical analysis and FES results, we conclude that

m⁰ ≈ 7, (30)

where the subscript “0” indicates that μ = 0. The constants on the right-hand side are universal, meaning that they are independent of any material parameters, such as the Young's modulus of the sheet.

Table 2 Properties of the elastic sheets used to investigate the onset of wrinkling instability. The table summarizes the method of investigation (experiment or FEM), material type, sheet radius a, sheet thickness h, and the coefficient of static friction μ. It also indicates the data markers used to represent each elastic sheet in Fig. 7–9 for both experimental and numerical results

Marker	Method	Material	a [mm]	h [mm]	μ
	FEM	Elite Double 8	100	0.3	0
	FEM	Elite Double 8	50–200	0.5	0
	FEM	Ecoflex 00-20	100	0.1–0.4	0
	FEM	Elite Double 8	100	0.1	0.01–2.0
	FEM	Elite Double 8	100	0.075–0.2	0.2
	FEM	Ecoflex 00-20	100	0.1	0.6–1.0
	FEM	Mold Star 15 SLOW	104	0.27	0.32
	FEM	HTV-4000	114	0.2	0.34–2.8
	Experiment	Elite Double 8	114	0.20 ± 0.02, 0.17 ± 0.01	0.34 ± 0.04
	Experiment	Ecoflex 00-20	114	0.34 ± 0.02, 0.64 ± 0.06	0.56 ± 0.03
	Experiment	Mold Star 15 SLOW	114	0.25 ± 0.02	0.32 ± 0.01
	Experiment	HTV-4000	114	0.20 ± 0.02	0.34 ± 0.04

---

Fig. 7 The number of wrinkles and the critical displacement at the onset of wrinkling instability in the frictionless case. Details of each data point are summarized in Table 2. (a) Number of wrinkles m as a function of a/ℓ_g. For a/ℓ_g ≳ 20, all data points are located near m = 7. (b) Relationship between the critical displacement d_w and sheet thickness h. Color indicates the value of a/ℓ_g. The dashed lines correspond to the expression given in eqn (30) and (31).

As a/ℓ_g decreased, the number of wrinkles tended to decrease to five and the deviation in Fig. 7(b) became more pronounced. This trend is likely due to the finite-size effects of the sheet, which were not taken into account in the theoretical analysis above. We will return to this point in Section 7.

6.2 Frictional case (μ > 0)

When frictional interactions occur between the sheet and the substrate, d_w/h and m are no longer constant and they increase monotonically with μ as shown in Fig. 8. Qualitatively, this can be understood as follows: friction suppresses |u_r(R)|, which in turn reduces |ε_θθ(R)| compared to the frictionless case at the same d; therefore, a larger d is required to destabilize the axisymmetric state.^16,44

Fig. 8 The number and critical displacement of wrinkles in the frictional cases, obtained from FES for varying friction coefficients μ, keeping all other parameters fixed. Both the critical number of wrinkles m [panel(a)] and the rescaled critical displacement d_w/h [panel(b)] increase monotonically with μ. The dotted lines in (a) and (b) indicate those without friction μ = 0, given in eqn (30) and (31). The inset images show the wrinkled shapes in FES for μ = 0.01 and μ = 2.0. The parameters used are those of “Elite double 8” given in Table 2.

As shown in eqn (26) and (27), the contribution of friction to the hoop stress can be estimated as ∼μρga, whereas its effect on the radial stress appears only at a higher orders. By normalizing the contraction force h|σ_θθ| with the characteristic stress scale , we find the relative frictional stress:

which is a dimensionless measure of the frictional contribution to the hoop stress. As we demonstrate below, this parameter τ plays a central role in characterizing the frictional effects in our problem.

We now refine the preceding scaling relations in eqn (30) and (31), by including frictional effects up to the first order in τ. The scaling relation d_w/h ∼ m², obtained by balancing the bending force and the radial tension, is expected to remain valid even for μ > 0, because σ_rr is insensitive to μ. Indeed, FES performed using the parameter sets listed in Table 2 confirms that the scaling prediction

agrees well with the FES results [Fig. 9(a)]. Here, the numerical coefficient 0.33 is determined within the scaling argument for μ = 0; that is, (d⁰_w/h)/(m⁰)² ≈ 0.33 based on eqn (30) and (31).

Fig. 9 Relationships among the number of wrinkles m, the normalized critical displacement d_w/h, and the normalized frictional stress τ defined in eqn (32). Details of the symbols are summarized in Table 2. (a) Relationship between d_w/h and m. The dashed line corresponds to the expression given in eqn (33). The color bar represents the magnitude of τ. (b) and (c) The number of wrinkles m and the normalized displacement d_w/h for μ > 0, plotted as functions of τ. For small τ, the rescaled data are captured by the master curves described by eqn (34) and (35), with the common fitting parameter c = 0.14. Data over a wide range of τ are also shown in the insets of these figures.

Based on the balance between radial and hoop stresses in eqn (29), we obtain m ∼ (σ_rr(R)/|σ_θθ(R)|)^1/2 ∼ 1 +O(τ), which is approximately equal to (d_w/h)^1/2 since eqn (33) holds. Therefore, we obtain the following:

m ≈ m⁰(1 + cτ), (34)

where c is a numerical constant common to both the expressions.

To test eqn (34) and (35), we performed additional experiments. In our experiments, visually detecting the presence or absence of wrinkles was challenging, and the force–displacement curve did not provide any information related to wrinkling. Thus, we determined the presence of wrinkles by carefully observing the point-cloud data of the sheet obtained via 3D scanning. From our measured data, we identified the maximum displacement for the axisymmetric state, d⁻_w and the minimum displacement for the wrinkled state, d⁺_w, respectively, and determined the critical displacement in our experiments as d_w ≡ (d⁺_w + d⁻_w)/2. The number of wrinkles m was determined from wave profiles: for several radii r we analyzed the angular dependence of w(r, θ) in the point-cloud data and computed a mean number of wrinkles from m = 2πr/λ(r) for each sample. In our experiments, the wave profiles were occasionally unevenly distributed along the circumferential direction. The irregular wavelengths λ may have originated from small defects, such as impurities introduced during sample fabrication, and/or from residual pre-stress induced when the sheet was placed on the substrate. Indeed, we did not observe such a irregularities in the wave profiles in our FES. We prepared four to five samples from the same material and with similar thicknesses (thickness variation within 10%, as shown in Table 2). We performed the above analysis for each sample, and calculated the mean and standard deviation for each group with similar thickness.

Our theoretical predictions, experimental data, and FES results are summarized in Fig. 9. The parameter sets investigated in the experiments and FES, together with the corresponding symbols, are listed in Table 2. The dashed lines represents the theoretical predictions of eqn (33)–(35), which consist of three dimensionless parameters m⁰, d⁰_w/h, and c. We fixed m⁰ and d⁰_w/h as the frictionless values given in eqn (30) and (31), respectively. The remaining parameter c was determined by fitting the FEM data in Fig. 9(b) to eqn (34) for τ < 2, yielding c ≈ 0.14. We use this value of c for the dashed curve in Fig. 9(c). The scaling plots of our data for a wide range of values of E, h, and μ in Fig. 9 supports our predictions, eqn (34) and (35), up to moderately large values of τ. The prominent deviations of some experimental results from the FES predictions are likely due to small fabrication defects in the samples, and possibly residual pre-tension. In addition, as τ increases, the deviation of each markers in our FEM data becomes prominent [see the insets of Fig. 9(b) and (c)]. This trend likely reflects the limitations of the present first-order-in-τ description (see SI Section IV).

6.3 Relative frictional stress

The dependence of m and d_w/h solely on τ reflects the fact that our problem is fully characterized by the two dimensionless parameters d/h and τ.⁴⁹ To make this explicit, we choose a non-dimensionalization that does not contain d:These scaling factors can also be obtained from eqn (10) or (18) with d = h. With this non-dimensionalization, we find that eqn (4) becomesIt can be seen that ℓ_g is the essential length scale in this problem, since normalizing all lengths by ℓ_g scales out the FvK equations.

In the absence of friction, since eqn (37) does not contain any parameters; only the boundary condition,

characterizes the physical behavior of the sheet. We thus expect d_w/h ∼ O(1) at the instability, and consequently m ∼ (d_w/h)² ∼ O(1).

In the frictional case, we focus the compressive hoop stress term in eqn (37), as it is most influenced by the friction. Using eqn (26) and (27), we find

which contains an additional parameter, τ. Accordingly, the problem is governed by two parameters, (d/h,τ).

The relative frictional stress τ is reminiscent of the parameter referred to as “mechanical bendability” in the previous studies.^41,49,58 Note, however, that while the mechanical bendability in these works is defined as the ratio of the radial tension to the bending force, τ instead quantifies the frictional contribution to the hoop stress. In the problem considered in the above references, the bendability is so large that the physically relevant regime is far from threshold, where the stress state is well described by tension-field theory.^53,58,59 In that limit, the hoop stress is asymptotically small compared with the radial stress, i.e., |σ_θθ|/σ_rr → 0. In contrast, our study focuses on a moderately bendable regime, in which σ_rr and σ_θθ remain of the same order, and the contribution of τ to σ_θθ in eqn (39) is subdominant. Therefore, our problem lies in a near-threshold regime, distinct from the highly bendable limit, and is more closely related to problems of wrinkling in the indentation of moderately thick floating sheets⁴⁹ and in liquid blister tests.³⁹

7 Global buckling

In this section, we focus on the instabilities induced by large indentation, in particular the morphological transitions leading to global buckling.

In the post-wrinkled state, the lifted region continues to grow with the indentation height d, as described in eqn (1). As shown in Fig. 10(a), the initially well-ordered m-fold wrinkles become distorted as d increases. Eventually, three to four dominant wrinkles grow preferentially, while the others disappear, resulting in a polygonal lifted shape with typically m = 3–4. The tips of these dominant wrinkles extend to the end of the sheet, as indicated by the orange arrow in Fig. 10(a-iii). When the hoop stress at r ∼ a, which is estimated as |σ_θθ(a)| ∼ E|ε_θθ(a)| ∼ E(d/a)², exceeds a critical threshold, localized buckling occurs near the sheet's periphery,³⁹ as shown in Fig. 10. This peripheral buckling is accompanied by discontinuous changes in force and stress,⁴⁰ as seen in Fig. 2 and 10. In both experiments and FES, we find that multiple wrinkle tips sometimes buckle simultaneously, whereas in other cases a single ruck forms first (Fig. 10).

Fig. 10 Evolution of the deformation from (i) the seven-fold symmetric wrinkle pattern to (ii) a distorted, less-ordered pattern, (iii) a triangular-like shape, and finally (iv) the post-buckling configuration. (a) Shape transition. (b) and (c) Radial and azimuthal components of the stress distribution. The color indicates the value of dimensionless stress. The dotted circle denotes the radius R calculated by eqn (1) and (23). All results are obtained from FES with parameters (E, h, a, μ) = (226 kPa, 0.5 mm, 150 mm, 0).

Here, we assume that the energy of each dominant wrinkle is nearly identical before global buckling, and estimate the energy per unit area of a single ruck formed at the periphery as:

where A and λ are the amplitude and width of the ruck, respectively. Here, we take the pre-buckling state as the reference, with e_pre = 0, corresponding to A = 0. Assuming conservation of the arc-length along the azimuthal direction at r = a, A is related to λ via A ∼ |ε_θθ|^1/2λ; therefore e is a function of the λ only. By minimizing e with respect to λ, the characteristic wave length of the ruck is obtained as follows:

λ ∼ ℓ_gb|ε_θθ|^1/6, (41)

where ℓ_gb ≡ (B/ρgh)^1/3 is the gravito-bending length,⁵⁰ which is distinct from ℓ_g = (B/ρg)^1/4 as noted in Section 3. The minimized energy for the case A≠0 is obtained substituting eqn (41) into eqn (40) asThe buckling threshold is determined by the requirement that the formation of a ruck reduces the energy: e_buckled < e_pre. This condition can be summarized as d_c ∼ a(ρgh/E)^1/4, or in dimensionless form:We plot d_c/h as obtained from the experiments and FES as a function of a/ℓ_g in Fig. 11. As predicted, d_c/h and a/ℓ_g exhibit an approximately proportional relationship, with a coefficient nearly equal to 1.

Fig. 11 Critical displacement of global buckling. The results obtained from FES and experiments are represented by open and filled symbols, respectively. The markers , *, and + indicate the results from FES without friction. The color bar indicates the magnitude of τ. The solid line represents d_c/h = a/ℓ_g, which corresponds to eqn (43). The dashed line denotes the value of a/ℓ_g = 16.3 given in eqn (44).

7.1 Few-mode wrinkle formation driven by global buckling

When various sheet parameters are fixed, and only the thickness is increased, or the radius is decreased, the critical displacement for global buckling d_c, eventually becomes smaller than the wrinkle-onset threshold d_w. The condition d_w > d_c appears to suggest that global buckling can occur directly from the axisymmetric state without the intermediate formation of wrinkles. Neglecting friction for simplicity, and using the numerical coefficients indicated by the dotted lines in Fig. 7 and 11, the condition d_w > d_c can be written as

a ≲ 16.3ℓ_g. (44)

The colormap in Fig. 12 shows the number of wrinkles obtained from FES at the earliest onset of instability prior to global buckling. Fig. 12(a) indicates that eqn (44) acts as a boundary that separates the phases with fewer (possibly only transient) wrinkles (typically m = 5) and those with statically formed wrinkles (m ≥ 7). We emphasize that this trend is consistent with the behavior observed in Fig. 7 for the range a/ℓ_g ≲ 20. We classify the wrinkling observed in the regime a/ℓ_g ≳ 16.3—which has been the focus of this study—as Type-I Wrinkling, and define the other regime as Type-II Wrinkling. Type-II wrinkles were observed only briefly, immediately before the onset of global buckling, and disappeared as soon as the periphery of the sheet buckled. Therefore, Type-II wrinkles are presumed to form as a transient feature during the direct transition from the axisymmetric state to the globally buckled state, and are likely governed by a mechanism distinct from that of the Type-I wrinkles discussed in Section 6.

Fig. 12 (a) Phase diagram classifying the observed wrinkling patterns into Type-I (with seven or more folds) and Type-II (with fewer folds). Details of the symbols are summarized in the inset of Fig. 11. The solid line indicates the equality condition in (44). The CG illustrations generated from FES show the transition of the shape from a wrinkled to a globally buckled state. (b) and (c) Force–displacement curves obtained from finite-element simulations with varying μ, plotted for two cases: (b) a transition from Type-I wrinkling to global buckling, and (c) a transition from Type-II wrinkling state to global buckling.

As the sheet becomes even thicker and/or its radius decreases, physically different behaviors may also emerge, including a simple case in which the sheet is trivially lifted off without any buckling at all. A detailed investigation of these phenomena, including the transient Type-II wrinkles, is beyond the scope of the present study, and is left for future works.

7.2 Dependence of friction

In the preceding argument, which leads to eqn (43), we neglected the contribution of friction. However, a weak dependence on the friction coefficient can be observed in Fig. 11. For each τ, the deviations from eqn (43) in both the experiments and FES exhibit qualitatively similar trends.

Fig. 12(b) and (c) show the force–displacement curves obtained for sheets in the Type-I and Type-II phases, respectively, for different values of μ. We clearly see that while d_c in Type I is significantly affected by friction, d_c in Type II is rather insensitive to μ. This might be reasonable, considering that, whereas Type-II wrinkle is the only transient pattern appearing just before the global buckling, Type-I wrinkle is a static stable structure formed at the preceding stage of the global buckling instability. During the emergence of the Type-I wrinkles, the contact area and geometry of the sheet on the frictional surface may become highly complicated, which could substantially influence the onset of the global buckling. Indeed, as shown in Fig. 11, the data exhibit scatter for a/ℓ_g > 16.3, where the deviation of d_c/h from the scaling prediction is larger for larger values of τ.

8 Hysteresis

Beyond the global buckling, the shape of the sheet undergoes large geometrically nonlinear deformations, self-contact, and folding. Fig. 13 presents the experimental and FES data for a protocol in which d is increased and then reduced to d = 0, in a quasi-static manner. In such a cyclic indentation test, we observe dramatic shape changes, including ruck formation and folding. Whether these remain permanent after the load is removed depends on the maximum indentation height, d_max. The morphological process of the buckled sheet is closely related to the magnitude of hysteresis observed in the force curves in Fig. 13. Specifically, we performed both experiments and FES to investigate the following three cases: (A) d_max < d_c; (B) d_max is slightly larger than d_c; and (C) d_max ≫ d_c. Below, we describe the three qualitatively distinct behaviors of the sheet observed in the experiment and FES.

Fig. 13 (a) Force–displacement curve for the one cyclic test. We changed the maximum indentation height, d_max. The blue, green, and red colors represent the indentation forces for d_max < d_c, d_max > d_c, and d_max ≫ d_c, respectively. The open symbols and solid lines represent experimental and numerical data, respectively. The illustrations of sheets are obtained from FES. (b) von Mises stress distribution in the inset (vii) of panel (a), calculated in Abaqus. Stress focusing is observed in some areas. All results are obtained from the experiment and FES with parameters (E, h, a, μ) = (477 kPa, 0.38 mm, 104 mm, 0.32).

For case (A), the global buckling of the sheet is absent because d_max < d_c, and no appreciable hysteresis in the force curve is observed, confirming that the cycle is reversible even in the presence of the Type-I wrinkles.

For case (B), a sheet is deformed beyond the global buckling point, and its corresponding force curve exhibits pronounced hysteresis. As d increases, significant stretching is stored in the pre-buckled state [see Fig. 13(a-i)], and this in-plane strain is released abruptly at the onset of global buckling [see Fig. 10 and 13(ii)].^35,36,39,40 Because the stress state in the sheet is largely different, the shape change behaviors can naturally differ between the ascending and descending processes (iii–iv). However, in this regime, the rucks subsequently slip and unfold as d is reduced²⁵ such that no rucks or folds remain once d becomes sufficiently small. The force–displacement curve then gets back to the same path as that in the increasing process (v),⁴⁰ and at d = 0 the sheet returns to its initial flat configuration. In this sense, the cycle process in (B) may be characterized by a single closed loop in the configurational space.

In contrast, for case (C), the sheet is lifted well beyond the global buckling (d_max ≫ d_c, and (vi) in Fig. 13(a)), we observe the rucks remain even when d decreases to zero.²⁵ The rucks formed are initially symmetric, but they often collapse to one side, either the left or right (by spontaneous symmetry breaking), as previously observed in heavy elastica,^50,60 leading to a partially folded sheet configuration. Such a self-folded ruck exhibits a stress focused structure³² as shown in Fig. 13(b). Moreover, these behaviors correspond to a seemingly mysterious discontinuous increase in indentation force F during the descent process (vii). Some of the residual rucks and folds prevent the sheet from returning to its initial flat configuration and remain permanently owing to the coupling among the intricate 3D sheet geometry, self-contacts, and frictional interactions. The full sequence is also shown in SI Movie, where the experiment and the FES are shown together for case (C) in Fig. 13.

Although these are the typical behaviors in cases (A)–(C), the sheet can behave in more complicated ways depending on the magnitude of friction. For example, unfolding events such as those observed in case (B) occur only when the surface friction is sufficiently weak. In addition, folds produced by collapsed rucks may either slide along the sheet surface while remaining in self-contact⁶¹ or lift back up into rucks during unloading. Furthermore, for large values of τ, irreversible shape change of the sheet can be observed even for d_max < d_c. Understanding and quantitatively classifying such morphological diversity is far from straightforward, and further investigations will be needed to uncover the underlying mechanisms.

9 Conclusion

In this study, we combined analytical theory, experiments, and numerical simulations to investigate the indentation response of a heavy elastic sheet resting on a frictional substrate. We uncovered the fundamental principles governing pattern formation under indentation, including nonlinear and discontinuous force responses, wrinkling instabilities, and global buckling. A characteristic length scale, ℓ_g, was identified, which not only controls the size of the uplifted region and the force response but also emerges consistently in the onset of wrinkling and the transition to global buckling. In the frictionless case, we discovered a remarkably simple rule: the number of wrinkles formed is constant (m ≈ 7) and the critical indentation height depends solely on the sheet thickness (d_w ≈ 16.3h). When friction is present, the wrinkling behavior becomes more complex; however, it can still be characterized in terms of another dimensionless parameter, the relative frictional stress τ. Together, these results provide a quantitative framework for the indentation-induced morphology of heavy sheets on frictional substrates.

One might suppose that the post-buckled shapes observed in our experiment resemble the developable configurations previously studied for suspended sheets under gravity.¹ However, these systems differ fundamentally in their internal stresses and relevant boundary conditions. Obviously, a draped sheet is force- and torque-free at its outmost edge. In contrast, because the globally buckled sheet studied here is still partially in contact with the substrate, it is generally subjected to more complicated boundary conditions at its edge owing to the residual in-plane strains within the sheet. Indeed, we have frequently observed a discontinuous transition to a developable surface in our experiment, which occurs when the outmost edge of a buckled sheet fully detaches from the substrate. The resulting developable configuration again contains wrinkles, but their number generally does not coincide with that in the globally buckled phase. This difference likely explains the discontinuity of the transition, although further investigations will be needed.

Thus far, we have discussed only the near-threshold behavior of the wrinkling instability. As shown in Fig. 10, well-developed wrinkles become distorted, with only a few dominant wrinkles growing preferentially. The linear stability analysis cannot explain this collapse or the process by which the number of dominant wrinkles is determined. Furthermore, our findings in eqn (34) and (35) apply only in the regime of relatively weak friction. As seen in Fig. 9, when τ ≳ 10, m and d_w/h are no longer characterized solely by τ. The materials and sheet sizes used in our experiments do not reach τ > 10 for typical friction coefficients. However, if we are to apply our study to extremely large systems, such as those in aerospace engineering⁷ or geophysical contexts,⁶² our theory should be extended to account for the high-τ regime, since τ is proportional to system size.

Finally, despite the extreme simplicity of the loading protocol, the sheet exhibits a remarkably rich set of configurations ranging from axisymmetric uplifts and wrinkle patterns to global buckling. Upon unloading from sufficiently large lifts, we observed a folded configurations stabilized by self-contact and friction, as shown in Fig. 13 and in the “ghost” in Fig. 1(a). The resulting structures can be viewed as gravity-driven self-organizing origami.^63–65 Our findings may contribute to the controlled design of well-ordered wrinkles or artistic complex folds through contact interactions under gravity, fluid drag, or high pressure.

Author contributions

Keisuke Yoshida: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, validation, visualization, writing – original draft, writing – review & editing. Hirofumi Wada: conceptualization, formal analysis, funding acquisition, project administration, resources, supervision, validation, writing – review & editing.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data that support the findings of this article are available within the article and within its supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6sm00062b.

Acknowledgements

We acknowledge financial support from JSPS KAKENHI (Grant No. 23K22463 to H. W.); a Grant-in-Aid for JSPS Research Fellows (DC1; Grant No. 21J22837 to K. Y.); the Early-career Researcher Development Program of Ritsumeikan University (to K. Y.); and a research grant of 2025 (during 2025-2026) from Amano Institute of Technology (K. Y.). ChatGPT (OpenAI) was used for language editing and minor assistance with data-processing scripts. The authors take full responsibility for the content and results.

References

1. E. Cerda, L. Mahadevan and J. M. Pasini, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 1806–1810.
2. J.-S. Chen and Y.-Y. Fang, Int. J. Solids Struct., 2010, 47, 2767–2774.
3. G. Boedec and J. Deschamps, Int. J. Solids Struct., 2021, 233, 111214.
4. H. Vandeparre, M. Piñeirua, F. Brau, B. Roman, J. Bico, C. Gay, W. Bao, C. N. Lau, P. M. Reis and P. Damman, Phys. Rev. Lett., 2011, 106, 224301.
5. M. Taffetani, F. Box, A. Neveu and D. Vella, Europhys. Lett., 2019, 127, 14001.
6. J. D. Johnston, J. R. Blandino and K. C. McEvoy, J. Spacecr. Rockets, 2006, 43, 762–770.
7. D. A. Spencer, L. Johnson and A. C. Long, Aerosp. Sci. Technol., 2019, 93, 105276.
8. G. Karner and A. Watts, J. Geophys. Res.: Solid Earth, 1983, 88, 10449–10477.
9. L. Mahadevan, R. Bendick and H. Liang, Tectonics, 2010, 29, TC6002.
10. S. Protière, C. Josserand, J. M. Aristoff, H. A. Stone and M. Abkarian, Phys. Rev. Lett., 2017, 118, 108001.
11. E. Cerda and L. Mahadevan, Phys.Rev. Lett., 2003, 90, 074302.
12. E. P. Chan, E. J. Smith, R. C. Hayward and A. J. Crosby, Adv. Mater., 2008, 20, 711–716.
13. J. Y. Chung, A. J. Nolte and C. M. Stafford, Adv. Mater., 2011, 23, 349–368.
14. T. Yamaguchi, S. Ohmata and M. Doi, J. Phys.: Condens. Matter, 2009, 21, 205105.
15. Y. Aoyanagi, J. Hure, J. Bico and B. Roman, Soft Matter, 2010, 6, 5720–5728.
16. A. Chawla and D. Kumar, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2320068121.
17. R. H. Plaut, S. Suherman, D. A. Dillard, B. E. Williams and L. T. Watson, Int. J. Solids Struct., 1999, 36, 1209–1229.
18. B. Roman and A. Pocheau, J. Mech. Phys. Solids, 2002, 50, 2379–2401.
19. N. Stoop, F. K. Wittel and H. J. Herrmann, Phys. Rev. Lett., 2008, 101, 094101.
20. C.-W. Liu and J.-S. Chen, Eur. J. Mech. A-Solid., 2013, 39, 50–59.
21. S. Alben, Proc. R. Soc. A, 2022, 478, 20210742.
22. S. Deboeuf, S. Protière and E. Katzav, Phys. Rev. Res., 2024, 6, 013100.
23. L. Mahadevan and J. B. Keller, Proc. R. Soc. A, 1996, 452, 1679–1694.
24. M. Habibi, N. Ribe and D. Bonn, Phys. Rev. Lett., 2007, 99, 154302.
25. D. Vella, A. Boudaoud and M. Adda-Bedia, Phys. Rev. Lett., 2009, 103, 174301.
26. J. M. Kolinski, P. Aussillous and L. Mahadevan, Phys. Rev. Lett., 2009, 103, 174302.
27. T. G. Sano, T. Yamaguchi and H. Wada, Phys. Rev. Lett., 2017, 118, 178001.
28. P. Grandgeorge, T. G. Sano and P. M. Reis, J. Mech. Phys. Solids, 2022, 164, 104885.
29. M. Tani and H. Wada, Phys. Rev. Lett., 2024, 132, 058204.
30. G. K. Curtis, I. M. Griffiths and D. Vella, Int. J. Solids Struct., 2026, 326, 113702.
31. B. Audoly and Y. Pomeau, Elasticity and geometry: from hair curls to the non-linear response of shells, Oxford University Press, 2010.
32. T. A. Witten, Rev. Mod. Phys., 2007, 79, 643.
33. J. Hure, B. Roman and J. Bico, Phys. Rev. Lett., 2011, 106, 174301.
34. J. Hure, B. Roman and J. Bico, Phys. Rev. Lett., 2012, 109, 054302.
35. T. Suzanne, J. Deschamps, M. Georgelin and G. Boedec, Phys. Rev. E, 2022, 106, 065002.
36. L. Stein-Montalvo, A. Guerra, K. Almeida, O. Kodio and D. P. Holmes, Phys. Rev. E, 2023, 108, 035002.
37. M. Krieg, G. Fläschner, D. Alsteens, B. M. Gaub, W. H. Roos, G. J. Wuite, H. E. Gaub, C. Gerber, Y. F. Dufrêne and D. J. Müller, Nat. Rev. Phys., 2019, 1, 41–57.
38. D. Akinwande, C. J. Brennan, J. S. Bunch, P. Egberts, J. R. Felts, H. Gao, R. Huang, J.-S. Kim, T. Li and Y. Li, et al. , Extreme Mech. Lett., 2017, 13, 42–77.
39. J. Chopin, D. Vella and A. Boudaoud, Proc. R. Soc. A, 2008, 464, 2887–2906.
40. D. P. Holmes and A. J. Crosby, Phys. Rev. Lett., 2010, 105, 038303.
41. D. Vella and B. Davidovitch, Phys. Rev. E, 2018, 98, 013003.
42. H. Dannenberg, J. Appl. Polym. Sci., 1961, 5, 125–134.
43. Z. Dai, Y. Hou, D. A. Sanchez, G. Wang, C. J. Brennan, Z. Zhang, L. Liu and N. Lu, Phys. Rev. Lett., 2018, 121, 266101.
44. Z. Dai, D. A. Sanchez, C. J. Brennan and N. Lu, J. Mech. Phys. Solids, 2020, 137, 103843.
45. V. L. Popov, Contact mechanics and friction, Springer, 2010.
46. K. L. Johnson, Contact mechanics, Cambridge University Press, 1987.
47. S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, McGraw-Hill, 1959.
48. D. A. Dillard, B. Mukherjee, P. Karnal, R. C. Batra and J. Frechette, Soft Matter, 2018, 14, 3669–3683.
49. F. Box, D. Vella, R. W. Style and J. A. Neufeld, Proc. R. Soc. A, 2017, 473, 20170335.
50. C. Wang, Int. J. Mech. Sci., 1986, 28, 549–559.
51. D. Vella and B. Davidovitch, Soft Matter, 2017, 13, 2264–2278.
52. R. C. Benson, J. Appl. Mech., 1991, 58, 484–492.
53. E. H. Mansfield, The bending and stretching of plates, Cambridge University Press, 2nd edn, 1989.
54. Z. Dai and N. Lu, J. Mech. Phys. Solids, 2021, 149, 104320.
55. D. Systèmes, Abaqus Version 6.6 Documentation, 2006, https://classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/v6.6/index.html.
56. T. Baumberger and C. Caroli, Adv. Phys., 2006, 55, 279–348.
57. S. Timoshenko and J. Goodier, Theory of Elasticity, McGraw-Hill, 1969.
58. B. Davidovitch, R. D. Schroll, D. Vella, M. Adda-Bedia and E. A. Cerda, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 18227–18232.
59. D. Vella, Nat. Rev. Phys., 2019, 1, 425–436.
60. G. Domokos, W. Fraser and I. Szeberényi, Phys. D, 2003, 185, 67–77.
61. V. Démery, B. Davidovitch and C. D. Santangelo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 042401.
62. D. L. Turcotte and G. Schubert, Geodynamics, Cambridge university press, 2002.
63. L. Mahadevan and S. Rica, Science, 2005, 307, 1740.
64. S. Felton, M. Tolley, E. Demaine, D. Rus and R. Wood, Science, 2014, 345, 644–646.
65. T.-H. Kim, D.-Y. Lee and J.-H. Han, Sci. Rep., 2025, 15, 19615.

---