A New Modal Theory for Wrinkling Analysis of Stretched Membranes

Journal Pre-proof

A New Modal Theory for Wrinkling Analysis of Stretched Membranes A.D. Martins , N. Silvestre , R. Bebiano PII: DOI: Reference:

S0020-7403(19)34439-X https://doi.org/10.1016/j.ijmecsci.2020.105519 MS 105519

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

22 November 2019 28 January 2020 8 February 2020

Please cite this article as: A.D. Martins , N. Silvestre , R. Bebiano , A New Modal Theory for Wrinkling Analysis of Stretched Membranes, International Journal of Mechanical Sciences (2020), doi: https://doi.org/10.1016/j.ijmecsci.2020.105519

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Highlights

 A new modal theory to analyze the wrinkling of stretched membranes is proposed.  Modal nature provides novel insight into the localized deformation of wrinkling phenomenon.  Rigorous pre-wrinkling analysis with axial, shear and transverse extension modes is required.  Duplicity and quadruplicity of symmetric/antisymmetric wrinkling modes is found.  Cosine longitudinal approximation for displacements is valid for a limited range of aspect ratio.

A New Modal Theory for Wrinkling Analysis of Stretched Membranes A.D. Martins a, N. Silvestrea*, R. Bebianob a

b

IDMEC, DEM, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal

CERIS, DECivil, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal * corresponding author: [email protected]

Abstract This paper presents a new modal theory for the wrinkling analysis of membranes and discusses the stretch-induced wrinkling behaviour of fixed-ended thin sheets. This modal theory corresponds to an extension of Generalised Beam Theory (GBT) to analyse the wrinkling behaviour of highly stretched sheets. The inherent GBT modal nature adopted in this investigation enables the acquisition of an indepth knowledge (when compared with the traditional shell finite element analysis) of the behaviour underlying this peculiar buckling phenomena. In order to conduct this study, the work begins by presenting the formulation, describing its main steps and concepts involved in the determination of the “deformation modes” (cross-section analysis), followed by the description of the procedures involved in the development and implementation of the aforementioned formulation that accounts for pre-wrinkling and wrinkling analysis (member analysis). Then, an illustrative example is presented and discussed in detail, which involves the characterisation of the sheet “signature curve” (tensile critical stress vs. length), respective critical wrinkling modes, modal participation diagrams and pre-wrinkling stress distribution patterns. For comparison and validation purposes, refined ABAQUS shell finite element wrinkling results are also reported. Finally, several conclusions are drawn regarding the wrinkling behaviour of stretched sheets, mostly based on the modal nature of this formulation. Keywords: Thin sheets; Wrinkling; Tensile buckling; Pre-wrinkling analysis; Generalised Beam Theory (GBT).

1. Introduction Ultra-thin structures are widely used in aerospace applications, such as solar sails, inflatable antennas and radars, due to their high strength-to-weight ratios and low space requirement [1]. However, their reduced thickness makes them invariably prone to instability phenomena that affects their structural performance and must be accounted in design stages. Although instability phenomena are typically associated with compressive applied forces/stresses, members subjected to tensile forces/stresses may also be affected. Indeed, this peculiar instability phenomenon occurs in fixed-ended ultra-thin structures 2

(membranes or thin sheets) under uniaxial tension, i.e., corresponds to wrinkling phenomenon of stretched thin sheets. However, it should be noted that wrinkling is not restricted to uniaxial loading but to a variety of loading conditions, e.g., shear and two pairs of equal and opposite diagonal forces at the corners [2] or tension and shear [3] have also been analysed. Generally speaking, an ultra-thin member increasingly stretched in one direction may wrinkle in the orthogonal direction if compressive normal transverse stresses develop – indeed, this constitutes a pre-requisite for wrinkling. The wrinkling phenomenon of stretched sheets has attracted the research community in the last decades, comprising (mostly) numerical simulations, analytical formulations and experimental investigations. Friedl et al. [4] (i) provided a mechanical explanation for the development of compressive transversal stresses in fixed-ended stretched sheets, (ii) proposed approximate analytical formulas, based on uniform pre-wrinkling stresses and a single parameter correlating the longitudinal tensile and transversal compressive stresses, to estimate the critical tensile loading causing instability including the respective minimum half-wave numbers in the transversal direction. Jacques and Potier-Ferry [5] also developed analytical expressions to determine the critical tensile load in fixed-ended sheets assuming (i) uniform longitudinal force and (ii) variation of transversal compressive force in the longitudinal direction according to a single function. In order to apply the developed formulation, the shape of the transverse stresses and transversal-to-longitudinal stress ratio was estimated by 1st order elastic shell finite element analysis. The approximate results obtained were in good agreement with finite element results for both critical tensile forces (differences below 10%) and mode shapes. Puntel et al. [6] developed an energy formulation to calculate the critical strain of clamped rectangular sheets under high tension. They adopted a cosine function to approximate the longitudinal variation of displacements of the wrinkling modes. Healey et al. [7] studied the wrinkling of highly stretched thin rectangular sheets and pointed out the potential dangers in modeling thin-sheet problems involving large mid-plane strains via the classical Föppl–von Kármán theory. They also proposed a simplified approach to calculate the wrinkling strain by adopting the cosine function to approximate the variation of displacements along the stretched direction. Nayyar et al. [1] reported a finite element investigation on stretch-induced wrinkling of hyperelastic thin sheets. The investigation begins with the determination of transverse normal stress distribution patterns of rectangular sheets with distinct aspect ratios and under tensile longitudinal strains – this was summarised in a phase diagram comprising four different distribution patterns. They also conducted eigenvalue and post-wrinkling analyses. The latter showed that the wrinkle amplitude first increases when the longitudinal strain varies between 5% to 10% and then the wrinkle amplitude drops with an increase of the longitudinal strain, after which all wrinkles become flattened at longitudinal strains approximately of 30% (restabilisation)1. A few years later Nayyar

1

Note that Zheng [8] was the first investigator that explored this restabilization phenomenon in highly stretched sheets.

3

et al. [9] conducted an experimental investigation on 19 stretch-induced wrinkling of polyethylene sheets using the 3D Digital Image Correlation (DIC) technique comprising two aspect ratios (2.0 and 2.5) and two different strain rates (0.0169/s and 0.00169/s). Contrarily to the numerical simulations conducted on hyperelastic thin sheets, the wrinkle amplitude of the specimens does not vanish at large strains (30%). Then they implemented in ABAQUS two models with distinct viscoelastic constitutive models and showed a significant improvement in predicting the experiments, concluding that the nonlinear viscoelasticity is essential to capture to structural response of these thin sheets (with these models the wrinkle amplitude of the specimens does not vanish at large strains as occurred with an hyperelastic constitutive model). Carvalho [10] also conducted and experimental and numerical investigation on thin sheets subjected to uniaxial tension. The experimental study involved 12 Kapton® thin sheets with distinct aspect ratios ranging between 1.0 and 4.0 also employing a 3D DIC technique to quantify wrinkles along the loading procedure. A shell finite element model was used to predict the wrinkling (linear stability) and post-wrinkling behaviour of the sheets and a qualitatively good agreement was found when comparing with the experiments. Silvestre [11] reported an analytical investigation on the approximation of stresses, strains and displacements of pre-wrinkling (1st order elastic analysis) state of stretched sheets. By adding sequentially three displacement fields (or “modes”) in the governing solution (axial extension, transverse and warping shear modes), it was concluded that it is not only the Poisson’s effect the sole cause of compressive transversal stresses but also the warping shear of the membrane (shear lag). Then, the analytical formulae developed were compared with results obtained by ABAQUS shell finite element analysis (SFEA) in a sheet characterised by an aspect ratio of 2.0. In general, a good correlation with SFEA was found and the (small) discrepancies can only be eliminated by increasing the number of (deformation) modes2. In 2016, Li and Healey [12] proposed a systematic global bifurcation/continuation approach to deal with wrinkling of stretched rectangular sheets using the more accurate elastomer models (neoHookean and Mooney–Rivlin materials) and provided an efficient strategy for the computation of stability boundaries (aspect ratio vs. macroscopic strain) for different models. Sipos and Feher [13] presented a Föppl–von Kármán formulation extended to finite in-plane strains for thin sheets under uniform tension and validated it by means of experimental results, namely the disappearance of the wrinkled pattern of clamped. Their work provided experimental proof about the second bifurcation point along the trivial equilibrium path, at which the planar solution regains stability, and confirmed that the aspect ratio of the rectangular sheet significantly affects the location of the bifurcation points. Since 2018, several works have been conducted on wrinkling phenomenon and towards its mechanical understanding. Iwasa [14] presented a simplified equation for estimating wavelength and maximum 2

However, the ensuing analytical solutions would be much more complex to obtain. An additional comment is added in the end of Section 5 regarding the additional number of deformation modes.

4

