GBT-Based Structural Analysis of Elastic-Plastic Thin-Walled Members

Structural systems made of high-strength and/or high-ductility metals are usually also rather slender, which means that their structural behavior and ultimate strength are often governed by a combination of plasticity and instability effects. Currently, the rigorous numerical analysis of such systems can only be achieved by resorting to complex and computationally costly shell finite element simulations. This work aims at supplying to designers/researchers an efficient and structurally clarifying alternative to assess the geometrically and/or materially non-linear behavior (up to and beyond the ultimate load) of prismatic thin-walled members, such as those built from cold-formed steel. The proposed approach is based on Generalized Beam Theory (GBT) and is suitable for members exhibiting arbitrary deformation patterns (e.g., global, local, distortional, shear) and made of non-linear isotropic materials (e.g., carbon/stainless steel grades or aluminum alloys). The paper begins by providing a critical overview of the physically and geometrically non-linear GBT formulation recently developed and validated by the authors (Abambres et al. 2012a), which is followed by the presentation and thorough discussion of several illustrative numerical results concerning the structural responses of 4 members (beams and columns) made of distinct (linear, bi-linear or highly non-linear) materials. The GBT results consist of equilibrium paths, modal participation diagrams and amplitude functions, stress contours, displacement profiles and collapse mechanisms  some of them are compared with values obtained from ABAQUS shell finite element analyses. It is shown that the GBT modal nature makes it possible (i) to acquire in-depth knowledge on the member behavioral mechanics at any given equilibrium state (elastic or elastic-plastic), as well as (ii) to provide evidence of the GBT computational efficiency, which is achieved by excluding from the analyses all the deformation modes that do not play any role in a particular member structural response.


