GBT-Based Elastic-Plastic Post-Buckling Analysis of Stainless Steel Thin-Walled Members

When compared with carbon steel, stainless steel exhibits a more pronounced non-linearity and no well-defined yield plateau, as well as appealing features such as aesthetics, higher corrosion resistance and lower life cycle cost. Due to its considerably high ductility/strength and cost, stainless steel structural solutions tend to be adopted mostly for slender/light structures, thus rendering the assessment of their structural behaviour rather complex, chiefly because of the high susceptibility to instability phenomena. The first objective of this paper is to present the main concepts and procedures involved in the development of a geometrically and physically non-linear Generalised Beam Theory (GBT) formulation and numerical implementation (code), intended to analyse the behaviour and collapse of thin-walled members made of materials with a highly non-linear stress-strain curve (e.g., stainless steel or aluminium). The second objective is to validate and illustrate the application of the proposed GBT formulation, by comparing its results (equilibrium paths, ultimate loads, deformed configurations, displacement profiles and stress distributions) with those provided by shell finite element analyses of two lean duplex square hollow section (SHS) columns previously investigated, both experimentally and numerically, by Theofanous and Gardner [1]. The stainless steel material behaviour is modelled as non-linear isotropic and the GBT analysis includes initial geometrical imperfections, but neglects corner strength enhancements and membrane residual stresses. It is shown that the GBT unique modal nature makes it possible to acquire in-depth knowledge concerning the mechanics of the column behaviour, by providing "structural x-rays" of the (elastic or elastic-plastic) equilibrium configurations: modal participation diagrams showing the quantitative contributions of the global, local, warping shear and transverse extension deformation modes  moreover, this feature makes it possible to exclude, from future similar GBT analyses, those deformation modes found to play a negligible role in the mechanics of the behaviour under scrutiny, thus further reducing the number of degrees of freedom involved in a GBT analysis, i.e., increasing its computational efficiency.


Introduction
Stainless steel has been used in the construction industry for over 70 years, even if its wide dissemination has been severely restricted by fairly large productions costs (e.g., much larger than for carbon steel). However, recent developments in material technology [2] are changing this situation quite rapidly, thus (i) making stainless steel nowadays one of the world's most profitably recycled material [3] and (ii) leading to a renewed interest on stainless steel structural members and systems -since the year 2000, there has been an increasing number of significant structural applications of stainless steel [4,5].
Stainless steel types are classified according to their main alloy constituents and the austenitic and duplex (or austenitic-ferritic) are, by far, the most frequently used alloys in building and construction.
Nonetheless, in spite of the several attractive features of stainless steel, when compared with carbon steel, such as better appearance, higher corrosion resistance and more cost-effective and longer life cycle, an increased use of stainless steel in common applications, such as office or residential buildings, requires (i) the development of efficient and safe design rules that can fully exploit the stainless steel structural potential and (ii) the dissemination of easy-to-use tools to perform the structural design. Concerning the first aspect, Eurocode 3, part 1.4 (EN 1993-1-4 [6]) was published in 2006 as a full European standard prescribing supplementary rules for the design of stainless steel structures. However, it is well known that some of its rules are mere extensions of similar rules for carbon steel design, included in Eurocode 3, part 1.1 (EN 1993-1-1 [7]). For instance, an aspect that might severely hamper the design of stainless steel elements is the assumption of an elastic-perfectly plastic constitutive relation, which is particularly punitive for stocky elements [8]. Indeed, there are great differences concerning the material behaviours of carbon and stainless steel alloys, since the latter are characterised 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 significant amount of strain-hardening.
Due to its considerably high ductility, strength and cost, stainless steel structural solutions tend to be adopted mostly for slender/light structures (e.g., formed by thin-walled members), thus achieving a sizeable weight economy that is often combined with a strong visual impact. Nevertheless, the high slenderness of thin-walled structural members makes them prone to instability (geometrically non-linear) phenomena, thus rendering the assessment of their buckling/collapse behaviour a very complex task. Since experimental investigations are invariably limited, due to their very high cost and time consumption (including the careful preparation of the test set-up and specimens), alternative complementary approaches must be sought. The most universally employed one is the performance of sophisticated shell finite element analyses (SFEA), using non-linear constitutive laws and incremental-iterative techniques. However, this approach has some drawbacks, namely (i) the still excessively high computational effort, (ii) the time-consuming and error-prone data input, and (iii) laborious output data processing and interpretation, particularly in the context of one-dimensional members (bars)  the results consist of nodal stresses, instead of cross-section stress resultants (axial force, bending moment, etc.), the traditional and more perceptible "language" usually adopted by the technical/scientific community. Despite its relatively narrow field of application (prismatic, straight and non-perforated thin-walled members) and fairly limited dissemination, Generalised Beam Theory (GBT) has been widely recognised as a powerful, versatile, elegant and efficient approach to analyse thin-walled members and structural systems. These elegance and efficiency arise mostly from its modal nature -the displacement field is expressed as a linear combination of cross-section deformation modes with amplitudes varying along the member length. 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 [9,10], 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 moments) and (iii) materials (steel, steel-concrete, FRP). With a few exceptions, the material models adopted in these works were always elastic, with no degradation (plasticity) involved. A physically non-linear GBT formulation was first reported by Gonçalves and Camotim [11] in the context of elastic-plastic bifurcation analyses  more recently, the same authors [12,13] proposed GBT beam finite elements 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. [14,15,16] developed alternative elastic-plastic GBT formulations, also based on the J 2 -flow plasticity theory, that differ from the previous ones in the fact that (i) the deformation modes are determined by means of the procedure proposed by Silva et al. [17] and (ii) an additional degree of freedom (warping rotation) is considered.
The aim of this paper is two-fold: (i) to present the main concepts and procedures involved in the development of a physically and geometrically non-linear GBT formulation and numerical implementation (code) intended to analyse the behaviour and collapse of thin-walled members made of highly non-linear materials, and also (ii) to illustrate its application and potential, by analysing the post-buckling behaviour of lean duplex stainless steel square hollow section (SHS) columns that were experimentally and numerically investigated by Theofanous and Gardner [1]. The GBT analyses include initial geometrical imperfections, exhibiting a local and/or global nature, but does not account for membrane residual stresses and corner strength enhancement effects. The stainless steel material behaviour is modelled as nonlinear isotropic and the three-stage stress-strain curve proposed by Quach et al. [18], involving only three parameters (Young modulus E, 0.2% proof stress σ 0.2 and strain-hardening power n), is adopted. The GBT results obtained (equilibrium paths, ultimate loads, deformed configurations, displacement profiles and stress distributions) are compared with the values provided by SFEA performed in the code ABAQUS [19]. Moreover, in order to assess how membrane residual stresses and corner strength enhancements affect the column structural response, the GBT results are compared with the experimental ones reported in [1].

