The broad wrinkling landscape of hyperelastic parallelogram-shaped membranes: from wrinkle migration to restabilization and their subsequent reappearance elsewhere

Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 The broad wrinkling landscape of hyperelastic parallelogram-shaped membranes: from wrinkle migration to restabilization and their subsequent reappearance elsewhere M.H. NejabatMeimandi1,2 and F. Dal Corso∗1 1DICAM, University of Trento, via Mesiano 77, I-38123 Trento, Italy 2Tensys Ltd, 122 Wells Rd, Bath, UK December 11, 2025 Abstract Wrinkling is a commonly observed out-of-plane instability in membrane structures due to their extremely low bending-to-stretching stiffness ratio. It has been extensively investigated for symmetric membrane geometries and boundary conditions that induce planar non-uniform stress states by preventing the lateral contraction at the edges, and is also known to potentially display self-restabilization. This study investigates an initially flat, parallelogram-shaped hyperelastic membrane, focusing on the effect of the inclination angle that defines its deviation from rectangular geometry. It is shown that wrinkling can occur either centrally or at the two opposite obtuse-angled corners—even for small inclination angles—during stretching with unconstrained lateral contraction, a condition under which the flat configuration for the rectangular counterpart remains always stable. Three distinct evolutions of the wrinkling pattern are numerically identified, all ultimately leading to cornerlocalized wrinkles. This final state may arise (i) directly, without a prior bifurcation, or after the appearance of central wrinkling that either (ii) restabilizes or (iii) separates and migrates toward the corners. A closed-form expression for the critical wrinkling condition is derived by combining a perturbation approach with an energybased method in the framework of linear elasticity. This provides an accurate estimate of the onset and pattern of central wrinkling. The present findings reveal new pathways in wrinkling pattern evolution and introduce a novel approach to unconventional boundary-value problems, with potential applications ranging from lightweight structural systems to flexible electronics. Keywords: Self-restabilization; quasi-rectangular membranes; wrinkling pattern morphing; perturbation approach. 1 Introduction Thin membranes and films are highly flexible bi-dimensional structures abundant in nature (cells, insect wings, and leaves) and widely manufactured for diverse engineering applications—including architecture, aerospace, electronics, and medicine—due to their lightweight and adaptable nature. These structures offer mechanically efficient and aesthetically appealing solutions for applications spanning from covering of large spaces [12] to hosting electronic components in flexible devices [33]. Their extremely low bending-to-stretching stiffness ratio however facilitates the onset of structural elastic instabilities, posing a significant challenge by compromising both structural integrity and visual appeal. Among the possible instabilities, membranes are particularly prone to wrinkling, exhibiting a sinusoidal-like out-of-plane displacement on the surface, with a short wavelength aligned to the direction of the principal compressive stress, creating localized curvature. Wrinkling is commonly experienced in everyday life with items like plastic wraps, clothing fabrics, curtains, and balloons. ∗Corresponding author: [email protected] 1 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 State of the art: Tension-field theory, initially developed in 1929 [57] and later extended [45], has been the foundational framework for analyzing wrinkling instability in thin-walled structures. This two-dimensional nonlinear model assumes that wrinkles have infinitesimal wavelengths, with the material being considered to have no compressive bearing capacity. However, despite several enhancements [50, 9] and generalizations [14, 49, 42], tension-field theory still faces limitations, particularly in its neglect of bending energy [13]. Consequently, researchers have increasingly turned attention to more advanced plate and shell models that account for both stretching and bending effects, such as the F¨oppl–von K´arm´an model and its extended version, the latter taking into account the large in-plane deformations of the membrane. Through these models, the analysis of thin sheets gained more precision also in terms of predicting the post-critical response and the corresponding wrinkling morphology. Deep insights into the wrinkling mechanisms have then been gained in several different setups, mainly restricted to axial or axis symmetry conditions. As expected, the position of wrinkling within a membrane changes with variations in the membrane geometry and loading conditions. Wrinkling may be spread over the whole membrane, as under shear loading [64], or emerges as an isola-center bifurcation in the two-axis symmetric problem of a stretched rectangular membrane with constrained lateral contraction at the edges [27, 52, 58]. Localized wrinkling can also occur near the membrane’s free edge, with wrinkles oriented orthogonally to the boundary, as observed in elastic plate stamping [26], elastically supported, prestressed incompressible isotropic plates [15], spinning elastic membranes [10], and membranes under pressure loading [11]. Alternatively, wrinkles may appear close to and align parallel to the free edge, particularly in twisted, pre-stretched membranes [58]. Edge wrinkling has also been investigated as instability patterns in growing curled petals and leaves [59]. In the context of inelastic membrane behaviour, the Mullins effect has been shown to induce a distinctive wrinkling response: while no wrinkles appear during the initial loading, they emerge during the first unloading and persist throughout all subsequent loading cycles [17]. Drawing an analogy with the restabilization of the trivial path observed in variable-length rods under compression [2, 3], both theoretical and experimental investigations of highly stretched, initially isotropic thin sheets have revealed a surprising restabilization phenomenon [20, 21, 23, 34, 41, 48, 53, 69, 70]. In this context, wrinkle amplitude initially increases with stretching but subsequently diminishes as stretching continues, eventually vanishing entirely—thereby demonstrating a recovery of stability of the flat configuration. Such a restabilization phenomenon has been also addressed in orthotropic [35, 47] and anisotropic [22] membranes, as well as in soft shells [62]. A recent analysis for anisotropic membrane has also disclosed the possibility of wrinkling reappearance at the same central location after its disappearance due to restabilization [8]. Wrinkling instability has been also studied in membranes attached to a substrate, showing wrinkling patterns in trapezoidal film/substrate bilayers [68], characterized by period doubling [67], with axisymmetric/diamond-like mode transition [66], with hexagonal geometry in spheres and toroids with an elastic core [4, 29, 51, 61], with smooth-wrinkle-ridge-sagging transitions [16], with crystallography on spherical surfaces [4], and in differentially growing bilayers [46]. Owing to the complexity of stress fields within membranes and the potential involvement of nonlinearities, only a limited number of studies have produced closed-form analytical solutions for predicting the critical conditions for wrinkling and the resulting patterns. The observation of stretch-induced wrinkling in rectangular sheets [18] thus motivated simplified buckling models for complex stress states lacking direct analytical solutions. Scaling laws for wrinkle wavelength and amplitude have been established through energy minimization [6, 7, 38, 55], and energybased models for wavelength selection have been proposed [28]. Furthermore, the dependence of critical wrinkling strain on aspect ratios has validated scaling relationships between applied stretch and wrinkling behaviour [7, 30, 44]. Finally, the transition between periodic wrinkling and global buckling has been analytically predicted for a thin elastic ring bound to an equally curved 2D substrate that contains an inner cavity [32]. The deep understanding of wrinkling mechanisms unlocked innovative applications across several technological fields. Recent advancements have leveraged the unique properties of wrinkling in auxetic membranes and nematic elastomer sheets to achieve on-demand, non-standard wrinkling patterns, suppressing unwanted instabilities through microstructure formation [43, 56]. These developments are critical for lightweight deployable space structures, such as solar sails and antennas, where precise control of membrane behaviour enhances performance 2 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 [27]. In terrestrial applications, wrinkling informs the design of complex fabric roofs and smart mechanical devices, including dielectric elastomeric layers and bilayer gel beams, where controlled instabilities enable structural morphing, sensing, and actuation in biomedical and soft robotics contexts [40]. For a further detailed overview of tension-induced film wrinkling, readers are directed to the review by Wang et al. [63]. Article contribution: The mechanics and wrinkling instability of an unconventional class of two-dimensional structures—specifically, initially flat, parallelogram-shaped hyperelastic membranes—are investigated under a stretching process that allows lateral contraction. In comparison with the usual setups analyzed in wrinkling problems, besides removing the edge constraint of zero lateral contraction, the investigated boundary-value problem is also less symmetric, as it has a rotational symmetry of order 2, and it allows refined tuning of the inhomogeneity level of the mechanical fields through the inclination angle θ0, which quantifies the deviation of the undeformed parallelogram shape from the rectangular geometry (θ0= 0). The inhomogeneous fields and wrinkling phenomena are investigated for varying initial inclination angle θ0and the in-plane aspect ratio β. By contrast to the limiting case of a rectangular membrane allowed to laterally contract under elongation—for which the stress field remains purely uniaxial under tension, and no compressive stresses arise during elongation, thereby precluding instability—it is shown that introducing an inclination angle generates a compressive principal stress component, which can induce wrinkling even for small deviations (|θ0|≪1) from the rectangular configuration. Finite Element analyses reveal that when wrinkling appears, it develops either in the central region of the membrane or close to the two obtuse-angled corners, Fig. 1. More specifically, three distinct evolutions of the wrinkling pattern are identified, all ultimately leading to the appearance of wrinkles at the corners. The key difference between these evolutions lies in the behaviour prior to this final state: central wrinkling may either not appear at all (case θ0= 2◦in Fig. 1) or manifest during an earlier stage of deformation. In the latter case, it either restabilizes (through a complete wrinkling disappearance, case θ0= 2.9◦in Fig. 1) or separates and migrates toward the two corners (cases θ0= 3◦and 4◦in Fig. 1) as the elongation strain ϵincreases. The numerical analysis is complemented by an analytical treatment within a linear elasticity framework. Using a perturbation approach based on small values of the initial inclination angle θ0, a second-order expansion of the stress field is derived. This analysis identifies the cause of potential central wrinkling formation as a second-order compressive principal stress within the parallelogram membrane. A relatively simple closed-form expression is finally obtained for evaluating the critical elongation strain and the central wrinkling pattern by considering the approximated stress fields within an energy principle. The analytical expression is shown to successfully predict the wrinkling condition disclosed through Finite Element simulations with varying of the parallelogram membrane geometry, Fig. 2. Article outline: The geometrically extended version of the F¨oppl–von K´arm´an model is recalled in Sect. 2, along with the derivation of the constitutive relations for specific choices of classical hyperelastic materials and a description of the membrane geometry and boundary conditions. Details of the Finite Element analyses, numerical results for the planar response, critical conditions, post-critical response, and wrinkling pattern evolution, along with the related discussion, are provided in Sect. 3. The perturbation approach is applied in Sect. 4 to obtain an approximated analytical description of the planar stress fields, which is in turn used into the energy approach to evaluate closed-form expressions for predicting the critical elongation strain of wrinkling and its pattern. The results in terms of both the planar stress field description and the critical deformation and mode are successfully validated through comparison with the numerical predictions. Concluding remarks are finally provided in Sect. 5. Article significance: This study introduces an unconventional class of boundary-value problems in which lateral contraction remains unconstrained, leading to wrinkling instabilities despite the stress fields being less inhomogeneous than those typically observed under classical clamped-edge conditions. Novel evolutions in wrinkling patterns are revealed, including the formation of a second set of wrinkles close to the two opposite obtuse-angled corners of the membrane after the stabilization of the initial central wrinkling. Another possible evolution involves 3 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 µ0[º] w/t 2 3 4 1 2.9 -0.2 -0.1 -0.05 -0.15 0 Figure 1: Maps of the out-of-plane displacement w, normalized through division by the membrane thickness t, showing the wrinkling pattern evolution at four stages of increasing elongation strain ϵ={1.5,7.8,14,30}%. The maps are obtained from Finite Element simulations for five compressible Neo-Hookean membranes with the same width-to-thickness ratio α= 1500 and (in-plane) aspect ratio β= 3, but differing in the initial inclination angle θ0defining the undeformed parallelogram shape. Homogeneous (green) color implies a planar state while inhomogeneous color shows the wrinkling pattern, which appears aligned parallel with the elongation direction. In addition to the case (θ0= 1◦) that shows no out-of-plane displacement at every deformation stage, three different types of evolutions are shown, all ending with wrinkle close to the two obtuse-angled corners. the separation of the central wrinkles and their subsequent migration toward the two opposite corners. These numerical findings are further supported by a closed-form analytical expression that predicts the critical elongation strain leading to central wrinkling, along with the corresponding wrinkling pattern. This research paves the way for designing and tuning the critical conditions and modes of wrinkling, offering valuable insights into how these instabilities can be controlled. It opens up new opportunities for applications across a wide range of fields, from enhancing the performance of lightweight structures to advancing wearable electronic devices. 2 Formulation 2.1 Mechanical model The equilibrium equations governing the large strains of a hyperelastic membrane are derived by following the geometrically extended version of the F¨oppl-von K´arm´an model (eFvK) [19, 21, 23]. The unloaded state of the 4 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 ²cr=1.860% ²cr=0.310% ²cr=0.414% ²cr=0.150% µ0=5 ²cr( =6.7)=n0.313% ²cr( =4.5)=n0.418% ²cr( =4.6)=n0.152% ²cr( =11.0)=n2.010% º µ0=3 ºµ0=3 º µ0=5 º 3, 5, Figure 2: Comparison of the out-of-plane displacement maps, disclosing the wrinkling pattern for β= 2 and θ0= 3◦and 5◦, for β= 3 and θ0= 5◦, and for β= 5 and θ0= 5◦. An excellent agreement is found between the patterns from finite element (FE) simulations (above) and from the analytical approach (bottom) in terms of either the critical elongation strain ϵcr and of the wrinkling pattern. The corresponding analytical prediction of the critical elongation strain ϵcr is also reported by highlighting the minimizing positive value of ndefining the wrinkling ansatz, Eq. (68). membrane is assumed as initially flat and with a geometry described by a prism of domain volume Vwith thickness tand a (flat) base domain Bx, of area Aand for which t≪√Adue to the thin membrane assumption. Describing the undeformed position xin a three-dimensional Cartesian reference system through its coordinates x1–x2–z, the undeformed membrane volume Vis given by V:= x{x1, x2}∈Bx, z ∈[−1,1] t 2.(1) The deformation function f(x)=x+u(x1, x2, z), mapping the reference configuration xinto the current one, is considered to be described by a displacement vector ufollowing the Kirchhoff hypothesis, for which the surface normal to the mid-plane remains normal during deformation, u(x1, x2, z) =      u1(x1, x2, z) u2(x1, x2, z) w(x1, x2)     =    eu1(x1, x2) eu2(x1, x2) w(x1, x2)     −z     w,1(x1, x2) w,2(x1, x2) 0     ,(2) where (the subscript , γ represents the partial derivative ∂/∂xγand) euγ(x1, x2) and w(x1, x2) are the primary kinematic fields, respectively representative of the in-plane and the out-of-plane displacement components of the mid-plane surface along xγ(γ= 1,2) and along z. By introducing the gradient operator ∇={∂/∂x1, ∂/∂x2, ∂/∂z} and the deformation gradient F=∇f=I+∇u, the (symmetric) Green-Lagrange strain tensor E=FTF−I/2 evaluated by considering the displacement vector u(2) finds the following decomposition E(x1, x2, z) = E[p](x1, x2)+zE[χ](x1, x2)−z2 4E[q](x1, x2),(3) where the tensors E[p],E[χ], and E[q](whose components are reported in Appendix A.1) are respectively representative of the constant, linear, and quadratic deformation contributions through the out-of-plane variable z. Considering this decomposition and a hyperelastic response, the (volume) strain energy density ψcan be integrated through the membrane thickness tto achieve the surface strain energy density Ψ as a function of E[p],E[χ], and E[q] Ψ = Zt 2 −t 2 ψ(E) dz= Ψ E[p],E[χ],E[q].(4) 5 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 By assuming an out-of-plane displacement wwith small firstand second-gradients (|∇sw|, t |∇2 sw|≪1), the surface strain energy density Ψ (4) can be approximated according to the extended F¨oppl-von K´arm´an theory as the sum of membrane Ψ[m]E[m]and bending Ψ[b]E[m], κenergy contributions as [48, 60] Ψ≈ΨE[m],κ=t ψ[m]E[m] | {z } Ψ[m](E[m]) +t3 24κ·∂2ψ[m]E[m] ∂E[m]∂E[m]E[m] κ | {z } Ψ[b](E[m],κ) ,(5) where E[m]and κare respectively the (symmetric) membrane Green-Lagrange strain and curvature (or bending strain) tensors, defined as the surface components of E[p]and E[χ] E[m]= E[p] 11 E[p] 12 E[p] 12 E[p] 22  =    eu1,1+(eu1,1)2+ (eu1,2)2+ (w,1)2 2eu1,2+eu2,1+eu1,1eu2,1+eu1,2eu2,2+w,1w,2 2 eu1,2+eu2,1+eu1,1eu2,1+eu1,2eu2,2+w,1w,2 2eu2,2+(eu2,1)2+ (eu2,2)2+ (w,2)2 2     , κ= E[χ] 11 E[χ] 12 E[χ] 12 E[χ] 22  =− w,11 w,12 w,12 w,22 , (6) while ψ[m]is the reduced (volume) strain energy density evaluated through variational dimensional reduction by the partial minimization over the out-of-plane strain components, ψ[m]E[m]= min {E[p] 13 ,E[p] 23 ,E[p] 33 } ψ       E[m]E[p] 13 E[p] 23 E[p] 13 E[p] 23 E[p] 33       .(7) Considering the second Piola-Kirchhoff stress, S=∂ψ(E)/∂E, the partial minimization (7) is mechanically equivalent to solving the through-thickness equilibrium that enforces plane stress conditions [24, 25], which is expressed by Sj3=∂ψ(E) ∂Ej3 = 0, j = 1,2,3.(8) The constitutive relations for the membrane force Nand bending moment M(symmetric) tensors per unit undeformed length1follow from the surface strain energy density Ψ(5), obtained under the extended F¨oppl-von K´arm´an approximation, as N=∂ΨE[m],κ ∂E[m],M=∂ΨE[m],κ ∂κ,(9) which, by considering the membrane and bending contributions (5), reduce to a nonlinear relation for the membrane force Nin the membrane Green-Lagrange strain E[m]and a linear relation for the bending moment Min the curvature tensor κ(with stiffness possibly depending nonlinearly on the membrane Green-Lagrange strain E[m]) NE[m]=t∂ψ[m]E[m] ∂E[m],ME[m],κ=t3 12 ∂2ψ[m]E[m] ∂E[m]∂E[m]κ,(10) where quadratic terms in the curvature κare neglected in the expression for the membrane force N(10)1since t|∇2 sw| ≪ 1. 1It is noted that the resultant membrane force Nand bending moment Mare associated with an undeformed unit length as these are evaluated through an integration of the second Piola-Kirchhoff stress distributions across the membrane thickness. 6 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 By introducing the 2D (membrane) identity tensor Is, the outward unit-normal nand the unit tangent t vectors along ∂Bx, the equilibrium of the extended F¨oppl-von K´arm´an membrane is governed by the indefinite equations (derivation details are deferred to Appendix A.2) (∇s·(∇s·M+N∇sw) = 0, ∇s·[(Is+∇se u)N] = 0.for x1, x2∈Bx,(11) These are complemented by the dual set of stress-free and kinematic boundary conditions on the boundary ∂Bx of the (flat) surface Bx,      n·Mn = 0,or (∇sw)·n=w,n, hN∇sw+∇s·M+t·(∇s(Mt))i·n= 0,or w=w, (Is+∇se u)Nn =0,or e u=u, for x1, x2∈∂Bx,(12) and the dual set of stress-free and kinematic boundary condition at the possible corner point Γjalong ∂Bx (j= 1, ..., Q) Jn·MtK|Γj= 0,or w|Γj=w, j = 1, ..., Q, (13) where (the symbol J·Kstands for the jump value of the relevant argument at the specific corner point, and) u,w, and w,n are the imposed in-plane displacement, out-of-plane displacement and its normal derivative, respectively. 2.2 Hyperelastic constitutive models Compressible Neo-Hookean (NH) model. The strain energy density ψNH associated to the compressible Neo–Hookean (NH) material is defined as [31] ψNH (E) = 1 2(µ"tr (I+ 2E) 3 pdet (I+ 2E)−3#+Khpdet (I+ 2E)−1i2),(14) where µand Kare the ground-state shear and bulk moduli, the latter connected to the ground Lam´e constant λ through K=λ+2µ/3. As the partial minimization (8) performed on the strain energy density ψNH (E) (14) leads to a set of nonlinear equations for the out-of-plane strains E[p] i3E[m](i= 1,2,3), the reduced strain energy density ψ[m] NH E[m]cannot be expressed in a closed form and must be evaluated numerically. Then, nonlinear relations for the membrane force NNH E[m]and bending moment MNH E[m],κfollows through the constitutive relation (10), the latter with a bending stiffness varying with the membrane strain E[m](except in the small-strain limit). Saint Venant–Kirchhoff (SVK) model. The strain energy density ψSVK is given for the Saint Venant– Kirchhoff (SVK) material by [54] ψSVK(E) = 1 2nλ[tr(E)]2+ 2µtrE2o.(15) Being the strain energy ψSVK (E) (15) a quadratic form in the Green-Lagrange strain E, the partial minimization (8) leads to a set of linear equations whose solution is given by the following out-of-plane strain components E13 =E23 = 0, E33 =−λ(E11 +E22) λ+µ,(16) and therefore the reduced strain energy density ψ[m] SVK E[m]follows as ψ[m] SVK E[m]=1 2nλps tr E[m]2+ 2µtr E[m]2o,(17) 7 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 where λps is the Lam´e constant under plane stress conditions, λps = 2λµ/(λ+ 2µ). By considering the relation connecting the Lam´e constants λand µto the Young’s modulus Eand Poisson’s ratio ν, and by introducing the membrane bending stiffness D, λ=E ν (1+ν)(1 −2ν), µ =E 2(1 + ν), D =E t3 12(1 −ν2),(18) the membrane force NSVK and bending moment MSVK reduce to linear functions of the membrane Green-Lagrange strain E[m]and of the curvature κ(with constant bending stiffness D), respectively, NSVKE[m]=E t 1−ν2hνtrE[m]Is+(1−ν)E[m]i,MSVK (κ)=DhνtrκIs+(1−ν)κi.(19) Linear elastic (LE) model. When the displacement gradient is sufficiently small (|∇u| ≪ 1), the higher-order (nonlinear) terms become negligible. It follows that the Green–Lagrange strain Eand its membrane counterpart E[m]reduce to the corresponding linearized strain tensors ϵand ϵ[m], namely E=ϵ=∇u+ (∇u)T/2 and E[m]=ϵ[m]=∇se u+ (∇se u)T/2, and therefore ϵ[m] 11 =eu1,1, ϵ[m] 12 =eu1,2+eu2,1 2, ϵ[m] 22 =eu2,2.(20) In this circumstance, both the compressible Neo-Hookean (14) and Saint Venant–Kirchhoff (17) models reduce to the linear elastic (LE) model by providing the celebrated linear elastic expressions of the F¨oppl-von K´arm´an plate [5, 54] for the membrane force NLE and bending moment MLE NLEϵ[m]=E t 1−ν2hνtrϵ[m]Is+(1−ν)ϵ[m]i,MLE (κ) = DhνtrκIs+(1−ν)κi.(21) 2.3 Geometry and boundary conditions In contrast to the rectangular and circular geometries with edge contraction prevented along portions of their boundaries, which are usually considered to investigate wrinkling, a parallelogram-shaped membrane that is free to contract laterally is considered here, Fig. 3 (left). Differently from the conventional setups, this ‘unconventional’ Figure 3: (Left) Undeformed parallelogram geometry of the hyperelastic membrane and boundary conditions. (Right) Reparameterization of the physical parallelogram domain Bxon the auxiliary unit square domain Bξ through the transformation matrix A(32). configuration allows refined tuning of the inhomogeneity level through the obliquity angle (π/2−θ0) of the parallelogram. In particular, it leads to less inhomogeneous mechanical fields during membrane elongation, as 8 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 it avoids the strong edge-decay effect typically associated with constrained lateral contraction. Moreover, the mechanical fields smoothly converge to the homogeneous solution in the right-angle limit (θ0= 0), where the configuration reduces to a rectangular membrane under uniform tension. The low level of inhomogeneity can be appreciated from Figure 4, where the deformed configuration for an elongation strain ϵ= 44% is reported as result of a planar Finite Element simulation (Sect. 3). Figure 4: Undeformed (gray) and deformed (green) planar configurations for the considered parallelogram-shaped membrane, the latter under boundary conditions (27)3, (28), and (29) described by an elongation strain ϵ. A deformed mesh is reported as a result from Finite Element simulation (Sect. 3) for ϵ= 44%. The current inclination angle θcorresponding to Eq. (37) is also shown as a measure of the (average) inclination of the stressfree edges. The undeformed base domain Bxof parallelogram shape is considered to have two sides parallel to the x2-axis, while the other two are inclined with respect to the x1-axis by the initial (undeformed) inclination angle θ0(Fig. 3, left), Bx:= x|− L 2≤x1≤L 2,−W 2+x1tan θ0≤x2≤W 2+x1tan θ0,(22) where Land Wdenote the undeformed length and width of the original parallelogram flat membrane. For the following analysis, it is instrumental to introduce two dimensionless parameters defining the width-tothickness ratio αand the (in-plane) aspect ratio βas α=W t≫1, β =L W.(23) The boundary ∂Bxof the domain Bxcan be described as the union of four portion boundaries ∂B[a] x,∂B[b] x,∂B[l] x, and ∂B[r] xdefined as ∂B[a] x ∂B[b] x):= x|− L 2≤x1≤L 2, x2=±W 2+x1tan θ0, ∂B[l] x ∂B[r] x):= x|x1=∓L 2,−W±Ltan θ0 2≤x2≤W∓Ltan θ0 2, (24) associated with the following corresponding outward unit normal vectors n[a]=−n[b]=(−sin θ0 cos θ0),n[r]=−n[l]=(1 0).(25) 9 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 2 3 4 1 µ0[º] 2.9 ¾min,c /E -0.2 -0.15 -0.1 -0.05 Figure 8: Maps of the compressive (minimum negative) principal in-plane stress σmin,c (35) associated with the five compressible Neo-Hookean membranes displayed in Fig. 1, evaluated from FE simulations under planar configuration. ϵin Fig. 9 for the membranes with θ0={2,2.9,3,4}◦. The different evolutions confirm linearity in θ0for the considered small values range and the nonlinear response in ϵ, which is displayed only for large strain ϵ > 10% at the quadrant centers (ξ1=−ξ2=±1/4) but also for intermediate small deformations (ϵ>3%) at the center (ξ1=ξ2= 0). The compressive stress exhibits a non-monotonic trend at both locations, although the reduction in its magnitude occurs at ϵ≈5% for the center and at ϵ≈30% for the quadrant centers. Moreover, while at small strains the compressive stress at the quadrant centers is around half of that at the center, at large strain the compressive stress becomes very small at the center and very large at the quadrant centers. The interplay of these nonlinear planar behaviours with the small bending stiffness of the membrane leads to the different scenarios for the wrinkling pattern evolution. 3.3.2 Inclined edges realignment with x1-axis during the stretching process During the elongation process, the parallelogram membrane modifies its shape towards a closer rectangular geometry, with its deformed free edges ∂B[a] xand ∂B[b] xreducing their inclination and becoming more aligned with the x1-axis. To quantitatively assess this phenomenon, the current inclination angle θis introduced as (Fig. 4) θ(θ0, β, ϵ) = arctan 1 1+ϵtan θ0+u2(ξ1= 1/2, ξ2= 0) −u2(ξ1=−1/2, ξ2= 0) L∈[0, θ0],(37) which measures the transition between the two limit cases, the undeformed parallelogram (θ=θ0) and the limit rectangular shape (θ= 0). The current inclination θhas been numerically evaluated from several simulations 16 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 Figure 9: Compressive (minimum negative) principal in-plane stress σmin,c at the center (ξ1=ξ2= 0) and at the quadrant centers (ξ1=−ξ2=±1/4) evaluated from FE simulations (under planar configuration) as a function of the elongation strain ϵfor Neo-Hookean membranes of Fig. 8 at different initial angle values θ0={2,2.9,3,4}◦. performed for six different moduli of the initial angle |θ0|={0.1,1,3,5,10,20}◦and for three different aspect ratios β={1,3,10}and is reported in Fig. 10(a) with varying of the elongation strain ϵ. For completeness, the curves of the current inclination θwith varying the strain ϵfor |θ0|= 0.1◦and β={1,2,3,4,5,10}are also reported in Fig. 10(b). From Fig. 10(a), it can be concluded that, although the evolution law for the Figure 10: Current inclination θ(37), normalized through division by its initial value θ0, evaluated from FE simulations (under planar configuration) as a function of the elongation strain ϵfor (a) different moduli of the initial inclination |θ0|={0.1,1,3,5,10,20}◦and aspect ratios β={1,3,10}and for (b) initial inclination |θ0|= 0.1◦ and different aspect ratios β={1,2,3,4,5,10}. current inclination θis highly nonlinear, the following linear expression holds for small initial inclination angles θ0(especially within the considered aspect ratio range β∈[2,5]) lim θ0→0θ(θ0, β, ϵ)=θ0R(β, ϵ),(38) where R(β, ϵ)∈[0,1] is a dimensionless function describing the reduction of the inclination angle magnitude due to both stretching and rigid-body motion. From Fig. 10(a), it is also noted that the function R(β, ϵ) is smooth 17 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 and therefore can be expressed through a polynomial expression in ϵtruncated at the M-th power (M∈N) R(β, ϵ) = 1 + M X j=1 (−1)jrj(β)ϵj.(39) From the theoretical point of view, the values for rj(β) can be obtained from the coefficients of the series expansion for R(β, ϵ) through its j-th derivative as rj(β) = (−1)j j! ∂jR(β, ϵ) ∂ϵj, j ∈N.(40) Nevertheless, the rj(β) coefficients obtained through this approach exhibit diverging magnitudes, resulting in poor convergence even for small ranges of elongation strain ϵ. For this reason, alternatively, the coefficients rj(β) can be defined from curve fitting of the numerical results reported in Fig. 10 to maintain, although in an approximate way, a good representation still with small Mvalues. As example, the coefficients rj(β) (j= 1,2,3) evaluated through curve fitting of the polynomial expression (39) for R(β, ϵ) limited to M=3 with the numerical data within the elongation range ϵ∈[0,10]% are reported in Tab. 2 for different aspect ratio values β={2,3,4,5}. Table 2: Values of the first three coefficients rj(β) (j= 1,2,3) for different aspect ratio values β={2,3,4,5} evaluated from cubic curve fitting of the numerical data in Fig. 10(b) within the elongation range ϵ∈[0,10]%. β2 3 4 5 r15.938 10.49 15.86 21.34 r227.96 82.24 166.54 266.42 r374.84 301.35 713.66 1245.4 It is finally highlighted that the magnitude of the current inclination θstrongly decreases even for limited small strain values, indeed for example at an elongation strain ϵ= 1% (10%) the current inclination θdiffers from the initial value θ0by 5.3% (38.4%) when β= 2 and by 19.1% (69%) when β= 5. 3.4 A fifth possible response: edge wrinkling parallel to the x2-axis displayed closely to the two acute-angled corners It is worth to mention that in addition to the four above-described behaviours (O,A1,B3,B1), a fifth one exists. This additional response falls outside the β-θ0domain reported in Fig. 7 and is associated to membrane geometries characterized by βθ0≳1/3. In this response the bifurcation from the flat state occurs through wrinkling parallel to the x2-axis displayed closely to the two acute-angled corners. The response transition from region B1to C is displayed in Fig. 11, where the maps of the out-of-plane displacement w(x) and of the compressive principal in-plane stress σmin,c(x) are reported for a membrane with an aspect ratio β= 5 at increasing initial angle θ0={3,4,5,7}◦and corresponding first critical elongation strain (ϵcr ={0.418,0.156,0.0024,0.018}%). While the out-of-displacement maps show that the central wrinkling is no longer the first critical mode when θ0= 5◦and 7◦, the compressive principal stress maps highlight how the central region becomes relatively unloaded and how the state changes close to the two acute-angled corners with the increase of the inclination angle θ0. An analytical explanation to this response is provided at the end of Sect. 4.1. 4 Analytical evaluation of the critical elongation strain for central wrinkling from perturbation approach and potential energy With reference to a small parameter a(with |a|≪1), the total potential energy Vassociated with a small outof-plane displacement abw(x1, x2) perturbing a flat (w= 0) equilibrated configuration described by the in-plane 18 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 ¾ ² min,c /E ¾ ² min,c /E Figure 11: Maps of the out-of-plane displacement w(x) (first line) and of the compressive principal in-plane stress σmin,c(x) (second line) from FE simulations for a membrane with an aspect ratio β= 5 and increasing initial angle θ0={3,4,5,7}◦(from left to right) at the corresponding first critical elongation strain ϵcr. The first bifurcation occurs as central wrinkling parallel to the x1-axis for θ0= 3◦and 4◦, and as edge wrinkling parallel to the x2-axis close to the two acute-angled corners for θ0= 5◦and 7◦. displacement e u(x1, x2) can be evaluated through the following second-order expansion in a V(e u, a bw)≈ V(e u,0)+aV′(e u, a bw)a=0 +a2 2V′′(e u, a bw)a=0 ,(41) where the prime symbol (′) stands for derivative in a. Because of equilibrium the first derivative vanishes, V′(e u, a bw)|a=0 = 0, while the second derivative V′′(e u, a bw)|a=0 follows from the generic expression of the total potential energy V(75) as (repeated indices imply summation) V′′(e u, a bw)a=0 =ZBxNγζ(e u)∂bw ∂xγ ∂bw ∂xζ +∂Mγζ ∂κρσ κ=0 ∂2bw ∂xγ∂xζ ∂2bw ∂xρ∂xσdx1dx2.(42) It follows that the stability of the flat configuration for the membrane is strictly connected to the sign of V′′(e u, a bw)|a=0, namely flat configuration (w= 0) is (stable if V′′(e u, a bw)|a=0 >0∀bw, unstable otherwise,(43) and therefore the critical elongation strain ϵcr is associated with the annihilation condition for the second derivative, critical condition: V′′(e u, a bw)a=0 = 0.(44) With reference to the linear response in bending, holding for the SVK (19)2and LE (21)2models, and reparameterizing the integral onto the unit square domain, the second derivative V′′ (42) can be rewritten as V′′(e u, a bw)|a=0 =βZBξ(N11(e u)1 β ∂bw ∂ξ1−tan θ0 ∂bw ∂ξ22 +N22(e u)∂bw ∂ξ22 + 2N12(e u)∂bw ∂ξ21 β ∂bw ∂ξ1−tan θ0 ∂bw ∂ξ2 +D β2W2"1 β ∂2bw ∂ξ2 1−2 tan θ0 ∂2bw ∂ξ1∂ξ2 +β1 + tan2θ0∂2bw ∂ξ2 22 −2(1 −ν) ∂2bw ∂ξ2 1 ∂2bw ∂ξ2 2−∂2bw ∂ξ1∂ξ22!#)dξ1dξ2. (45) 19 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 As the in-plane displacement e uand the membrane force Nover the hyperelastic parallelogram and the outof-plane displacement bwat wrinkling depend on the inclination angle θ0, the following second-order expansions in θ0can be considered e u(ξ) = e u(0)(ξ)+θ0e u(1)(ξ)+θ2 0e u(2)(ξ),N(ξ)=N(0)+θ0N(1)(ξ)+θ2 0N(2)(ξ),bw(ξ) = bw(0)(ξ)+θ0bw(1)(ξ)+θ2 0bw(2)(ξ), (46) from which the second-order expansion of the second derivative V′′(e u, a bw)|a=0 follows as V′′ =V′′(0) +θ0V′′(1) +θ2 0V′′(2),(47) where the explicit expressions for V′′(j)(j= 0,1,2) as functions of N(k)(ξ) and bw(k)(ξ) (k= 0,1,2) are deferred to Appendix C. To further proceed with the identification of the critical elongation strain ϵcr, an approximation for N(ξ) has be disclosed and a class for bw(ξ) describing possible wrinkling patterns has to be introduced. As only the magnitude of the inclination angle θ0, and not its sign, can affect the critical elongation strain ϵcr value, the first-order expansion of the second variation has to vanish V′′(1) = 0,(48) and the critical condition (44) reduces to V′′(0) +θ2 0V′′(2) = 0.(49) The evaluation of the membrane force N(ξ) and the critical condition for the potential energy second derivative V′′ are respectively addressed in the next two Subsections. 4.1 Approximate in-plane stress field for quasi-rectangular geometries through perturbation approach Under the assumption of small-strain (F≈I, det [F]≈1, E[m]≈ϵ[m]) and planar state (w= 0, M=0) the nonlinear relation (36) between membrane force Nand the (planar) Cauchy stress σreduces to σ=N/t and therefore from the linearized constitutive response (21) the Cauchy stress σis dependent on the membrane strain ϵ[m]by the linear elastic relation σ=E 1−ν2hνtrϵ[m]Is+(1−ν)ϵ[m]i.(50) The Cauchy stress σsatisfies the linearized version of the equilibrium equations (11) and of the static boundary conditions (26) and (27) ∇s·σ=0,for x1, x2∈Bx,   σn=0,for x1, x2∈∂B[a] x∪B[b] x, σ12 = 0,for x1, x2∈∂B[l] x∪B[r] x, (51) while the in-plane displacement e uremains subject to the kinematic boundary conditions (28) and (29). By assuming small values for θ0, the Cauchy stress σand the membrane strain ϵ[m]over the hyperelastic parallelogram, as well as the outward unit normals n(25), can be expanded at the second-order in the initial inclination angle θ0as σ=σ(0) +θ0σ(1) +θ2 0σ(2),ϵ[m]=ϵ[m](0) +θ0ϵ[m](1) +θ2 0ϵ[m](2),n=n(0) +θ0n(1) +θ2 0n(2),(52) where the strain-displacement (20), the stress-strain (50) and the stress-membrane force relations hold at the different orders. The described expansion considers the parallelogram domain Bxas a perturbation of the corresponding rectangular domain by inclining its two sides parallel to the x1axis by the small (initial) inclination angle θ0. Therefore, as the domain Bxdepends on θ0, it is instrumental to perform the perturbation analysis by referring the mechanical fields to the unit square domain Bξ, where the two domains are related through each 20 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 other via the linear transformation (30), which can be approximated at the second order in θ0under the condition β≪1/|θ0|as ξ=1 W   1 β0 0 1  +θ0"0 0 −1 0 # x.(53) Performing the analysis by referring to the auxiliary (normalized) coordinate ξdoes not only allow to define the differential problem on a non-varying domain with respect to the perturbation parameter θ0, but also to have a simple boundary description (namely, each boundary portion is defined only one variable coordinate, while the other one is constant), in a similar vein of the recently introduced approach for treating tapered beams [36] and [37]. The indefinite equilibrium equations (51)1can be expressed at the different orders in θ0in the unit square domain Bξas σ(0) 1j,1+βσ(0) 2j,2= 0, σ(1) 1j,1+βσ(1) 2j,2=β >0 σ(0) 1j,2, σ(2) 1j,1+βσ(2) 2j,2=βσ(1) 1j,2, j = 1,2,(54) and the (kinematic and static) boundary conditions (28) and (51)2,3as                  σ(0) j2ξ1, ξ2=±1 2= 0, σ(0) 12 ξ1=±1 2, ξ2= 0, eu(0) 1ξ1=−1 2, ξ2= 0, eu(0) 1ξ1=1 2, ξ2=ϵL,                      σ(1) j2ξ1, ξ2=±1 2=: {Eϵ,0} σ(0) 1jξ1, ξ2=±1 2, σ(2) j2ξ1, ξ2=±1 2=σ(1) j1ξ1, ξ2=±1 2+ : 0 1 2σ(0) j2ξ1, ξ2=±1 2, σ(k) 12 ξ1=±1 2, ξ2= 0, eu(k) 1ξ1=±1 2, ξ2= 0. j, k = 1,2.(55) While solving the differential problem (54) complemented by the boundary conditions (55) is straightforward at the 0th order, as the solution is given by the uniform tension along ξ1(corresponding to a rectangular domain under uniform stretching with lateral contraction allowed), achieving the solution at the 1st and 2nd order becomes awkward. Indeed, the respective boundary conditions enforce a discontinuity of the tangential stress σ(k) 12 (k≥1) at each corner of the unit square (ξ1=±ξ2=±1/2), similarly to the Timoshenko paradox problem [1], and therefore an analytical solution could be only obtained through series expression. Nevertheless, an approximate evaluation of the higher-order solution can still be pursued by replacing the local annihilation of the tangential stress and of normal displacement, Eqs. (55)7and (55)8, with corresponding global conditions along the interested edges Z1/2 −1/2 σ(k) 12 (ξ1=±1/2, ξ2) dξ2=Z1/2 −1/2eu(k) 1(ξ1=±1/2, ξ2) dξ2= 0, k = 1,2,(56) providing a null value for either the resultant shear force and the average normal displacement on the two parallelogram edges parallel to B[l] ξand B[r] ξ. Under this assumption, the second-order expansion (with rotational symmetry of order 2) for the stress field is found to be provided by the following closed-form over the parallelogram domain Bxthrough the expanded linear transformation (53), leading to (see Appendix D for details)      σ11(x) σ12(x) σ22(x)     =E ϵ            1 0 0     +θ0           −12x1x2 W2 −1 2+6x2 2 W2 0            +θ2 0              1 + β2−3ν 2+ 12x2 1−2x2 2 W2−C −24x1x2 W2 −2 + 12x2 2 W2                     ,x∈Bx, (57) where Cis a constant completing the description of the first-order in-plane displacement field e u(1) and remains arbitrary in the present analysis as it is based on the integral boundary conditions (56). 21 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 The minimum principal stress component σmin and its inclination ϕwith respect to the x1-axis can be evaluated through expression (35) from the stress field σ(x) (57), which, as β≪1/|θ0|, result σmin =−9 41−4x2 2 W22 θ2 0Eϵ < 0 (if ϵ>0), ϕ =π 2−θ01 2−6x2 2 W2−θ2 05−12x2 2 W26x1x2 W2.(58) Eq. (58) discloses a second-order compressive stress in the initial inclination angle θ0within a parallelogram under elongation (ϵ>0), with maximum magnitude for x2≈0 and acting approximately parallel to the x2-axis. With regards to the Cparameter, which would remain arbitrary from the present analytical method and defines only a second-order term in the stress component σ11, a comparison of the numerical value of σ11 at x1=x2= 0 from FE simulations with the corresponding value from the analytical expression (57)1for different β∈[1,10] and θ0∈ {0.2,10}◦(and disregarding the possible effect of the Poisson’s ratio ν) provides the following expression C=−3 41−2β2.(59) For this value of C, the displacement component u2that can be obtained from integration of the stress components (57) leads to the evaluation of the current inclination θthrough expression (37) as θ≈θ01−ϵ23 12 +β2,(60) which is in excellent agreement with the linear trend of the numerical curves reported in Fig. 10 associated to small values of the elongation strain ϵ. It is also noted that the obtained analytical expression (60) for the current inclination θimplies the normalized derivative 1/θ0(dθ/dϵ)|ϵ=0 =−23/12+β2, confirming that the current inclination angle θdramatically decrease for β≳2. As the closed-form expression (57) is obtained by disregarding the boundary layer effects associated with the edges ∂B[l] xand ∂B[r] x, the corresponding prediction is expected to be mostly reliable far from such boundaries. As a comparison, the stresses σ11,σ12,σ22, and σmin (normalized through division by Eϵ) are reported with varying ξ2∈[−1/2,1/2] for constant values of ξ1=−{3/8,1/4,1/8,0}as provided by Eqs. (57) and (58) and by the numerical simulations in Fig. 12. The results are shown for a parallelogram with aspect ratio β= 3 and an initial inclination angle θ0= 2.5◦at two different elongation strains, ϵ= 1% (left) and 5% (right). The analytical representation of the stress components is fully confirmed by the FE results at 1% except for σ22 and σmin evaluated at ξ1=−3/8, as this is closest to the boundary ∂B[l] x(ξ1=−1/2). Differently, although still within the small strain regime, at 5% the analytical predictions loses fidelity for σ12,σ22, and σmin even for central coordinates (ξ1=−1/4,−1/8 and 0). This loss of reliability for the analytical solution is associated to the dramatic change in the current inclination angle θ(discussed in Sect. 3.3.2), which introduces strong nonlinearities in the mechanical response even under limited values of the elongation strain ϵ. A further comment is also provided with reference to membranes with large aspect ratios β. Despite the obtained solution (57) is formally valid only for β≪1/|θ0|as it is based on the transformation expansion (53), it is interesting to observe that the stress component σ11 truncated at the first-order in θ0evaluated at the two acute-angled corners is σ11(ξ1= sign[θ0]ξ2=±1/2) = E ϵ (1 −3β|θ0|),(61) providing an hint about the possibility for a compressive stress parallel to the x1axis when β > 1/(3|θ0|). In a rough sense, this is in agreement with the stress maps in Fig. 11 for the membrane with β= 5 and θ0= 5◦and 7◦ (for which βθ0= 0.436 and 0.611, respectively), associated with edge wrinkling parallel to the x2-axis displayed closely to the two acute-angled corners (Sect. 3.4). 4.2 Analytical evaluation of the critical elongation strain ϵcr As the in-plane membrane force has been determined in the previous Subsection, the critical condition (49) can be fully elucidated by selecting a specific out-of-plane displacement bw. In general terms, the out-of-plane displacement 22 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 Figure 12: Stress components σ11,σ12,σ22,σmin as functions of ξ2, normalized through division by Eϵ and evaluated at ξ1=−{3/8,1/4,1/8,0}. Continuous and dotted curves represent results respectively evaluated from the analytical expression (57) and the FE simulations for β= 3 and θ0= 2.5◦. Two elongation strains are considered, ϵ= 1% (left columnn) and 5% (right column). bwcan be introduced as the following linear superposition bw(x)=H·h(x),(62) 23 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 where Hand h(x) are the vectors respectively collecting the amplitude Hjand the functions hj(x) (j=1, ..., N, being Nthe number of functions considered as basis for bw), the latter such that the kinematic boundary condition (27)2is satisfied and the first-order term V′′(1) of the potential energy vanishes, (48). Under representation (68), the second order expansion in θ0of the second derivative V′′, Eq. (47), is described by the following quadratic expression in H V′′ =E W α3β3H·nK(0) e+α2β2K(0) gϵ+θ2 0hK(2) e+α2β2K(2) gϵioH,(63) where K(j) eand K(j) gare respectively the (normalized) elastic and geometric stiffness (symmetric N×N) matrices at the j-th order (j= 0 and 2). The critical condition (49) for the second derivative of the potential energy (63) defines a generalized eigenproblem where His the eigenvector and the critical elongation strain ϵcr is the eigenvalue obtained by imposing the determinant annihilation det nK(0) e+α2β2K(0) gϵcr +θ2 0hK(2) e+α2β2K(2) gϵcrio= 0.(64) Therefore, the critical condition (64) defines ϵcr as the root(s) of a polynomium of order N cNϵN cr +cN−1ϵN−1 cr +... +c2ϵ2 cr +c1ϵcr +c0= 0,(65) with cj(j= 0, ..., N) being coefficients depending on the elastic (K(0) eand K(2) e) and geometric (K(0) gand K(2) g) stiffness matrices. In order to maintain simplicity in the treatment, the number Nof functions hj(x) is considered as N= 1, in which case the matrices reduce to a single component, K(j) e=K(j) eand K(j) g=K(j) g(j= 0,2), and the polynomial equation (65) reduces to a linear equation α2β2hK(0) g+θ2 0K(2) giϵcr +K(0) e+θ2 0K(2) e= 0,(66) for which the critical elongation strain ϵcr follows as ϵcr =−K(0) e+θ2 0K(2) e α2β2hK(0) g+θ2 0K(2) gi.(67) Considering that the closed-form expression (57) for the in-plane stress state is accurate only in the central region and for small elongation strains ϵ(Fig. 12), the critical condition is therefore evaluated exclusively for the central wrinkling patterns numerically identified in Sect. 3.2. These appear as wrinkles oriented parallel to the x1axis in the undeformed configuration, with an amplitude that decays away from the origin (x1=x2= 0) along both the x1and x2axes. Therefore, the following expression for the function h1(x) is considered h1(x) = "1−2x1 βW 2#2 cos4πx2 Wcos nπ 1 2+x2 W,(68) being na positive number (not restricted to be natural) defining the wrinkling wavelength along x2. Once the function h1is selected, the four elastic and geometric stiffness coefficients K(0) e(n), K(0) g(n), K(2) e(n), and K(2) g(n) can be evaluated as functions of n(in addition to the dimensionless membrane properties α,β,ν) and the critical elongation ϵcr,n(n) assessed from expression (67). While by definition K(0) e(n)>0, from the integration over the domain it is found that K(2) e(n)=0∀nand the expression (67) for the critical elongation strain ϵcr reduces to ϵcr,n =−K(0) e(n) α2β2hK(0) g(n)+θ2 0K(2) g(n)i,(69) 24 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 which, considering that from the domain integration it is found that K(0) g(n)>0∀n, implies that the critical elongation strain is never positive for rectangular membranes, ϵcr,n(θ0=0)<0∀n, and therefore confirming that no tension instability occurs for this limit geometry. However, the question remains if a non-vanishing inclination angle θ0may trigger or not the central wrinkling under elongation strain, ϵcr >0. With this regard, it can be noted that whenever K(2) g(n)<0 then a corresponding critical elongation strain exists, ϵcr,n >0, for a sufficiently large inclination angle θ0(although within the small range). It follows that for a specific membrane material (ν) and geometry (α,βand θ0), the critical elongation strain ϵcr can be evaluated through a minimization over the positive parameter nof the positive critical elongation strains ϵcr,n >0, namely ϵcr(θ0, α, β, ν) = min n∈R+[ϵcr,n(θ0, α, β, ν, n)] ,restricted to ϵcr,n(θ0, β, α, ν, n)>0.(70) Leaving the minimization over the whole set of positive values for nto specific membrane geometry and material, in the case of natural values of nthe following closed-form expressions for the critical elongation strain ϵcr,n can be obtained as ϵcr,n=1 =−50π2 (1 −ν2)α2β2 679π4β4+ 1320π2β2+ 3528 50400π2−θ2 0(175π2(1549β2+ 432ν+ 72) + 33000π4β2−643272), ϵcr,n=2 =−100π2 (1 −ν2)α2β2 386π4β4+ 732π2β2+ 3087 88200π2−θ2 0(7π2(52267β2+ 3150(6ν+ 1)) + 36600π4β2−1627206), ϵcr,n=3 =−2450π2 (1 −ν2)α2β2 7663π4β4+ 9384π2β2+ 15624 10936800π2−θ2 0(49π2(2064019β2+ 55800(6ν+ 1)) + 11495400π4β2−168194184), ϵcr,n=4 =−9800π2 (1 −ν2)α2β2 4288π4β4+ 3840π2β2+ 4473 12524400π2−θ2 0(20π2(8318078β2+ 156555(6ν+ 1)) + 18816000π4β2−220379757), ϵcr,n≥5=−2π2n2 5 (1 −ν2)α2β2 17640 + 120π2β27n2+ 16+π4β435n4+ 480n2+ 512 2016π2n2−θ2 0Gn(β, ν), (71) where Gn(β, ν) is defined for n∈N+and with n≥5 by Gn(β, ν) = −382205952 2π2β2−21+ 7π2205 + 24π2β2n20 −4n18 π220489 + 2424π2β2 −126(6ν+ 1)) + 8610) + 18n16 35π22899β2−288ν−48+ 12216π4β2+ 114800 −4n14 π24966619 + 618024π2β2−182196(6ν+ 1)+ 12450060 +n12 π2104529355 + 14060328π2β2−9082080(6ν+ 1)+ 620608800 −72n10 3π2805483 + 153096π2β2−291389(6ν+ 1)+ 59734745 −16n8π237922381 + 1212312π2β2+ 15191820(6ν+ 1)−1045945908 +128n6π2212034225 + 860424π2β2+ 3885903(6ν+ 1)−288398145 −36864n4π2436065 + 2082π2β2+ 12915(6ν+ 1)−1398495 +2654208n2π22865 + 24π2β2+ 63(6ν+ 1)−17220/n2−1n2−4n2−9n2−162. (72) The critical elongation strain ϵcr,n obtained from the closed-form expressions (71) is reported at varying natural values of nas a function of the absolute value of the inclination angle |θ0|for α= 1500 and ν= 0.33 in Fig. 13 for different β={2,3,4,5}. Each curve associated to a different nis reported as dashed except for the continuous portion representative of the minimum within the set of considered n. The natural value nproviding the minimization is highlighted below the horizontal axis label for the corresponding ranges of θ0. This information shows how the wrinkling wavelength (which is inversely proportional to n) decreases with the decrease of the aspect ratio βand of the initial inclination angle |θ0|. The union of the minimum curve portions associated with natural values of napproximately represents the envelope ϵcr that would be obtained by performing the minimization over the positive set values for nas described by Eq. (70). The critical conditions numerically evaluated from the FE model are also included as circles in 25 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 References [1] J.R. Barber. Elasticity. Springer, 2002. [2] D. Bigoni, F. Bosi, F. Dal Corso, and D. Misseroni. “Instability of a penetrating blade”. In: Journal of the Mechanics and Physics of Solids 64 (2014), pp. 411–425. [3] F. Bosi, D. Misseroni, F. Dal Corso, S. Neukirch, and D. Bigoni. “Asymptotic self-restabilization of a continuous elastic structure”. In: Phys. Rev. E 94.6 (2016), p. 063005. [4] M. Brojan, D. Terwagne, R. Lagrange, and P.M. Reis. “Wrinkling crystallography on spherical surfaces”. In: Proceedings of the National Academy of Sciences 112.1 (2015), pp. 14–19. [5] C.R. Calladine. Theory of shell structures. Cambridge University Press, 1983. [6] E. Cerda and L. Mahadevan. “Geometry and physics of wrinkling”. In: Physical Review Letters 90.7 (2003), p. 074302. [7] E. Cerda, K. Ravi-Chandar, and L. Mahadevan. “Wrinkling of an elastic sheet under tension”. In: Nature 419.6907 (2002), pp. 579–580. [8] P.-P. Chai, Y. Liu, and F.-F. Wang. “Stretch-induced wrinkling of anisotropic hyperelastic thin films”. In: Thin-Walled Structures 200 (2024), p. 111961. [9] C.D. Coman. “On the applicability of tension field theory to a wrinkling instability problem”. In: Acta Mechanica 190.1 (2007), pp. 57–72. [10] C.D. Coman. “Wrinkling of a normally loaded, spinning, elastic membrane: An asymptotic approximation”. In: International Journal of Non-Linear Mechanics 156 (2023), p. 104482. [11] C.D. Coman and A.P. Bassom. “On the nonlinear membrane approximation and edge-wrinkling”. In: International Journal of Solids and Structures 82 (2016), pp. 85–94. [12] A. Comitti, H. Vijayakumaran, M.H. Nejabatmeimandi, L. Seixas, A. Cabello, D. Misseroni, M. Penasa, Ch. Paech, M. Bessa, A.C. Bown, F. Dal Corso, and F. Bosi. “Ultralight Membrane Structures Toward a Sustainable Environment”. In: Sustainable Structures and Buildings. Ed. by A. Bahrami. Cham: Springer International Publishing, 2024, pp. 17–37. isbn: 978-3-031-46688-5. [13] N. Damil, M. Potier-Ferry, and H. Hu. “New nonlinear multi-scale models for wrinkled membranes”. In: Comptes Rendus M´ecanique 341.8 (2013), pp. 616–624. [14] D.A. Danielson and S. Natarajan. “Tension field theory and the stress in stretched skin”. In: Journal of Biomechanics 8.2 (1975), pp. 135–142. [15] M. Destrade, Y. Fu, and A. Nobili. “Edge wrinkling in elastically supported pre-stressed incompressible isotropic plates”. In: Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472.2193 (2016), p. 20160410. [16] M. Ding, F. Xu, T. Wang, and C. Fu. “Nanosleeves: Morphology transitions of infilled carbon nanotubes”. In: Journal of the Mechanics and Physics of Solids 152 (2021), p. 104398. [17] E. Feh´er, T.J. Healey, and A.A. Sipos. “The Mullins effect in the wrinkling behavior of highly stretched thin films”. In: Journal of the Mechanics and Physics of Solids 119 (2018), pp. 417–427. [18] N. Friedl, F.G. Rammerstorfer, and F.D. Fischer. “Buckling of stretched strips”. In: Computers & Structures 78.1-3 (2000), pp. 185–190. [19] G. Friesecke, R. James, and S. M¨uller. “A Hierarchy of Plate Models Derived from Nonlinear Elasticity by Gamma-Convergence.” In: Arch. Rational Mech. Anal. 180 (2006), pp. 183–236. [20] C. Fu, H.-H. Dai, and F. Xu. “Computing wrinkling and restabilization of stretched sheets based on a consistent finite-strain plate theory”. In: Computer Methods in Applied Mechanics and Engineering 384 (2021), p. 113986. 32 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 [21] C. Fu, T. Wang, F. Xu, Y. Huo, and M. Potier-Ferry. “A modeling and resolution framework for wrinkling in hyperelastic sheets at finite membrane strain”. In: Journal of the Mechanics and Physics of Solids 124 (2019), pp. 446–470. [22] C. Fu, Y. Yang, T. Wang, and F. Xu. “A consistent finite-strain plate model for wrinkling of stretched anisotropic hyperelastic films”. In: Thin-Walled Structures 179 (2022), p. 109643. [23] T.J. Healey, Q. Li, and R.B. Cheng. “Wrinkling behavior of highly stretched rectangular elastic films via parametric global bifurcation”. In: Journal of Nonlinear Science 23 (2013), pp. 777–805. [24] M.G. Hilgers and A.C. Pipkin. “Bending energy of highly elastic membranes”. In: Quarterly of Applied Mathematics 50(2) (1992), pp. 389–400. [25] M.G. Hilgers and A.C. Pipkin. “Bending energy of highly elastic membranes II”. In: Quarterly of Applied Mathematics 54(2) (1996), pp. 307–316. [26] J. Hure, B. Roman, and J. Bico. “Stamping and wrinkling of elastic plates”. In: Physical Review Letters 109.5 (2012), p. 054302. [27] T. Ishida, S. Matsubara, S. Nagashima, and D. Okumura. “Deformation in the wrinkle–crease transformation”. In: International Journal of Solids and Structures 298 (2024), p. 112876. [28] N. Jacques and M. Potier-Ferry. “On mode localisation in tensile plate buckling”. In: Comptes rendus. M´ecanique 333.11 (2005), pp. 804–809. [29] F.L. Jim´enez, N. Stoop, R. Lagrange, J. Dunkel, and P.M. Reis. “Curvature-controlled defect localization in elastic surface crystals”. In: Physical Review Letters 116.10 (2016), p. 104301. [30] T.-Y. Kim, E. Puntel, and E. Fried. “Numerical study of the wrinkling of a stretched thin sheet”. In: International Journal of Solids and Structures 49.5 (2012), pp. 771–782. [31] A. Kossa, M.T. Valentine, and R.M. McMeeking. “Analysis of the compressible, isotropic, neo-Hookean hyperelastic model”. In: Meccanica 58 (2023), pp. 217–232. [32] R. Lagrange, F. L. Jim´enez, D. Terwagne, M. Brojan, and P.M. Reis. “From wrinkling to global buckling of a ring on a curved substrate”. In: Journal of the Mechanics and Physics of Solids 89 (2016), pp. 77–95. [33] C.M. Landis, R. Huang, and J.W. Hutchinson. “Formation of surface wrinkles and creases in constrained dielectric elastomers subject to electromechanical loading”. In: Journal of the Mechanics and Physics of Solids 167 (2022), p. 105023. [34] Q. Li and T.J. Healey. “Stability boundaries for wrinkling in highly stretched elastic sheets”. In: Journal of the Mechanics and Physics of Solids 97 (2016), pp. 260–274. [35] F. Liu, F. Xu, and C. Fu. “Orientable wrinkles in stretched orthotropic films”. In: Extreme Mechanics Letters 33 (2019), p. 100579. [36] G. Migliaccio and F. D’Annibale. “On the inadequacy of a stepped-beam approach in predicting shear stresses in tapered slender solids”. In: European Journal of Mechanics-A/Solids 111 (2025), p. 105590. [37] G. Migliaccio and G. Ruta. “Rotor blades as curved, twisted and tapered beam-like structures subjected to large deflections”. In: Engineering Structures 222 (2020), p. 111089. [38] A. Mirandola, A. Cutolo, A.R. Carotenuto, N. Nguyen, L. Pocivavsek, M. Fraldi, and L. Deseri. “Toward new scaling laws for wrinkling in biologically relevant fiber-reinforced bilayers”. In: Journal of Applied Physics 134 (2023), p. 154702. [39] M. Misseroni, P.P. Pratapa, K. Liu, and G.H. Paulino. “Experimental realization of tunable Poisson’s ratio in deployable origami metamaterials”. In: Extreme Mechanics Letters 53 (2022), p. 101685. [40] P. Nardinocchi and E. Puntel. “Swelling-induced wrinkling in layered gel beams”. In: Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473.2207 (2017), p. 20170454. 33 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 [41] V. Nayyar, K. Ravi-Chandar, and R. Huang. “Stretch-induced stress patterns and wrinkles in hyperelastic thin sheets”. In: International Journal of Solids and Structures 48.25-26 (2011), pp. 3471–3483. [42] A.C. Pipkin. “The relaxed energy density for isotropic elastic membranes”. In: IMA journal of applied mathematics 36.1 (1986), pp. 85–99. [43] P. Plucinsky and K. Bhattacharya. “Microstructure-enabled control of wrinkling in nematic elastomer sheets”. In: Journal of the Mechanics and Physics of Solids 102 (2017), pp. 125–150. [44] E. Puntel, L. Deseri, and E. Fried. “Wrinkling of a stretched thin sheet”. In: Journal of Elasticity 105 (2011), pp. 137–170. [45] E. Reissner. “On tension field theory”. In: Proc. of the 5th Int. Congr. for Applied Mechanics Harvard Univ. & MIT (1938), pp. 88–92. [46] J. Shen, Y. Fu, A. Pirrera, and R.M.J. Groh. “Wrinkling of differentially growing bilayers with similar film and substrate moduli”. In: Journal of the Mechanics and Physics of Solids 193 (2024), p. 105900. [47] A.A. Sipos and E. Feh´er. “Disappearance of stretch-induced wrinkles of thin sheets: a study of orthotropic films”. In: International Journal of Solids and Structures 97 (2016), pp. 275–283. [48] D.J. Steigmann. “Koiter’s shell theory from the perspective of three-dimensional nonlinear elasticity”. In: Journal of Elasticity 111 (2013), pp. 91–107. [49] D.J. Steigmann. “Tension-field theory”. In: Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 429.1876 (1990), pp. 141–173. [50] M. Stein and J.M. Hedgepeth. Analysis of partly wrinkled membranes. National Aeronautics and Space Administration, 1961. [51] N. Stoop, R. Lagrange, D. Terwagne, P.M. Reis, and J. Dunkel. “Curvature-induced symmetry breaking determines elastic surface patterns”. In: Nature materials 14.3 (2015), pp. 337–342. [52] M. Su˜n´e, C. Arratia, A.F. Bonfils, D. Vella, and J.S. Wettlaufer. “Wrinkling composite sheets”. In: Soft Matter 19.45 (2023), pp. 8729–8743. [53] M. Taylor, K. Bertoldi, and D.J. Steigmann. “Spatial resolution of wrinkle patterns in thin elastic sheets at finite strain”. In: Journal of the Mechanics and Physics of Solids 62 (2014), pp. 163–180. [54] S. Timoshenko and S. Woinowsky-Krieger. Theory of plates and shells. Mc Graw-Hill, 1959. [55] H. Vandeparre, M. Pi˜neirua, F. Brau, B. Roman, J. Bico, C. Gay, W. Bao, CH.N. Lau, P.M. Reis, and P. Damman. “Wrinkling hierarchy in constrained thin sheets from suspended graphene to curtains”. In: Physical Review Letters 106.22 (2011), p. 224301. [56] S.P. Venkata, V. Balbi, M. Destrade, D. Accoto, and G. Zurlo. “Programmable wrinkling for functionallygraded auxetic circular membranes”. In: Extreme Mechanics Letters 63 (2023), p. 102045. [57] H. Wagner. “Flat sheet metal girders with very thin metal web”. In: Z. flugtechn. motorluftschiffahrt 20 (1929), pp. 200–314. [58] F.-F. Wang, T. Wang, X. Zhang, Y. Huang, I. Giorgio, and F. Xu. “Wrinkling of twisted thin films”. In: International Journal of Solids and Structures 262 (2023), p. 112075. [59] T. Wang, C. Fu, M. Potier-Ferry, and F. Xu. “Morphomechanics of growing curled petals and leaves”. In: Journal of the Mechanics and Physics of Solids 184 (2024), p. 105534. [60] T. Wang, C. Fu, F. Xu, Y. Huo, and M. Potier-Ferry. “On the wrinkling and restabilization of highly stretched sheets”. In: International Journal of Engineering Science 136 (2019), pp. 1–16. [61] T. Wang, M. Potier-Ferry, and F. Xu. “A nonlinear toroidal shell model for surface morphologies and morphogenesis”. In: Journal of the Mechanics and Physics of Solids 200 (2025), p. 106135. [62] T. Wang, Y. Yang, C. Fu, F. Liu, K. Wang, and F. Xu. “Wrinkling and smoothing of a soft shell”. In: Journal of the Mechanics and Physics of Solids 134 (2020), p. 103738. 34 Accepted for publication in Journal of the Mechanics and Physics of Solids (2026) doi: https://doi.org/10.1016/j.jmps.2025.106461 [63] T. Wang, Y. Yang, and F. Xu. “Mechanics of tension-induced film wrinkling and restabilization: a review”. In: Proceedings of the Royal Society A 478.2263 (2022), p. 20220149. [64] W. Wong and S. Pellegrino. “Wrinkled membranes I: experiments”. In: Journal of Mechanics of Materials and Structures 1.1 (2006), pp. 3–25. [65] W. Wong and S. Pellegrino. “Wrinkled membranes III: numerical simulations”. In: Journal of Mechanics of Materials and Structures 1.1 (2006), pp. 63–95. [66] F. Xu and M. Potier-Ferry. “On axisymmetric/diamond-like mode transitions in axially compressed core– shell cylinders”. In: Journal of the Mechanics and Physics of Solids 94 (2016), pp. 68–87. [67] F. Xu, M. Potier-Ferry, S. Belouettar, and H. Hu. “Multiple bifurcations in wrinkling analysis of thin films on compliant substrates”. In: International Journal of Non-Linear Mechanics 76 (2015), pp. 203–222. [68] F. Xu, T. Wang, C. Fu, Y. Cong, Y. Huo, and M. Potier-Ferry. “Post-buckling evolution of wavy patterns in trapezoidal film/substrate bilayers”. In: International Journal of Non-Linear Mechanics 96 (2017), pp. 46– 55. [69] E. Yang, M. Zhang, J. Zeng, and F. Tian. “Wrinkling and restabilization of a hyperelastic PDMS membrane at finite strain”. In: Soft Matter 18.29 (2022), pp. 5465–5473. [70] L. Zheng. Wrinkling of dielectric elastomer membranes. California Institute of Technology, 2009. 35