Introduction
The use of high-strength and/or high-ductility metals, like steel, aluminium or titanium alloys, as well as recent advances in manufacturing technology, has made it possible to build structural systems (often made of thin-walled members) that (i) are highly efficient (large strength-to-weight ratios), thus leading to fairly economic solutions (reduced construction time and considerable material/transportation savings) and (ii) very often achieve a remarkable visual/aesthetic impact. Therefore, it is not surprising that thinwalled structures are widely employed in several areas of Aeronautical, Automotive, Civil, Mechanical, Chemical and Offshore Engineering (Loughlan 2004). Nowadays, the most commonly used thin-walled in the particular behavior under scrutiny, thus further reducing the number of d.o.f. involved in a GBT analysis (i.e., increasing its computational efficiency). GBT has attracted the interest of several researchers worldwide, leading to the development of new formulations and applications. In particular, GBT has been extensively upgraded at the Technical University of Lisbon (Camotim et al. 2010a(Camotim et al. , 2010b, where it has been applied to different (i) types of analysis (first-order, buckling, vibration, postbuckling, dynamic), (ii) boundary and loading conditions (e.g., localized supports and non-uniform internal forces) and (iii) materials (steel, steel-concrete, FRP). With a few exceptions, the material models adopted in all these works were always linear elastic. A physically non-linear GBT formulation was first reported by Gonçalves and Camotim (2004) in the context of elastic-plastic bifurcation analyses  more recently, the same authors (Gonçalves & Camotim 2011 proposed GBT beam finite elements (BFE) based on the J 2 -flow plasticity theory and aimed at performing member first-order and second-order elastic-plastic analyses. In parallel, Abambres et al. (2012aAbambres et al. ( -c, 2013a developed and validated alternative elastic-plastic GBT formulations, also based on the J 2 -flow plasticity theory but differing from the previous ones in the fact that (i) the deformation modes are determined by means of the procedure proposed by Silva et al. (2008) and Silva (2013), and (ii) a novel degree of freedom (warping rotation) is considered. Very recently, GBT was first employed (and effectively validated) to perform post-buckling analyses of thin-walled members made of elastoplastic materials exhibiting strain-hardening, namely stainless steel tubular columns (Abambres et al. 2012c). When compared with carbon steel, stainless steel (and also aluminum) alloys are characterized by (i) the absence of a well-defined yield plateau, and (ii) a pronounced non-linearity beyond the proportional limit, generally associated with the presence of a large amount of strain-hardening. This work begins by providing a critical overview of the main concepts and procedures involved in the recent development and numerical implementation of the aforementioned geometrically and materially (elastic-plastic) non-linear GBT formulation. Then, the paper is essentially devoted to illustrate the application and potential of the developed GBT approach  this is done by presenting and discussing in great detail results concerning the structural response (mostly second-order, including the post-collapse stages) of 4 thin-walled members (beam and columns) that (i) display different cross-section types (open/closed, branched/unbranched), (ii) exhibit distinct material behaviors (linear, bi-linear or highly non-linear) and (iii) experience several deformation patterns (global, local, distortional, transverse and shear). Concerning initial imperfections, only the geometrical ones are incorporated in the non-linear (post-buckling) elastic-plastic analyses  both residual stresses and cross-section corner effects are neglected. The most relevant GBT features, stemming from its unique modal nature, are highlighted  in particular, (i) the evolution of the member deformed configuration is analyzed through GBT modal participation diagrams and amplitude functions, and (ii) equilibrium paths and deformed configurations obtained with different deformation mode sets are compared for each analysis. Moreover, for validation purposes, some GBT-based results are compared with those determined through ABAQUS SFEA.

GBT Non-Linear Elastic-Plastic Formulation -Overview
Only a critical overview of the main concepts and procedures involved in the development and numerical implementation of the proposed GBT formulation is presented in this section. After presenting some aspects related to (i) the kinematics and cross-section analysis adopted and (ii) the establishment of the member equilibrium equations, the development of the non-linear BFE incremental equilibrium equations is addressed (the plasticity model employed is described ahead). For detailed information on these issues, the interested reader is referred to a very recent publication by the authors (Abambres et al. 2012c).

Kinematics and Equilibrium Equations
Consider the local coordinate system (x, s, z) shown in Figure 1(a) at each wall mid-surface of a thinwalled bar, where x, s and z are the longitudinal (0 ≤ x ≤ L  L is the member length), transverse (0 ≤ s ≤ b  b is the wall width) and through-thickness (-t/2 ≤ z ≤ t/2  t is the wall thickness) coordinates, respectively  the corresponding displacements are u (axial or warping), v (transverse) and w (flexural). The GBT analysis of a thin-walled member consists of two main steps, namely (i) a cross-section analysis and (ii) a member analysis. The former leads to the determination of the deformation modes, i.e., their displacement profiles u k (s), v k (s) and w k (s) (k stands for a generic mode) along the cross-section mid-line (z=0), and to the evaluation of the associated modal mechanical properties. As for the member analysis, it consists of determining the deformation mode amplitude functions  k (x), providing the longitudinal variations of the corresponding displacement profiles. Then, the GBT displacement field at the member mid-surface is given by a linear combination of products between each cross-section modal displacement profile and the respective amplitude function, where subscript k, denoting each deformation mode, satisfies Einstein's summation convention. As mentioned earlier, the cross-section analysis adopted follows the approach developed by Silva et al. (2008) and Silva (2013), which considers 4 deformation mode families: conventional, cell shear flow, warping shear and transverse extension  information about this procedure can also be found in Silvestre et al. (2011) and Abambres (2013). The existing GBT cross-section analyses, which consider either three (u, v, w) or four (u, v, w,  x  rotation about the x-axis) d.o.f. per section node, require fine discretizations to capture the highly non-linear variation of the warping displacements, along the cross-section mid-line, due to shear deformation and/or spread of yielding. This fact brings about the convenience of enhancing the representation of these displacements, which is done here by considering additional nodal d.o.f., denoted "warping rotations" and consisting of rotations about the z-axis ( z )  one per intermediate node and more than one per natural node (there is an independent rotation  z per converging wall direction). This novel feature makes it possible to approximate each warping profile u k (s) by means of piecewise (within each wall segment) cubic polynomials, instead of the piecewise linear functions employed in the previous approaches  see the warping shear deformation mode profiles shown in subsection 4.1.
It should be noted that the proposed formulation retains the GBT fundamental plane-stress assumption, which means that all stress (σ xz , σ sz , σ zz ) and strain (γ xz , γ sz , ε zz ) components concerning the through-thickness Figure 1: (a) Local coordinate system at each wall mid-surface, (b) general external distributed load q(x,s) and (c) non-null stress components (plane-stress assumption). direction (z) are deemed null everywhere, regardless of the problem under consideration. Accordingly, and making use of the Principle of Virtual Work, the GBT equilibrium equation reads where (i) σ xx , σ ss , σ xs are the longitudinal normal, transverse normal and shear 2 nd Piola-Kirchhoff stress components (see Fig. 1(c)), (ii) q x , q s , q z are the local components of a general external distributed force applied at the member mid-surface (see Fig. 1(b)), and (iii)  xx ,  ss , γ xs are the longitudinal, transverse and shear virtual Green-Saint-Venant strain components (non-linear functions of displacements u, v, w), expressed in terms of the modal amplitude functions  k (x) and their derivatives, through relations (1). Moreover, the formulation accounts for (i) residual stresses and initial geometrical imperfections, and (ii) arbitrary loadings dependent on a single load parameter , including concentrated forces and/or moments.

Non-Linear Beam Finite Element
The rigorous determination of equilibrium configurations in a non-linear analysis requires the use of an incremental-iterative strategy. The cylindrical arc-length method (Clarke & Hancock 1990, De Borst et al. 2012) is adopted in this work and its implementation involves (i) establishing incremental equilibrium equations, based on the tangent stiffness matrix (K tan ), and (ii) evaluating internal force vectors (f int ) -a path-independent iterative strategy (Powell & Simons 1981) is adopted: the strain increments at any Gauss integration point are evaluated with respect to the last (converged) equilibrium configuration.
The numerical results presented and discussed in section 4 were obtained by adopting Hermite cubic polynomials to approximate the GBT modal amplitude functions in each FE, i.e., where (i) the 1×4 H Ψ vector stores the Hermite cubic polynomials, (ii) d k is the 4×1 displacement vector related to the approximation of mode k, (iii) the first expression concerns only the axial extension and warping shear modes (no in-plane displacements: v=w=0) 3 , and (iv) the second expression applies to all other deformation modes (v≠0 and/or w≠0). After writing all amplitude functions according to their FE approximations, including the non-linear virtual strain components in the first member of Eq. (2) yields where f is the external force vector corresponding to a unit load parameter The internal force vector i th component (4×1 sub-vector associated with deformation mode i) can be obtained by means of where   int i mn F is a 4×1 vector. Note that all displacement vectors (d k in Eq. (3)) and stress components that appear in Eq. (5) concern a generic "equilibrium configuration" (during the iteration procedure and when the structural response is non-linear, this configuration does not satisfy equilibrium for the applied loads under consideration). Once the internal force vector is defined, the incremental equilibrium equation for an arbitrary member deformed configuration (j) can be established as where (i) d is the displacement vector and (ii) the tangent stiffness matrix (i, p) th component, K ip , tan , is a 4x4 sub-matrix concerning deformation modes i and p, reading (Silvestre & Camotim 2003) which corresponds to the Jacobian of the internal force vector i th component, with respect to the mode p displacement vector  its components are d pl (l=1,…,4). The Jacobian columns are defined by the second expression in Eq. (7), where each term is expressed as (recall Eq. (5)) with the stresses (recall Fig. 1(c)), stress gradients and displacement vectors computed at "equilibrium configuration" j. For an arbitrary elastic-plastic material, the stress gradients are given at every point by where the deformation gradients ∂ε/∂d pl are obtained by writing all strain-displacement relations according to the FE approximations in Eq.
(3). The stress components (Eqs. (5) and (8)) and their gradients ∂σ/∂ε (Eq. (9)) are obtained from the J 2 -flow plasticity model, resorting to implicit and/or explicit numerical integration schemes wherever plastic deformation takes place  this issue will be addressed in section 3.
Finally, a few words devoted to the imposition of boundary conditions. In first-order analyses of members with cross-sections exhibiting null primary warping and non-null shear stresses σ xs , it has been concluded that shear locking effects (due to the beam theory inadequacy to describe rigorously static and kinematic measures) are best minimized by ensuring the satisfaction of kinematic boundary conditions:  k,x =0 in every section with null warping. Besides refining the FE mesh near cross-sections with null warping and non-null shear stresses (e.g., a clamped support), it was also noted that imposing (at those sections)  k,xx =0 for the warping shear modes also helps improving the accuracy of the analysis. However, several validation examples provided evidence that, in second-order analyses, the above procedure should be kept, but without imposing  k,xx =0 for the warping shear modes  due to the non-linearity of the straindisplacement relations, this condition is no longer essential to obtain accurate results.

Overview of the Plasticity Model
Since a detailed description of the (i) J 2 -flow theory (with the associated flow rule) and (ii) implicit and explicit numerical integration schemes considered in this work, can be found in the literature (e.g., Wu 2005, De Borst et al. 2012 and also in recent publication by the authors (Abambres et al 2012b,c), only an overview of the main concepts and procedures involved in their development is reported herein.

J 2 -flow Theory
For an isotropic material with strain-hardening, the von Mises yield criterion reads where (i) f and F(σ) are the von Mises yield function and stress, respectively, (ii) σ y is the yield stress (or hardening function), (iii) κ is the hardening parameter and (iv) J 2 is the 2 nd invariant of the deviatoric stress. If a point (σ, κ), located on the yield surface, experiences an infinitesimal plastic flow, Prager´s consistency equation states that where (i) dη ≥ 0 is the plastic proportionality factor (used to define the flow rule dε ij p = dη ∂f / ∂σ ij ) and (ii) h ≥ 0 is the hardening modulus (null for perfectly-plastic models). This modulus plays a crucial role in defining the elastoplastic constitutive matrix and, therefore, also in the development of implicit and explicit integration techniques used to update σ and κ at any material point experiencing plastic flow.
Since von Mises yield criterion is generally recognized as appropriate to model metals, it was decided to derive dκ on the basis of the modified work-hardening relation (De Borst et al. 2012) which leads to dκ=dη for the plasticity model under consideration, and can also be used to define the hardening modulus of any elastoplastic material exhibiting linear elasticity as where (i) E is the material Young's modulus, (ii) the gradients dσ y /dε p and dσ /dε refer to the uniaxial stress-strain curve, (iii) ε 0 y is the initial yield strain and (iv) the total strain ε ≥ ε 0 y must be evaluated where its plastic component (ε p ) equals κ. For materials exhibiting linear strain-hardening, it turns out that (i) h is constant, and given by E sh /(1 -E sh / E) at every material point (E sh is the hardening slope), and (ii) the yield stress reads σ y (ε p =κ) = σ 0 y + h κ. For all other materials, one has ε = σ y /E + ε p and, therefore, the relation ε(ε p ) must be obtained numerically. In order to avoid the use of computationally demanding numerical techniques, like Newton-Raphson's method (N-R), a database was created to provide the (i) yield stress (σ y ), (ii) plastic strain (ε p = εσ y /E ≥ 0) and (iii) hardening modulus (h), for equally spaced (Δε=10 -6 ) total strains such that y 0      u (ε u is the ultimate strain) -values of σ y (ε p ) and h(ε p ) not specified in the database are obtained by linear interpolation.
In order to compute the internal force vector and tangent stiffness matrix mentioned in subsection 2.2, it is mandatory to update, at every material point for each "equilibrium configuration", (i) first the stresses and hardening parameter κ, and (ii) then the stress gradients ∂σ/∂ε. For that purpose, at every Gauss integration point and after each iteration of the arc-length procedure, one must (i) compute the elastic stress variation, due to the imposed strain variation (Δε) w.r.t. the last equilibrium configuration i (stress vector σ i and hardening parameter κ i ), and then (ii) check whether the corresponding final stresses (σ e =σ i + D e Δε) fall outside the yield surface f (σ, κ i )=0 (see Eq. (10)) or not 4 . If they do, numerical integration of the incremental constitutive equations is performed, in order to update the stresses and the hardening parameter κ  several implicit and explicit Euler-type methods were implemented to perform this task (an overview is presented in the next subsection). Once the stresses and the hardening parameter are updated, the stress gradients ∂σ/∂ε can be computed straightforwardly, through either the conventional or the consistent elastoplastic constitutive matrix  the latter was used to obtain the results presented in section 4, since it leads to a faster convergence of the arc-length method.

Implicit and Explicit Integration Schemes
A key step in physically non-linear FE analysis concerns the numerical integration of the elastoplastic constitutive relations, so that the unknown increments of both the stresses and the hardening parameter can be obtained wherever plastic deformation occurs. The integration methods available are usually classified as explicit or implicit 5 . The explicit ones can be applied to a wide range of models, although this should be done with some care, since the solution is not enforced to satisfy the yield criterion to a specified tolerance, unlike in the implicit schemes. According to Sloan et al. (2001), the accuracy and efficiency of the explicit methods can be significantly enhanced by adopting strain sub-incrementation and solution correction aimed at bringing the final stresses and hardening parameter to the yield surface  both these recommendations were implemented in this work.
The Backward Euler and Forward Euler schemes are the most popular implicit and explicit methods used in elastoplastic problems. They belong to a family of methods aimed at computing, for each material point, the solution at configuration t + Δt, on the basis of the solution at configuration t + αΔt (0 ≤α ≤1), assuming that (i) t corresponds to the onset of yielding (initial or subsequent) and (ii) Δt reflects the imposition of the elastoplastic strain increment Δε ep (Δε=Δε e + Δε ep ) 6 -the Backward and Forward Euler methods correspond to α=1 and α=0, respectively. However, there are other schemes proposed in the literature that adopt intermediate α values  e.g., α=1/2 has been used to define the method known as Mean Normal. Although the above three methods were implemented, the Backward Euler is addressed here in more detail, since it was the one employed to obtain the results discussed in the next section. Taking (σ j , κ j ) as the solution at a generic material point at configuration t, the Backward Euler method, formulated on the basis of the conventional elastoplastic constitutive matrix, states that the final state (σ j+1 , κ j+1 ) is such that where (i) Δε is the elastoplastic strain increment, (ii) D e is the GBT elastic constitutive matrix and (iii) η j+1 and n j+1 are the plastic proportionality factor and gradient vector evaluated at the (yet unknown) final 4 If f (σ e , κ i ) ≤ 0, the strain increment occurs in the elastic range and all variables are updated under that condition (κ unchanged). If f (σ e , κ i ) > 0, plastic flow occurs and, thus, the updated final stresses (not σ e ) are located on the (updated) yield surface. 5 Explicit integration is based on known variables and, strictly speaking, requires no iteration to predict the unknown increments. Implicit integration, on the other hand, is based on unknown variables and, therefore, involves the iterative solution of a system of non-linear equations at each Gauss point. 6 It was implemented a closed-form expression to determine the initial elastic strain variation Δε e , required to drive a stress state located inside the yield surface to reach it (or to move a stress state between two yield surface points). state, respectively. Since the unknowns in Eq. (14) outnumber the equations by one unit (η j+1 is the extra unknown), an additional equation must ensure that the final stress point σ j+1 is located on the yield surface, i.e., f (σ j+1 , η j+1 )=0 (recall that η j+1 =κ j+1 ). Since this constitutes a non-linear problem, the N-R method is chosen to solve the set of equations where (i) the trial solution (σ 1 j+1 , η 1 j+1 ) for the first iteration is obtained by means of the Forward Euler method, implemented with three strain sub-increments and correction of the final solution (Abambres et al. 2012c), and (ii) the convergence criterion adopted to end the iterative process reads where (σ k+1 j+1 , η k+1 j+1 ) is the solution at the end of iteration k.

Illustrative Examples
The application of the GBT formulation for elastoplastic post-buckling analysis, which was implemented in MATLAB R2010b (MathWorks 2010), is now illustrated. The numerical results presented and discussed in subsections 4.1-4.5 concern four analyses involving (i) LiteSteel (hollow-flange channel) and I-section beams, and (ii) lipped channel and equal-leg angle columns. Note that not all GBT results are shown  only those deemed more significant to (i) highlight the advantages of the GBT analyses and (ii) provide validation against SFEA values from ABAQUS (DS Simulia 2008). The material models considered exhibit a linear elastic regime and aim at simulating the constitutive behaviors of (i) carbon steel (perfectly-plastic and bi-linear models  no hardening and linear hardening) and (ii) austenitic, ferritic and duplex stainless steel alloys (three-stage stress-strain relation due to Quach et al. 2008). Regarding the latter model, it is defined by the three basic Ramberg-Osgood parameters (E, σ 0.2 and n) and was shown to (i) provide accurate predictions over a wide range of tensile and compressive strains, and (ii) outperform the full-range model proposed by Rasmussen (2003)  detailed information about its implementation can be found in Abambres et al. (2012c). Table 1 shows the material properties considered in each analysis, where (i) E, ν, σ 0 y are the Young's modulus, Poisson's ratio and initial yield stress, (ii) E sh is the strain-hardening slope (bi-linear model) and (iii) σ 0.2 , n are the 0.2% proof stress and hardening power (non-linear model 7 ). In the ABAQUS models, which are also based on J 2 -flow theory, the compressive stress-strain input data must be expressed in terms of the true stress (σ t ) and true plastic extension (ε t(p) ), obtained from the nominal stress (σ n ) and extension (ε n ) absolute values as 7 In this model, (i) the 0.01% proof stress (σ 0.01 ) is taken as the initial yield stress (proportional limit) and (ii) superscripts c and t denote "compressive" and "tensile" (see Table 1  angle column row). Since only one hardening function σ y (κ) can be used with von Mises yield criterion, the compressive (instead of the tensile) behavior was modeled -indeed, most structural engineering problems involve members under compression (those prone to instability phenomena).
Concerning the FE meshes adopted, Table 2 shows the number of (i) cross-section wall segments (see 4.1) and BFEs used in GBT, (ii) longitudinal (along x) and transverse (along s) SFEs used in ABAQUS, and (iii) d.o.f. involved in the BFE and SFE models (to assess the efficiency of the GBT analyses). Lastly, Table 3 indicates the number of Gauss integration points used per BFE wall segment direction (s, z, x in Fig. 1(a)) in each GBT analysis -the SFEA involved fairly uniform fine meshes of 4-node isoparametric SFEs with full integration (S4 elements  ABAQUS nomenclature), through 2 integration points per surface direction and the same number as in the GBT analysis for the through-thickness direction (z-axis).   The beams and columns analyzed are depicted in Figs. 2 to 5  after providing information about their cross-section deformation modes, in 4.1, the results of each analysis are addressed individually in 4.2-4.5.

Angle column
The first analysis concerns a fixed-pinned 8 LiteSteel beam (LSB) acted by vertical forces F=1000 λ N applied at the x=L/4 and x=3L/4 cross-section web mid-points. It is worth noting that this cross-section is aimed primarily for flexural members and its unique geometry leads to a high structural efficiency, when compared with equivalent hot-rolled profiles (Anapayan et al. 2011, Anapayan & Mahendran 2012. A first-order analysis is performed for a beam made of a perfectly-plastic material, with length L=4000 mm and the cross-section centreline dimensions given in Fig. 2(b). The second problem deals with a fixedpinned I-section beam acted by a vertical patch load (p=2λ/3 N/mm 2  resultant F=7125 λ N) applied along the top flange and centred in the member length  see Fig. 3(a). A post-buckling (non-linear) analysis is performed for a beam exhibiting a bi-linear material behavior, with length L=2850 mm and the cross-section centreline dimensions given in Fig. 3(b). The beam contains critical-mode (lateral-torsional  λ cr =3.10) initial geometrical imperfection with a 3.02 mm amplitude (maximum in-plane displacement).  The first column analyzed is a fixed-ended lipped channel under a compressive end force (F=1000 λ N) uniformly distributed along the section contour (q=2.78 λ N/mm). A post-buckling analysis is performed for a column made of a perfectly-plastic material, with length L=2200 mm and the crosssection centreline dimensions given in Fig. 4(b). The column exhibits a critical-mode (predominantly distortional with three half-waves, even if a small contribution from the first symmetric GBT local mode is perceptible  λ cr =125.14) initial geometrical imperfection with a 2.06 mm amplitude (b f /50  b f is the flange width). Lastly, the final set of results concerns the fixed-ended equal-leg angle column depicted in Fig. 5(a), which (i) is made of a ferritic stainless steel alloy (material properties, reported by Becque & Rasmussen (2009), given in Table 1 and adopted stress-strain curve shown in Figs. 6(a)-(b)), (ii) is acted by a compressive end force (F=10 5 λ N) uniformly distributed along the section contour (q=333.33 λ N/mm), and (iii) has length L=1500 mm and the cross-section centreline dimensions given in Fig. 5(b). In order to assess the influence of the material and/or geometrically non-linear effects on the column response, three analyses are performed: first-order elastoplastic, second-order elastic and second-order elastoplastic. In the second and third analyses the column incorporates a critical-mode (flexural-torsional, but predominantly torsional  λ cr =11.14) initial geometrical imperfection with a 5.40 mm amplitude.
In each analysis, the results presented and discussed (all displacements/stresses concern the member midsurface) correspond to three particular equilibrium states (points on the equilibrium path), denoted as BP (ascending branch  before peak), P (peak or beginning of plateau  ultimate strength) and AP (descending branch after the peak or well within the plateau)  the GBT and ABAQUS values are selected   to ensure a comparison between equilibrium states as close as possible (similar displacements). Finally, two last words to mention that (i) no residual stresses are included in the analyses and (ii) the GBT modal participation diagrams are obtained as explained in Abambres et. al. (2013b) 9 .

GBT Deformation Modes
Information about the GBT deformation modes included in each analysis is now provided. First, Table 4 provides the total numbers of modes per type: global (G), distortional (D), local (L), warping shear (WS) and transverse extension (TE). Since such modes are always a fraction of the full set obtained from the cross-section analysis, they were renumbered sequentially to make the presentation of the results clearer. The in-plane and out-of-plane (warping) configurations of the most relevant GBT modes in each analysis are displayed in Figs. 7-10  the cross-section discretisation adopted is shown together with the transverse extension mode shapes 10 . Regarding the mode selection, a few remarks are appropriate: (i) On the basis of the expected primary member deformed configuration, which is heavily dependent on the characteristics of the loading and initial geometrical imperfections, it is fairly straightforward matter to select which global (mostly), distortional and local deformation modes to include in the GBT analysis. For instance, it is clear that the LSB analysis must necessarily include modes 1 (majoraxis bending), 2 (torsion) and 3 (web lateral distortion)  see Fig. 7. Similarly, (i 1 ) the I-section beam analysis will certainly involve modes 2 (major-axis bending), 3 (minor-axis bending), 4 (torsion) and 8+9 (flange-driven local modes) (see Fig. 8), (i 2 ) the lipped channel column analysis must include modes 1 (axial extension), 2 (minor-axis bending 11 ), 3 (symmetric distortion) and 4+5 (web-driven symmetric local modes) (see Fig. 9), and (i 3 ) the angle column analysis must involve modes 1 (axial extension), 2 (major-axis bending), 3 (minor-axis bending 12 ) and 4 (torsion) (see Fig. 10). (ii) The selection of the warping shear and transverse extension (and also some local) modes, mostly intended to capture "higher-order effects" stemming from pronounced geometrical and/or material non-linearity, is much less obvious. It must be made by a combination of (ii 1 ) common sense, to exclude deformation modes exhibiting excessively high strain energies (e.g., local modes having multiple curvature within a single wall), (ii 2 ) symmetry considerations, which must comply with the member loading, and (ii 3 ) past experience from the user -e.g., lower-order warping shear modes arising from global bending are easily identified. For instance, (ii 1 ) the LSB analysis only needs to include modes anti-symmetric w.r.t. the centroidal horizontal axis, (ii 2 ) the I-section beam analysis must include almost all the modes, as no clear selection criterion can be found a priori (only the higher-order local modes can be readily excluded), (ii 3 ) the lipped channel column analysis must To account for the horizontal effective centroid shift stemming from the normal stress redistribution associated with the distortional and local deformations (e.g., Young & Rasmussen 1999). 12 Again to account for the effective centroid shift stemming from (i) the normal stress redistribution associated with torsion (Dinis et al. 2012) and also from (ii) the spread of plasticity. include at least all the symmetric (w.r.t. the centroidal horizontal axis) warping shear and transverse extension modes, essential to capture the "higher-order effects" associated with the distortional deformations, and (ii 4 ) the angle column analysis involves only a few lower-order local modes, but all transverse extension modes and most of the shear modes symmetric w.r.t. the cross-section major axis. (iii) Of course, it is always possible to adopt a "safe approach" to the selection of the deformation modes, which consists of (iii 1 ) including "all but the clearly irrelevant modes" in the analysis and, then, (iii 2 ) remove from subsequent similar analyses all those modes found to have no (or minute) contribution to the member deformed configuration evolution.

LiteSteel Beam (LSB)
The equilibrium paths λ(δ z ) shown in Fig. 11(a), where δ z is the vertical displacement of the mid-span top flange node shown in Fig. 2(b), were obtained by means of an ABAQUS SFEA and three GBT analyses, which include either (i) the initially selected 43 modes indicated in Table 4, corresponding to 6170 d.o.f., (ii) the most relevant 17 modes (2 global, 2 distortional, 3 local, 5 warping shear and 5 transverse extension), including all those shown in Fig. 7 and involving 2436 d.o.f., or (iii) 12 modes (transverse extension modes removed from the previous set), corresponding to 1721 d.o.f. 13 . As for Fig. 11(b), it shows the GBT modal participation diagram concerning the evolution (as δ z increases) of the beam most deformed cross-section configuration, located at the vicinity of x=3000 mm 14 . The observation of the results presented in Figs (ii) The modal participation diagram shows that the beam behavior is governed by modes 1 (major axis bending  38.4-47.4%), 2 (torsion  15.9-38.9%) and 3 (web distortion  19.1-31.1%). The remaining modes, among which mode 5 (web local) is the most visible (its participation grows up to 2.5%), only emerge beyond the elastic regime. (iii) When only 17 deformation modes are included in the GBT analysis, which corresponds roughly to a 60% d.o.f. reduction, the results obtained are still reasonably good -the plastic load at the AP state only overestimates its ABAQUS counterpart by 4.3%. (iv) However, the further removal of the five transverse extension modes considerably reduces the accuracy of the results  e.g., the ABAQUS λ value at the AP state is overestimated by 14.6%. This provides evidence of the key role played by the transverse extension modes, even if their joint modal participation is fairly small  see Fig. 11(b).

Figs. 12(a)-(c)
show the evolution of the beam deformed configuration, expressed in terms of the modal amplitude functions ζ (x) concerning the 4 deformation modes appearing in the modal participation diagram of Fig. 11(b), namely modes 1 (major-axis bending), 2 (torsion), 3 (web lateral distortion) and 5 (web local) (see Fig. 7)  these amplitude functions are evaluated at the BP, P and AP equilibrium states. It can be observed that: (i) At equilibrium states BP and P only participations of modes 1, 2 and 3 are visible. The participation of mode 5 is only perceptible at the AP equilibrium state and its impact is mostly restricted to the vicinity of the x=3000 mm cross-section. Figure 13: LSB collapse mechanism (AP state) representations obtained from the ABAQUS and GBT analyses.
(ii) As the load increases, the maximum modal in-plane displacement contributions, proportional to the transverse extensions ε ss , approach the cross-section located at x=3000 mm, where the force closer to the pinned support acts. On the other hand, the first and second derivatives of the above amplitude functions, respectively proportional to the shear strains γ xs and longitudinal extensions ε xx , exhibit the higher values near the cross-sections located at x=0 (fixed support) and, again, x=3000 mm  the locations of the plastic hinges leading to the collapse mechanism shown in Fig. 13. It is worth noting the remarkable similarity between the GBT and ABAQUS collapse mechanism representations  both show clearly the web yield-line and distortion at the x=3000 mm cross-section.
Finally, Figs. 14(a)-(b) concern the P equilibrium state and show the top flange node (see Fig. 2(b)) axial and transverse displacement profiles  all GBT profiles correspond to λ GBT =104.1. It is worth noting that: (i) The three GBT profiles follow the same trend and, naturally, the displacements approach the SFEA solution as the number of deformation modes included in the analysis increases. When the 43 selected modes are included, there is an excellent agreement between the GBT and ABAQUS results. (ii) The comparison between the displacement profiles provided by the GBT analyses including 17 and 12 modes makes it possible to assess the relevant role played by the transverse extension modes.
(a) (b) Figure 14: Top flange node (see Fig. 2(b)) (a) axial (δ x ) and (b) transverse (δ z ) displacement profiles at the P state (λ GBT = 104.1). Fig. 15(a) shows the GBT and ABAQUS equilibrium paths λ(|δ y |), where |δ y | is the absolute value of the mid-span top flange node lateral displacement, as indicated in Fig. 3(b). As for Fig. 15(b), it shows the GBT modal participation diagram concerning the evolution, as | δ y | increases, of the beam most deformed GBT ABAQUS cross-section configuration, located in the mid-span region. Finally, Fig. 16 depicts the beam collapse mechanism (AP equilibrium state). The observation of these results leads to the following comments: (i) The GBT and ABAQUS SFEA equilibrium paths are in very good agreement -the maximum difference is 4.8% and occurs well into the descending branch (AP equilibrium state). (ii) Global deformation clearly governs the beam behavior, namely through contributions from modes 2 (major-axis bending -15.4-60.7%), 3 (minor-axis bending -12.0-28.7%) and 4 (torsion -21.6-44.9%). Concerning the transverse extension modes 99 (web) and 100 (top flange), their combined participations rise gradually along the equilibrium path, reaching 13% in the descending branch  excluding these two deformation modes from the analysis would considerably "stiffen" the beam global behavior, by severely restraining the significant wall in-plane motions clearly observed in Fig. 16. Finally, modes 8 and 9, associated with the top/loaded flange transverse bending due to the patch loading (see Fig. 8   at the pinned support -the warping contributions from modes 3 and 4 reinforce (counteract) each other at this support (see Fig. 8). (v) Since the pinned support only allows for axial displacements, the top flange in-plane rotation leads to high tensile transverse stresses, which may "artificially" spread along the beam length if the BFE mesh is too coarse near that support. Although the BFE discretisation adopted, which is finer at the fixed support and patch load locations (see Table 2), leads to an accurate prediction of the beam displacement field evolution, similarly accurate stress results require a further BFE mesh refinement.

I-Section Beam
Concerning the vertical displacement of the mid-span top flange node (δ z see Fig. 3(b)), Fig. 17 shows equilibrium paths provided by the ABAQUS SFEA and three GBT analyses, which adopt either (i) the initially selected 137 modes (see Table 4 The "switch" (increase/decrease) between the participations of modes 4 (torsion) and 2 (majoraxis bending) that occurs prior to the peak load (see Fig. 15(b)) is responsible for the observed stiffness increase and subsequent (localised) displacement reversal (loop). (iii) The relevance of the transverse extension modes can be attested by the clearly improved (although still erroneous) GBT non-linear equilibrium path determined with 50 deformation modes, when compared with that obtained with 43 modes -note the "incipient" loop formation.
Finally, Fig. 18 shows the top flange node (see Fig. 3(b)) axial and vertical displacement profiles at the (i) BP, P ( Fig. 18(a)) and (ii) AP ( Fig. 18(b)) equilibrium states (the GBT P profiles concern λ GBT =2.40). The following conclusions can be drawn: (i) There is an excellent correlation between all ABAQUS and GBT (137 modes) results.

GBT ABAQUS
(ii) The largest vertical displacements occur at the mid-span region and are associated with abruptly high gradients, leading to the formation of a yield-line at the top flange (see Fig. 16). (iii) When evolving between equilibrium states BP and AP, the transverse displacements in the member left quarter-span change from downwards to upwards, stemming from the prevalence of the torsional deformations over their major-axis bending counterparts. (iv) Naturally, increasing the number of deformation modes included in the GBT analysis improves the results obtained. Note also that, in spite of the significant quantitative differences between the displacement profiles provided by the "exact" (137 modes) and approximate (43 or 50 modes) GBT analyses, they are qualitatively very similar. associated with localized lip deformations. Fig. 19(b) shows the GBT modal participation diagram obtained from the 79 mode GBT analysis and concerning the evolution, as Δ increases, of the beam midspan cross-section deformed configuration. The analysis of these results prompts the following remarks: (i) The GBT equilibrium path obtained with 128 modes correlates quite well with its ABAQUS SFEA counterpart. For instance, the corresponding ultimate loads only differ by 3.7% (λ u , GBT =152.37 and λ u ,ABQ=147.00) and the {Δ, λ} values concerning the BP, P and AP equilibrium states read ( GBT and {4.24,112.00} ABQ , respectively. (ii) The GBT equilibrium paths obtained with 128 and 79 modes virtually coincide over practically the whole deformation range  slight differences only occur well into the descending branch. This means that the 49 additional deformation modes contribute next to nothing to the column response and, thus, can be removed without sacrificing the accuracy of the analysis. (iii) Fig. 19(b) shows clearly that modes 1 (axial extension -3.8-29.4%), 2 (minor-axis bending -2.3-18.9%), 3 (symmetric distortional -47.9-74.7%), 4 (local -3.1-13.8%) and 5 (local -0.7-5.8%) govern the column behavior in the whole deformation range. The dominance of the distortional deformation stems from the fact that the critical buckling load is distortional and considerably lower than the first local buckling counterpart (λ cr,dist =125.14 vs λ b,loc =158.61). (iv) The participation of minor-axis flexure (mode 2) increases visibly a bit before the peak load and becomes even more noticeable in the descending branch. This is due to the load eccentricity arising from the effective centroid shift towards the web caused the longitudinal normal stress redistribution associated with the occurrence of distortional and/or local deformations. Moreover, yielding at the flange-lip corner and lip regions also contribute to the load eccentricity (the onset of yielding at the member mid-surface corresponds approximately to the dashed line in Fig. 19(b)). (v) Beyond the onset of yielding (dashed line), increasing contributions from modes 2 (minor-axis flexure) and 4-9 (local) gradually replace the participation of mode 3 (symmetric distortional) -the column collapse mechanism, shown in Fig. 20, clearly shows web local deformations at mid-span. show the evolution of the column deformed configuration, expressed by the amplitude functions of modes 2, 3 and 4-5, evaluated at the BP, P and AP equilibrium states. It can be observed that: (i) Since the local/distortional deformations and severe yielding erode the mid-span region bending stiffness, forming a kind of "plastic hinge", localized bending emerges there in the post-peak stages. (ii) The local half-wave number grows as the column evolves from the BP to the AP equilibrium statesafter failure, this is mostly due to increasing compressive forces developing at the mid-span region. Note also that, close to peak load ( Fig. 21(b)), the participation of the local mode 4 is equally relevant at the column mid-span and end regions  localization at mid-span only occurs afterwards.

Lipped Channel Column
Lastly, Figs. 22(a)-(b) and 23(a)-(b) provide vertical, lateral and axial displacement profiles of the top flange node (Fig. 4(b)) at the BP, P and AP equilibrium states  the GBT profiles correspond to the load parameter values indicated in the first item of this subsection. The following remarks are appropriate: (i) The GBT displacement profiles obtained with 128 and 79 modes are practically coincident with one minor exception: the axial displacements at the AP equilibrium state ( Fig. 23 (b)), which may differ up to 2%. Thus, the GBT analysis with 79 modes may be viewed as "exact" for this particular case. (ii) The GBT and ABAQUS SFEA axial displacement profiles compare quite well, slightly better than the remaining ones. In the latter case, the column response predicted by GBT is visibly stiffer for the P and AP states, particularly in the mid-span (most deformed) region -a BFE mesh refinement and/or the inclusion of additional transverse extension modes should eliminate this discrepancy 15 . (iii) The vertical displacement profiles indicate that, during collapse (between P and AP), there is a decrease (increase) of the positive (negative) displacements and respective curvatures. This feature stems from the simultaneous (iii 1 ) formation of the yield-line mechanism depicted in Fig. 20 and (iii 2 ) elastic unloading at the positive amplitude regions  the axial displacement profiles in Fig. 23(b) are also affected by the formation of the yield-line mechanism, through the considerable axial stiffness loss at mid-span and slope reduction at remaining (stiffer) regions.

Equal-Leg Angle Column
The GBT and ABAQUS equilibrium paths λ(Δ) shown in Fig. 24, where Δ is the column end-shortening, were obtained by means of three types of analyses: materially non-linear (MNA), geometrically nonlinear imperfect (GNIA), and geometrically and materially non-linear imperfect (GMNIA) analyses.  Fig. 10. Fig. 25(a) compares the thee GMNIA equilibrium paths and Fig. 25(b) provides the GBT modal participation diagram (obtained with 72 modes) concerning the evolution of the column most deformed cross-section configuration (located near mid-span). Lastly, Fig. 26 displays the column collapse mechanisms yielded by the GBT and ABAQUS GMNIA. The observation of the results presented in Figs. 24-26 prompts the following remarks: (i) The GBT MNA, GNIA and GMNIA curves compare very well with their ABAQUS counterpartsindeed, the GMNIA differences never exceed 1.5%. As expected, the structural response yielded by the MNA resembles the material compressive stress-strain curve shown in Fig. 6(b). (ii) The GNIA results evidence a significant post-critical strength, which is a bit surprising in view of the flexural-torsional nature of the critical buckling mode (λ cr =11.14) and initial geometrical imperfection  this peculiar behavioral feature, recently unveiled by Dinis et al. (2012), stems from (ii 1 ) the mechanical similarity between the torsional and local deformations (for this particular cross-section) and (ii 2 ) the fact that the end section (secondary) warping is prevented. (iii) After the separation between the GMNIA and GNIA equilibrium paths (the onset of yielding occurs for λ y ≈6.15), the spread of plasticity erodes the column strength quite rapidly (λ u =6.79). (iv) The {Δ, λ} values provided by the GBT (72 modes) and ABAQUS GMNIA for the BP, P and AP equilibrium states read (iv 1 ) {0.92, 3.74} GBT and {0.93, 3.79} ABQ , (iv 2 ) {2.14, 6.88} GBT (λ u,GBT =6.89) and {2.11, 6.79} ABQ , and (iv 3 ) {14.18, 6.08} GBT and {14.00, 6.08} ABQ , respectively. (v) The modal participation diagram shows that the angle column behavior is governed by contributions from modes 1 (axial extension  4.0-19.9%), 3 (minor-axis bending  1.8-17.7%) and 4 (torsion  55.8-81.1%)  the presence of minor-axis bending is due to a combination of effective centroid shift and interaction effects (Dinis et al. 2012). As for the participation of mode 2 (major-axis bending), which appears in the column critical buckling mode (and, therefore, also in the initial geometrical imperfection), it remains fairly small in the whole deformation range. As explained earlier (see 4.3), large in-plane rigid-body motions, such as those exhibited by the column collapse mechanism central region (see Fig. 26), always entail relevant transverse extensions, reflected by the contributions from modes 45 (0.6-7.4%) and 46 (0.6-7.5%) -see Fig. 10.  (vi) During collapse, the dominant participation of torsion (mode 4) is partially and gradually replaced by increasing contributions from minor-axis bending and transverse extension (modes 3 and 45-46). (vii) The GBT GMNIA with 47 deformation modes yields an equilibrium path that practically coincides with the ABAQUS one up to the peak load, which is overestimated by about 3.9%  this overestimation remains fairly throughout the whole descending branch (the maximum difference is 5.7%).

GBT ABAQUS
Figs. 27(a)-(b) show the GMNIA axial and lateral displacement profiles concerning the leg free node indicated in Fig. 5(b) and the BP, P, AP equilibrium states  the GBT P profiles concern λ GBT =6.88. The observation of these displacement profiles makes it possible to conclude that: (i) All the 72 mode GBT and ABAQUS profiles compare very well, with one sole exception: the AP axial displacements, even if they are qualitatively similar and show good agreement at the end and mid-span regions. It is still worth mentioning that the GBT analysis reveals important contributions from the warping shear modes 20-21 (associated with the major/minor-axis bending modes 2-3). (ii) The lateral displacement profiles exhibit an inner half-wave and two outer quarter-waves (to comply with the fixed boundary conditions and the imposed geometrical imperfection). Unlike in the BP and P equilibrium states, where the inner half-wave exhibits a fairly sinusoidal shape, the AP equilibrium state is characterized by a "centrally flattened" inner half-wave  this is due to the development of a torsional collapse mechanism occurring in the column central region, which is triggered by the formation of yield-lines roughly at the 1/3 and 2/3-span cross-sections (see Fig. 26). (iii) Although the λ(Δ) equilibrium path provided by the 47 mode GBT GMNIA agrees reasonably well with the ABAQUS one up to the peak load (see Fig. 25(a)), Fig. 27(a) shows that the corresponding P displacement profiles exhibit considerable discrepancies, even if the correct trends are captured. It was found that these discrepancies are due to the removal of several warping shear modes  in spite of the global nature of the column deformation, these modes play a relevant role. Finally, Fig. 28 shows the Mises stress distributions obtained with the 72 mode GBT and ABAQUS GMNIA at the P equilibrium state. Besides the remarkable similarity exhibited by the two stress contours, it is also interesting to observe that the highest stresses, already above the initial yield stress (220 N/mm 2 ) at this stage, occur in those regions that were found to exhibit lower axial stiffness (see Fig. 27). Figure 28: Angle column Mises stress (N/mm 2 ) contour at the P equilibrium state.

Concluding Remarks
This work provided an overview of the main concepts and procedures involved in the development of an elastic-plastic geometrically non-linear GBT formulation, which is based on the J 2 -flow plasticity theory and applicable to isotropic materials exhibiting arbitrary strain-hardening. Moreover, the paper also includes a vast illustration of the application and potential of the developed GBT formulation: numerical results concerning the structural response (up to and beyond the ultimate load) of 4 thin-walled members (i) made of distinct materials (linear, bi-linear or highly non-linear) and (ii) exhibiting several deformation patterns (global, local, distortional, transverse and shear) were presented and discussed in detail. These results were obtained by means of (i) a LiteSteel beam first-order elastoplastic analysis, (ii) an I-section beam second-order (post-buckling) elastoplastic analysis, (iii) a lipped channel column second-order elastoplastic analysis and (iv) equal-leg angle column first-order elastoplastic, second-order elastic and second-order elastoplastic analyses. The GBT unique modal features were highlighted (i) by analyzing the evolution of the member deformed configurations through GBT modal participation diagrams and amplitude functions, and (ii) by comparing equilibrium paths and deformed configurations obtained with different deformation mode sets. It was shown that the GBT modal analyses make it possible to acquire in-depth knowledge on the member behavioral mechanics in the elastic and elastic-plastic regimes. For validation and numerical efficiency assessment purposes, most GBT results (equilibrium paths, collapse mechanisms, displacements profiles and stress contours) were compared with values yielded by shell finite element analyses performed in the code ABAQUS. It was shown that similarly accurate results are provided by GBT analyses including judiciously selected deformation mode sets that involve only 25% (on average) of the number of d.o.f. required by the corresponding ABAQUS analyses. Moreover, it was also shown that GBT analyses including fairly small deformation mode sets ( numbers) are able to provide an excellent qualitative assessment of the member structural response  they are ideally suited for the performance of preliminary/exploratory analyses. Finally, a few closing words to mention that this work demonstrated that, in the context of the assessment of the structural behavior, ultimate strength and collapse of prismatic thin-walled members, GBT constitutes a very promising alternative to shell finite element analysis or other traditional methods (e.g., yield-line analysis).