amplitude of wrinkles using a tension-field solution for various types of membrane structures. The effectiveness of the formula was discussed for the two classical membrane models and later experimentally confirmed [15] on the basis of two membrane models (in-plane shear and a corner tension load) and photogrammetry and laser displacement sensor. Zhu et al [16] presented a study with an interesting experimental observation of stretching a highly orthotropic thin film: no wrinkles were observed in the sample when stretched in the high-stiffness direction, but wrinkles occurred when stretched in the low-stiffness direction. This counterintuitive phenomenon triggered the rethinking of the existing theories. A comprehensive analysis of this problem was carried out using energy approach by decomposing the total energy into several components with different physical meanings. In 2019, Wang et al. [17] investigated the wrinkling behaviour of highly stretched thin sheets and provided a 3D phase diagram on the stability boundary. They showed that Poisson’s ratio makes later onset of wrinkling, lower amplitude and earlier disappearance of wrinkles and when Poisson’s ratio is below a threshold, no wrinkles occur. Fu et al. [18] developed a model for the nonlinear phenomenon of highly stretch-induced wrinkling of hyperelastic membranes and showed that reducing the Poisson’s ratio leads to a decrease of transverse compressive stresses and lower probability of wrinkling. Huang et al. [19] presented a multiscale model to study membrane instabilities, by combining the Fourier analysis and the bridging domain method. Their model was firstly established based on the Föppl-Von Karman plate equations and the numerical results showed that this approach was able to simulate the membrane instabilities with high efficiency and accuracy. Dadgar-Rad and Imani [20] presented and validated a formulation of gradient-elastic membranes at finite strains, by introducing the internal energy density with first and second deformation gradient tensors and constitutive laws for isotropic and anisotropic materials. Liu et al. [21] applied an extended Föppl–von Kármán nonlinear plate model for orthotropic materials and used an asymptotic-numerical method to investigate wrinkling in uniaxially stretched orthotropic films. They revealed that the degree of orthotropy and shear modulus significantly influence the appearance and disappearance of wrinkles. These very recent investigations show the relevance of studying wrinkling phenomena and understanding it. Typically, three ways have been explored: (i) the use of shell finite element simulations, which provide accurate outputs but are not easy to apprehend/interpret due to their lack of mechanical reasoning, (ii) the solution of nonlinear plate equations, which give good results but fail to provide a general framework of the problem, and (iii) the development of approximate analytical formulas, which offer straightforward calculation of critical stresses and strains but only applicable for rather specific cases. The aim of the present paper is to show an alternative and very promising approach to investigate wrinkling phenomena, which grounds on a modal theory that is mechanically sound because not only it separates the contribution of different mechanical contributions (extension, shear, warping, local) but also is more computationally efficient (similar 5

results to shell finite analysis are obtained with much less degrees-of-freedom). The problems unveiled by recent works may be further explained with this modal approach, such as the influence of Poisson’s ratio on the wrinkling mode shape and its restabilisation. This modal approach is based on a “higher-order” beam theory that has emerged in the last 20 years – Generalised Beam Theory (GBT) [22-25]  to obtain equally accurate results in a more efficient (lower d.o.f.’s required) and mostly clarifying manner (modal approach instead of nodal approach). The GBT-based analysis requires two consecutive and independent main tasks, namely: (i) cross-section analysis, which is devoted to the determination of structurally meaningful d.o.f.’s, i.e., the so-called GBT “deformation modes”, and (ii) member analysis, which consists of solving the appropriate differential equilibrium equation system (e.g., 1st order elastic and linear buckling analyses in the present work) enabling to unveil and quantify the contributions of each deformation mode thus making it possible to acquire a much deeper insight on that structural response. GBT has been used to analyse several buckling phenomena occurring in thin-walled structures but it has never been applied to analyse the wrinkling phenomenon, despite stress concentration and nonuniform longitudinal pre-buckling stresses have been studied in the past, e.g., [26-29]. The present paper presents a new modal theory to investigate the wrinkling behaviour of stretched fixed-ended elastic thin sheets, based on GBT concepts. The work begins in Section 2 by presenting the GBT fundamental concepts, followed by describing the main steps involved in the determination of its “deformation modes” in Section 3. Section 4 describes the procedures involved in the development and implementation of the modal formulation that accounts for prewrinkling (Section 4.1) and wrinkling analysis (Section 4.2). Then, an illustrative example is presented and discussed in detail in Section 5, comprising the characterisation of the (i) sheet “signature curve” (tensile critical stress 3 vs. length), (ii) critical wrinkling modes, (iii) modal participation diagrams and also (iv) pre-wrinkling stress distribution patterns – for comparison and validation purposes, the eigenvalues provided by refined ABAQUS [30] shell finite element results are also reported. Lastly, the paper closes highlighting the most relevant conclusion drawn from this investigation.

2. GBT fundamental concepts 2.1 Kinematics of thin sheets Consider the thin sheet of length 𝐿 depicted in Fig. 1(a), with a rectangular cross-section of width 𝑊 and height ℎ = 𝑡 (i.e. the sheet thickness 𝑡). The coordinate system (𝑥, 𝑠, 𝑧) is shown in Fig. 1(b), where 0 ≤ 𝑥 ≤ 𝐿, 0 ≤ 𝑠 ≤ 𝑊 and − 𝑡⁄2 ≤ 𝑧 ≤ 𝑡⁄2 and the displacement components in 3

The tensile critical stress and tensile critical strain correlate proportionally according to the elastic law. In this work, we use both in some figures. 6

the system are (𝑢, 𝑣, 𝑤 ). According to Love-Kirchhoff’s hypothesis ( 𝜀𝑧𝑧 = 𝛾𝑥𝑧 = 𝛾𝑧𝑥 = 𝛾𝑠𝑧 = 𝛾𝑧𝑠 = 0), the displacement vector of an arbitrary point is given by 𝑢(𝑥, 𝑠, 𝑧) = 𝑢 − 𝑧𝑤,𝑥 𝑣 (𝑥, 𝑠, 𝑧) = 𝑣 − 𝑧𝑤,𝑠

(1)

𝑤(𝑥, 𝑠, 𝑧) = 𝑤 where (i) 𝑢, 𝑣 and 𝑤 are the mid-plane (𝑧 = 0) displacement field components, and (ii) (∙),𝑥 = 𝜕⁄𝜕𝑥 and (∙),𝑠 = 𝜕⁄𝜕𝑠. 𝑠 (𝑣)

plan view

3D view

𝑠 (𝑣)

𝑡

𝑧 (𝑤)

𝑊 𝑥 (𝑢)

𝑥 (𝑢)

𝑧 (𝑤) 𝐿

𝑧 (𝑤)

cross-section

𝑡 𝑊

𝑠 (𝑣)

(a) (b) Figure 1. Geometrical configuration of thin sheet (dimensions, coordinate system, displacement components)

Moreover, in accordance with the classical thin-walled bar theory [31], the mid-plane displacement field can be conveniently expressed as (summation convention applies to subscript 𝑘) 𝑢(𝑥, 𝑠) = 𝑢𝑘 (𝑠)𝜑𝑘,𝑥 (𝑥 )

(2a)

𝑣 (𝑥, 𝑠) = 𝑣𝑘 (𝑠)𝜑𝑘 (𝑥 )

(2b)

𝑤(𝑥, 𝑠) = 𝑤𝑘 (𝑠)𝜑𝑘 (𝑥 )

(2c)

where (i) 𝑢𝑘 (𝑠) , 𝑣𝑘 (𝑠) and 𝑤𝑘 (𝑠) are the mid-line functions defining the deformation mode 𝑘 (or “GBT mode 𝑘 ”) of the sheet cross-section (unknown at this stage), (ii) 𝜑𝑘 (𝑥) or 𝜑𝑘,𝑥 (𝑥 ) are the amplitude functions4 describing their variation along the sheet length and (iii) 1 ≤ 𝑘 ≤ 𝑁, where 𝑁 is the total number of deformation modes employed. This way, for a given value of 𝑘, Eqs. (2a)-(2c) define the contribution of mode 𝑘 to the mid-plane displacement field. Therefore, the sheet deformed configuration can be expressed as a sum of contributions from the 𝑁 deformation modes  Eqs. (2a)-(2c) can also be written in matrix form as

4

The reason why Eq. (2a) involves one order of differentiation above Eqs. (2b)-(2c) is because that is a necessary (though not 𝑀 sufficient) condition for the validity of Vlasov’s hypothesis (i.e., 𝛾𝑥𝑠 = 0).

7

𝑢 = 𝒖𝑻 𝝋,𝑥

𝑣 = 𝒗𝑻 𝝋

𝑤 = 𝒘𝑻 𝝋

(3)

where (i) u, v and w are vectors containing the 𝑢𝑘 (𝑠), 𝑣𝑘 (𝑠) and 𝑤𝑘 (𝑠) functions, respectively, and (ii) 𝝋 is a vector containing the corresponding amplitude functions 𝜑𝑘 (𝑥 ). The linear strains associated with this displacement field are provided by 𝑀 𝐹 −𝑧𝑤,𝑥𝑥 𝑢,𝑥 𝜀𝑥𝑥 𝜀𝑥𝑥 𝑀 𝐹 𝑀 𝐹 𝛆 = 𝛆 + 𝛆 = { 𝜀𝑠𝑠 } + { 𝜀𝑠𝑠 } = { 𝑣,𝑠 } + { −𝑧𝑤,𝑠𝑠 } 𝑀 𝐹 𝑢,𝑠 + 𝑣,𝑥 −2𝑧𝑤,𝑥𝑠 𝛾𝑥𝑠 𝛾𝑥𝑠

(4)

