First order elastoplastic GBT analysis of tubular beams

: This paper aims to present an original formulation of Generalised Beam Theory (GBT) in-tended to perform first order elastoplastic analysis of thin-walled members, made of isotropic non-linear material and subjected to arbitrary deformation. The J 2 -flow theory is used to plasticity in conjunc-tion with the Euler-Backward return-mapping algorithm. After presenting the formulation, its application is illustrated by means of the first order analysis of beams with (i) rectangular hollow section (RHS) and (ii) LiteSteel section, made of an elastic-perfectly plastic material and subjected to distributed and point loading, respectively. The GBT results, which include equilibrium paths, displacement profiles, stress diagrams, 3D stress/displacement contours and deformed shapes, are compared with the ones obtained by A BAQUS code using a shell finite element model. GBT and A BAQUS results display a very good agreement .


INTRODUCTION
Materials with non-linear behaviour (e.g. metals) are widely used in construction industries, like civil, mechanical and aeronautical engineering. This use is due to the fact that these materials are able to simultaneously assure outstanding properties such as strength, stiffness, ductility and tenacity. These high performance properties often lead to the design of structural elements displaying thin walls and rather slender cross-sections and exhibiting high susceptibility to local deformations. One very interesting thinwalled steel member is the LiteSteel beam (LSB), which has a channel section combining two hollow flanges and a slender web. This shape has proven to be very competitive in comparison with hot rolled profiles.
The numerical determination of accurate collapse loads of structural members has been only possible by means of shell finite element (FE) analyses using non-linear constitutive laws and incremental-iterative techniques. Since this task requires time-consuming procedures, a very promising alternative to this approach is the use of one-dimensional models (beam FE) based on Generalised Beam Theory (GBT). Currently, GBT is widely recognized by the scientific community as a powerful, versatile, elegant and efficient thin-walled beam theory. It is valid for prismatic bars only and its elegance arises from the fact that the displacements field is obtained as a linear combination of cross-section deformation modes, which can be carefully determined and whose amplitudes along the longitudinal axis are the problem un-knowns. Until the end of last century, GBT was developed to perform first order and linear stability analysis. Since 2002, GBT has attracted the interest of several researchers, which led to the development of new formulations and applications. In this regard, Camotim et al. (2010) successively developed, validated and illustrated the application of GBT to (i) stability and vibration analyses of orthotropic members, and (ii) post-buckling analyses of isotropic members. The theory has been extensively upgraded at the Technical University of Lisbon (Camotim et al. 2010) and has been applied to distinct materials (steel, steel-concrete, FRP). In these contributions, the material was always assumed elastic, with no degradation (plasticity) involved. The first physically non-linear GBT formulation was presented by Gonçalves and Camotim (2004), in the context of elastoplastic bifurcation analysis. Recently, Gonçalves and Camotim (2011a, b) proposed a finite element based on the J 2 -flow plasticity theory in the context of GBT. Abambres et al. (2011a) developed a GBT first order formulation for the elastic-plastic analysis of members undergoing only global (axial, bending, torsion) deformation. In order to overcome this limitation, the formulation presented in this paper (i) is valid for a non-linear isotropic material subjected to arbitrary (local and/or global) deformation, and (ii) is based on the J 2 -flow plasticity theory. Although this plasticity theory was already considered in the aforementioned work of Gonçalves and Camotim (2011a, b), the GBT formulation presented herein is based on different types of deformation modes (see section 2). In order to illustrate the application of this formulation, analyses of RHS and LSB fixed-pinned beams made of elastic-perfectly plastic material were performed. In section 3, GBT equilibrium paths and several displacement/stress results are compared with the ones yielded by ABAQUS software (Simulia 2004) using a shell FE model.