Brief Overview of the GBT Kinematics
Consider the local coordinate system (x, s, z) at each wall mid-surface of a thin-walled bar ( Fig. 1(a)), where x, s and z are, respectively, the longitudinal coordinate (0≤x≤L, L is the member length), the transverse coordinate (0≤s≤b, b is the wall width), and the through-thickness coordinate (-t/2≤z≤t/2, t is the wall thickness). The corresponding displacements are u (axial or warping displacement), v (transverse displacement) and w (flexural displacement), respectively. The GBT analysis of a thin-walled member consists of two main steps: (i) a cross-section analysis and (ii) a member analysis. The cross-section analysis comprises the determination of the deformation modes, i.e., their displacement profiles u k (s), v k (s) and w k (s) along the cross-section mid-line, and the evaluation of the associated modal mechanical properties. The member analysis involves the determination of the each deformation mode amplitude function  k (x), providing the longitudinal variation of the corresponding displacement profile. The GBT displacement field at the member wall mid-surface is then a linear combination of products between the cross-section modal displacement profiles and the corresponding amplitude functions, where subscript k denotes the deformation mode and satisfies the summation (Einstein) convention.  As mentioned earlier, the cross-section analysis adopted in this work is based on the approach developed by Silva et al. [17], which considers four deformation mode families: conventional, warping shear, transverse extension and cell shear flow  additional information about this cross-section analysis procedure can also be found in [20,21]. Although the existing GBT cross-section analyses consider either three (u, v, w) or four (u, v, w,  x  rotation about the x-axis) degrees of freedom (d.o.f.) per crosssection node, the fact that shear deformation and/or spread of yielding may cause a highly non-linear variation of the warping displacements along the cross-section mid-line brings about the convenience of enhancing the warping displacement representation. This is done in this work by considering a fifth d.o.f. per node, denoted as "warping rotation" and consisting of a rotation  z about the z-axis. Instead of the piecewise linear approximation of the warping displacement profiles u k (s), associated with the existing approaches, the consideration of the warping rotation  z leads to an approximation by means of piecewise cubic polynomials (as illustrated in [14]).
Finally, note that the formulation developed retains the GBT fundamental plane-stress assumption, which means that the (i) stress components σ xz , σ sz , σ zz , and (ii) strain components γ xz , γ sz and ε zz (see Fig. 1(a)) are deemed null everywhere, regardless of the material behaviour under consideration.

