Geometrically and Physically Non-Linear GBT-Based Analysis of Thin-Walled Steel Members

This paper presents and illustrates the application of an elastic-plastic Generalised Beam Theory (GBT) formulation, based on J 2 -flow plasticity theory, that makes it possible to perform physically and geometrically non-linear (post-buckling) analyses of prismatic thin-walled members (i) with arbitrary cross-section shapes, (ii) exhibiting any type of deformation pattern (global, local, distortional, warping, shear), (iii) made from non-linear materials with isotropic strain-hardening and (iv) containing initial imperfections, namely residual stresses and/or geometric imperfections, having generic distributions. After providing a brief overview of the main GBT assumptions, kinematical relations and equilibrium equations, the development of a novel non-linear beam finite element (FE) is addressed in some detail. Moreover, its application is illustrated through the presentation and discussion of numerical results concerning the post-buckling behaviour of a fixed-ended I-section steel column exhibiting local initial geometrical imperfections, namely (i) non-linear equilibrium paths, (ii) displacement profiles, (iii) stress diagrams/distributions and (iv) deformed configurations. For validation purposes, the GBT results are also compared with values yielded by A BAQUS rigorous shell FE analyses.


Introduction
The use of high performance steels leads to considerable weight savings in built structures. In fact, steel structural elements, namely cold-formed ones, often exhibit slender thin walls, which makes them highly susceptible to global (flexural, torsional and flexural-torsional) and local instability phenomena, which play a key role in their ultimate strength and design. Up to now, the numerical determination of accurate collapse loads for such elements (e.g., columns or beams) has only been possible by resorting to complex non-linear shell finite element analyses (SFEA). Since this task is very timeconsuming, due to the large amount of degrees of freedom (d.o.f.) involved and the laborious/difficult data input and result interpretation, a very promising alternative to the above approach is the use of one-dimensional (beam) finite element formulations based on Generalised Beam Theory (GBT).
Due to its unique modal features, GBT is widely recognized as an elegant, computationally efficient and structurally clear approach to analyse prismatic thin-walled members and structural systems. The displacement field is expressed as a linear combination of cross-section deformation modes whose amplitudes vary along the member length. In the last decade, GBT formulations have been developed to cover different (i) analyses (first-order, buckling, vibration, post-buckling) and (ii) materials (steel and composite − FRP and steel-concrete) [1]. Most of these works assume a linear elastic material behaviour, with no plasticity (or any other degradation source) involved. The first physically non-linear GBT formulation was developed by Gonçalves and Camotim [2], in the context of elastic-plastic bifurcation, and these authors [3,4] recently proposed GBT-based FE formulations based on the J2-flow plasticity theory and intended to perform physically and geometrically non-linear analyses. In parallel, Abambres et al. [5,6] developed alternative GBT first-order elastic-plastic formulations, also based on the J2-flow plasticity theory, to analyse members with (i) global (axial, flexural, torsional) and (ii) an arbitrary deformation pattern. They differ from those reported

Brief Overview of the GBT Assumptions and Kinematics
At each wall mid-surface, consider the local coordinate system (x, s, z) illustrated in Fig 1(a), where x, s e z are 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. The corresponding local displacements are (i) u (along x − warping), (ii) v (along s − transverse) and (iii) w (along z − flexural). A GBT analysis consists of (i) a cross-section analysis and (ii) a member analysis. While the former comprises the determination of the cross-section deformation modes (uk(s), vk(s), wk(s)) and associated mechanical properties, the latter leads to the corresponding modal amplitude functions k(x), characterising the variation of the deformation mode amplitude along the member length. Each member mid-surface displacement field component is expressed as a linear combination of products involving cross-section deformation modes and their amplitude functions, where subscript k follows the Einstein (summation) convention.
The cross-section analysis procedure proposed in [7] is adopted: four deformation mode families (conventional, warping shear, transverse extension and cell shear flow) are defined and obtained by solving a sequence of auxiliary eigenvalue problems − detailed information in [7,9,10]. However, this cross-section analysis also includes an innovation recently proposed by the authors [11]: the inclusion of the warping rotation (about z) as a fifth d.o.f. in each node. Finally, it is worth noting that the formulation developed retains the GBT fundamental plane-stress assumption (after all, GBT is a theory for thin-walled members), which means that (i) the stress components σ xz , σ sz , σ zz , and (ii) the 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, stating that [11] ( ) ( ) ( ) ,         dz ds dx q u q v q w ds dx , (2) where (i) σ xx , σ ss , σ xs are the axial normal, transverse normal and shear stress components, (ii) qx, qs, qz are the local components of an external distributed force applied at the member mid-surface (see Fig.   1(b)), and (iii)  xx ,  ss , γ xs are the virtual Green-Saint-Venant strain components: axial extension, transverse extension and shear strain, expressed in terms of the displacement components u, v, w 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 . Although this choice has consistently led to fairly accurate results [1,4], the formulation proposed retains the terms dependent on z (not those dependent on z 2 ), in order to assess their relevance. Moreover, all stress and strain components are evaluated taking into account the possible presence of initial imperfections, namely residual stresses and/or geometrical imperfections.