where superscripts (∙)𝑀 and (∙)𝐹 denote membrane and flexural terms, respectively. After incorporating (2a)-(2c) into (4), it becomes possible to express the strain components in terms of the deformation modes (1 ≤ 𝑘 ≤ 𝑁𝑑 ), as 𝑢𝑘 𝜑𝑘,𝑥𝑥 𝑣𝑘,𝑠 𝜑𝑘 } 𝛆𝑀 = { (𝑢𝑘,𝑠 + 𝑣𝑘 )𝜑𝑘,𝑥

𝑤𝑘 𝜑𝑘,𝑥𝑥 𝛆𝐹 = −𝑧 { 𝑤𝑘,𝑠𝑠 𝜑𝑘 } 2𝑤𝑘,𝑠 𝜑𝑘,𝑥

(5)

2.2 Elastic strain energy – linear stiffness matrices For an isotropic linear elastic material under plane stress state, the stress field components are obtained from their strain field counterparts by means of the relations 𝑀 𝐹 1 𝜈 𝜎𝑥𝑥 𝜎𝑥𝑥 𝐸 𝑀} + { 𝐹 } = 𝛔 = 𝛔𝑀 + 𝛔𝐹 = { 𝜎𝑠𝑠 [𝜈 1 𝜎𝑠𝑠 1−𝜈2 𝑀 𝐹 0 0 𝜏𝑥𝑠 𝜏𝑥𝑠

0 0 ] (𝛆𝑀 + 𝛆𝐹 ) 1−𝜈

(6)

2

where E and ν are the material Young’s modulus and Poisson’s ratio, respectively. The member elastic strain energy 𝑈, which plays a pivotal role in several variational principles and energy methods employed to solve a wide variety of structural problems, is given by (𝑉 is the member volume) 1

𝑈 = 2 ∫𝑉 (𝛔𝑀 ∙ 𝛆𝑀 + 𝛔𝐹 ∙ 𝛆𝐹 )𝑑𝑉

(7)

and, after the incorporation of Eqs. (5)-(6), can be expressed in terms of the deformation modes as 1

𝐿

𝑈 = 2 ∫0 (𝛗𝑇,𝑥𝑥 𝐂𝛗,𝑥𝑥 + 𝛗𝑇,𝑥 𝐃𝛗,𝑥 + 𝛗𝑇 𝐁𝛗 + 𝛗𝑇,𝑥𝑥 𝐄𝛗 + 𝛗𝑇 𝐄 𝑇 𝛗,𝑥𝑥 )𝑑𝑥

(8)

in which 𝑪 , 𝑩 , 𝑫 and 𝑬 are 𝑁 × 𝑁 linear stiffness matrices, associated with several cross-section mechanical properties, namely (i) primary/secondary warping, (ii) transverse extension/flexure, (iii) sheet

8

shear distortion/torsion and (iv) membrane/flexural Poisson effects5 – their components are given by 𝑊 𝐸𝑡

𝑀 𝐹 𝐶𝑖𝑘 = 𝐶𝑖𝑘 + 𝐶𝑖𝑘 = ∫0

𝑊 𝐸𝑡

𝑀 𝐹 𝐵𝑖𝑘 = 𝐵𝑖𝑘 + 𝐵𝑖𝑘 = ∫0

1−𝜈2

𝑊

𝑀 𝐹 𝐷𝑖𝑘 = 𝐷𝑖𝑘 + 𝐷𝑖𝑘 = ∫0

𝑤𝑖 𝑤𝑘 𝑑𝑠

𝐸𝑡 3

𝑊

12(1−𝜈2 )

(9a)

𝑤𝑖,𝑠𝑠 𝑤𝑘,𝑠𝑠 𝑑𝑠 𝑊

2(1+𝜈) 1−𝜈2

12(1−𝜈2 )

𝑣𝑖,𝑠 𝑣𝑘,𝑠 𝑑𝑠 + ∫0

𝐸𝑡

𝑊 𝜈𝐸𝑡

𝑀 𝐹 𝐸𝑖𝑘 = 𝐸𝑖𝑘 + 𝐸𝑖𝑘 = ∫0

𝐸𝑡 3

𝑊

𝑢 𝑢 𝑑𝑠 + ∫0 1−𝜈2 𝑖 𝑘

(𝑢𝑖,𝑠 + 𝑣𝑖 )(𝑢𝑘,𝑠 + 𝑣𝑘 )𝑑𝑠 + ∫0 𝜈𝐸𝑡 3

𝑊

𝑢𝑖 𝑣𝑘,𝑠 𝑑𝑠 + ∫0

12(1−𝜈2 )

𝑤𝑖 𝑤𝑘,𝑠𝑠 𝑑𝑠

𝐸𝑡 3 6(1+𝜈)

(9b) 𝑤𝑖,𝑠 𝑤𝑘,𝑠 𝑑𝑠

(9c) (9d)

The tensor components appearing in Eqs. (9a)-(9d) are determined on the basis of the cross-section deformation modes, i.e., functions 𝑢𝑘 (𝑠), 𝑣𝑘 (𝑠) and 𝑤𝑘 (𝑠), which are still unknown at this point. The next section concerns their determination.

3. Cross-section analysis As mentioned before, the first step of a GBT structural analysis is the determination of the cross-section deformation modes (𝑢𝑘 (𝑠), 𝑣𝑘 (𝑠) and 𝑤𝑘 (𝑠) functions) of the sheet and their associated mechanical properties (𝐶𝑖𝑘 , 𝐵𝑖𝑘 , 𝐷𝑖𝑘 and 𝐸𝑖𝑘 components), which is done through a systematic procedure termed cross-section analysis. This work is based on a recently developed version of this procedure, which is applicable to arbitrary flat-walled members and is described in detail in [32,33]. Due to space limitations and the complexity of the procedure, only a very brief overview is provided here. It starts with the specification of a discretisation of the sheet cross-section6 with 𝑛 nodes along the its width b ( 0 ≤ 𝑠 ≤ 𝑊 ), where three elementary modes are associated with each node: they are associated with the imposition of three unit displacements, along 𝑥, 𝑠 and 𝑧 directions, which yields a total of 𝑁 = 3 × 𝑛 such modes. Matrices given by tensorial components (9a)-(9b) are then calculated on the basis of these elementary modes and become fully populated. Then, these matrices are involved in a series of simultaneous diagonalisation operations – which make it possible to (i) identify/separate modes from different mechanical families and (ii) arrange them within a given family, according to some hierarchy – leading to a new set of “coordinates” which are the 𝑁 deformation modes of GBT. The three main “mechanical families” involved are:  Vlasov (or conventional) modes (for which 𝛾𝑥𝑠 = 𝜀𝑠𝑠 = 0),  Shear modes (for which 𝛾𝑥𝑠 ≠ 0; 𝜀𝑠𝑠 = 0),  Transverse extension modes (for which 𝜀𝑠𝑠 ≠ 0).

5

6

Matrix 𝐄 is associated with the work done by the normal stresses stemming from Poisson effects (see the constitutive relation defined by Eq. (7)) on the strain field 𝛆. In Schardt’s original notation [22] and many current papers, the flexural component 𝐄 𝐹 is referred to as “𝐃2 ” and included in matrix 𝐃, given by 𝐃 = 𝐃1 + (𝐃2 + 𝐃𝑇2 ), where 𝐃1 ≡ 𝐃𝐹 . The number, nature and quality of the deformation modes obtained depend on this nodal discretisation.

9

For the case of single-wall sections (plates), these 3 main families can be further sub-divided into 2 subfamilies [30]: global and local modes, which represent global or localised deformation and stem from natural and intermediate nodes, respectively. The 6 (3 × 2) global modes of a plate section are: (i) the four classical rigid-body modes (axial extension, major and minor-axis bending and torsion – Fig. 2, modes 1-4), (ii) a linear symmetrical warping distribution (Fig. 2, mode 8) and (iii) an isotropic extension stemming from the centroid (Fig. 2, mode 12). For illustrative purposes, consider the sheet cross-section discretised into 4 elements of equal width (each with 𝑊/4). In this case, we have 5 nodes (𝑛 = 5), being 2 edge (natural) nodes and 3 equally-spaced internal nodes. Fig. 2 displays the in and/or out-of-plane (warping) shapes of the 𝑁 = 5 deformation modes, divided into the main families. The first 4 modes are the classical rigid-body modes: axial extension (1), major and minor-axis bending (2-3) and torsion (4) – for major-axis bending, the in-plane displacements are accompanied by the associated primary warping counterpart. Local Vlasov modes (57) involve transverse flexural curvatures with increasing number of half-waves. Shear modes (8-11) comprise only warping displacements – the global shear mode (8) is the warping component of majoraxis bending (2) taken alone. Finally, the isotropic extension mode (12) and the local extension modes (13-15), involving extension/contraction of the wall segments while preserving the overall width. At this point, depending on the particular problem under consideration, a selection can be made involving any sub-set of 𝑁𝑑 (1 ≤ 𝑁𝑑 ≤ 𝑁) deformation modes to be included in the analysis (i.e., to be part of the problem solution). This modal selection capability makes it possible to reduce the number of degrees of freedom involved in the solution of the problem and specify the nature (mode family) of the deformation pattern(s) to be considered. 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

Vlasov

In-Plane Warping

Transv. Extension

Shear

In-Plane

Warping

In-Plane

Figure 2. GBT deformation modes for a single plate section: Vlasov modes – global (1-4), and local (5-7) , Shear modes – global (8) and local (9-11) , and Transverse extension modes – isotropic (12) and local (13-15).

10

4. Member analysis Once determined the cross-section deformation modes and respective modal mechanical properties, the next step is to perform the member analysis, leading to the solution of the sheet wrinkling problem: wrinkling loads (eigenvalues) and associated mode shapes (eigenvectors), defined in terms of the 𝜑𝑘 (𝑥 ) functions (1 ≤ 𝑘 ≤ 𝑁𝑑 ). The wrinkling of thin sheets under uniaxial tension is characterized by a non-uniform stress state, which is in the origin of wrinkle formation. This requires the performance of a preliminary linear analysis (pre-wrinkling) to identify the nonuniform stress state. Then, using this stress state as input, a wrinkling analysis is done by considering the geometrical stiffness degradation associated with compressive stresses. 4.1 Pre-wrinkling analysis The first-order equilibrium equations may be established by means of the Principle of Stationary Potential Energy, which reads 𝛿 (𝑈 + 𝛱 ) = 0

(10)

where 𝛱 is the potential of the applied loads. In order to define this potential, let 𝑄𝑥 be a load distributed over the sheet mid-plane and acting along the 𝑥 direction. It is further assumed that 𝑄𝑥 may be expressed by 𝑄𝑥 (𝑥, 𝑠) = 𝑞𝑥 (𝑠)𝜙(𝑥)

(11)

where 𝑞𝑠 (𝑠) describes the variation of loading profile along the sheet width (0 ≤ 𝑠 ≤ 𝑊) and 𝜙(𝑥 ) is a function that describes the variation of 𝑞𝑠 along the sheet length (0 ≤ 𝑥 ≤ 𝐿). Then, the potential of this applied load can be expressed as 𝐿

𝑊

𝐿

𝑊

𝛱 = − ∫0 ∫0 𝑄𝑥 𝑢𝑑𝑠𝑑𝑥 = − ∫0 (∫0 𝑞𝑥 𝑢𝑑𝑠) 𝜙𝑑𝑥

(12)

Introducing (2) into (12) and integrating along the sheet width (0 ≤ 𝑠 ≤ 𝑊), it is obtained 𝐿

𝛱 = − ∫0 𝛗𝑻 𝐪𝜙𝑑𝑥

(13)

where the components of vector 𝒒 are generalised modal forces given by 𝑊

𝑞𝑖 = ∫0 𝑞𝑥 𝑢𝑖 𝑑𝑠

(14)

Fig. 3 shows an illustrative example of a load 𝑄𝑥 (𝑥, 𝑠) that is uniform along the sheet width (𝜕𝑞𝑥 /𝜕𝑠 = 0) and varies nonlinearly along the sheet length (according to a prescribed function 𝜙(𝑥 )).

11

Figure 3. Illustrative example of a load 𝑄𝑥 (𝑥, 𝑠) corresponding to uniform 𝑞𝑥 (𝑠) = 𝑞 and nonlinear 𝜙(𝑥).

Finally, introducing (8) and (12) into (10), one is led to the weak form of the member first-order equilibrium equations, which reads ∫𝐿 (𝛿𝛗𝑇,𝑥𝑥 𝐂𝛗,𝑥𝑥 + 𝛿𝛗𝑇,𝑥 𝐃𝛗,𝑥 + 𝛿𝛗𝑇 𝐁𝛗 + 𝛿𝛗𝑇,𝑥𝑥 𝐄𝛗 + 𝛿𝛗𝑇 𝐄 𝑇 𝛗,𝑥𝑥 − 𝛿𝛗𝑻 𝐪𝜙)𝑑𝑥 = 0

(15)

Then, a general pre-wrinkling analysis can be performed by means of (15). Its solution is the vector 𝛗𝟎 , whose components are the pre-wrinkling modal amplitude functions 𝜑𝑗0 (𝑥 ). On its basis, the pre-wrinkling deformed configuration is given by (3) and reads 𝑢0 = 𝐮𝑇 𝛗𝟎,𝑥

𝑣 0 = 𝐯 𝑇 𝛗𝟎

𝑤 0 = 𝐰 𝑇 𝛗𝟎

(16)

As for the membrane components of the pre-wrinkling stresses, they are provided by the expressions 𝐸

0 𝜎𝑥𝑥 = 1−𝜈2 (𝐮𝑇 𝛗𝟎,𝒙 + 𝜈𝐯,𝑠𝑇 𝛗𝟎 ) 𝐸

0 𝜎𝑠𝑠 = 1−𝜈2 (𝐯,𝑠𝑇 𝛗𝟎 + 𝜈𝐮𝑇 𝛗𝟎,𝒙 ) 𝐸

0 𝜏𝑥𝑠 = 2(1+𝜈) (𝐮𝑇,𝑠 𝛗𝟎 + 𝐯 𝑇 𝛗𝟎,𝒙 )

(17a) (17b) (17c)

where the vectors generally contain the 𝑁𝑑 GBT deformation modes obtained from the crosssection analysis – while this is obviously the more accurate way to represent the stress field, it also possible to achieve quite good results by considering only a given subset of them (𝑛𝑑0 ≤ 𝑁𝑑 ), resulting in a pre-wrinkling stress field combining only some specific patterns. 4.2 Wrinkling analysis In wrinkling analyses, the equilibrium equations may also be established by means of the Principle of Stationary Potential Energy, which now reads 𝛿 (𝑈 + 𝜆𝑈 0 ) = 0

12

(18)

where 𝑈 is the linear strain energy given by (8). Additionally, 𝑈 0 is strain energy standing for the work performed by the pre-wrinkling stresses, defined in (17a)-(17c), on the non-linear (quadratic) strain components, given by (𝑣𝑖 𝑣𝑘 + 𝑤𝑖 𝑤𝑘 )𝜑𝑖,𝑥 𝜑𝑘,𝑥 1 𝑤𝑖,𝑠 𝑤𝑘,𝑠 𝜑𝑖 𝜑𝑘 } 𝛆𝑁𝐿.𝑀 = 2 { (19) 2(𝑣𝑖 𝑣𝑘,𝑠 + 𝑤𝑖 𝑤𝑘,𝑠 )𝜑𝑖,𝑥 𝜑𝑘 and 𝜆 is a load parameter. The final form of 𝑈 0 is given by, 0 𝛗𝑇,𝑥 [𝐗 𝛔−𝐱 𝜑𝑗,𝑥𝑥 + 𝐗 𝛔−𝐱𝐏 𝜑𝑗0 ]𝛗,𝑥 + 𝐣 𝐣 𝐿

1 𝑁𝑑 0 𝜑𝑗,𝑥𝑥 + 𝐗 𝛔−𝐬 𝜑𝑗0 ]𝛗 + } 𝑑𝑥 𝑈 0 = 2 ∫0 ∑𝑗=1 { +𝛗𝑇 [𝐗 𝛔−𝐬𝐏 𝐣 𝐣 0 ]𝛗,𝑥 + +𝛗𝑇 [𝐗 𝛕𝐣 𝜑𝑗,𝑥

(20)

0 ]𝛗 𝛗𝑇,𝑥 [𝐗 𝛕𝐣 𝑻 𝜑𝑗,𝑥

where 1 ≤ 𝑗 ≤ 𝑁𝑑0 and 𝑿𝝈−𝒙 , 𝑿𝝈−𝒙𝑷 , 𝑿𝝈−𝒔𝑷 , 𝑿𝝈−𝒔 and 𝑿𝝉𝒋 are the cross-section geometric 𝒋 𝒋 𝒋 𝒋 stiffness matrices associated with (i) normal longitudinal (𝑿𝝈−𝒙 and 𝑿𝝈−𝒙𝑷 ), (ii) normal transverse 𝒋 𝒋 (𝑿𝝈−𝒔𝑷 and 𝑿𝝈−𝒔 ) and (iii) shear (𝑿𝝉𝒋 ) pre-wrinkling stresses7, which read 𝒋 𝒋 𝑊 𝐸𝑡

𝜎−𝑥 𝑋𝑗𝑖𝑘 = ∫0

𝑢 (𝑣 𝑣 1−𝜈2 𝑗 𝑖 𝑘

+ 𝑤𝑖 𝑤𝑘 )𝑑𝑠

𝑊 𝜈𝐸𝑡

𝜎−𝑥𝑃 𝑋𝑗𝑖𝑘 = ∫0

𝑣 (𝑣 𝑣 1−𝜈2 𝑗,𝑠 𝑖 𝑘

(21a)

+ 𝑤𝑖 𝑤𝑘 )𝑑𝑠

(21b)

𝑊 𝐸𝑡

𝜎−𝑠 𝑋𝑗𝑖𝑘 = ∫0

𝑣 𝑤 𝑤 𝑑𝑠 1−𝜈2 𝑗,𝑠 𝑖,𝑠 𝑘,𝑠

(21c)

𝑊 𝜈𝐸𝑡

𝜎−𝑠𝑃 𝑋𝑗𝑖𝑘 = ∫0 𝑊

𝜏 𝑋𝑗𝑖𝑘 = ∫0

𝑢 𝑤 𝑤 𝑑𝑠 1−𝜈2 𝑗 𝑖,𝑠 𝑘,𝑠

𝐸𝑡 2(1+𝜈)

(21d)

(𝑢𝑗,𝑠 + 𝑣𝑗 )(𝑣𝑖,𝑠 𝑣𝑘 + 𝑤𝑖,𝑠 𝑤𝑘 )𝑑𝑠

(21e)

Note that the 𝜑𝑗0 (𝑥 ) functions appearing in Eq. (20), which define the pre-wrinkling state, are already known at this point, after having solved Eq. (15). Then, adding 𝑈 (given by (8)) and 𝑈 0 (given by (20)), and using (18), the equilibrium condition of the eigenfunction problem is obtained, in which (i) the eigenvalues of parameter 𝜆 correspond to the wrinkling loads and (ii) the eigenfunctions 𝛗(𝑥 ) correspond to the wrinkling modes: 𝛿𝛗𝑇,𝑥𝑥 𝐂𝛗,𝑥𝑥 + 𝛿𝛗𝑇,𝑥 𝐃𝛗,𝑥 + 𝛿𝛗𝑇 𝐁𝛗 + 𝛿𝛗𝑇,𝑥𝑥 𝐄𝛗 + 𝛿𝛗𝑇 𝐄 𝑇 𝛗,𝑥𝑥 + 𝑁

𝑑 0 {𝛿𝛗𝑇,𝑥 [𝐗𝑗σ−x 𝜑𝑗,𝑥𝑥 +𝜆 ∑𝑗=1 + 𝐗𝑗σ−xP 𝜑𝑗0 ]𝛗,𝑥 +

∫𝐿

𝑑𝑥 = 0

0 +𝛿𝛗𝑇 [𝐗𝑗σ−sP 𝜑𝑗,𝑥𝑥 + 𝐗𝑗σ−s 𝜑𝑗0 ]𝛗 + 𝑇

(

7

0 0 ]𝛗,𝑥 + 𝛗𝑇,𝑥 [(𝐗𝑗τ ) 𝜑𝑗,𝑥 ] 𝛗} +𝛿𝛗𝑇 [𝐗𝑗τ 𝜑𝑗,𝑥

Note that the 𝐗 (∙)𝐏 terms stem from Poisson effects.

13

)

(22)

4.3 Implementation of a 1D finite element In a previous work, Silvestre [11] used a modal theory to determine analytically an approximate pattern of stresses, strains and displacements in the pre-wrinkling state of stretched sheets. Despite giving a fair approximation to the nonlinear distribution of pre-wrinkling stresses, the analytical approach proposed by [11] is restricted to three modes and an exact distribution usually requires a higher number of modes. In order to obtain such exact stress patterns, both the pre-wrinkling (Eq. (15)) and the wrinkling (Eq. (22)) analyses of the sheet should be performed numerically using (i) several deformation modes to approximate the deformed shape of sheet cross-sections and (ii) several one-dimensional (1D) finite elements. The degrees of freedom involved in such a one-dimensional (1D) finite element are the nodal values and derivatives of the 𝜑𝑘0 (𝑥 ) or 𝜑𝑘 (𝑥 ) functions, approximated within the element by means of cubic Hermite8 polynomials, ℎ𝑖 (𝜉 ), i.e., 𝛗(𝑥 ) = ℎ1 (𝜉 )𝒅1𝑒 + ℎ2 (𝜉 )𝒅𝑒2 + ℎ3 (𝜉 )𝒅𝑒3 + ℎ4 (𝜉 )𝒅𝑒4

(23)

where 𝜉 = 𝑥⁄𝐿𝑒 (𝐿𝑒 is the element length) and the cubic polynomials are defined as ℎ1 (𝜉 ) = 𝐿𝑒 (𝜉 3 − 2𝜉 2 + 𝜉 ) ℎ3 (𝜉 ) = 𝐿𝑒 (𝜉 3 − 𝜉 2 )

ℎ2 (𝜉 ) = 2𝜉 3 − 3𝜉 2 + 1 ℎ4 (𝜉 ) = −2𝜉 3 + 3𝜉 2

(24a) (24b)

and vectors 𝒅𝑒𝑖 contain the degrees of freedom: 𝒅1𝑒 ≡ 𝛗,𝑥 (0), 𝒅𝑒2 ≡ 𝛗(0), 𝒅𝑒3 ≡ 𝛗,𝑥 (1) and 𝒅𝑒4 ≡ 𝛗(1). Incorporating (23)-(24) into (15) and performing the integrations (over 𝐿𝑒 ), the finite element linear stiffness matrix and force vector for pre-wrinkling analysis are obtained: 𝑒 𝑲𝑒 = 𝑪 ⊗ 𝒌𝑒22 + 𝑫 ⊗ 𝒌11 + 𝑩 ⊗ 𝒌𝑒00 + 𝑬 ⊗ 𝒌𝑒20 + 𝑬𝑇 ⊗ 𝒌𝑒02

(25)

𝒇𝑒 = 𝒒𝑥𝑋 ⊗ 𝒇1𝑒

(26)

where 4 × 4 matrices 𝒌𝑒𝑝𝑞 and 4 × 1vector 𝒇1𝑒 , which components are integrals of products of Hermite polynomials or derivatives, are given in Appendix. Performing the assembly of the global stiffness matrix (𝑲) and force vector (𝐟), and taking into account the member support conditions, the global pre-wrinkling equilibrium system is obtained, 𝑲𝒂 = 𝒇

(27)

where 𝒂 is a vector containing the global generalised displacements, which leads to the definition of 𝝋𝟎 (𝑥 ). Likewise, the introduction of (23)-(24) into (22) leads to the determination of the finite element linear (Eq. (25)) and geometrical stiffness matrices for the wrinkling analysis. The geometrical stiffness matrix is expressed in the form

8

While Hermite polynomials could be used to approximate all the 𝜑𝑘 (𝑥) functions, it is usual to use cubic Lagrange polynomials for those modes involving only warping displacements, namely (i) axial extension (mode 1) and (ii) global and (iii) local shear modes.

14

𝑒 𝑒 ) + 𝑿𝝈−𝒙𝑷 )+ 𝑿𝝈−𝒙 ⊗ (𝑐𝑗2.𝑖 𝒈11.𝑖 ⊗ (𝑐𝑗0.𝑖 𝒈11.𝑖 𝒋 𝒋 𝑁𝑑 ⊗ (𝑐𝑗2.𝑖 𝒈𝑒00.𝑖 ) + 𝑿𝝈−𝒔 ⊗ (𝑐𝑗0.𝑖 𝒈𝑒00.𝑖 ) +} ∑3𝑖=0 {+𝑿𝝈−𝒔𝑷 𝑮𝑒 = ∑𝑗=1 𝒋 𝒋

+𝑿𝝉𝒋

⊗

(𝑐𝑗1.𝑖 𝒈𝑒01.𝑖 )

+

𝑿𝝉𝒋 𝑻

⊗

(28)

𝑒 (𝑐𝑗1.𝑖 𝒈10.𝑖 )

in which (i) 𝑐𝑗𝑘.𝑖 is the coefficient of 𝑖-th order (𝑖 = 0, … ,3) term of function (𝜑𝑗0 )

(𝑘)

(𝑘 is the order of

derivation, 𝑘 = 0, … ,2) – e.g., 𝜑𝑗0 (𝜉 ) = ∑3𝑖=0 𝑐𝑗0.𝑖 𝜉 𝑖 – and (ii) 𝐠 𝑒𝑝𝑞.𝑖 are 4 × 4 matrices. The assembly of the linear (𝑲) and geometrical (𝑮) stiffness matrices leads to the eigensystem, (𝑲 − 𝜆𝑮)𝒃 = 𝟎

(29)

that governs the wrinkling problem, where 𝜆 (eigenvalue) in the wrinkling load/factor and 𝒃 (eigenvector) leads to the definition of the wrinkling mode, 𝛗(𝑥 ). The sheet mid-surface is thus (i) longitudinally discretised into GBT-based finite elements and (ii) transversally discretized into “wall segments” (limited by consecutive nodes), in a manner similar to a shell finite element mesh. Finally, one word to mention that because different sets of GBT modes may be used in both pre-wrinkling and wrinkling analyses (totalling 𝑛𝑑0 and 𝑛𝑑 modes, respectively), the dimensions of systems (27) and (29) need not be the same – upper bound estimates on those dimensions are 2(𝑛𝑓𝑒 + 1)𝑛𝑑0 and 2(𝑛𝑓𝑒 + 1)𝑛𝑑 , respectively, where 𝑛𝑓𝑒 is the number of finite elements.

5. Illustrative examples This section presents illustrative examples of the formulation. In this study, the main aim is to investigate the effect of the (i) sheet geometry (i.e., the aspect ratio L/W) and (ii) the Poisson’s ratio on the critical tensile stress cr – the sheets are made of an elastic isotropic material with Young’s modulus E=1.8GPa. Two values of the Poisson’s ratio were considered: =0.30 or 0.509. With this purpose, a sheet with width W=100mm, thickness t=0.1mm and length L was considered, which varies between L=100mm (short sheet, L/W=1.0) and L=1000mm (long sheet, L/W=10.0) – in total 30 distinct lengths were considered to provide a smooth cr vs. L curve. Due to distinct nature of the wrinkling modes (discussed later), it was necessary to consider different nodal discretisations, which involved 30 (in 25 L/W ratios), 40 (in 3 L/W ratios) or 50 (in 2 L/W ratios) intermediate nodes (equally spaced). Since only conventional modes participate in the wrinkling mode shapes, independently on the L/W ratio, Fig. 4(a)+(b) shows the inplane displacements of all these modes yielded by two distinct nodal discretisations10: 30 (Fig. 4(a)) and

9

10

The Poisson’s ratio =0.5 corresponds to an infinite value of the elastic bulk modulus and, therefore, to an incompressible material. In the adopted shell finite element S4R (Abaqus/Standard) under plane stress state, the material is free to deform in the thickness direction. In this case, special treatment of the volumetric behaviour is not necessary as the use of regular stress/displacement elements has been proven satisfactory. The GBT conventional deformation modes concerning a nodal discretisation with 50 intermediate nodes are omitted in this document.

15

40 (Fig. 4(b)) intermediate nodes – note that, the local modes 5 to 34 (Fig. 4(a)) or 44 (Fig. 4(b)) are the most relevant in the GBT wrinkling modes. 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

(a) 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

(b) Figure 4: GBT conventional deformation modes for a single plate discretised with (a) 30 and (b) 40 intermediate equally spaced nodes.

Concerning the GBT boundary conditions, the d.o.f.’s at the end cross-sections concerning the axial extension and shear modes are free while those associated with the remaining conventional and transverse extension modes are prevented, i.e., the end cross-sections are locally/globally fixed and free to warp. In the ABAQUS SFEA counterparts, all the d.o.f.s at the cross-section ends are prevented with exception of the axial displacements that are free to deform. In order to preclude longitudinal rigid-body motions, the axial extensional d.o.f. at mid-span and the central node were prevented in the GBT and 16

SFEA models, respectively – in both cases the longitudinal edges are free to deform. Since the load is applied at both ends, the models (GBT and SFEA) are symmetric at mid-span. Regarding the numerical determination of the eigenvalue problem (Eq. (29)) it is worth mentioning that Lanczos method was employed in ABAQUS SFEA by setting the lower bound to zero, i.e., only positive eigenvalues are determined. In GBT analysis, all the pairs of eigenvalues and eigenvectors were determined in some computational inexpensive cases (reduced number of d.o.f.s), e.g., in sheets with aspect ratios of 1.4≤L/W<4.0, while in the remaining cases a similar strategy akin to ABAQUS analyses was conducted. Fig. 5(a) shows (i) two GBT-based “signature curves” (cr or cr vs. L/W curve) of the fixed-ended stretch-induced thin sheet corresponding to the two Poisson’s ratios (0.30 and 0.50) with the corresponding ABAQUS SFEA results (for validation purposes only), and (ii) representative 2D configurations of relevant critical (symmetric) wrinkling modes while Fig. 5(b) presents the corresponding GBT modal participation diagrams (pk vs. L/W) for the sheet with v=0.50. For L/W >1, the obtained critical modes appear in pairs (independently of the Poisson’s ratio), i.e., duplicate critical modes (one symmetric and one antisymmetric, with practically identical eigenvalues). This fact was also observed by Healey et al. [7], in case of a non-linear bifurcation shell finite element analysis of a stretched sheet with L/W=2 (L=200mm, W =100mm, =0.45). They specified that coincidence of bifurcation points corresponding to symmetric and antisymmetric modes as “striking” and remarked that “there is no generic explanation for the occurrence of double eigenvalue”. Therefore, two distinct modal participation diagrams are presented in Fig. 5(b): Fig. 5(b1) and Fig. 5(b2) correspond to the symmetric and antisymmetric mode shapes, respectively. The observation of the GBT results presented in Fig. 5 enables the following remarks: (i) First of all, there is an excellent agreement between the GBT and ABAQUS SFEA results (wrinkling values and modes). The GBT discretisation (nodal and/or longitudinal11) was increased until the difference in the critical wrinkling values between the very (very) refined SFEA results does not exceed 5% despite the significantly lower d.o.f.’s required – generally speaking, the GBT-based analysis requires at least one order of magnitude lower. For instance, the SFEA model in the sheet with L/W=1.0 and L/W=10.0 comprised about 250000 and 1700000 d.o.f.s while the GBT counterparts required approximately 16000. (ii) Fig. 5(a) shows that the “signature curves” of sheets with v=0.30 and 0.50 are qualitatively similar. In fact, the comparison between the ensuring critical wrinkling modes, for identical L/W ratios, showed that these shapes are practically indistinguishable.

11

All the GBT results were obtained by employing equally-spaced finite elements, this number varying between 50 and 180. The exceptions were the sheets with L/W=10.0, which were modelled using 180 finite elements, following a stepwise distribution along the x-axis: 72 in 0≤x≤ L/3, 36 in 0
17

cr (%)

=0.30 =0.50

GBT

cr (MPa)

=0.30 =0.50

SFEA

20

375

16

300

12

225

8

150

4

75

0

1

2

3

4

5

6

7

0

8 9 10

L/W (a) 34+36+38+40+42+44+46+48+50+52+54

33+35+37+39+41+43+45+47+49+51+53 29+31

21

25+27

17 15

30+32

11

23 19

18

28 26 24

13

16

14 12

20

10

22

9

8

7 1

1

2

pk (%)

17

100

15

19

14 80

11

60

16

100

13

80

2

pk (%) 12

60 10

40

40

9

8

20

20 7 0 1

2

5 3

4

5

6

6

0 1

7 8 9 10

2

3

4

5

6

7 8 9 10

L/W

L/W

(b1)

(b2)

Figure 5. Stretch-induced thin sheet (a) cr (or cr) vs. L/W curves (GBT and SFEA results with v=0.30 and 0.50) and (b) GBT modal participation diagrams (pk vs. L/W) for the (1) symmetric and (2) antisymmetric wrinkling modes (v=0.50).

(iii) The sheet “signature curves” show that, as the aspect ratio increases, (iii1) cr (or cr) decreases abruptly at L/W=1.00 until it reaches a minimum of 8.0MPa (v=0.50) or 16.2MPa (v=0.30) at L/W1.70, and (iii2) then increases moderately to 22.0MPa (v=0.50) or 45.7MPa (v=0.30) at L/W3.75. For higher aspect ratios (long sheets) the critical tensile stress maintains approximately that constant value, i.e., it is independent of the length. It is worth mentioning that there are no positive eigenvalues for L/W<1.0 or, in other words, y(x,s) is always positive along the sheet width and length, i.e., the “signature curve” has a vertical asymptote in the vicinity of L/W1.0. This means that rectangular sheets stretched along their lower dimension (L/W<1.0) do not wrinkle. It should be 18

outlined that Friedl et al. [4] also obtained a qualitatively similar “signature curve” for a different sheet (E=70GPa, =0.30, W=200mm, h=0.05mm). In this study, it is stated that transversal compressive stresses (and wrinkling) vanish for L/W<1.1, which evidences the current observation. In 2011, Nayyar et al. [1] also shown that stretched hyperelastic sheets with L/W<1.0 do not wrinkle due to the absence of compressive stresses. (iv) The modal participation diagrams of Fig. 5(b1) (critical symmetric modes) and Fig. 5(b2) (critical antisymmetric modes) are qualitatively very similar. Moreover, a closer inspection of these diagrams and the different mode participations pk shows that there are three mechanically distinct wrinkling modes associated to different L/W ratios:  equal/close to 1.0 – Type I  between 1.05 and 1.30 – Type II  higher than 1.40 – Type III Indeed, the wrinkling modes of shorter sheets involve a higher number of local modes and higher strain energy modes (thus, a refined nodal discretisation is needed) than those associated with larger aspect ratios. For instance, the most relevant participations in the critical wrinkling mode of sheet with L/W=1.0 (Type I) are: p15=2.6%, p17=4.1%, p19=3.6%, p21=16.4%, p23=8.3%, p25=17.4%, p27=14.4%, p29=10.5%, p31=10.3%, p33=4.0% and p35=3.9% (symmetric mode) or p14=1.3%, p16=4.3%, p18=2.6%, p20=12.0%, p22=18.1%, p24=6.0%, p26=23.6%, p28=1.3%, p30=16.1%, p32=2.8%, p34=6.3%, p36=1.2% and p38=1.7% (antisymmetric mode). As the aspect ratio increases to 1.15 (Type II), it can be observed that the critical symmetric and antisymmetric wrinkling modes may be obtained by employing solely modes 15+17+19+21 (p15+17+19+21=23.6+37.6+25.3+8.4=94.9%) or 14+16+18+20+22 (p14+16+18+20+22=10.4+33.3+33.1+15.6+3.5=95.9%), respectively. For L/W≥1.4 (Type III), Fig. 5(b) shows a drastic change in the mode nature since the corresponding wrinkling modes now involve mainly the contributions of modes 7+9+11+13+15 (symmetric mode) or 6+8+10+12+14+16 (antisymmetric mode) with significantly different relative participations between L/W=1.3 and 5.0. For L/W>5.0 the relative participation of the aforementioned modes remains virtually uniform: p7=15.6%, p9=30.0%, p11=28.8%, p13=16.7%, p15=6.1% (symmetric mode) and p6=6.3%, p8=24.3%, p10=31.5%, p12=23.3%, p14=10.7%, p16=3.2% (antisymmetric mode). Now, let us turn our attention to the longitudinal variation of deformation modes that participate in the wrinkling modes. Figs. 6 and 7 show GBT modal amplitude functions, 3D mode shape configuration and (at least) one representative 2D configuration, concerning the symmetric and antisymmetric wrinkling modes of the sheet with v=0.50 belonging to Type I (L/W=1.00, Fig. 6) and Type II (L/W=1.10, Fig. 7). However, note that an exception occurs for the quadrangular sheet (L/W =1) which exhibits 4 identical eigenvalues, i.e. quadruplicate critical modes as seen in Fig. 6. Regarding the Type III (L/W≥1.4) wrinkling modes, it was observed that the longitudinal 19

variation of these varied too much with the aspect ratio L/W. 0.1

0.1 15 17 19 21 23 25 27 29 31 33 35 37 39

0.05

0

0.05

0

-0.05

-0.1 0

16 18 20 22 24 26 30 32 34 38

-0.05

0.2

0.4

0.6

0.8

-0.1 0

1

0.2

0.4

(a1)

0.8

1

(a2)

0.1

0.1 14 18 20 22 24 26 28 30 32 34 36

0.05

0

-0.05

-0.05

0.2

0.4

0.6

0.8

17 19 21 23 25 27 29 31 33 35 37

0.05

0

-0.1 0

0.6

-0.1 0

1

0.2

0.4

(a3)

0.6

1

(a4)

x

x

x=0.5L

x

0.8

x=0.5L

x

(b1)

x=0.5L

(b2)

x=0.5L

(b3) (b4) Figure 6. L=100mm (L/W=1.0): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the quadruplicate critical mode shapes (1) mode 1 (symmetric), (2) mode 2 (anti-symmetric), (3) mode 3 (antisymmetric) and (4) mode 4 (symmetric) with v=0.50.

20

0.15

0.15 11 13 15 17 19 21 23 25 27

0.1

0.05

12 14 16 18 20 22 24 26 28

0.1

0.05

0

0

-0.05

-0.05

-0.1

-0.1

-0.15

-0.15 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

(a1)

0.6

0.8

1

(a2)

x

x

x=0.5L

x=0.5L

(b1) (b2) Figure 7. L=110mm (L/W=1.1): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the critical (1) symmetric and (2) antisymmetric mode shapes (v=0.50). 1.0

1.0 7 9 11 13 15 17

0.5

6 8 10 12 14 16

0.5

0

0

-0.5

-0.50

-1.0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

(a1)

0.6

0.8

1

(a2)

x

x

x=0.5L

x=0.5L

(b1) (b2) Figure 8. L=170mm (L/W=1.7): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the critical (1) symmetric and (2) antisymmetric mode shapes (v=0.50). 21

1

1.5

(a1)

(a2)

0.6

1

0.2

0.5

-0.2

6 8 10 12 14 16

0 5 7 9 11 13 15

-0.6

-0.5

-1

-1 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

x

x

x=0.50L

x=0.50L

(b1) (b2) Figure 9. L=260mm (L/W=2.6): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the critical (1) symmetric and (2) antisymmetric mode shapes (v=0.50). 1

1

(a1)

(a2)

0.5

0.5

0

0 5 7 9 11 13 15 17

-0.5

-0.5

-1 0

0.2

0.4

6 8 10 12 14 16 18

0.6

0.8

-1

1

0

x

0.2

0.4

0.6

0.8

1

x

x=0.5L

x=0.5L

(b1) (b2) Figure 10. L=350mm (L/W=3.5): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the critical (1) symmetric and (2) antisymmetric mode shapes (v=0.50).

22

1

1

0.5

0.5

0

0

5 7 9 11 13 15 17 19

-0.5

5 7 9 11 13 15 17 19

-0.5

-1

-1 0

0.2

0.4

0.6

0.8

1

0

0.2

(a1)

0.4

0.6

0.8

1

(a2)

x

x x=0.50L

x=0.50L

x=0.15L

x=0.15L

(b1) (b2) Figure 11. L=550mm (L/W=5.5): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the critical (1) symmetric and (2) antisymmetric mode shapes (v=0.50).

Therefore, five different Type III wrinkling modes of sheets with distinct L/W values are depicted in Figs. 8 to 12: 1.70 (Fig. 8), 2.60 (Fig. 9), 3.50 (Fig. 10), 5.50 (Fig. 11) and 10.00 (Fig. 12). Noting that the sheet wrinkling modes shown in Figs. 6 to 12 are not shell finite element illustrations, but rather tridimensional deformed sheet representations obtained from a one-dimensional modal theory, which is GBT. The observation of Figs. 6 to 12 leads to the following remarks: (i) As it was mentioned previously, the critical wrinkling modes of the sheet with L/W=1.0 are mechanically distinct from all the remaining cases. Indeed, this sheet exhibits four positive equal eigenvalues (quadruplicate modes), represented in Fig. 6(b1)-(b4), contrarily to the sheets with larger aspect ratios that exhibit only two (duplicate modes). On the other hand, the deformed configurations of all four cases are located near the free longitudinal edges (close to s=W/4 and s=3W/4 – see Fig. 1) with negligible deformations in the central region (s=W/2) – once again distinct from the other sheet configurations – compare Fig. 6(b) with Figs. 7-12(b)). In fact, the configuration of the critical modes exhibits 7 or 8 (unequal) half-waves near each longitudinal edge that are independent, i.e., the critical mode 1 and 2 exhibit both 8+8 transverse half-waves but the former is symmetric at mid-span while the latter is antisymmetric (compare Fig. 6(b1)+(b2)). Similarly, the critical mode 3 and 4 exhibit 7+7 transverse half-waves with symmetry (mode 4) or antisymmetry (mode 3) along the mid-span. 23

1

1 5 7 9 11 13 15 17 19

0.5

0.5

0

0

-0.5

-0.5

-1 0

6 8 10 12 14 16 18

0.2

0.4

0.6

0.8

-1

1

0

0.2

(a1)

0.4

0.6

0.8

1

(a2)

x

x x=0.50L

x=0.50L

x=0.095L x=0.095L

(b1) (b2) Figure 12. L=1000mm (L/W=10): (a) GBT modal amplitude functions and (b) bifurcation mode shapes concerning the critical (1) symmetric and (2) antisymmetric mode shapes (v=0.50).

(ii) As the aspect ratio slightly increases from 1.0 to 1.30, the critical wrinkling modes change hugely their shape/configuration (sheets with L/W=1.10, 1.15, 1.20 and 1.30 compared with the previous case L/W=1.00). In fact, the high number of transversal half-waves exhibited by these modes decreases with increasing L/W ratios. The sheets with L/W=1.10, 1.15, 1.20 and 1.30 exhibit either 11 (see Fig. 7(b1)), 9, 7 and 7 half-waves in symmetric modes, respectively, or 12 (see Fig. 7(b2)), 10, 8 and 6 half-waves in antisymmetric modes, respectively. Lastly, the critical wrinkling mode shape of the sheet with L/W=1.05 provides the “transition” between the L/W=1.00 and 1.10 – see a representative 2D mode configuration of the critical symmetric mode in Fig. 5(a). (iii) Although the critical wrinkling modes of sheets with L/W=1.7 (Fig. 8), 2.6 (Fig. 9), 3.5 (Fig. 10), 5.5 (Fig. 11) and 10.0 (Fig. 12) involve the same local deformation modes 7+9+11+13+15 (see Fig. 4(a)), their relative contributions (see/recall Fig. 5(b)) and longitudinal evolution are remarkably distinct – see Figs. 8(a)-12(a). Nonetheless, the number of transversal half-waves is identical in all the sheets with aspect ratios higher than 1.40, either involving 5 (symmetric critical mode) or 6 (antisymmetric critical mode) unequal half-waves. In addition, note that the experimental

24

investigation of polyethylene sheets [6] showed also 5 transverse half-waves characterised with an aspect ratio of 2.5. Now, the longitudinal variation of displacement w(x) at the sheet half-line (s=W/2) for the symmetric modes is shown in Fig. 13 – this function is herein divided to the maximum displacement so that the maximum value is scaled to 1.0 (recall this is a mode). Moreover, the curve of the cosine function 1

1

𝑤 (𝑥) = 2 − 2 cos (

2𝜋𝑥 𝐿

) ,0 ≤ 𝑥 ≤ 𝐿

(30)

was also added to Fig. 13. Note that similar cosine expressions were adopted in some analytical formulations to calculate the critical wrinkling stresses and strains [6,16].

Figure 13. GBT longitudinal profile w(x,s=W/2) for the critical symmetric mode shape of sheets with aspect ratios of 1.0, 1.1, 1.7, 2.6, 3.5, 5.5 and 10.0 (=0.5).

From the observation of Fig. 13 it becomes clear that the longitudinal variation of wrinkling modes is similar to the longitudinal variation provided by Eqs. (30) only if the aspect ratio is close to 1.7 (minimum critical stress). As L/W increases to 2.6, the longitudinal variation of wrinkling modes ceases to exhibit a cosine shape and becomes flattened for 0.4 < 𝑥⁄𝐿 < 0.6. Increasing further the aspect ratio, the longitudinal variation of wrinkling modes is no longer adequately approximated by a cosine function since the maximum amplitude becomes gradually closer to the ends while marginal amplitude appears at the sheet center (particularly for long sheets, e.g., L/W=10.0). Contrarily to the previous cases, the maximum amplitude occurs at the sheet center for values of L/W ratio between 1.7 to 1.0. However, the half-wave width reduces significantly – indeed, for L/W=1.0, the longitudinal profile is concentrated mostly at 0.35 < 𝑥⁄𝐿 < 0.65. Although the expression (30) satisfies the clamped boundary conditions of the sheet, they provide a good approximation of the longitudinal variation of wrinkling modes for a very limited range of aspect ratios. Therefore, the analytical expressions developed in some works [6,16] should be used with caution. 25

Now let us turn our attention to the pre-wrinkling transversal membrane normal stress distributions (𝜎𝑦𝑀 (𝑥, 𝑠)) of the above seven sheets (with v=0.50), shown in Fig. 14(a)-(g).

(a)

(b)

(c)

(d)

(e)

(f)

(g) Figure 14. GBT pre-wrinkling transverse membrane normal stress distribution (𝜎𝑦𝑀 (𝑥, 𝑠)) for sheet with v=0.50 and an aspect ratio of (a) 1.0, (b) 1.1, (c) 1.7, (d) 2.6, (e) 3.5, (f) 5.5 and (g) 10.0 (stresses in MPa).

As expected, the observation of the pre-wrinkling membrane transversal normal stress distributions represented in Figs. 14(a)-(g) shows compressive stresses in the regions where wrinkling occurs. Indeed, 𝑀 the wrinkles of the stretched sheet depends on the magnitude and distribution of 𝜎𝑦𝑦 (𝑥, 𝑠). For L/W=1.1

26

and 1.7, the nearly constant maximum compressive stresses at mid-span region of Figs. 14(b)+(c) coincides with the maximum amplitude of the longitudinal profiles displayed in Figs. 7(a)+8(a). For 𝑀 L/W=2.6, 3.5, 5.5 and 10.0, the peaks of 𝜎𝑦𝑦 in Figs. 14(d)-(g) coincide with the maximum amplitude of

the winkles represented in Figs. 9(a)-12(a). On the other hand, it is confirmed that most part of the longer sheets (e.g., L/W>3.5) are practically under a uniaxial tension (𝜎𝑥𝑥 > 0) since 𝜎𝑦𝑦 0 (except in the regions near the fixed-ends). In addition, note that, all the pre-wrinkling stress distributions in Fig. 14 were obtained by employing the complete set of deformations modes12, i.e., conventional, shear and transverse deformation modes. As it was anticipated in [11], the non-linear variation of the stress field cannot be achieved by resorting only to the axial deformation mode, one shear mode and one transverse extension mode (local modes are not necessary in this case) – indeed, a significant number of shear and transverse extension modes are deemed required, which varies on the aspect ratio of the sheet and makes the derivation of accurate analytical solutions a (more) complex task. Finally, Fig. 15 shows the variation of critical strain cr (and stress cr) with the Poisson’s ratio from 0.3 to 0.5, for three different aspect ratios (L/W=1.0, 1.7 and 5.5). It becomes clear that the increase of Poisson’s ratio always decreases the critical strain, but this drop is distinct for different aspect ratios. For sheets with moderate aspect ratio, L/W=1.7 (green curve), the value of v does not play a major role in the variation of the critical strain. For long and very long sheets, L/W=5.5 (red curve), the Poisson’s ratio is slightly more influential. However, it is for square sheets, L/W=1.0 (blue curve), that the Poisson’s ratio is most felt: the critical strain decreases from about 0.18 to roughly 0.08, less than half value. cr (%)

cr (MPa)

20

375 L/W=1.0 L/W=1.7

16

L/W=5.5

300

12

225

8

150

4

75 0

0

0.30

0.35

0.40 v

0.45

0.50

Figure 15. Fixed-ended stretch-induced thin sheet cr (or cr) vs. v curves for L/W=1.0, 1.7 and 5.5.

12

As it was mentioned in the end of Section 4, different sets of GBT modes may be used in pre-wrinkling and wrinkling analysis, however, usually, the number of GBT modes in the former case is much higher in order to compute an accurate pre-buckling stress field (𝜎𝑥𝑥 , 𝜎𝑦𝑦 , 𝜎𝑥𝑦 ).

27

6. Conclusion This paper presented a GBT-based modal theory to investigate the stretch-induced wrinkling of fixedended rectangular thin sheets. The following aspects deserve to be highlighted: (i)

The results obtained by the modal theory showed an excellent agreement with the ones yielded by shell finite elements.

(ii)

The “signature curve” shape is characterised by a vertical asymptote in close vicinity of L/W=1.0. As L/W increases, the critical tensile stress decreases abruptly until it reaches a minimum, then it increases mildly for moderate L/W and tends to a value, being lengthindependent for longer sheets.

(iii) The wrinkling modes of sheets with identical L/W ratios but distinct Poisson’s ratios are practically indistinguishable. (iv)

The critical wrinkling modes can be obtained by employing only local deformation modes. However, the pre-wrinkling stress fields must be computed with a higher number of modes: axial extension, shear and transverse extension modes.

(v)

The modal analysis showed three mechanically distinct critical wrinkling mode types that are associated with (v1) L/W1.0, characterised by a high number of transversal half-waves, located near the longitudinal edges (v2) 1.0
(vi)

The longitudinal variation of transverse displacement only displays an approximate cosine variation for sheets with L/W closed to 1.7 (minimum critical stress). This hypothesis, adopted in [6,16], is no longer satisfactory for other aspect ratios.

(vii) All the wrinkling modes exhibit duplicate critical modes characterised by symmetric/antisymmetric configurations at mid-span – the exception is the square sheet (quadruplicate critical modes). (viii) Lastly, the increase of Poisson’s ratio always decreases the critical strain, but this drop is distinct for different aspect ratios, being marginal/reasonable/high for moderate/long/short aspect ratios.

Acknowledgements This work was supported by FCT, through IDMEC, under LAETA, project UIDB/50022/2020. 28

Appendix

𝒌𝑒00

4𝐿2𝑒 22𝐿𝑒 𝐿 = 𝑒 420 −3𝐿2 𝑒 [ 13𝐿𝑒

𝑒 𝒌11

4𝐿2𝑒 3𝐿𝑒 1 = 30𝐿 2 𝑒 −𝐿𝑒 [−3𝐿𝑒

𝒌𝑒22

2𝐿2𝑒 3𝐿𝑒 2 = 3 𝐿𝑒 𝐿2𝑒 [−3𝐿𝑒

𝒌𝑒20

−4𝐿2𝑒 −3𝐿𝑒 1 = 30𝐿 𝑒 𝐿2𝑒 [ 3𝐿𝑒

22𝐿𝑒 156 −13𝐿𝑒 54

−3𝐿2𝑒 −13𝐿𝑒 4𝐿2𝑒 −22𝐿𝑒

3𝐿𝑒 −𝐿2𝑒 36 3𝐿𝑒 3𝐿𝑒 4𝐿2𝑒 −36 −3𝐿𝑒 3𝐿𝑒 6 3𝐿𝑒 6

𝐿2𝑒 3𝐿𝑒 2𝐿2𝑒 −3𝐿𝑒

−33𝐿𝑒 −36 −3𝐿𝑒 36

13𝐿𝑒 54 −22𝐿𝑒 156 ]

(A1)

−3𝐿𝑒 −36 −3𝐿𝑒 36 ]

(A2)

−3𝐿𝑒 −6 −3𝐿𝑒 6 ]

𝐿2𝑒 −3𝐿𝑒 −4𝐿2𝑒 3𝐿𝑒

3𝐿𝑒 36 33𝐿𝑒 −36 ]

(A3)

𝒌𝑒02 = (𝒌𝑒20 )𝑇

André D. Martins: Investigation, Software, Formal analysis, Validation, Data curation, Visualization, Writing- Original draft preparation Nuno Silvestre: Supervision, Conceptualization, Methodology, Investigation, WritingReviewing and Editing Rui Bebiano: Software

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Graphical Abstract

29

(A4)

1

1 3D view of sheet wrinkling obtained from the 1D modal theory

0.5

0.5

x=0.50 L 0

x

0

5 7 9 11 13 15 17 19

-0.5

x=0.15 L

5 7 9 11 13 15 17 19

-0.5

-1

-1 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Variation of modes along the sheet length

References [1] Nayyar V, Ravi-Chandar K, Huang R (2011). Stretch-induced stress patterns and wrinkles in hyperelastic thin sheets, International Journal of Solids and Structures, 48(25-26), 15 December 2011, 3471-3483. [2] Wong YW, Pellegrino S (2006). Wrinkled membranes. Part I: experiments, Journal Mechanics of Materials and Structures, 1(1), 1-23. [3] Jenkins CH, Haugen F, Spicher WH (1998). Experimental measurement of wrinkling in membranes undergoing planar deformations, Experimental Mechanics, 38(2), 147-152. [4] Friedl N, Rammerstorfer FG, Fischer FD (2000). Buckling of stretched strips, Computers & Structures, 78(1-3), 185-190. [5] Jacques J, Potier-Ferry M (2005). On mode localisation in tensile plate buckling, Comptes Rendus Mécanique, 333(11), 804-809. [6] Puntel E, Deseri L, Fried E (2011). Wrinkling of a stretched thin sheet, Journal of Elasticity, 105, 137–170. [7] Healey TJ, Li Q, Cheng RB (2013). Wrinkling behavior of highly stretched rectangular elastic films via parametric global bifurcation, Journal of Nonlinear Science, 23, 777–805. [8] Zheng L (2009). Wrinkling of Dielectric Elastomer Membranes, PhD Thesis, California Institute of Technology, Pasadena, California (USA). [9] Nayyar V, Ravi-Chandar K, Huang R (2014). Stretch-induced wrinkling of polyethylene thin sheets: experiments and modeling, International Journal of Solids and Structures, 51(9), 1 May 2014, 1847-1858.

30

[10] Carvalho RV (2015). Wrinkling of Thin Sheets Under Tension, MSc Thesis in Aerospace Engineering, University of Lisbon, Portugal. [11] Silvestre N (2016). Wrinkling of stretched thin sheets: is restrained Poisson’s effects the sole cause?, Engineering Structures, 106, 195-208. [12] Li Q, Healey TJ (2016). Stability boundaries for wrinkling in highly stretched elastic sheets, Journal of the Mechanics and Physics of Solids, 97, 260-274. [13] Sipos AA, Fehér E (2016). Disappearance of stretch-induced wrinkles of thin sheets: A study of orthotropic films, International Journal of Solids and Structures, 97–98, 275-283. [14] Iwasa T (2017). Approximate estimation of wrinkle wavelength and maximum amplitude using a tension-field solution, International Journal of Solids and Structures, 121, 201-211. [15] Iwasa T (2018). Experimental verification on simplified estimation method for envelope curve of wrinkled membrane surface distortions, Thin-Walled Structures, 122, 622-634. [16] Zhu J, Zhang X, Wierzbicki T (2018). Stretch-induced wrinkling of highly orthotropic thin films, International Journal of Solids and Structures, 139–140, 238-249. [17] Wang T, Fu C, Xu F, Huo Y, Potier-Ferry M (2019). On the wrinkling and restabilization of highly stretched sheets, International Journal of Engineering Science, 136, 1-16. [18] Fu C, Wang T, Xu F, Huo Y, Potier-Ferry M (2019). A modeling and resolution framework for wrinkling in hyperelastic sheets at finite membrane strain, Journal of the Mechanics and Physics of Solids, 124, 446-470. [19] Huang Q, Kuang Z, Hu H, Potier-Ferry M (2019). Multiscale analysis of membrane instability by using the Arlequin method, International Journal of Solids and Structures, 162, 60-75. [20] Dadgar-Rad F, Imani A (2019). Theory of gradient-elastic membranes and its application in the wrinkling analysis of stretched thin sheets, Journal of the Mechanics and Physics of Solids, 132, 103679. [21] Liu F, Xu F, Fu C (2019). Orientable wrinkles in stretched orthotropic films, Extreme Mechanics Letters, 33, 100579. [22] Schardt R (1989). Verallgemeinerte Technische Biegetheorie, Berlin, Springer-Verlag. (German) [23] Camotim D, Silvestre N, Gonçalves R, Dinis PB (2004). GBT analysis of thin-walled members: new formulations and applications, Thin-Walled Structures: Recent Advances and Future Trends in Thin-Walled Structures Technology, J Loughlan (ed.), Canopus Publishing Ltd. (Bath), 137-168. [24] Camotim D, Basaglia C, Bebiano R, Gonçalves R, Silvestre N (2010). Latest developments in the GBT analysis of thin-walled steel structures, Proceedings of International Colloquium on Stability and Ductility of Steel Structures (SDSS’Rio 2010 – Rio de Janeiro, 8-10/9), E. Batista, P. Vellasco, L. Lima (eds.), 33-58 (Vol. 1). 31

[25] Camotim D, Basaglia C, Silva NF, Silvestre N (2010). Numerical analysis of thin-walled structures using Generalised Beam Theory (GBT): recent and future developments, Computational Technology Reviews, vol. 1, B. Topping et al. (ed.), Saxe-Coburg (Stirlingshire), 315-354. [26] Silva NMF, Camotim D, Silvestre N (2008). “Generalized beam theory formulation able to capture load application and localized web buckling effects”, Proceedings of 19th International Specialty Conference on Cold-Formed Steel Structures, (CCFSS 2008 – St. Louis, 14-15/10), 33-59. [27] Natário P, Silvestre N, Camotim D (2012). Localised web buckling analysis of beams subjected to concentrated loads using GBT, Thin-Walled Structures, 61(December), 27-41. [28] Basaglia C, Camotim D (2013). Enhanced Generalised Beam Theory buckling formulation to handle transverse load application effects, International Journal of Solids and Structures, 50(3-4), 531-547. [29] Bebiano R, Basaglia C, Camotim D, Gonçalves R (2018). “GBT buckling analysis of generally loaded thin-walled members with arbitrary flat-walled cross-sections”, Thin-Walled Structures, 123(February), 11-24. [30] Simulia Inc. (2009). ABAQUS Standard (version 6.9-3). [31] Vlasov VZ (1959). Thin-Walled Elastic Bars, Fizgmatiz (Moscow). (Russian  English translation: Israel Program for Scientific Translation, Jerusalem, 1961) [32] Bebiano R, Gonçalves R, Camotim D (2015). A cross-section analysis procedure to rationalise and automate the performance of GBT-based structural analyses, Thin-Walled Structures, 92(July), 29-47. [33] Gonçalves R, Bebiano R, Camotim D (2014). On the shear deformation modes in the framework of Generalised Beam Theory, Thin-Walled Structures, 84(November), 325-334.

32