Equilibrium Equations
The GBT equilibrium equation can be obtained from the Principle of Virtual Work [14], where (i) σ xx , σ ss , σ xs are the longitudinal normal, transverse normal and shear stress components (2 nd Piola-Kirchhoff tensor), (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 strain components  Green-Saint-Venant strain tensor, whose components are expressed in terms of the components u, v, w (see Fig. 1(a)) as where u, v and w are given in Eq. (1) and α is a factor equal to 0 or 1, depending on whether the non-linear bending terms dependent on z are neglected or not. Note that all existing geometrically non-linear GBT formulations neglect the non-linear bending terms dependent on z and z 2 , and, although this choice has consistently led to fairly accurate results [9,13], the proposed formulation accounts for the terms dependent on z (not those dependent on z 2 ), making it possible to assess their relevance. Furthermore, the formulation contemplates (i) initial imperfections, namely residual stresses and geometrical imperfections, and (ii) arbitrary loadings, including concentrated forces and/or moments, as long as all the loads depend on a single (load) parameter .

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 [22,23] is adopted in this work and its implementation involves (i) evaluating internal force vectors f int and (ii) establishing incremental equilibrium equations, based on the tangent stiffness matrix K tan . A path-independent iterative strategy [24] was considered, meaning that the strain increments at any Gauss integration point are evaluated with respect to the last (converged) equilibrium configuration. Information on the implementation of the arc-length procedure can be found in references [14,21].
After introducing the strain components (Eqs. (1) and (3)) in the first member of Eq. (2), all amplitude functions are replaced by their finite element (FE) approximations, yielding where f is the external force vector corresponding to a unit load parameter  The numerical results presented and discussed further ahead were obtained by adopting Hermite cubic polynomials to approximate the GBT modal amplitude functions in each FE, i.e.
where (i) H Ψ is a 1x4 vector storing the Hermite cubic polynomials, (ii) d k is the corresponding displacement vector (4×1) concerning 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), and (iv) the second expression applies to all other deformation modes (v≠0 and/or w≠0).
According to the above procedure, the i th component (4x1 sub-vector) of the internal force vector f int can be obtained by means of (i) PΨ H is the primitive of Ψ H , (ii) subscripts i and k identify GBT deformation modes (i is a free index and k satisfies Einstein's convention), and (iii) the displacement vector d k and stress components σ xx , σ ss and σ xs are associated with 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). After determining the internal force vector f int , 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 [25] Paper Presented by Prof. Dinar Camotim -dcamotim@civil.ist.utl.pt (12) and corresponding to the internal force vector i th component (see (6)) Jacobian with respect to the displacement vector concerning deformation mode p  its components are d pl (l=1,…,4). The Jacobian columns are defined by the second expression in Eq. (12), where (i) each term is given by in Eqs. (7)- (9) and (iv) the stresses, stress gradients and displacement vectors are computed at "equilibrium configuration" j. Finally, for an arbitrary elastic-plastic material, the stress gradients in Eq. (13) are given at every point by where the deformation gradients ∂ε/∂d pl are defined after rewriting all deformation components in Eq. (3) according to the FE approximation in Eq. (5). The stress components and their gradients ∂σ/∂ε are obtained according to the J 2 -flow plasticity model and resorting to implicit and/or explicit numerical integration techniques, to be described in the following section.
Finally, the imposition of the boundary conditions is highlighted. For the first-order analysis of members with cross-sections having null warping and non-null shear stresses σ xs , the authors [14] concluded that the best way to reduce the effect of a shear locking phenomenon (stemming from the inadequacy of beam theories to describe rigorously both static and kinematic measures) was to guarantee that kinematic boundary conditions are satisfied, meaning that one has  k,x =0 in every section with null warping. Besides the FE mesh refinement in the vicinity of cross-sections with null warping and non-null shear stresses (e.g., a clamped support), it was also shown that the imposition of  k,xx =0 only for the warping shear modes helps in improving the accuracy of the analysis. In the case of second-order analysis, and after validating the formulation presented herein (e.g., [16]), the authors concluded that the aforementioned procedure should be kept, but without the need to impose  k,xx =0  this condition is no longer essential to guarantee accurate results because the strain-displacements relations are now non-linear.

The Prandtl-Reuss model
A detailed description of the Prandtl-Reuss model, combining the (von Mises) J 2 -flow theory with its associated flow rule, can be found in the literature (e.g., [26]) and also in recent work by the authors [14]. Therefore, only a brief overview of the main concepts and procedures involved in its development is presented herein, together with some remarks regarding the implementation of this model for non-linear materials exhibiting pronounced strain-hardening, such as stainless steel. For a material with isotropic hardening, 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 deviatoric stress invariant. If a point (σ, κ), located on the yield surface, experiences an infinitesimal plastic flow, Prager´s consistency equation [22] states that where dη ≥ 0 is the plastic proportionality factor and h is the hardening modulus (null for perfectlyplastic models). Since von Mises's yield criterion is generally appropriate to model metals, it was decided to derive dκ on the basis the modified work-hardening relation suggested by Borst [22]  which leads to dκ=dη for the plasticity model under consideration and can be used to obtain, in a straightforward way, the hardening modulus h, given by [14] y 0 where the gradient function dσ y /dε P is taken from the uniaxial stress-strain curve.
The hardening modulus plays a crucial role in defining the elastoplastic constitutive matrix [14] and, therefore, in the development of implicit and explicit integration schemes to update the stresses and the hardening parameter at any material point that has experienced plastic flow (see section 4.2). In order to simplify the application of this plasticity model to an arbitrary isotropic material, a more straightforward expression for h is derived next. For a general non-linear stress-strain curve σ(ε) with a linear elastic region, the yield stress increment (dσ y ) at any strain ε ≥ ε y 0 is given by where (i) ε y 0 is the total strain at the onset of yielding, (ii) the current total strain ε ≥ ε y 0 is divided into its elastic (ε e ) and plastic (ε p ) parts, and (iii) E is the material Young's modulus. Then, taking into account Eqs. (18) and (19), the hardening modulus can be expressed as in which the total strain ε ≥ ε y 0 appearing in the second member must be evaluated at the corresponding plastic strain ε p =κ. However, in view the relation between the total strain and its elastic and plastic parts (see Eq. (19)), the relation ε(ε p ) can only be obtained numerically. In order to avoid resorting to computationally costly numerical methods, like Newton-Raphson's (N-R), and also to enable the determination of σ y (ε p ), a database was created to provide the (i) yield stress (σ y ), (ii) plastic strain (ε p =εσ y /E ≥ 0) and (iii) hardening modulus (h), for each total strain inside the range 0 During the analysis, σ y (ε p =κ) and h(ε p =κ) are obtained from the database, using linear interpolation.
In order to compute the internal force vector and tangent stiffness matrix mentioned in section 3.2, it is mandatory to update (i) first the stresses and hardening parameter κ, and (ii) then the stress gradients (∂σ/∂ε) at each material point and for every "equilibrium configuration". For that purpose, at every Gauss integration point and after each iteration of the arc-length procedure mentioned in section 3.2, one must (i) compute the elastic stress variation, due to the imposed strain variation (Δε) w.r.t. the last equilibrium configuration i (stresses σ 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. (15)) or not 1 .
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/or explicit Euler-type methods were implemented to perform this task (an overview is presented in the next section). Once the stresses and the hardening parameter are updated, the stress gradients (∂σ/∂ε) at each integration point can be computed straightforwardly, through either the conventional or the consistent elastoplastic constitutive matrices (derived in [14]). The consistent matrix was used to obtain the results presented in this work, in order to reduce the computational cost of the incremental-iterative procedure [22].

Implicit and Explicit Integration Schemes
A key step in physically non-linear FE analysis concerns the numerical integration of the constitutive relations, so that the unknown increments of both the stresses and the hardening parameter can be obtained wherever plastic deformation takes place. The integration methods available are usually classified as explicit or implicit 2  concerning the latter, one should be aware that, when using the N-R algorithm, divergence or non-convergence may occur, particularly in the case of yield surfaces exhibiting vertices and/or high curvature gradients (in such case, some form of strain subincrementation may be necessary [27]). As for the explicit schemes, they 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 methods. According to Sloan et al. [27], the accuracy and efficiency of explicit methods can be significantly enhanced by adopting strain sub-incrementation and error control, in conjunction with a correction aimed at bringing the final stresses and hardening parameter to the yield surface  some of these recommendations were followed in this work and are discussed in this section.
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 aiming to compute, 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 plastic flow and (ii) Δt reflects the imposition of the elastoplastic strain increment (Δε=Δε e + Δε ep ) 3 -backward and forward Euler methods correspond to α=1 and α=0, respectively. However, there are other schemes proposed in the literature that adopt intermediate α values  for instance, α=1/2 has been used to define the method known as Mean Normal [22], which is unconditionally stable (Ortiz and Popov showed that this is so for methods adopting α≥1/2 [22]). Details about the definition and implementation of the forward Euler, mean normal and backward Euler methods (without any complementary algorithm) can be found in [14] and, therefore, are not discussed herein. In previous work [14,15,16,28], several GBT-based results concerning first and second-order elastoplastic analyses adopting elastic perfectly-plastic and bi-linear stress-strain curves were validated against ABAQUS SFEA results. As described in [14], either the backward Euler or the mean normal algorithms (initial solution given by the backward Euler scheme) can be used to integrate the constitutive relations  nevertheless, after performing the analyses whose results are presented in section 6, for a typical stainless steel stress-strain curve (pronounced non-linear strainhardening), the authors realised that the backward Euler method has convergence problems, which could not be eliminated even after trying (i) different arc length sizes, (ii) better approximations of the trial solution used in the N-R procedure, and/or (iii) different reliable convergence criteria [14]. Indeed, as mentioned before, convergence issues are likely to occur in the presence of high yield surface curvature gradients, which is the case of stainless steel. In view of this, it was decided to implement improved forward Euler and mean normal algorithms to efficiently perform analyses involving highly non-linear materials (avoiding convergence problems). The accuracy and robustness of these explicit methods can be improved if complemented with (i) sub-incrementation and (ii) final solution correction. Concerning the latter, the Consistent and Normal iterative schemes described in [27] are adopted if the final stresses (σ f ) and hardening parameter (κ f ) obtained by the explicit method violate the yield condition within a specified tolerance (TOL), i.e., if | f (σ f , κ f ) | > TOL. The Consistent scheme, which should be employed first, is more reliable, because it ensures that the strain increment remains unchanged (consistent with the FE method). It prescribes that   , , , , where (i) δσ and δκ are small iterative corrections that must be applied until | f (σ f , κ f ) | ≤ TOL, (ii) D e is the GBT elastic constitutive matrix and (iii) h f and n f are the hardening modulus and gradient vector, evaluated at (σ f , κ f ) [14]. If there is lack of convergence, indicated by the fact that the "corrected" stress state is further away from the yield surface than the uncorrected one, the Consistent correction scheme should be abandoned and replaced by the Normal method, which assumes that the hardening parameter remains unchanged  then, the final stresses are iteratively obtained from Regarding the convergence criterion used in both methods, it reads where F and σ y are the von Mises and yield stresses. According to Sloan et al. [27], both the algorithms presented before are efficient and robust, typically converging within 5 to 10 iterations  this assertion was confirmed by the performed GBT analyses.
The other improvement to the explicit methods that was implemented in this work is the strain subincrementation, which is only performed for the forward Euler method and always in combination with the iterative correction of the solution (see (21)- (23)). It consists of (i) dividing the total strain increment at any Gauss point into N equal sub-increments (N=3 was used), (ii) computing the corrected solution after imposing each strain sub-increment, and (iii) then using this corrected solution as the initial state for the next sub-increment, until the full total strain increment is imposed. The mean normal scheme was used to obtain the GBT results presented in section 6, but considering an initial solution at configuration t + Δt/2, provided by the forward Euler method (with sub-incrementation and correction). Moreover, since the final step of the mean normal scheme is always explicit, a correction is also applied to it.

Material Modelling of Stainless Steel Alloys
Several experimental investigations ( [29][30][31][32][33][34]) showed that the mechanical behaviour of stainless steel is asymmetric (distinct in tension and compression) and anisotropic (varies with the direction -e.g., distinct for the rolling/longitudinal and transverse directions). The relevance of these features depends on several factors, such as the stainless steel grade, manufacturing process (e.g., cold or hot-rolling, cold-forming) and strain history. In the GBT formulation proposed in this work, the stainless steel elastoplastic behaviour is modelled as non-linear isotropic, according to the J 2 -flow theory presented in section 4. This model was selected due to its frequent use in FE simulations [1,35,36], leading to a fairly good agreement between numerical and experimental results. Some researchers [31,32] have investigated the influence of the constitutive model (isotropic and anisotropic) on the structural response of stainless steel (austenitic, ferritic or duplex) thin-walled elements  they concluded that (i) the isotropic strainhardening model provides fairly accurate results and that (ii) the influence of anisotropy is negligible and may be neglected in numerical analyses involving monotonic loading.

Full-Range Stress-Strain Curve
In order to model the stainless steel alloy material behaviour, the uniaxial stress-strain curve is obtained by the full-range three-stage stress-strain relation proposed by Quach et al. [18], which (i) is valid for the tensile and compressive behaviours of any alloy (austenitic, ferritic or duplex), (ii) was shown to compare quite well with experimental curves, up to the ultimate stress, and (iii) outperforms the alternative relations proposed in the literature (e.g., Gardner and Nethercot [34] or Rasmussen [37]).
Since only one hardening function σ y (κ) can be used with von Mises's yield criterion (see section 4), it was decided to model the compressive behaviour, instead of the tensile  this choice is justified by the fact that the overwhelming majority of problems to be tackled involve compressed thin-walled members (those whose behaviour is often governed by instability phenomena Regarding this stress-strain relation, note that (i) it involves the stress and strain absolute values, (ii) the upper and lower signs in the third-stage expression of (24) concern tension and compression, respectively, and (iii) σ 0.2 , σ 1 and σ 2 are the 0.2%, 1% and 2% proof stresses  ε 0.2 , ε 1 and ε 2 are the respective total strains. The hardening power (n) definition ensures that Eq. (24) passes through the experimental point associated with the 0.01% proof stress (σ 0.01 ).
Using linear regression models, applied to a large number of uniaxial test results concerning different stainless steel alloys (either in their virgin state or cut from flat portions of cold-formed sections), Quach et al. [18] 1 , where (i) "t" and "c" stand for tension and compression and (ii) the σ t u expressions (a) and (b) are valid for (ii 1 ) ferritic and (ii 2 ) austenitic and duplex alloys, respectively. The approximation for stress σ 2 can be obtained from Eqs. (24) and (25), using an iterative scheme like N-R to solve a non-linear equation. However, in order not to increase the computational cost of the analysis, the approximation proposed in [18] Regarding the computation of the hardening modulus h(κ) (recall section 4.1 and Eq. (20)) and the yield stress σ y (κ), the database to allow for their determination by linear interpolation, for any hardening parameter κ, was created on the basis of a range of equally spaced (Δε=1x10 -6 ) total strains 0      y u . Since relation (24) does not allow a direct stress computation from a range of pre-defined strains, its accurate inverse, proposed by Abdella et al. [38], was adopted to obtain the stress range corresponding to the above equally spaced strains  the incorporation of those stresses into Eq. (24) leads to the final total strain database. Note that these calculations have a low computational cost, since they are performed only once in each analysis  before the beginning of the arc-length procedure.
where "ln" stands for the natural logarithm.

Illustrative Examples
In order to (i) validate the proposed GBT formulation, and the corresponding numerical implementation (using MATLAB R2010b [41]), and (ii) illustrate its application and potential, the geometrically and physically non-linear behaviours of two SHS cold-formed stainless steel columns, previously investigated by Theofanous and Gardner [1], are presented and discussed next: one fixed-ended stub column 5 and one long pin-ended column, both cold-rolled and seam-welded. The stainless steel alloy considered is the EN 1.4162, known as "lean duplex", which is beginning to make inroads in the building and construction industry, because (i) it is generally considerably less expensive than the traditional austenitic grades (due to the much lower nickel content), (ii) its mechanical resistance is higher than those of the most common ferritic and austenitic grades, and (iii) it still retains adequate weldability and resistance to several types of corrosion [1]. Its material behaviour is modelled by means of the Prandtl-Reuss model, presented in section 4.1, and considering the stress-strain relation given in Eq. (24). The Ramberg-Osgood material parameters adopted were those reported in [1], obtained from tensile and compressive tests of flat coupons cut from the centre of the face opposite to the weld -the uniaxial stress-strain relations used for each column are shown in Figs. 2(a)-(b). Although it is well known that cold-forming induces strength enhancements in the cross-section corner regions, their effects were not taken into account in this study (and the corners were modelled with right angles).
Residual stresses (longitudinal and transverse) may arise in cold-formed tubular sections as a result of (i) the plastic deformations taking place during sheet coiling and/or forming, and (ii) the seam-welding operation required to "close" the section. Generally speaking, there are two types of residual stresses, characterised by non-null (membrane stresses) and null (bending stresses) resultants in the throughthickness direction  rigorous measurements [42] in stainless steel hollow sections showed that bending residual stresses are clearly dominant (membrane residual stresses are rather minute). After cutting the coupons to be tested, they exhibit a longitudinal curvature due to bending residual stress release. Nevertheless, because they are essentially reintroduced when the coupons are straightened, during the test, their effect is implicitly taken into account by the material stress-strain properties [42,43]. For the above reasons, residual stresses were not explicitly modelled in this work, even if it is acknowledged that they exist and affect the column behaviour. Regarding the initial geometrical imperfections included in the analyses, which are specified in the next sub-sections, they are based on those that yielded the best results (ultimate load/displacement) in the SFE simulations reported in [1].
In the following sub-sections, the GBT results are validated through the comparison with ABAQUS SFEA values, obtained with column discretisations involving fine meshes of 4-node isoparametric SFEs with full integration ("S4 elements", in the ABAQUS nomenclature). Each SFE contains 2 integration points per surface direction and adopts the same number considered in the GBT analysis for the throughthickness direction (z-axis) -the meshes and integration point numbers considered will be given later.
In ABAQUS, the material behaviour is also based on the Prandtl-Reuss model and the compressive stressstrain input data must be expressed in terms of the true stress (σ t ) and true plastic extension (ε t(p) ) absolute values  they are obtained, from their nominal counterparts, as where σ n and ε n are the nominal stress and extension absolute values. A meaningful comparison between the GBT and ABAQUS results requires that all GBT normal stresses be transformed into true stresses at a post-processing stage. However, the GBT von Mises stresses are based on nominal stress components  those adopted to implement the J 2 -flow plasticity theory. The results presented and discussed comprise equilibrium paths, GBT modal participation diagrams, displacement profiles, stress diagrams, stress contours and column deformed configurations. The comparison between GBT and ABAQUS results concerns three equilibrium states, termed (i) BP (before peak), (ii) P (peak  ultimate load) and (iii) AP (after peak)  moreover, for each column, the experimental [1] and GBT-based equilibrium paths (including the peak load and associated displacement values) are also compared. Finally, note that (i) the cross-section dimensions are referred to the wall mid-lines, (ii) all displacements and stresses concern the column mid-surface (z=0  membrane type) and (iii) the stress diagrams are plotted along the mid-line coordinate s for the three cross-section walls labelled W1, W2, W3 in Figs. 3 and 11.

Long Pin-Ended Column
Consider the pin-ended lean duplex SHS column (labelled as 80x80x4-2000 in [1]) shown in Fig. 3(a), (i) with length L=1999 mm and the cross-section dimensions (width B, height H, and thickness t) given in Fig. 3(b), and (ii) submitted to a compressive load were employed to materialise the pin-ended boundary conditions in the experimental test performed by Theofanous and Gardner [1], only (i) the minor-axis rotation at both supports and (ii) the global axial displacement at the loaded end cross-section, were not fully restrained in both the GBT and ABAQUS analyses. In the latter, these boundary conditions were imposed at the reference nodes of the rigid plates attached to the column end cross-sections.
In the GBT analysis, the SHS was discretised into 24 wall segments (5 intermediate nodes per wall  see Concerning the longitudinal discretisation, the GBT analysis considers 48 symmetric FEs with the following distribution in half of the column (see Fig. 3(a)): (i) 3 FEs for x ≤ 0.03L, (ii) 18 FEs for 0.03L ≤ x ≤ 0.47L and (iii) 3 FEs for 0.47L ≤ x ≤ 0.5L, corresponding to a total of 3909 d.o.f.. As for the ABAQUS model, it involves 150 SFEs along the column length (uniform mesh) and 24 SFEs along the cross-section mid-line (also uniform mesh), which corresponds to 21747 d.o.f.. In each FE, the GBT numerical integration involves three Gauss points per wall segment direction (s, z and x  see Fig. 1(a)).
Regarding the column initial geometrical imperfection, they are a linear combination of the first local (symmetric) and the first global (flexural) buckling modes (λ b,glob =1.43, λ b,loc =5.76), with amplitudes    Fig. 7(a), representing the horizontal displacement of the cross-section point (node) indicated in Fig. 3(b). Fig. 5(a) displays the GBT and ABAQUS equilibrium paths λ(δ), where δ is the horizontal displacement of the mid-span cross-section point indicated in Fig. 3(b). The load parameter values obtained for the equilibrium states BP, P and AP (indicated in Fig. 5(a)) are: (i) λ GBT =0.79 and λ ABQ =0.81 (BP), (ii) λ GBT =1.04 and λ ABQ =1.04 (P), and (iii) λ GBT =0.72 and λ ABQ =0.71 (AP). Figure 5(b) shows the GBT modal participation diagram concerning the evolution of the column mid-span cross-section deformed configuration as the loading progresses (and δ increases). The observation of the results presented in Figs. 5(a)-(b) prompts the following remarks: (i) There is an excellent agreement between the GBT and ABAQUS equilibrium paths, including the descending (post-peak) branch  the maximum difference is 1.5%. (ii) The GBT column ultimate load and associated δ value read F u.GBT =376376 N (λ u =1.04) and δ u.GBT =16.45 mm, values that compare fairly well with the experimental results reported in [1]: F u.test =361900 N, λ u.test =1.00 and δ u.test =20 mm -see the dashed curve in Fig. 5(a). The closeness between the GBT and test values makes it logical to infer that neglecting both the rounded corner strength enhancement and membrane residual stresses did not affect too much the quality of the GBT analysis (e.g., the ultimate load is almost exactly predicted). (iii) The absence of post-buckling strength (λ u < λ b.glob =1.43) stems from the fact that the onset of yielding occurs considerably below the global buckling load (λ y =1.11, leading to F =A σ y 0 ).
(iv) Theofanous and Gardner [1] observed that the column failed in a global mode, which agrees with the peak load deformed configurations provided by the GBT and ABAQUS analyses -see Fig. 6. (v) The GBT modal participation diagram in Fig. 5(b), obtained as explained in [28], clearly shows that global deformation governs the column behaviour, as reflected in the dominant participations of modes 1 (axial extension) and 3 (minor-axis bending 6 ) -see Fig. 4(a). Moreover, the observation of the modal diagram leads to the following comments: (v.1) The participation of mode 1 reaches a maximum of 23.5% close to the BP equilibrium state and decreases subsequently, tending to 7.0% in the post-peak stage. As for the participation of mode 3, it reaches a minimum value of 75.2% (practically coincident with the maximum of mode 1) and then gradually increases until about 92%, value exhibited in the post-peak stage. (v.2) In spite of appearing in the initial geometrical imperfection, mode 6 (local  see Fig. 4(b)) only shows some (minute) relevance for δ >40 mm, in agreement with the fact that visible local deformations only develop (at column mid-span region) beyond the AP state.
GBT ABAQUS Figure 6. GBT and ABAQUS column post-collapse deformed configurations (AP equilibrium state). Fig. 7(b) depicts the longitudinal profiles of the column mid-height lateral displacement (see Fig. 3(b)) at the equilibrium states P and AP  the GBT and ABAQUS analyses provide almost coincident results. Local bending is not visible, because it is completely overshadowed by global bending (as shown in Fig. 5(b))  however, note that it clearly appears in the column initial (imperfect) configuration, as attested by the corresponding longitudinal profile shown in Fig. 7(a) (see the column mid-span region).

Figs. 8(a)-(b)
show the GBT and ABAQUS longitudinal normal ( xx ) and von Mises ( Mises ) stress diagrams at the mid-span cross-section (x=999.5 mm), for the equilibrium states BP, P and AP. Their observation makes it possible to draw the following conclusions: (i) There is a very good agreement between all the GBT and ABAQUS stress diagrams  as expected, they are symmetric with respect to the cross-section horizontal principal axis (see Fig 3(b)). to the absolute values of the corresponding longitudinal stress diagrams, which means that the shear and transverse normal stresses are very small (negligible) at the mid-span cross-section. This assertion also holds for equilibrium state AP, but only for points where the  xx values are reasonably high.
(iii) Fig. 8(a)  have already yielded in tension (see Fig. 8(b) and recall that σ 0 y =347.34 N/mm 2 ). (iv) At equilibrium state BP, the whole mid-span cross-section is in the initial elastic regime (σ Mises ≤ σ 0 y ) and, therefore, the stress diagrams are linear or constant in each wall. At equilibrium states P and AP, on the other hand, the stress diagram non-linearity is clearly visible in the W1 and W3 wall yielded areas (σ Mises >σ 0 y )  it reflects the non-linear hardening exhibited by the lean duplex stainless steel.
Figs. 9(a)-(b) show the GBT and ABAQUS  Mises contours concerning the equilibrium state AP and one notices a remarkable resemblance in both pairs of column side views displayed, showing either the compressed or the tensioned region. As mentioned earlier (see also Fig. 8(b)), the column side under tension, depicted in Fig. 9(b)), has also yielded, but to a lesser extent than its compressed counterpart. Note that the column regions exhibiting lower stress values (coloured in blue) correspond to the vicinity of the column neutral surface  these "blue regions" tend to move away from the principal (centroidal) plane as one travels from the end section to mid-span, which is a direct consequence of the load eccentricity increase (null at the end section and highest at mid-span).   since the stresses are much lower in the remaining column regions, thus leading to a homogeneous contour, only this specific region is dealt with. It is observed that the GBT transverse normal stresses do not vary as smoothly as their ABAQUS counterparts  nevertheless, there is clear good qualitative agreement between the two contours. It should be mentioned that the current GBT formulation may lead to visible transverse normal stress discontinuities between adjacent cross-section wall segments, which is due to the linear approximation adopted for the membrane transverse displacements v(s).
Nevertheless, previous work from the authors [14,16,28] (i) included GBT  ss diagrams that compared much better with their SFEA counterparts, provided that the mean  ss value in each wall segment was considered (this fact explains the good correlation between the GBT and ABAQUS stress contours), and (ii) showed that the  ss discontinuities have very little impact on the overall (high) quality/accuracy of the GBT results.
GBT ABAQUS Figure 10. Transverse normal stress ( ss , N/mm 2 ) contours at collapse (vicinity of column end section  along its height).

Fixed-Ended Stub Column
Consider now the fixed-ended lean duplex SHS stub column shown in Fig. 11(a), with length L=400 mm and the cross-section dimensions (width B, height H, and thickness t) given in Fig. 11(b)  its is labelled as 100x100x4 -SC1 in reference [1]. The column is subjected to a compressive force F=1022000 λ N (λ is the load parameter) and, as before, the load is uniformly distributed along the end section mid-line in both the GBT and ABAQUS models (q=2618.63 N/mm). The stainless steel (EN 1.4162) is modelled according to the Ramberg-Osgood basic parameters (i) 0. corresponds to a total of 865 d.o.f.. The ABAQUS model involved 50 SFEs along the column length (uniform mesh) and 48 SFEs along the cross-section mid-line (also uniform mesh), corresponding to a total of 14401 d.o.f.. Finally, the GBT numerical integration adopted three Gauss points per wall segment direction in each FE, and the initial geometrical imperfection included in the analyses exhibits the shape of the first symmetric local buckling mode (λ b =1.8368) and amplitude t / 100=0.0393 mm -the corresponding longitudinal profile, concerning the vertical displacement δ z (see Fig. 11(b)), is depicted in Fig. 15(a). Fig. 13(a) displays the GBT and ABAQUS equilibrium paths λ(Δ), where Δ is the axial shortening (see Fig. 11(a)). The load parameter values corresponding to the equilibrium states BP, P and AP (indicated   Fig. 13(a)) are: (i) λ GBT =0.701 and λ ABQ =0.701 (BP), (ii) λ GBT =0.893 and λ ABQ =0.893 (P), and (iii) λ GBT =0.590 and λ ABQ =0.577 (AP). As for Fig. 13(b), it provides the GBT modal participation diagram concerning the evolution of the column mid-span cross-section deformed configuration as the loading progresses (and  increases). The observation of these results prompts the following remarks: (i) The agreement between the GBT and ABAQUS equilibrium paths is excellent, up to the peak load, and very good in the post-collapse stage  the maximum difference is 4.9% and occurs at the beginning of the descending branch, characterised by a very pronounced non-linearity. (ii) The GBT column ultimate load and associated  value read F u.GBT =912646 N (λ u =0.893) and  u.GBT =3.13 mm, values that compare fairly well with the experimental results reported in [1]: F u.test =102000 N, λ u.test =1.000 and  u.test =3.63 mm -see the dashed curve in Fig. 13(a)  the differences between the GBT and test values point to the need to improve the developed GBT code to be able to simulate numerically this particular stub column behaviour. The SFEA results reported by Theofanous and Gardner [1] suggest that including in the model the corner strength enhancement effects (obtained from coupon tests) will lead to a much better prediction (even if the influence of membrane residual stresses is neglected). (iii) The lack of post-buckling strength (λ u < λ b. ) is just a logical consequence of the column low slenderness, which implies that premature yielding (λ y =0.586, leading to F=A σ y 0 ) precludes its strength to approach even the local buckling load level. (iv) Theofanous and Gardner [1] reported that the stub column failed (iv 1 ) in a mode exhibiting local deformations and (iv 2 ) after considerable plastic deformation, which fully complies with the postcollapse deformed configurations obtained from the GBT and ABAQUS analyses -see Fig. 14. (v) The GBT modal participation diagram in Fig. 13(b) clearly shows the presence of two markedly different behaviours: (i) a first one that is predominantly global and involves basically the (prebuckling) axial deformations associated with the axial extension mode 1 7 , and (ii) a second one that corresponds to the emergence and gradual growth of local deformations, mostly due to the contribution of mode 6 (see Figs. 12(a)-(b))  it is worth noting that this mode becomes much more dominant in the post-collapse stage (it literally "sweeps aside" mode 1). Moreover, the observation of Fig. 13(b) also leads to the following remarks: (v.1) In the ascending part of the column equilibrium path ( ≤ Δ u =3.13 mm), there are contributions from four modes, even if the clear dominance belongs naturally to mode 1, whose contribution rises up to 77.3% (recall the content of footnote 7) before coming down a little. Mode 6, whose growing emergence is partly (together with plasticity) responsible for the non-linearity of the column equilibrium path, has a participation that rises up to 28% at collapse (and keeps rising in the descending branch). As for the transverse extension modes 166 and 168, whose presence is intrinsically linked to that of mode 6, they exhibit increasing contributions that go up to about 6% each  it should be pointed out that, in spite of their fairly small contributions, the removal of these two transverse extension modes from the GBT analysis would considerably affect the accuracy of the results obtained  the simulated column behaviour would be stiffer, because the local deformation would be "artificially restrained".
features of the column behaviour and, by removing the non relevant deformation modes from the analysis, (ii) to further reduce the number of degrees of freedom involved and, therefore, increase the computational efficiency. For these two particular examples, the GBT analyses involved about 18% (long column) and 6% (stub column) of the number of d.o.f. required by the SFEA, in order to provide similarly accurate results  for instance, the ultimate loads obtained from the two analyses were found to be almost coincident for both columns.
In order to assess the influence of the (i) corner strength enhancement and (ii) membrane residual stresses on each column load-carrying capacity, the GBT ultimate loads were compared with the experimental ones reported in [1]  it was concluded that the above effects have (i) a small influence on the long pin-ended column ultimate strength (λ u.GBT =1.04 vs. λ u.test =1.00), but (ii) a moderate impact on the fixed-ended stub column ultimate strength (λ u.GBT =0.893 vs. λ u.test =1.000). In order to overcome this shortcoming and be able to predict the behaviour and ultimate strength of thin-walled stainless steel members more accurately, the corner strength enhancement and membrane residual stress effects will be soon included in GBT-based analyses.
Finally, just two words to mention that the modal nature of the GBT analysis, which was only very partially explored in this work, is expected to provide fresh insight on the mechanics underlying the behaviour and collapse of thin-walled stainless steel members under various loadings. The authors plan to investigate this issue in the not too distant future, now that the developed geometrically and physically non-linear GBT formulation, and respective beam finite element code, may be deemed properly validated.