BRIEF OVERVIEW OF GBT FORMULATION
Due to space limitations, the GBT formulation is briefly presented in this section. The interested reader is referred to detailed information in Abambres et al. (2011a, b). A GBT analysis of a structural element consists of two main stages: (i) a cross-section analysis followed by (ii) a member analysis. One of the features that better characterize this theory is its modal nature, since the displacement field is obtained as a linear combination of cross-section deformation modes. These modes approximate the displacements along the cross-section mid-line (s coordinate). Fig. 1 illustrates the local coordinates and displacements used in GBT. The cross section analysis considered in this work is based on the approach developed by Silva et al. (2008), which consists of (i) 4 deformation mode families obtained via different eigenvalue problems and (ii) 4 dof per cross-section node (three displacements and one transverse flexural rotation). Although these mode families are considered in the present formulation, the initial (elementar) deformation modes are based on a set of 5 degrees-of-freedom (dof) per crosssection node 1 . Besides the usual 4 dof per node, the present formulation also includes a rotation corresponding to a vector contained in the cross-section plane, which is an original contribution of the present work. This new dof, called warping rotation, is used in each node and wall direction, making it possible to approach the warping profile u(s) by means of piecewise Hermite cubic polynomials, instead of usual piecewise linear functions (Abambres et al. 2011a). The GBT equilibrium equation states that It is well known that the achievement of a nonlinear equilibrium path requires the use of an incre-1 Every node defined in the cross-section (natural or intermediate) will be designated in this paper as "section node". mental-iterative strategy. Consider a generic bar FE and its incremental relation of an arbitrary point of the equilibrium path int tan where f int is the internal force vector, λ is the load parameter and K tan is the tangent stiffness matrix at the equilibrium point. Due to use of a non-linear constitutive law (), the stress-displacements relations (d) are non-linear. Thus, the relation expressed by (2) can be determined based on the "linearization" of relations where  mn is a generic 2 nd Piolla-Kirchhoff stress com-ponentSince (i) these components depend on the deformation components (ii) depend on the longitudinal deformations Abambres et al. 2011a), and (iii)  depend on the generalized displacements d (unknowns of the FE method), relation (3) can be rewritten as where k and i are indexes related with the deformation mode k and the component i of displacements vector d (satisfying the summation convention). The gradients / are given by the sectional deformations E Abambres et al. 2011a) and /d depend on the type of approximation chosen for the FE. In this work, it is based on the Hermite cubic polynomials (Ψ H ) and is given by where k denotes the deformation mode k, Eq. (5) is related to all the modes k with in-plane displacements (v, w≠0), Eq. (6) is related to axial and warping shear modes only (v, w=0) and d k is the elementar displacements vector (4×1). Once defined the FE approximation, the insertion of Eqs. (4)-(6) in Eq. (1) gives rise to the incremental equilibrium relation where the component i-k of the tangent stiffness matrix is a 4×4 sub-matrix related to modes i and k and is defined in Abambres et al. (2011b) for the equilibrium point j. However, it is important pointing out that expressions presented in that work are not valid for axial and warping shear modes. The way to obtain the components i-k due to those modes consists in replacing the correspondent vectors Ψ and their derivatives by vectors with an order of differentiation reduced by 1 (check Eqs. (5)- (6)).
In order to complete the definition of all the variables needed to perform an incremental-iterative procedure, it is also necessary to define the vector of internal forces for a generic point at the equilibrium path, which arises from the first member of Eq. (1), if Eqs. (5)-(6) are taken into account. For all modes displaying in-plane displacements, the component k (4×1 sub-vector) of the elementar (FE with length L e ) internal force vector at the j th equilibrium point is obtained by The way of defining the components related to axial and warping shear modes is identical to the aforementioned procedure regarding tangent stiffness matrix. Concerning the plasticity model, it was considered a J 2 -flow theory with associated flow rule (Borst and Sluys 2007). In order to compute the tangent stiffness matrix presented in Eq. (7), the consistent elastoplastic constitutive matrix D ep* was used (Abambres et al. 2011b), which is known to reduce the computational cost of the incremental-iterative procedure when compared to the conventional constitutive matrix (Borst and Sluys 2007).
In order to compute the tangent stiffness matrix and the internal forces vector defined in Eqs. (7)-(8), one must be able to update stresses and the plastic proportional factor at any point and any equilibrium instant. For that purpose, a robust return-mapping scheme should be used for those points that undergo plastic flow during a deformation increment. In this work, the Euler-backward method based on a Newton-Raphson scheme was considered (Borst and Sluys 2007). The GBT formulation was implemented using MATLAB (R2010a) software (MathWorks 1998).