Non-Linear Beam Finite Element
In non-linear analysis, reaching an equilibrium configuration requires using an incremental-iterative strategy. This work adopts the cylindrical arc-length method [12,13], whose implementation involves evaluating internal force vectors f int and subsequently establishing incremental equilibrium equations, based on the tangent stiffness matrix Ktan. After introducing the strain components considered in the first member of Eq. (2), all amplitude functions in Eq. (1) are replaced by their 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 longitudinal variation of the GBT deformation mode amplitudes, i.e., where (i) k H Ψ is a 1x4 vector storing the Hermite polynomials adopted to approximate the mode k amplitude function, (ii) dk is the corresponding displacement vector (4×1), (iii) the first expression concerns only the axial 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 where , , , (i) subscripts i and k identify GBT deformation modes (i is a free subscript and k satisfies Einstein convention), and (ii) the displacement vectors dk and stress components σ xx , σ ss , σ xs are associated with a generic "equilibrium configuration" (during the iteration procedure and wherever 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 displacements vector and (ii) the tangent stiffness matrix (i, p) th component, Kip,tan, is a 4x4 sub-matrix concerning deformation modes i and p, which is given by [ (9), and (iii) 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 defined at every point by where the deformation gradients ∂ε/∂dpl 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 from (i) the elastic stress-strain law routinely used in GBT [1] and (ii) the J2-flow model (associated flow rule) − the consistent elastic-plastic constitutive matrix is used to obtain ∂σ/∂ε and the mean normal return-mapping algorithm is adopted to evaluate σ if plastic flow occurs [12]. For the plasticity model implemented, after each iteration one must sequentially compute, at each point, (i) the strain variation Δε w.r.t. the last equilibrium configuration (pathindependent strategy), (ii) the elastic stress variation Δσ (due to Δε), added to the previous stress σi (last equilibrium configuration) to provide the predicted stress σpred = σi + Δσ, and, if these stresses fall outside the previous yield surface (dependent on the deformation history if there is strainhardening), (iii) the stress correction Δσcor (to σpred), by using the return-mapping algorithm, thus leading to the final stress σf = σpred + Δσcor, "located" on the current yield surface. Due to space limitations, the implemented (J2-flow) plasticity model and arc-length procedure are not presented here − they can be found in [6,11].

Illustrative Example
The developed physically and geometrically non-linear GBT formulation was implemented using MATLAB R2010a [14] and, for validation purposes, the results obtained are compared with values yielded by an ABAQUS [8] SFEA based on the J2-flow theory − since ABAQUS provides true stress outputs, the GBT normal stresses results are converted to true stresses through σ t =σ n (1+ε n ), where t and n denote true and nominal (the GBT von Mises stresses are based on nominal stresses).
The numerical results presented and discussed concern the fixed-ended I-section steel column shown in Fig. 2(b), (i) with length L=1500 mm and the mid-line cross-section dimensions (flange width bf, web height bw and thickness t) given in Fig. 2(a), and (ii) under uniform compression F=10000 λ N (λ is the load parameter) − the load is deemed uniformly distributed along the end section mid-lines in both the GBT and ABAQUS models (q=30.30 N/mm). The steel constitutive behaviour is assumed linearelastic/perfectly plastic, and characterised by (i) Young's modulus E=200000 N/mm 2 , (ii) Poisson's ratio ν=0.3 and (iii) nominal uni-axial yield stress σ y =460 N/mm 2 (shear modulus G=E/2(1+ν)).The GBT cross-section discretisation, shown in Fig. 3  points involved in the numerical integration are (i) 4 Gauss per FE (x direction) and (ii) 4 (s direction) and 3 (z direction) per wall segment (see Fig. 1(a)).
In ABAQUS, four-node isoparametric shell FE with full integration (S4) are adopted and 3 Gauss points are involved in the through-thickness integrations. The FE mesh comprises 100 FE along the column length and 22 along the cross-section mid-line contour − 15x15 mm dimensions. In order to ensure fixed end supports, "layers" of 22 additional "rigid" (E=2x10 8 N/mm 2 ) FE, with t=100 mm, were attached to each column end cross-section − only all their axial d.o.f. are not fully restrained.
Regarding the initial imperfections, no residual stresses are considered but the column contains criticalmode (local) geometrical imperfections with amplitude 0.25 t=0.5 mm − in both the GBT and ABAQUS analyses, this shape was obtained by performing a buckling analysis. Note that the results presented include (i) stresses and displacements at z=0 (membrane) and (ii) stress diagrams concerning the mid-span cross-section (x=750 mm), plotted along the walls indicated in Fig. 2(a) (dashed arrow). Fig. 4(a) shows the GBT and ABAQUS elastic-plastic equilibrium paths λ(δ), as well as the ABAQUS elastic path (for comparison purposes), where δ is the horizontal displacement of the web mid-point at mid-span (see Fig. 2). It is observed that the two elastic-plastic curves are in very good agreement, as (i) they fully coincide in the elastic range and (ii) do not differ by more than 1.8% after the onset of yielding (λ≈15). Moreover, note that the elastic-plastic equilibrium paths strongly diverge from their ABAQUS elastic counterpart -δ keeps increasing monotonically at a decaying rate, unlike in the elastic case, where a steep/stiffer increase is followed by a clear decrease, which is due to emergence of minor-axis bending, involving δ values in the opposite direction. For the three equilibrium states identified by white circles in Fig. 4(a) Fig. 4(a). Minor-axis bending also emerged in the elastic-plastic behaviour, but its contributions only became meaningful after the failure load (limit point C) was reached − this is because, due to the spread of plasticity, the cross-section stiffness centre (effective centroid) shifts, creating an applied load eccentricity. Fig. 4(b) shows the GBT and ABAQUS deformed configurations at the PC equilibrium state (amplified 5 times) and it is possible to notice their remarkable similarity. Note also that the higher deformations occur at the mid-span region (not visible in Fig. 4(b)), due to minor-axis bending − this feature can be easily quantified by means of the corresponding GBT modal participation [11]: 0.019%, 0.057% and 8.56%, respectively for equilibrium states E, C and PC.
Finally, the substantial difference between the failure (λu=17.5) and critical (λcr=9.9) load parameter values, corresponding to λu /λcr=1.77, evidences the significant column post-buckling strength, as it would be logical to expect in view of the collapse mode local nature. The dominance of the secondorder effects can be assessed by the fact that the failure is only 57.64% of the squash load Py=A σ y . Figs. 5(a)-(b) concern configuration E, C and PC and compare the δ longitudinal profiles obtained with the GBT and ABAQUS analyses, evidencing either (i) an excellent agreement (E) or (ii) a fairly good correlation (C and PC). Moreover, note the difference between the positive and negative amplitudes in the PC equilibrium state, reflecting the presence of minor-axis bending (causing positive δ values). Figs. 6(a)-(b) concern configurations E and C and depict the mid-span axial (σxx) and Mises (σMises) stress diagrams. Even if the GBT and ABAQUS values agree fairly well, there are large discontinuities in the GBT C diagram, which are due to an approximation in the transverse stresses (σss) shown in Fig. 7(a). Indeed, since the v displacements are approximated by piecewise linear functions in each wall segment, there is no transverse extension (v,s) continuity between adjacent wall segments, which affects the σxx and σMises values (all stress components "interact" in the plastic range) − nevertheless, note that the ABAQUS and GBT σss values are quite close at wall segment mid-points. Moreover, the ABAQUS and GBT web σss profiles at state C are qualitatively similar (see Fig. 7(b)).

Conclusion
A second-order elastic-plastic GBT formulation was proposed and its application was illustrated by means of numerical results concerning an initially imperfect linear-elastic-perfectly plastic fixedended I-section steel column. The GBT results agree very well with ABAQUS SFEA values, in spite of the huge disparity between the d.o.f. numbers involved (2537 vs. 13984), which confirms that the GBT approach constitutes a viable and efficient alternative to perform geometrically and physically non-linear analyses of thin-walled members and structural systems.