ILLUSTRATIVE EXAMPLES
Next, two illustrative examples are presented and the GBT results are compared with the ones yielded by ABAQUS (Simulia 2004) using S4 shell FE model. The cross-section dimensions are referred to section mid-line. The material behavior is linear elastic-perfectly plastic and isotropic (Young modulus E=200000N/mm 2 , Poisson ratio ν=0.3, uniaxial nominal yield stress σ y =450N/mm 2 ). In order to perform an effective comparison between GBT and ABAQUS results, 3 points of the equilibrium path ( curve) were selected, namely: (A) elastic, (B) elastoplastic and (C) plastic. In order to make a correct comparison between GBT and ABAQUS stress results, it is worth pointing out that all GBT normal stress outputs were transformed to true stresses in post-processing. However, GBT Von Mises stresses were based on nominal stress components as those were the ones used when implementing J 2 -flow theory. The stresses distributions and displacements profiles are obtained at z=0 and correspond to membrane type (see Fig. 1).

RHS beam
Consider the RHS fixed-pinned beam depicted in Fig. 2. The beam has a length L=4000 mm and is acted by an uniformly distributed load p = 0.01λ N/mm 2 (Fig. 2(b)) applied at the top flange ( Fig.  2(a)). Only 17 (out of 74) modes were used in the GBT analysis. The RHS node discretization used in GBT, as well as some of the most relevant modes, are represented in Fig. 3. Regarding FE models, (i) 16 beam FE's were used in GBT model (6 FE's in x ≤ 0.15L and 10 FE's in x > 0.15L) and (ii) 160 "transverse rows" of 36 shell FE's were used in ABAQUS model. For the numerical integration purposes, (i) 5 Gauss points were used in all directions (s, z, x) of the GBT model and (ii) 5 Gauss points were used in through-thickness direction of the ABAQUS FE model. All stress diagrams presented in this example are referred to the section located at x = 1996 mm (see Fig. 2(b)) and to the section path depicted by the arrowed line in Fig.2(a). In Fig. 4, the equilibrium paths λ(δ) obtained by GBT and ABAQUS are presented, where δ is the monitoring displacement represented in Fig. 2 and located  at x=2000 mm. It is observed that the λ(δ) curves yielded by GBT and ABAQUS compare very well, with the maximum difference never exceeding 1.6%. The 3 equilibrium points chosen to compare the results correspond to the following load parameters (see Fig.4 Fig. 2(a) and show a very satisfactory agreement between ABAQUS and GBT. Regarding point C, the difference at the maximum absolute displacement is about 4.5%. In Fig. 6. the axial stress diagrams for points A and C are presented. Whereas for point A there is a quite satisfactory agreement between GBT and ABAQUS for the whole domain of the cross-section, for point C there are notorious similarities along the top flange (150≤s≤450 mm) but relevant differences at the webs (0≤s≤150 mm and 450≤s≤600 mm). However, one should bear in mind that point C is the one where load parameter  has greater difference between ABAQUS and GBT, and the results regarding the plastic response are more susceptible to numerical errors.  , but there are huge discrepancies along the webs (0≤s≤150 mm and 450≤s≤600 mm). Despite these dissimilarities, it can be noticed that the mid-points of GBT lines are very close to ABAQUS diagram. The existence of these notorious discrepancies may be attributed to the use of linear functions to approximate the transverse (v) displacements between section nodes. Additionally, the GBT analysis adopted only two nodes in the web (see Fig. 3) while the ABAQUS analysis used five.   Fig. 8 shows quite satisfactory results regarding the shear stress diagrams for points A and B, although some important differences are seen in the top flange for point B. However, one should bear in mind that GBT and ABAQUS load parameters coincide for point A but do not for point B (see Fig. 4). Concerning Mises stress diagrams, and taken into account the differences between GBT and ABAQUS load parameters, Fig. 9 shows that results from both analyses are very close. Because point C is the one where discrepancies are bigger, Von Mises stress contours for that point are shown in Fig. 10. It can be noted that the spread of plasticity in ABAQUS and GBT exhibit a very good agreement. Lastly, Fig. 11 displays the collapse deformed shape for both analyses: it is evident that local deformation takes place where plastic deformation is bigger. Moreover, it is also notorious the existence of shear deformation at clamp's vicinity.

LiteSteel beam
Now, consider the fixed-pinned LSB with length L=4000 mm depicted in Fig. 12. The beam is loaded at L/4 (x=1000 mm) and 3L/4 (x=3000 mm) by two vertical concentrated loads F=1000 λ N located in the middle of the web, as illustrated in Fig. 12. In this case, 70 (out of 93) GBT modes were used in the analysis. The cross-section discretization and some of the most relevant modes used in GBT are represented in Fig. 13. Regarding the FE models, (i) 40 beam FE's were used in GBT model (6 in 0<x<0.2L, 12 in 0.2L≤x≤0.3L, 4 in 0.3L≤x≤0.7L, 12 in 0.7L≤x≤0.8L and 6 in x > 0.8L) and (ii) 184 "transverse rows" of 28 shell FE's were used in ABAQUS model. For the numerical integration, (i) 5 Gauss points in s and x directions were adopted in the GBT model and (ii) 7 Gauss points in throughthickness direction were used in both ABAQUS and GBT models. All the stress diagrams presented in this example are referred to the section located at x = 2000 mm (see Fig.12(b)) and to the section path depicted by the arrowed line in Fig.12(a). In Fig. 14, the equilibrium paths λ(δ z ) obtained by means of GBT and ABAQUS are presented, where δ z is the monitoring displacement represented in Fig.12 and located at x=2000 mm. The curves yielded by both analyses compare very well, with a maximum difference of 2%. The 3 equilibrium points chosen to the comparison of results correspond to the following load parameters (Fig. 14): (A) λ GBT =17.25, λ ABQ =17.20, (B) λ GBT =98.90, λ ABQ =97.00 and (C) λ GBT =112.65, λ ABQ =111.00. In the bottom part of Fig. 14, the modal participation diagram is displayed. In accordance with the modal participation definition used in this work 2 , one can conclude that modes 2 (major axis bending - Fig. 13(a)), 4 (torsion - Fig. 13(b)) and 5 (local - Fig. 13(c 1 )) are the most relevant ones (participation ≥ 5%). In the spread of plasticity phase (after yielding and before the horizontal plat-eau), it is seen that modes 2 and 4 decrease their participation while the contribution of mode 5 clearly increases. Within the horizontal plateau, the participation of modes 2, 4 and 5, increase, decrease and remain uniform, respectively.   The displacement profiles δ z (x) and δ y (x) of the node indicated in Fig. 12(a) are depicted in Fig. 15 and the comparison between GBT and ABAQUS results display an excellent agreement for points A and B. In the case of point C, GBT and ABAQUS profiles show a maximum difference in the section where the maximum absolute displacement takes place, with differences about 4.5 and 8.5%, respectively for δ z (x) and δ y (x). However, one should bear in mind that point C doesn´t correspond to the same displacement in GBT and ABAQUS (see Fig. 14).
Figs. 16(a) and 16(b) show axial and Von Mises stresses diagrams, respectively. Once again, the similarities between GBT and ABAQUS results are notorious. Fig. 17(a) and 17(b) depicts the variations of shear stresses for points A and C, respectively. It is seen that the shear stresses in the web evolve from quadratic variation (point A -elastic) to an almost uniform profile (point C -plastic). There is also a good resemblance between both GBT and ABAQUS results.
The Von Mises stress and axial displacement contours for point B are presented in Figs.18 and 19, respectively, where the results are qualitatively and quantitatively fairly good. Lastly, similar collapse deformed shapes yielded by GBT and ABAQUS are depicted in Fig. 20 and it can be seen that most of the deformation occurs close to the load at x=3000 mm.

CONCLUSIONS
A GBT formulation to analyse the first-order elastoplastic behaviour of thin-walled members was presented and applied to illustrate the behaviour of RHS and LSB fixed-pinned beams subjected to distributed and concentrated loading, respectively. The GBT results were compared with the ones obtained from ABAQUS (Simulia 2004) code using shell FE's.
Equilibrium paths, displacement profiles, stress diagrams, stress / displacement 3D contours, as well as 3D deformed shapes were compared and a very good agreement was found for all results. Finally, it is important to highlight that the number of dof used in each ABAQUS model was 65.1 (RHS beam) and 5.3 (LSB beam) times higher than in GBT.