Asymmetric thin-walled sections such as steel angles and tees are widely used in truss structures and transmission towers. To address extreme limit states that these structures encounter due to extreme events such as hurricanes and earthquakes, it is important to capture their response due to large deformations caused by static or dynamic loading. In the nonlinear large deformation regime, these members have coupled axial -flexural -torsional deformation due to the so-called Wagner effect and the noncoincident shear center and centroid. A three -dimensi onal corotational total Lagrangian beam element is formulated and implemented in the OpenSees corotational framework to account for these coupling effects by invoking Green-Lagrang e strain s referenced to a basic system. In the basic system, shear forces an d torque are defined with respect to the shear center, axial force is referred to the centroid, and flexure is defined around the section principle axes but in the planes containing the shear center. The element tangent stiffness matrix is derived through linearization of the governing equation obtained from the principle of virtual work. Cubic Hermitian functions for the transverse displacements and a linear shape function for the axial and torsional deformation are adopted in the development. Before conducting the corotational transformation, all element end forces and displacements are transformed to act about the shear center. In order to remedy membrane locking in the inextensional bending mode, the high order bending terms in the axial strain are replaced by a constant effective strain. Cyclic material nonlinearity is considered by discretizing the cross section into a grid of fibers, tracking the steel uniaxial stress-strain constitutive at each fiber, and performing numerical integration over the cros s section to obtain the section stiffness matrix. The formulation is compared against a set of experimental and numerical results to validate that the element can model geometric and material nonlinearities accurately and efficiently. Keywords: Asymmetric section, Corotational transformation, Displacement -based element, Green -Lagrange strain, Axial-flexural -torsional interaction , Membrane locking, Fiber section Introduction The conventional approaches developed to address the mechanics of structural beam elements are often based on the assumption that different deformation modes (axial, bending and torsion) are independent . However, for a three-dimensional structural member, this assumption is only acceptable for small deformation s and only if t he member has a doubly symmetric cross section. For nonlinear large deformatio ns, this assumption breaks down even for members with doubly symmetric sections, because different modes of deformation may be coupl ed. To consider these coupling effects accurat ely, several phenomena need to be taken into account , such as the so -called Wagner effect, and the coupling caused by the noncoincident shear center and centroid of asymmetric cross section s. In finite element analysis, these coupling effects usually are modeled through including of the geometric stiffness matrix, which is usually obtained using three approaches: a total Lagrangian formulation, an updated Lagrangian formulation and a corotational formulation [1-5]. In the total Lagrangian formulation, the initial undeformed configuration is selected as the reference state, whereas in the updated Lagrangian formulation, the last calculated configuration is adopted for the same purpose [6]. In the corotational formulation, a basic coordinate system is lined up with each element chord and continuously translates and rotates with the element as the deformations proceed [7]. The element is formulated in the basic system then transformed to the g lobal system through the corotational transformation. Kitipornchai and Chan [8] and Chan and Kitipornchai [9] derived element geometric stiffness matrices for angle and tee beam -columns, as well as for generic asymmetric thin -walled beam - columns. A n updated Lagrangian formulation, coupled with the arc -length technique for iterative solution at each step, was adopted to trace the nonlinear load -deformation path. The validity and accuracy of the formulation was demonstrated on different problems, for which the buckling loads and the load- deformation relationships were derived. The material was assumed to be linear elastic. In later work, they introduced a lumped plasticity approach, coupled with the concept of a yield surface in stress-resultant space, to model the material nonlinearity of thin-walled structures [10] . Lee and McClure [11] developed a three-dimensional L -section beam finite element for elastoplastic large deformation analysis. They proposed a generalized interpolation scheme for the isoparametric formulation of a three -dimensional thin -walled beam element. The updated Lagrangian formulation was adopted considering large deformation behavior. The axial, bending and shearing actions were included and the ‘mixed interpolation of tensor ial components’ (MITC) technique was used for the locking removal. In order to facilitate consideration of loading and displacement eccentricities at the finite element level, the element longitudinal reference line can be arbitrarily positioned on the bea m section depending on the location of the nodal degrees of freedom. However, the eccentricities were implemented in the element derivation, which means the users cannot specify the eccentricities by themselves. Recently, Liu et al. [12] developed an efficient beam element implementation within the educational structural analysis software MASTAN2 [13], which is capable of doing large - deformation analysis of thin-walled members with asymmetric sections and modeling the Wagner effects. The updated Lagrangian method was adopted for the large deflection analysis. Then, they improved derivations and validations for the element linear and geometric stiffness matrices [14] . This work also showed the calculation details of a refined cross -section analysis algo rithm for arbitrarily -shaped open sections. Linear elastic material is assumed throughout their work. There are other researchers who developed beam -column elements with asymmetric sections under the corotational framework . Hsiao and Lin [15] and Chen et al. [16] derived a corotational total Lagrangian element formulation for geometric nonlinear analysis of thin-walled members with monosymmetric and asymmetric cross sections and linear elastic materials. Battini and Pacoste [17, 18] proposed the formulation of three -dimensional corotational beam elements for buckling and post-buckling analysis of frame structures, where the centroid and shear center of the cross section are not necessarily coincident . More recently, Rinchen et al. [19, 20] derived a di splacement -based beam element formulation for asymmetric thin -walled members, which was implemented in the OpenSees corotational framework [21] . In the work of Rinchen et al. [19, 20] , membrane locking occurs because linear and cubic shape functions are used for axial and lateral displacements, respectively, which may cause the element to overestimate the bending stiffness when the beam section is thin or the element geometry becomes curved [4, 22] . In this paper, a displacement -based beam element is developed following the corotational total Lagrangian formulation. Rigid body motions are considered in the corotational transformation, while element deformations with geometric nonlinear effects are modeled in the basic system through the total Lagrangia n approach. The basic system and beam section kinematics used in the original OpenSees corotational framework is modified to consider the noncoincident centroid and shear center. The fiber section method with uniaxial constitutive laws for each fiber is ut ilized to model the distributed inelastic behavior. Compared to the work of Battini and Pacoste [17, 18] , this research adopts a different kinematic model, which decouples displacements and rotations in the basic system [8, 9] and follows the classical beam theory where bending is referred to the centroid and torsion is defined with respect to the shear center. Due to this kinematic model, a cross-section transformation matrix is utilized in advance of the corotational transfor mation to move all degrees - of-freedom (DOFs) to act about the member reference axis. Compared to the work of Rinchen et al. [20] , in this research membrane locking is remedied through a special treatment of the higher order terms in the Green -Lagrange stra in; in addition, the cross -section transformation matrix and the shape functions are improved. Finally, a number of members and systems are analyzed to validate the accuracy of the results from the implementation within the OpenSees software. As shown in S ection 4, fewer of the elements developed in this work are needed to obtain comparable accuracy with the element developed by others in the literature. Geometric Nonlinear Analysis Formulation The corotational transformation decomposes the motion of the el ement into rigid body movements and pure deformations, through the use of a basic coordinate system attached to the current deformed configuration as the reference system, which continuously rotates and translates with the element. Within the basic system, a geometric linear formulation, a total Lagrangian formulation or an updated Lagrangian formulation can be adopted [15, 16, 23] . OpenSees employs the corotational transformation with a geometric linear formulation in the basic system for developing nonlinear displacement -based beam elements. Consequently, the high order terms in the strain displacement relationship are neglected in the basic system, which means that the axial, flexural and torsional deformations are independent of each other in the basic s ystem. The axial -flexural interaction can be captured in the global system by the corotational transformation; however, the axial -torsional and flexural -torsional interaction cannot be recovered in the global system. Therefore, the torsional behavior is independent to axial and flexural deformations for the original displacement -based beam element in OpenSees, which leads to the inability of this element to simulate torsional and flexural -torsional buckling behavior of members. In this research, the high order terms in the strain displacement compatibility equation are included through the total Lagrangian formulation in the basic system of the corotational transformation formulation. A brief introduction of this total Lagrangian corotational formulation is presented in this section. Total Lagrangian Formulation Fig. 1 shows the concept of the total Lagrangian method applied to a 3D beam element. In the total Lagrangian f ormulation, when describing strains and displacements, the reference configuration is the initial configuration at time 0. The element is formulated in the 12 -DOF local system of the initial configuration, then transformed to the global system using the tr ansformation matrix of the initial configuration. The transformation matrix is constant for each element. There are 12 DOFs in the local system because element ends have lateral displacements and rotations. Fig. 1. Illustration of the total Lagrangian formulation Theoretically, the principle of virtual work at time t+Δt is � 𝜏𝜏𝑖𝑖𝑖𝑖 𝑡𝑡+∆𝑡𝑡𝛿𝛿𝑒𝑒𝑡𝑡+∆𝑡𝑡 𝑖𝑖𝑖𝑖 𝑉𝑉 𝑡𝑡+∆𝑡𝑡𝑑𝑑𝑑𝑑 𝑡𝑡+∆𝑡𝑡=ℛ 𝑡𝑡+∆𝑡𝑡 (1) where 𝑑𝑑 𝑡𝑡+∆𝑡𝑡 is the volume at time t+Δt; ℛ 𝑡𝑡+∆𝑡𝑡 is the virtual work done by external loads; 𝜏𝜏𝑖𝑖𝑖𝑖 𝑡𝑡+∆𝑡𝑡 is the Cauchy stress tensor at time t+Δt, and 𝛿𝛿𝑒𝑒𝑡𝑡+∆𝑡𝑡 𝑖𝑖𝑖𝑖 is the variation in the small strains referred to the configuration at time t+Δt [1, 2] . Since the configuration at time t+Δt is unknow n and directly working with increments of Cauchy stresses is unreasonable, this equation cannot be used to derive a beam element. By employing the Green -Lagrange strains and 2nd Piola -Kirchhoff stresses, the above principle of virtual work at time t+Δt can be replaced by �𝜎𝜎𝑖𝑖𝑖𝑖0𝑡𝑡+∆𝑡𝑡𝛿𝛿 𝑉𝑉 0𝜀𝜀𝑖𝑖𝑖𝑖0𝑡𝑡+∆𝑡𝑡𝑑𝑑𝑑𝑑 0=ℛ 𝑡𝑡+∆𝑡𝑡 (2) where 𝑑𝑑 0 is the volume at time 0; 𝜎𝜎𝑖𝑖𝑖𝑖0𝑡𝑡+∆𝑡𝑡 is the 2nd Piola -Kirchhoff stress tensor and 𝜀𝜀𝑖𝑖𝑖𝑖0𝑡𝑡+∆𝑡𝑡 is the Green -Lagrange strain tensor [6]. Equation (2 ) takes the initial configuration at time 0 as the reference configuration , from which the element formulations can be developed. Corotational Transformation In the corotational transformation, when describing displacement s and strains at any time point , the reference system is the basic system at th e same time point ( Fig. 2). The element is formulated in the 6-DOF basic system at this time, then transformed to the global system using the corotational transformation matrix at this time. The t ransformation matrix is different at different time points because the basic system continuously rotates and translates with the element. The deformational response is considered at the level of the basic coordinate system, whereas the rigid body motion is captured by the transformation matrices relating the basic and global systems. There are 6 DOFs in the basic system because rigid body motion is separated from the deformation and only rotations and an elongation remain in the basic system. In the corotational transformation, the 6 x6 element stiffness matrix is derived in the basic system at t+Δt, then transformed to the global system using the corotational transformation matrix of the current configuration. Fig. 2. Illustration of the corotational total Lagrangian formulation The displacements, forces and stiffnesses transformations between the basic and global systems that were implemented in OpenSees are introduced here briefl y. The tangential relation between the d isplacements 𝑫𝑫 𝑏𝑏 in the basic system and the displacements 𝑫𝑫� in the global systems can be defined as 𝛿𝛿𝑫𝑫𝑏𝑏=𝑻𝑻𝛿𝛿𝑫𝑫� (3) where T is a transformation matrix connecting the global and basic systems. In geometric linear analysis T is a constant matrix, while in geometric nonlinear analysis T is a function of displacements 𝑫𝑫�. Equating the internal virtual work in both the basic and global systems, the relationship between element end forces in the basic and global systems can be ob tained as 𝑷𝑷�=𝑻𝑻𝑇𝑇𝑷𝑷𝑏𝑏 (4) where 𝑷𝑷� and 𝑷𝑷𝑏𝑏 are the element end forces at the global system and basic system, respectively. The element tangent stiffness in the global system is obtained from the linearization of Eq. (4 ), such that 𝛿𝛿𝑷𝑷�=𝛿𝛿(𝑻𝑻𝑇𝑇𝑷𝑷𝑏𝑏)=𝑻𝑻𝑇𝑇𝛿𝛿𝑷𝑷𝑏𝑏+𝛿𝛿𝑻𝑻𝑇𝑇𝑷𝑷𝑏𝑏=𝑻𝑻𝑇𝑇𝑲𝑲𝑏𝑏𝛿𝛿𝑫𝑫𝑏𝑏+𝛿𝛿𝑻𝑻𝑇𝑇𝑷𝑷𝑏𝑏=(𝑻𝑻𝑇𝑇𝑲𝑲𝑏𝑏𝑻𝑻+𝑲𝑲𝐺𝐺)𝛿𝛿𝑫𝑫� =𝑲𝑲�𝛿𝛿𝑫𝑫� (5) where 𝑲𝑲𝑏𝑏 is the tangent stiffness matrix in the basic system, and 𝛿𝛿𝑷𝑷𝑏𝑏=𝑲𝑲𝑏𝑏𝛿𝛿𝑫𝑫𝑏𝑏 is used in the above equation. Here, 𝑲𝑲�=𝑻𝑻𝑇𝑇𝑲𝑲𝑏𝑏𝑻𝑻+𝑲𝑲𝐺𝐺 (6) is the tangent stiffness matrix in the global system. On the right hand side of Eq. (6 ), the second term 𝑲𝑲𝐺𝐺 is called the external geometric stiffness matrix. The detailed derivation of T and 𝑲𝑲𝐺𝐺 is described by Crisfield [5] and de Souza [24] . The element formulation in the basic system is independent to the corotational transformation; therefore, this research uses the current corotational transformation in OpenSees with the element formulation in the basic system modified. Corotational Total Lagrangian Formulation The corotational total Lagrangian formulation is utilized in the present research t o consider geometric nonlinearity in the basic system so that the torsional DOF will be coupled with the axial and flexural DOFs. In the corotational total Lagrangian formulati on, the reference configuration is taken as the initial undeformed configuration but translates and rotates in accordance with the motion of the corotating basic system (see Fig. 2). When deriving the stiffness matrix and force recovery in the basic system, we can assume that the basic system is fixed and apply the total Lagrangian formulation in the basic system as shown in Fig. 3. This means that the Green -Lagrange strain tensor, the 2 nd Piola -Kirchhoff stress tensor and Eq. (2 ) will be adopted to describe the element responses with respect to the reference configuration illustrated in Fig. 3. Fig. 3. Total Lagrangian method in the basic system Displacement -Based Beam Element in the Basic System This section presents the development of a new d isplacement -based beam element in OpenSees, which can be used to simulate the nonlinear behavior of members with asymmetric cross secti ons. The new displacement -based element can address inelastic behavior through the use of fiber -based cross-section formulations. Coordinate Systems of the Original OpenSees elements OpenSees has two coordinate systems: the basic system with 6 DOFs per element, and the global system with 12 DOFs per element. The basic system and global system are represented in Fig. 4. The element ends node I and node J in Fig. 4 are centroids of the cross section. The element stiffness matrix is formulated in the basic system which has the 6 DOFs: one relative axial displacement 𝑢𝑢𝐽𝐽, two rotat ions relative to the chord 𝜃𝜃𝐼𝐼𝐼𝐼 and 𝜃𝜃𝐽𝐽𝐼𝐼, about the z axis, two rotations relative to the chord 𝜃𝜃𝐼𝐼𝐼𝐼 and 𝜃𝜃𝐽𝐽𝐼𝐼, about the y axis, and one relative angle of twist 𝜙𝜙𝐽𝐽. Note that OpenSees assumes the cross section is doubly symmetric and all the 6 DOFs are defined with respect to the centroid of the cross section. These relative displacement DOFs in the basic system are the minimum number of geometric variables necessary to describe the deformation modes of the element in space, which also make sure t here are no rigid body modes in the basic system. The six statically independent end forces corresponding to these displacement DOFs are the axial force 𝑁𝑁𝐽𝐽, the two bending moments in the xy plane, 𝑀𝑀𝐼𝐼𝐼𝐼 and 𝑀𝑀𝐽𝐽𝐼𝐼, the two bending moments in the xz plane, 𝑀𝑀𝐼𝐼𝐼𝐼 and 𝑀𝑀𝐽𝐽𝐼𝐼, and the torque 𝑇𝑇𝐽𝐽. These element forces and displacements are grouped in vectors 𝑫𝑫𝑏𝑏=[𝑢𝑢𝐽𝐽𝜃𝜃𝐼𝐼𝐼𝐼𝜃𝜃𝐽𝐽𝐼𝐼𝜃𝜃𝐼𝐼𝐼𝐼𝜃𝜃𝐽𝐽𝐼𝐼𝜙𝜙𝐽𝐽]𝑇𝑇 (7) and 𝑷𝑷𝑏𝑏=[𝑁𝑁𝐽𝐽𝑀𝑀𝐼𝐼𝐼𝐼𝑀𝑀𝐽𝐽𝐼𝐼𝑀𝑀𝐼𝐼𝐼𝐼𝑀𝑀𝐽𝐽𝐼𝐼𝑇𝑇𝐽𝐽]𝑇𝑇 (8) Fig. 4. Basic and global coordinate system s in space The OpenSees element has 12 DOFs in the global system, i.e., three translational displacements and three rotations at each node. These DOFs are grouped into a vector 𝑫𝑫�=�𝑼𝑼𝐼𝐼𝑇𝑇𝜸𝜸𝐼𝐼𝑇𝑇𝑼𝑼𝐽𝐽𝑇𝑇𝜸𝜸𝐽𝐽𝑇𝑇�𝑇𝑇 (9) where 𝑼𝑼𝐼𝐼𝑇𝑇 and 𝑼𝑼𝐽𝐽𝑇𝑇 are 3x1 vectors with the translational displacements, and 𝜸𝜸𝐼𝐼𝑇𝑇 and 𝜸𝜸𝐽𝐽𝑇𝑇 are 3x1 pseudo-vectors that define the rotations of nodes I and J, respectively. In the absolute sense, these rotations can be arbitrarily large, although it is assumed that the relative rotation of two element ends is small. As such, the OpenSees element can experience finite displacements and rigid body rotations; however, the deformations along the element are considered to be moderate. The corresponding work conjugate forces in the global system are 𝑷𝑷�=[𝑃𝑃�1𝑃𝑃�2𝑃𝑃�3𝑃𝑃�4𝑃𝑃�5𝑃𝑃�6𝑃𝑃�7𝑃𝑃�8𝑃𝑃�9𝑃𝑃�10𝑃𝑃�11𝑃𝑃�12]𝑇𝑇 (10) Coordinate Systems for Asymmetric Sections Since the new element is implemented in the OpenSees corotational framework, the original OpenSees coordinate systems described in Section 3. 1 are adopted with some modifications for the basic coordinate system. The same global system with the original OpenSees element is used here; however, due to the noncoincident shear center and centroid of the beam element with asymmetric sections, its basic system is slightly different as shown in Fig. 5. The modified basic system is defined by two set of coordinates: x , y, z and x, y, z. The coordinate system x , y, z is chosen such that x passes through the end cross-section centroids C and C’, and y and z are the section principal axes. A parallel set of coordinates x , y, z is chosen such that x passes through the end cross-section shear centers S and S’, and y and z are parallel to the p rincipal y and z axes of the cross section. Let v and w denote the displacements of the shear center in the y and z directions, u the axial displacement along the centroidal axis CC’, and ϕ the angle of twist about the shear center axis SS’. Such a coordinate system will uncouple displacements and rotations in the sense of the first order effect [8, 9] . Fig. 5. Basic and global coordinate systems in space for asymmetric sections The element stiffness matrix is formulated in the basic system with the 6 DOFs presented in Eqs. (7) and (8). However, the definitions of the DOFs are different as follows: one relative axial displacement 𝑢𝑢𝐽𝐽 of the centroids , two rotations relative to the chord 𝜃𝜃𝐼𝐼𝐼𝐼 and 𝜃𝜃𝐽𝐽𝐼𝐼, about the z axis, two rotations relative to the chord 𝜃𝜃𝐼𝐼𝐼𝐼 and 𝜃𝜃𝐽𝐽𝐼𝐼, about the y axis, and one relative angle of twist 𝜙𝜙𝐽𝐽 about the x axis. The six statically independent end forces related to these displacements are: the axial force 𝑁𝑁𝐽𝐽 acting along the centroidal axis CC’ ; the two bending moments acting about the principle axis z and in the xy plane, 𝑀𝑀𝐼𝐼𝐼𝐼 and 𝑀𝑀𝐽𝐽𝐼𝐼; the two bending moments acting about the principle ax is y and in the xz plane, 𝑀𝑀𝐼𝐼𝐼𝐼 and 𝑀𝑀𝐽𝐽𝐼𝐼; and the torsional moment 𝑇𝑇𝐽𝐽 about the shear center axis SS’ [25] . Following the corotational transformation procedure, the forces, displacements and stiffness matrix need to be transformed to the global system. However, in the current definition of the basic system, some DOFs are defined with respect to the centroi d, while others are defined with respect to the shear center, which makes the corotational transformation cannot be applied directly [20]. Consequently, all DOFs (end forces and displacements) need to be transformed to one reference point in advance of the corotational transformation. Note that in the basic system, the axial force is referred to the centroid, the shear forces and torque are defined with respect to the shear center, and the moments act in the planes containing the shear center. Thus, it is straightforward to transform all DOFs to the shear center, because only the axial force needs to be transformed. If other points (e.g., the centroid) are selected as the reference point, the shear forces need to be transformed, which is more d ifficult because the shear forces are only recovered after the corotational transformation. Therefore, it is proper to choose the shear center as the reference point and the shear center axis as the member reference axis before the corotational transformat ion. The coordinate system x , y, z is referred as the element basic reference system with the element end forces 𝑷𝑷 𝑟𝑟 and displacements 𝑫𝑫𝑟𝑟 act through the shear center. The following equations are used to transform all DOFs to the shear center 𝑷𝑷𝑟𝑟=𝑻𝑻𝑟𝑟𝑇𝑇𝑷𝑷𝑏𝑏 (11) 𝑫𝑫𝑏𝑏=𝑻𝑻𝑟𝑟𝑫𝑫𝑟𝑟 (12) where the cross-section transformation matrix is 𝑻𝑻𝑟𝑟𝑇𝑇= ⎣⎢⎢⎢⎢⎡1 00000 𝑦𝑦𝑠𝑠 10000 −𝑦𝑦𝑠𝑠01000 −𝑧𝑧𝑠𝑠00100 𝑧𝑧𝑠𝑠 00010 0 00001⎦⎥⎥⎥⎥⎤ (13) and ys and zs are coordinates of the shear center relative to the centroid. Consequently, the stiffness matrix in the element basic reference system is 𝑲𝑲𝑟𝑟=𝑻𝑻𝑟𝑟𝑇𝑇𝑲𝑲𝑏𝑏𝑻𝑻𝑟𝑟 (14) With this consideration, the coordinate system x , y, z in Fig. 5 is the same with the basic system in Fig. 4 in the sense of corotational transformation, but the displacement and force DOFs should be replaced by 𝑫𝑫𝑟𝑟 and 𝑷𝑷𝑟𝑟. The corotational transformation will remain the same as in the original OpenSees with the Eqs. (3 ), (4) and (6) being modified as 𝛿𝛿𝑫𝑫𝑟𝑟=𝑻𝑻𝛿𝛿𝑫𝑫� (15) 𝑷𝑷�=𝑻𝑻𝑇𝑇𝑷𝑷𝑟𝑟 (16) 𝑲𝑲�=𝑻𝑻𝑇𝑇𝑲𝑲𝑟𝑟𝑻𝑻+𝑲𝑲𝐺𝐺 (17) Beam Section Kinematics With the kinematic assumption of the Euler -Bernoulli beam theory and neglecting the in -plane distortion of the section, the motion of a material point P ( x, y, z) (P is expressed in the coordinate system x, y, z shown in Fig. 5) in the beam section is described in terms of the displacement components [26, 27] 𝑢𝑢𝑝𝑝=𝑢𝑢−𝑦𝑦𝑣𝑣′−𝑧𝑧𝑤𝑤′+𝑧𝑧𝜙𝜙𝑣𝑣′−𝑦𝑦𝜙𝜙𝑤𝑤′ 𝑣𝑣𝑝𝑝=𝑣𝑣−𝜙𝜙(𝑧𝑧−𝑧𝑧𝑠𝑠) 𝑤𝑤𝑝𝑝=𝑤𝑤+𝜙𝜙(𝑦𝑦−𝑦𝑦𝑠𝑠) (18) where ys and zs are defined in Eq. (13); in addition, u, v, w and ϕ are defined in Section 3.2. The longitudinal/axial component of the Green-Lagrange strain is expressed as 𝜀𝜀=𝜕𝜕𝑢𝑢𝑝𝑝 𝜕𝜕𝜕𝜕+1 2�𝜕𝜕𝑢𝑢𝑝𝑝 𝜕𝜕𝜕𝜕�2 +12�𝜕𝜕𝑣𝑣𝑝𝑝 𝜕𝜕𝜕𝜕�2 +12�𝜕𝜕𝑤𝑤𝑝𝑝 𝜕𝜕𝜕𝜕�2 (19) The high-order term 1 2�𝜕𝜕𝑢𝑢𝑝𝑝/𝜕𝜕𝜕𝜕�2 is negligible compared to 𝜕𝜕𝑢𝑢𝑝𝑝/𝜕𝜕𝜕𝜕 because it is assumed that the term 𝜕𝜕𝑢𝑢𝑝𝑝/𝜕𝜕𝜕𝜕 is small compared to unity. Hence, the Green -Lagrange strain becomes 𝜀𝜀=𝜕𝜕𝑢𝑢𝑝𝑝 𝜕𝜕𝜕𝜕+1 2�𝜕𝜕𝑣𝑣𝑝𝑝 𝜕𝜕𝜕𝜕�2 +12�𝜕𝜕𝑤𝑤𝑝𝑝 𝜕𝜕𝜕𝜕�2 (20) Taking derivatives of the displacement field in Eq. (18) with respect to x and substituting the results into Eq. (20) gives 𝜀𝜀=𝑢𝑢′−𝑦𝑦𝑣𝑣′′−𝑧𝑧𝑤𝑤′′+1 2[(𝑣𝑣′)2+(𝑤𝑤′)2]+1 2[(𝑦𝑦−𝑦𝑦𝑠𝑠)2+(𝑧𝑧−𝑧𝑧𝑠𝑠)2](𝜙𝜙′)2 +(𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝜙𝜙′+(𝑧𝑧𝑣𝑣′′−𝑦𝑦𝑤𝑤′′)𝜙𝜙 (21) Compared with the strain term 𝜀𝜀=𝑢𝑢′−𝑦𝑦𝑣𝑣′′−𝑧𝑧𝑤𝑤′′ used for the original OpenSees displacement - based element, Eq. (21 ) has the following extra terms 1 2[(𝑣𝑣′)2+(𝑤𝑤′)2] : geometric nonlinear term (coupling between the axial tension/compression and bending) 1 2[(𝑦𝑦−𝑦𝑦𝑠𝑠)2+(𝑧𝑧−𝑧𝑧𝑠𝑠)2](𝜙𝜙′)2 : Wagner term (coupling between the axial strain and torsion) (𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝜙𝜙′+(𝑧𝑧𝑣𝑣′′−𝑦𝑦𝑤𝑤′′)𝜙𝜙 : coupling term between bending and torsion (effect of this coupling on the axial strain) The shear strain at P resulting from uniform twisting 𝜙𝜙′ of a thin -walled open section member is approximated by [28] 𝛾𝛾=2𝑛𝑛𝜙𝜙′ (22) in which n is the perpendicular distance of point P ( x, y, z) from the mid -thickness line of the cross section. Shear strains due to bending are neglected. Displacement -Based Element Formulation In this section, the element governing equilibrium equation is derived from the principle of virtual work with the aid of the kinematic assumptions and the displacement shape functions. Linearization of the governing equation leads to the element tangent stiffness and the equation for force recovery. To start, the axial strai n and shear strain can be put together in vector form as 𝝐𝝐=�𝜀𝜀 𝛾𝛾� (23) The terms in the strain vector can be expressed in matrix form as 𝝐𝝐=𝒀𝒀𝒀𝒀 (24) where 𝒀𝒀=�1−𝑦𝑦𝑧𝑧 0 0 0 (𝑦𝑦−𝑦𝑦𝑠𝑠)2+(𝑧𝑧−𝑧𝑧𝑠𝑠)20 0 2𝑛𝑛� (25) and 𝒀𝒀= ⎣⎢⎢⎢⎢⎢⎡𝑢𝑢′+1 2[(𝑣𝑣′)2+(𝑤𝑤′)2]+(𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝜙𝜙′ 𝑣𝑣′′+𝑤𝑤′′𝜙𝜙 −𝑤𝑤′′+𝑣𝑣′′𝜙𝜙 1 2(𝜙𝜙′)2 𝜙𝜙′⎦⎥⎥⎥⎥⎥⎤ (26) Vector d is the section deformation vector. Matrix Y relates strains on a material point and the section deformations on the corresponding section. The variation of the strain vector is obtained as 𝛿𝛿𝝐𝝐=𝒀𝒀𝛿𝛿𝒀𝒀 (27) in which 𝛿𝛿𝒀𝒀= ⎣⎢⎢⎢⎡𝛿𝛿𝑢𝑢′+𝑣𝑣′𝛿𝛿𝑣𝑣′+𝑤𝑤′𝛿𝛿𝑤𝑤′+(𝑧𝑧𝑠𝑠𝛿𝛿𝑣𝑣′−𝑦𝑦𝑠𝑠𝛿𝛿𝑤𝑤′)𝜙𝜙′+(𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝛿𝛿𝜙𝜙′ 𝛿𝛿𝑣𝑣′′+𝜙𝜙𝛿𝛿𝑤𝑤′′+𝑤𝑤′′𝛿𝛿𝜙𝜙 −𝛿𝛿𝑤𝑤′′+𝜙𝜙𝛿𝛿𝑣𝑣′′+𝑣𝑣′′𝛿𝛿𝜙𝜙 𝜙𝜙′𝛿𝛿𝜙𝜙′ 𝛿𝛿𝜙𝜙′⎦⎥⎥⎥⎤ (28) And 𝛿𝛿𝒀𝒀 can be expressed as 𝛿𝛿𝒀𝒀=𝑵𝑵𝛿𝛿𝛿𝛿1𝛿𝛿𝒗𝒗 (29) where 𝑵𝑵𝛿𝛿𝛿𝛿1= ⎣⎢⎢⎢⎡1𝑣𝑣′+𝑧𝑧𝑠𝑠𝜙𝜙′𝑤𝑤′−𝑦𝑦𝑠𝑠𝜙𝜙′ 0 0 0 0 0 00 0 0 1𝜙𝜙𝑤𝑤′′ 𝜙𝜙−1𝑣𝑣′′𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′ 0 0 0 0 0 0 0 00 0 0 0 0 0 𝜙𝜙′ 1⎦⎥⎥⎥⎤ (30) 𝛿𝛿𝒗𝒗=[𝛿𝛿𝑢𝑢′𝛿𝛿𝑣𝑣′𝛿𝛿𝑤𝑤′𝛿𝛿𝑣𝑣′′𝛿𝛿𝑤𝑤′′𝛿𝛿𝜙𝜙𝛿𝛿𝜙𝜙′]𝑻𝑻 (31) Thus, Eq. ( 27) becomes 𝛿𝛿𝝐𝝐=𝒀𝒀𝑵𝑵𝛿𝛿𝛿𝛿1𝛿𝛿𝒗𝒗 (32) In the basic system, the axial elongation field, the transverse displacement field and the twist displacement field along the element can be interpolated between the element DOFs by the help of shape functions 𝑢𝑢=𝑵𝑵𝑢𝑢𝑇𝑇𝑫𝑫𝑏𝑏=[𝑁𝑁𝑢𝑢100000]𝑫𝑫𝑏𝑏 𝑣𝑣=𝑵𝑵𝑣𝑣𝑇𝑇𝑫𝑫𝑏𝑏=[0𝑁𝑁𝑣𝑣1𝑁𝑁𝑣𝑣2000]𝑫𝑫𝑏𝑏 𝑤𝑤=𝑵𝑵𝑤𝑤𝑇𝑇𝑫𝑫𝑏𝑏=[000𝑁𝑁𝑤𝑤1𝑁𝑁𝑤𝑤20]𝑫𝑫𝑏𝑏 𝜙𝜙=𝑵𝑵𝜙𝜙𝑇𝑇𝑫𝑫𝑏𝑏=[00000𝑁𝑁𝜙𝜙1]𝑫𝑫𝑏𝑏 (33) The shape functions in Eq. (33) are chosen as follows: cubic Hermitian functions for the transverse displacements, and a linear function for the axial deformation and twist. Specifically, the shape functions utilized in the current development are 𝑁𝑁𝑢𝑢1=𝑁𝑁𝜙𝜙1=𝜕𝜕 𝑙𝑙0 𝑁𝑁𝑣𝑣1=−𝑁𝑁𝑤𝑤1=𝜕𝜕�1−𝜕𝜕 𝑙𝑙0�2 𝑁𝑁𝑣𝑣2=−𝑁𝑁𝑤𝑤2=𝜕𝜕�𝜕𝜕 𝑙𝑙0��𝜕𝜕 𝑙𝑙0−1� (34) where 𝑙𝑙0 is the length of the undeformed element. With the help of shape functions, 𝛿𝛿𝒗𝒗 can be written in terms of the element DOFs, as shown 𝛿𝛿𝒗𝒗=𝑵𝑵𝛿𝛿𝛿𝛿2𝛿𝛿𝑫𝑫𝑏𝑏 (35) in which 𝑵𝑵𝛿𝛿𝛿𝛿2= ⎣⎢⎢⎢⎢⎢⎢⎡𝑁𝑁𝑢𝑢1′0 0 0 0 0 0𝑁𝑁𝑣𝑣1′𝑁𝑁𝑣𝑣2′0 0 0 0 0 0𝑁𝑁𝑤𝑤1′𝑁𝑁𝑤𝑤2′0 0𝑁𝑁𝑣𝑣1′′𝑁𝑁𝑣𝑣2′′0 0 0 0 0 0𝑁𝑁𝑤𝑤1′′𝑁𝑁𝑤𝑤2′′0 0 0 0 0 0𝑁𝑁𝜙𝜙1 0 0 0 0 0𝑁𝑁𝜙𝜙1′⎦⎥⎥⎥⎥⎥⎥⎤ (36) Therefore, Eq. (32) becomes 𝛿𝛿𝝐𝝐=𝒀𝒀𝑵𝑵𝛿𝛿𝛿𝛿1𝑵𝑵𝛿𝛿𝛿𝛿2𝛿𝛿𝑫𝑫𝑏𝑏 (37) Here, the Total Lagrangian formulation is used for the virtual work principle. In the basic system, the stresses and strains are referred to the undeformed configuration of the beam element. The principle of virtual work can be written in a general form as �𝛿𝛿𝝐𝝐𝑇𝑇𝝈𝝈 𝑉𝑉0𝑑𝑑𝑑𝑑−𝛿𝛿𝑫𝑫𝑏𝑏𝑇𝑇𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡=0 (38) where 𝝈𝝈=[𝜎𝜎𝜏𝜏]𝑇𝑇 is the corresponding stress vector, 𝑑𝑑0 is the volume of the undeformed element, and 𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡 is a vector of external forces. After substituting Eq. (37) into Eq. (38) and rearranging, we have 𝛿𝛿𝑫𝑫𝑏𝑏𝑇𝑇��𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇��𝒀𝒀𝑇𝑇𝝈𝝈𝑑𝑑𝑑𝑑 𝐴𝐴0� 𝑙𝑙0𝑑𝑑𝜕𝜕−𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡�=0 (39) and 𝛿𝛿𝑫𝑫𝑏𝑏𝑇𝑇��𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑺𝑺 𝑙𝑙0𝑑𝑑𝜕𝜕−𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡�=0 (40) where 𝑑𝑑0 is the cross-section area of the undeformed element; and S is the vector of section forces expressed as 𝑺𝑺=�𝒀𝒀𝑇𝑇𝝈𝝈𝑑𝑑𝑑𝑑 𝐴𝐴0= ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ �𝜎𝜎𝑑𝑑𝑑𝑑 𝐴𝐴0 −�𝑦𝑦𝜎𝜎𝑑𝑑𝑑𝑑 𝐴𝐴0 �𝑧𝑧𝜎𝜎𝑑𝑑𝑑𝑑 𝐴𝐴0 �[(𝑦𝑦−𝑦𝑦𝑠𝑠)2+(𝑧𝑧−𝑧𝑧𝑠𝑠)2]𝜎𝜎𝑑𝑑𝑑𝑑 𝐴𝐴0 �2𝑛𝑛𝜏𝜏𝑑𝑑𝑑𝑑 𝐴𝐴0 ⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ = ⎣⎢⎢⎢⎡𝑁𝑁 𝑀𝑀𝐼𝐼 𝑀𝑀𝐼𝐼 𝑊𝑊 𝑇𝑇⎦⎥⎥⎥⎤ (41) where N is the axial force, 𝑀𝑀𝐼𝐼 is the bending moment about the z -axis, 𝑀𝑀𝐼𝐼 is the bending moment about the y -axis, W is the Wagner stress resultant, and T is the St. Venant torque. The governing equilibrium equation can be rewritten as 𝒈𝒈=�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑺𝑺 𝑙𝑙0𝑑𝑑𝜕𝜕−𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡=𝟎𝟎 (42) The governing equation is in general a nonlinear function of the unknown displacement increment. It should be linearized to obtain the element tangent stiffness matrix. Such a linearization process can be achieved by exp anding the governing equation at ( i+1)th iteration based on previous configuration at ith iteration. The details of the linearization process are shown in Appendix A, from which the following iteration equation can be obtained 𝑲𝑲𝑏𝑏∆𝑫𝑫𝑏𝑏=𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖+1−𝑷𝑷𝑖𝑖𝑖𝑖𝑡𝑡𝑖𝑖 (43) where 𝑲𝑲𝑏𝑏=�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑮𝑮𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕+�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑲𝑲𝑠𝑠𝑵𝑵𝛿𝛿𝛿𝛿1 𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕 (44) 𝑷𝑷𝑖𝑖𝑖𝑖𝑡𝑡𝑖𝑖=�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑺𝑺𝑖𝑖 𝑙𝑙0𝑑𝑑𝜕𝜕 (45) 𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖+1=𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖+∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡 (46) 𝑲𝑲𝑏𝑏 is the element tangent stiffness matrix in the basic system, the first integral in which is the element internal geometric stiffness matrix and the second integral is the element material stiffness matrix. The definition of th e auxiliary matrix G and the section stiffness matrix 𝑲𝑲𝑠𝑠 can be found in Appendix A. In addition, 𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖+1 and 𝑷𝑷𝑖𝑖𝑖𝑖𝑡𝑡𝑖𝑖 are the external load vector in ( i+1)th iteration and the element internal resisting forces in ith iteration, respectively. The vector 𝑺𝑺𝑖𝑖 represents section forces in ith iteration. Fiber Section and Material Nonlinearity In order to consider material nonlinearity, the beam section is discretized into a set of fibers so that the section stiffness matrix 𝑲𝑲 𝑠𝑠 can be obtained through numerical integration. In this element, it is assumed that the uniform torsion behavior is linear elastic so that shear strain is always elastic, and GJ is used as the rigidity of uniform torsion. Thus, uniaxial constitutive law is used to check the yielding of the material. Note that the tangent Young’s modulus E can be updated for each fiber as it reaches plastic stage. The calculation of 𝑲𝑲 𝑠𝑠 through numerical integration is shown as 𝑲𝑲𝑠𝑠= ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡�𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1−�𝑦𝑦𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑧𝑧𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑝𝑝𝑖𝑖2𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=10 −�𝑦𝑦𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑦𝑦𝑖𝑖2𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1−�𝑦𝑦𝑖𝑖𝑧𝑧𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1−�𝑝𝑝𝑖𝑖2𝑦𝑦𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=10 �𝑧𝑧𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1−�𝑦𝑦𝑖𝑖𝑧𝑧𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑧𝑧𝑖𝑖2𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑝𝑝𝑖𝑖2𝑧𝑧𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=10 �𝑝𝑝𝑖𝑖2𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1−�𝑝𝑝𝑖𝑖2𝑦𝑦𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑝𝑝𝑖𝑖2𝑧𝑧𝑖𝑖𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1�𝑝𝑝𝑖𝑖4𝐸𝐸𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=10 0 0 0 0 𝐺𝐺𝐺𝐺⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ (47) where 𝐸𝐸𝑖𝑖, 𝑑𝑑𝑖𝑖, 𝑦𝑦𝑖𝑖, and 𝑧𝑧𝑖𝑖 and 𝑝𝑝𝑖𝑖 are the tangent Young’s modulus, area, y coordinate, z coordinate and p value ( see Eq. (A.8) in Appendix A for definition ) of the jth fiber, respectively; 𝑛𝑛1 is the number of fibers on the section. Similarly, the section forces can be calculated numerically as shown in Eq. (48). Here the torque due to uniform torsion is achieved directly from the first derivative of the twist rotation 𝜙𝜙′ and the equivalent section torsional rigidity GJ because of the assumption of linear elastic St. Venant torsional behavior. Note that the nonli near torsional behavior is considered through including the Wanger stress resultant W [29]. 𝑺𝑺= ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡�𝜎𝜎𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1 −�𝑦𝑦𝑖𝑖𝜎𝜎𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1 �𝑧𝑧𝑖𝑖𝜎𝜎𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1 ���𝑦𝑦𝑖𝑖−𝑦𝑦𝑠𝑠�2+�𝑧𝑧𝑖𝑖−𝑧𝑧𝑠𝑠�2�𝜎𝜎𝑖𝑖𝑑𝑑𝑖𝑖𝑖𝑖1 𝑖𝑖=1 𝐺𝐺𝐺𝐺𝜙𝜙′⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ (48) Thus, the matrix G in Eq. ( 44) can be obtained from results of the numerical calculation of section forces S. Numerical integration methods like Gaussian quadrature are used to calculate 𝑲𝑲𝑏𝑏 and 𝑷𝑷𝑖𝑖𝑖𝑖𝑡𝑡𝑖𝑖. Membrane Locking Membrane locking results from the inability of beam elements to represent inextensional modes of deformation [22]. For the formulation of displacement -based element in the previous section, membrane strain can appear in an inextensional bending mode, which results in underprediction of displacements and strains. As recommended by Crisfield [4, 5, 26] , in order to remedy the problem of membrane locking, the high order term 1 2[(𝑣𝑣′)2+(𝑤𝑤′)2] in Eq. (21) can be replaced by an effective membrane strain 1 𝑙𝑙0∫12[(𝑣𝑣′)2+(𝑤𝑤′)2]𝑙𝑙0 0𝑑𝑑𝜕𝜕. This integral can be obtained by assuming cubic shape functions for v and w; thus, Eq. (21) should be substituted by 𝜀𝜀=𝑢𝑢′−𝑦𝑦𝑣𝑣′′−𝑧𝑧𝑤𝑤′′+1 60𝜽𝜽𝑏𝑏𝑇𝑇𝑿𝑿𝜽𝜽𝑏𝑏+1 2[(𝑦𝑦−𝑦𝑦𝑠𝑠)2+(𝑧𝑧−𝑧𝑧𝑠𝑠)2](𝜙𝜙′)2 +(𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝜙𝜙′+(𝑧𝑧𝑣𝑣′′−𝑦𝑦𝑤𝑤′′)𝜙𝜙 (49) in which the vector 𝜽𝜽𝑏𝑏=[𝜃𝜃𝐼𝐼𝐼𝐼𝜃𝜃𝐼𝐼𝐼𝐼𝜃𝜃𝐽𝐽𝐼𝐼𝜃𝜃𝐽𝐽𝐼𝐼]𝑇𝑇 holds the local slopes in the basic system, and the matrix X is defined as 𝑿𝑿=�4 0−1 0 0 4 0−1 −1 0 4 0 0−1 0 4� (50) Thus, the modified d is 𝒀𝒀= ⎣⎢⎢⎢⎢⎢⎡𝑢𝑢′+1 60𝜽𝜽𝑏𝑏𝑇𝑇𝑿𝑿𝜽𝜽𝑏𝑏+(𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝜙𝜙′ 𝑣𝑣′′+𝑤𝑤′′𝜙𝜙 −𝑤𝑤′′+𝑣𝑣′′𝜙𝜙 1 2(𝜙𝜙′)2 𝜙𝜙′⎦⎥⎥⎥⎥⎥⎤ (51) The variation of strain vector is obtained as 𝛿𝛿𝝐𝝐=𝒀𝒀𝛿𝛿𝒀𝒀 (52) in which 𝛿𝛿𝒀𝒀= ⎣⎢⎢⎢⎢⎡𝛿𝛿𝑢𝑢′+1 30𝜽𝜽𝑏𝑏𝑇𝑇𝑿𝑿𝛿𝛿𝜽𝜽𝑏𝑏+(𝑧𝑧𝑠𝑠𝛿𝛿𝑣𝑣′−𝑦𝑦𝑠𝑠𝛿𝛿𝑤𝑤′)𝜙𝜙′+(𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′)𝛿𝛿𝜙𝜙′ 𝛿𝛿𝑣𝑣′′+𝜙𝜙𝛿𝛿𝑤𝑤′′+𝑤𝑤′′𝛿𝛿𝜙𝜙 −𝛿𝛿𝑤𝑤′′+𝜙𝜙𝛿𝛿𝑣𝑣′′+𝑣𝑣′′𝛿𝛿𝜙𝜙 𝜙𝜙′𝛿𝛿𝜙𝜙′ 𝛿𝛿𝜙𝜙′⎦⎥⎥⎥⎥⎤ (53) And 𝛿𝛿𝒀𝒀 can be expressed as 𝛿𝛿𝒀𝒀=𝑵𝑵𝛿𝛿𝛿𝛿1𝛿𝛿𝒗𝒗 (54) where 𝑵𝑵𝛿𝛿𝛿𝛿1= ⎣⎢⎢⎢⎢⎡14𝜃𝜃𝐼𝐼𝐼𝐼−𝜃𝜃𝐽𝐽𝐼𝐼 304𝜃𝜃𝐼𝐼𝐼𝐼−𝜃𝜃𝐽𝐽𝐼𝐼 30−𝜃𝜃𝐼𝐼𝐼𝐼+4𝜃𝜃𝐽𝐽𝐼𝐼 30−𝜃𝜃𝐼𝐼𝐼𝐼+4𝜃𝜃𝐽𝐽𝐼𝐼 30𝑧𝑧𝑠𝑠𝜙𝜙′−𝑦𝑦𝑠𝑠𝜙𝜙′0 0 0𝑧𝑧𝑠𝑠𝑣𝑣′−𝑦𝑦𝑠𝑠𝑤𝑤′ 0 0 0 0 0 0 0 1𝜙𝜙𝑤𝑤′′0 0 0 0 0 0 0 0𝜙𝜙−1𝑣𝑣′′0 0 0 0 0 0 0 0 0 0 0 𝜙𝜙′ 0 0 0 0 0 0 0 0 0 0 1⎦⎥⎥⎥⎥⎤ (55) 𝛿𝛿𝒗𝒗=[𝛿𝛿𝑢𝑢′𝛿𝛿𝜃𝜃𝐼𝐼𝐼𝐼𝛿𝛿𝜃𝜃𝐼𝐼𝐼𝐼𝛿𝛿𝜃𝜃𝐽𝐽𝐼𝐼𝛿𝛿𝜃𝜃𝐽𝐽𝐼𝐼𝛿𝛿𝑣𝑣′𝛿𝛿𝑤𝑤′𝛿𝛿𝑣𝑣′′𝛿𝛿𝑤𝑤′′𝛿𝛿𝜙𝜙𝛿𝛿𝜙𝜙′]𝑻𝑻 (56) So 𝛿𝛿𝝐𝝐=𝒀𝒀𝑵𝑵𝛿𝛿𝛿𝛿1𝛿𝛿𝒗𝒗 (57) with the aid of shape functions, 𝛿𝛿𝒗𝒗 can be written in terms of the element DOFs 𝛿𝛿𝒗𝒗=𝑵𝑵𝛿𝛿𝛿𝛿2𝛿𝛿𝑫𝑫𝑏𝑏 (58) in which 𝑵𝑵𝛿𝛿𝛿𝛿2 is rederived as 𝑵𝑵𝛿𝛿𝛿𝛿2= ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡𝑁𝑁𝑢𝑢1′0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0𝑁𝑁𝑣𝑣1′𝑁𝑁𝑣𝑣2′0 0 0 0 0 0𝑁𝑁𝑤𝑤1′𝑁𝑁𝑤𝑤2′0 0𝑁𝑁𝑣𝑣1′′𝑁𝑁𝑣𝑣2′′0 0 0 0 0 0𝑁𝑁𝑤𝑤1′′𝑁𝑁𝑤𝑤2′′0 0 0 0 0 0𝑁𝑁𝜙𝜙1 0 0 0 0 0𝑁𝑁𝜙𝜙1′⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ (59) Therefore, Eq. (57) becomes 𝛿𝛿𝝐𝝐=𝒀𝒀𝑵𝑵𝛿𝛿𝛿𝛿1𝑵𝑵𝛿𝛿𝛿𝛿2𝛿𝛿𝑫𝑫𝑏𝑏 (60) The same linearization process for Eq. (42) is adopted as in Appendix A with the new matrices 𝑵𝑵𝛿𝛿𝛿𝛿1 and 𝑵𝑵𝛿𝛿𝛿𝛿2 in the present section. The same element formulation (Eqs. (43 )-(46)) is obtained and the only change is for the matrix G that takes the following form 𝑮𝑮= ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡0 0 0 0 0 0 0 0 0 0 0 04 30𝑁𝑁 0−1 30𝑁𝑁 0 0 0 0 0 0 0 0 04 30𝑁𝑁 0−1 30𝑁𝑁 0 0 0 0 0 0 0−1 30𝑁𝑁 04 30𝑁𝑁 0 0 0 0 0 0 0 0 0−1 30𝑁𝑁 04 30𝑁𝑁 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0𝑁𝑁𝑧𝑧𝑠𝑠 0 0 0 0 0 0 0 0 0 0−𝑁𝑁𝑦𝑦𝑠𝑠 0 0 0 0 0 0 0 0 0𝑀𝑀𝐼𝐼 0 0 0 0 0 0 0 0 0 0𝑀𝑀𝐼𝐼 0 0 0 0 0 0 0 0𝑀𝑀𝐼𝐼𝑀𝑀𝐼𝐼 0 0 0 0 0 0 0𝑁𝑁𝑧𝑧𝑠𝑠−𝑁𝑁𝑦𝑦𝑠𝑠 0 0 0𝑊𝑊⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ (61) Rigid Offsets Generally, thin -walled asymmetric members are connected eccentrically and sometimes structural joints should be considered as being of a finite size. Therefore, the element reference ends should be offset to the conne cting point from the shear center, as shown in Fig. 6. Nodes I and J are shear centers, while nodes 1 and 2 are the connecting points of the element (element reference ends). Here, dXI, dYI and dZI are the offsets specified with respect to the global coordinate system for element end node I, and they are designated to be the coordinates of node I minus the coordinates of node 1. Similarly, d XJ, dYJ and dZJ are the coordi nates of node J minus the coordinates of node 2. Fig. 6. Element with rigid offsets This phenomenon can be considered systematically by utilizing a transformation matrix to model the rigid offsets [12, 14, 30] . The relationship between end forces of the element reference ends (1 and 2 in Fig. 6) and the shear centers ( I and J in Fig. 6) can be derived with the help of equilibriums of rigid offsets. The details of this transformation are shown in Appendix B. Element State Determination The following algorithm is utilized for element state determination of the displacement -based element. Once the incremental displacements and rotations are retrieved from the glob al solution, the element state determination process is as follows 1. Update the nodal translational displacements and rotations in the global system. 2. Compute the rotations and the axial displacement in the basic system through corotational transformation. 3. Compute axial strain at each fiber according to Eq. ( 24), and then update the stresses and fiber stiffnesses through constitutive models. 4. Compute section stiffnesses and forces according to Eqs. (47) and (48) using the updated stresses and fiber stiffnesses through numerical integration across each section. Note that the torsional rigidity and the St. Venant torque are calculate separately from the section numerical integration. 5. Compute the element tangent stiffness matr ix and the internal force vector in the basic system according to Eqs. ( 44) and (45). 6. Compute the global tangent stiffness matrix and global internal force vector through corotational transformation according to Eqs. ( 4) and (6). The transformation in Appendix B should be applied if rigid offsets need to be included. 7. Conduct the convergence test. If the residual is smaller than the threshold, then commence the next incr ement. Otherwise go back to step 1 and continue the iteration. Validation of Formulation Ten examples are shown here to validate the accuracy of this new element. In examples 1 -4 and 6- 8, the material is assumed to be linear elastic. Examples 5, 9 and 10 address both elastic and inelastic materials. These validations show that the new element can simulate both geometric and material nonlinearity accurately. In the following, “DBxx” means using xx number of the new displacement -based elements with membrane l ocking remedied; while “DBxx -with locking” means using xx number of the new displacement -based element elements with membrane locking not remedied . Inextensional Bending This example compares the simulations of inextensional bending of a cantilever (see Fig. 7) using the original OpenSees displacement -based element (which is appropriate only for symmetric cross sections) [21], the new element with membrane locking, and th e new element without membrane locking. Through this example, the accuracy of the element without membrane locking is validated. When a beam element cannot bend without stretching, the energy is incorrectly shifted to membrane energy, resulting in underpre diction of displacements and strains. In this structure, a concentrated lateral force P in the global Z direction is applied to the tip of the cantilever. Fig. 7. Inextensional bending of a cantilever Fig . 8 compares the force-displacement relationship using the original OpenSees element and the new element with membrane locking. With the increase of elements, both results converge to the same answer. However, the new element with membrane locking converges too slowly compared with the original OpenSees element. As seen in Fig. 8, 5 original OpenSees elements and 20 new elements with membrane locking can have comparable accuracy. Fig. 9 shows the results of the new element with membrane locking remedied. It can be seen that when membrane locking is alleviated, the new element can give more accurate answers compared with the o riginal OpenSees element. Five new elements without membrane locking are enough to simulate this phenomenon. Fig. 8. Comparison of original OpenSees element and new element with membrane locking Fig. 9. Com parison of original OpenSees element and new element with membrane locking remedied Nonlinear Torsion This example considers a cantilever beam of rectangular cross section (see Fig. 10) subjected to end torque T . The size of the cross section is 200 mm × 10 mm. The beam properties are Young’s modulus E = 200,000 MPa and G = 80,000 MPa, and length L = 1,000 mm. The results obtained using 1 and 20 new elements are shown in Fig. 11 and compared with the theoretical solution of Trahair [29]. It is shown that the new element can consider the Wagner effect, which induced the nonlinear phenome non. The cantilever has axial shortening deformation due to the applied torque. Since the results of using 20 displacement -based elements coincide with the results of using 1 element, it can be accepted that only 1 element is enough to model the nonlinear torsion behavior by adopting the newly developed displacement -based element. In addition, Battini and Pacoste [17] also presented a similar example of nonlinear torsion of a rectangular section, where 20 elements were used to get the results in their work. Fig. 10. Nonlinear torsion of a cantilever Fig. 11. Moment end- rotation relationship for nonlinear torsion Torsion of a C antilever Beam with an Angle Section This example considers a cantilever beam with an angle section subjected to an end torque acting on its shear center. As shown in Fig. 12, the geometric and material pr operties of the cantilever are: L = 177.8 mm, α = 90°, b = 14.605 mm, t = 0.9601 mm, Young’s modulus E = 89,632 MPa and shear modulus G = 33,445 MPa. The example was experimentally and numerically studied by different researchers [15, 31] . Fig. 13 compares the path of shear center B of the end cross section obtained by the present study, the experimental results from Gregory [31] , and the numerical results from Hsiao and Lin [15] using 20 elastic beam elements. This work uses 20 displacement -based elements, whose results match with the experimental and numerical results very well. Fig. 12. Torsion of a cantilever beam of angle section Fig. 13. Path of point B during torsional behavior Flexural -Torsional Interaction This example investigates the influence of bending moments on the torsional stiffness of a simply supported beam. The experiments were conducted by Engel and Goodier [32] on simply supported angle bars (see Fig. 14). The material was 61 S -T aluminum with Young’s modulus E = 10,000 ksi and shear modulus G = 3,750 ksi . The properties of the monosymmetric angle section are: L = 36.8 in., α = 90°, b = 1 in., and t = 0.030 in. Bending moments acting in the plane of symmetry were applied to the ends of the angle bar to produce constant bending moment between supports. A torque is applied to the mid -span of the member. Fig. 15 shows the comparison of torque -twist lines from experiments and simulations using 20 elements in the present study. The results indicate that the new element can successfully model the flexural -torsional interaction of beams. Fig. 14. Simply supported angle bar subjected to end moments and mid -span torque Fig. 15. Twis t per unit length versus applied Flexural -Torsional Buckling of Concentrically Loaded Angle Struts This example, analyzed theoretically and numerically by Kitipornchai and Chan [8] and Kitipornchai and Lee [33] , studies the elastic flexural -torsional buckl ing load of a concentrically loaded, pin-ended unequal -leg angle bar (L76x51x5 mm). The material properties are: Young’s modulus of elasticity E = 200,000 MPa, Poisson’s ratio ν = 0.3, and yielding stress F y = 312 MPa. Fig. 16 shows the elastic buckling factors for members with different modified slenderness. The definition of modified slenderness is 𝜆𝜆=�𝐹𝐹𝐼𝐼/(𝜋𝜋2𝐸𝐸)𝐿𝐿/𝑟𝑟𝑚𝑚𝑖𝑖𝑖𝑖, where L is the length of the angle bar and 𝑟𝑟𝑚𝑚𝑖𝑖𝑖𝑖 is the minimum radius of gyration. In the figure, Pc is the applied axial load and P y is the section yielding load (squash load). Twenty new displacement -based elements are used to simulate this flexural -torsio nal buckling behavior. It is shown that the buckling loads computed form the newly developed element coincide well with the theoretical values. Since the default reference line of the new element is the shear center axis, the rigid offsets developed in Section 3.7 are utilized to move the element reference line to centroidal axis so that the concentric loading condition is satisfied. Fig. 16. Flexural -torsional buckling of concentrically loaded angle struts (elastic) For the inelastic buckling case, the idealized residual stress distribution shown in Fig. 17 is used both in the research of Kitipornchai and Lee [33] and in the present simulation. An elastic-perfectly plastic stress- strain curve is employed, while in the plastic region the shear modulus is assumed to remain unchanged. In the calculation of inelastic buckling loads, Kitipornchai and Lee [33] utilized a simplified mathematical model similar to that proposed by Trahair and Kitipornchai [34] for inelastic buckling of steel I -beams. Using the new displacement -based and mixed elements, the inelastic buckling loads of the angle struts with different slenderness are calculated and then compared with the results of Kitipornchai and Lee [33] in Fig. 18. Fig. 17. Assumed residual stress distribution for angle sections Fig. 18. Flexural -torsional buckling of concentrically loaded angle struts ( inelastic) Cantilever with A ngle Section Subjected to an E ccentric Axial Force This example, studied numerically by several researchers [9, 16, 35] , addresses a cantilever beam with asymmetric cross section subjected to an eccentric axial force P . The cantilever is fixed at the left end in Fig. 19 . The axial force P is applied to the shear center of the free end. As shown in Fig. 19. the geometric and material properties of the cantilever are: L = 1,400 mm, a = 76 mm, b = 51 mm, t = 6.5 mm, Young’s modulus E = 193.05 MPa and Poisson’s ratio ν = 0.3. The load -deflection curves of the present study together with the results of Alsafadie et al. [35] are plotted in Fig. 20, in which the axial displacement u and lateral displacement v in Y direction are presented. As seen in this figure, three new displacement -based elements in the present study are enough to obt ain relatively accurate results, while Alsafadie et al. [35] utilized 30 displacement -based elements to get similar results. Fig. 19. Cantilever with asymmetric angle section subjected to an eccentric axial force Fig. 20. Cantilever with asymmetric angle section: load -deflection u and v curves Lateral -Torsional Buckling of a C antilever with T S ection This example is part of the buckling studies on cantilever beams subjected to an end force conducted by sever al researchers [15, 36-39] . Two cases, 4Ba65 and 4Ba50, in their research are simulated using the element developed in the present study. Here, 4Ba65 means that T section is used, larger flange is at the top, the load acts at the upper face, and the length of the cantilever is L = 65 in; while 4Ba50 mean the same with 4Ba65 except for the length is L = 50 in. The loading and geometry properties are shown in Fig. 21 with b = 1.239 in., t f = 0.1236 in., h = 2.8292 in., and t w = 0.0863 in. The buckling loads obtained using 5 elements developed in the present study are shown in Table 1 together with the results from other researchers for comparison, which proves the accuracy of the new element in the prediction of buckling loads. Fig. 21. Cantilever beam subjected to an end force Table 1 Buckling loads of canti lever beams subjected to an end force Case Buckling loads (lb) Exp. [39] Theory [39] FEM [38] FEM [36, 37] Present 4Ba65 21.2 21.4 21.1 21.6 21.2 4Ba50 33.6 32.7 32.1 33.2 32.4 Buckling of a R ight-Angled Frame under End Load The right -angled frame depicted in Fig. 22 was analyzed by many authors [3, 24, 26, 40] . The load P is applied in the X -direction at the member tip with a very small per turbation load of 0.0002 N in the Z-direction at the tip to induce buckling artificially. At first, the problem is solved using the new element with membrane locking, the computed response of the applied load P and the tip deflection in Z-direction is plotted in Fig. 23 together with other results of this work. It is shown that due to membrane locking, more than 20 elements in total are needed for the two members to get comparable results with other work. If fewer elements are used, the post -buckling stiffness is overest imated dramatically. Fig. 22. Right -angled frame subjected to end load Fig. 23. Load/tip Z -displacement for right -angled frame (with membrane locking) However, if the new element with membrane locking remedied is used, only 10 elements in total are enough to obtain an accurate response (see Fig. 24). Thus, it may be accepted that the adopted methodology can remove m embrane locking and improve the results significantly. Fig. 24. Load/tip Z -displacement for right -angled frame (without membrane locking) Flexural -Torsional Buckling of a T Beam This example addresses flexural -torsional buckling of a T beam under axial load considering both elastic and inelastic material. Battini and Pacoste [18] conducted numerical study of this example as shown in Fig. 25 with the following geometric material properties: L = 1,800, h = b = 60, t = 6, Young’s modulus E =70,000, Poisson’s ratio ν = 0.33, yielding stress Fy = 20, and post -yield strain hardening modulus E t = E/5. For the inelastic case, a bilinear plastic constitutive relation is assumed. The left end of the beam is fixed and only the axial displacement is allowed at the right end. A compressive axial force P is applied at the shear center of the right end. In addition, a small perturbation force F = P/1000 is introduced at the midspan point in the Z -direction. In the present research, 4 and 8 elements with 304 fibers in the cross section are used to simulate the buckling and yielding behavior. Note that Battini and Pacoste [18] used 40 displacement -based elements in their simulation. Fig. 25. Flexural -torsional buckling of a T beam The load/Z -displacement curves of point O in both elastic and inelastic cases are shown in Fig. 26 and Fig. 27, respectively. The agreement between the present simulation results an d the results of Battini and Pacoste [18] is very good, but fewer elements are needed for the element in this work . Fig. 26. Load/Z -displacement of point O curve (elastic case) Fig. 27. Load/Z -displacement of point O curve (inelastic case) Lee’s Frame Lee et al. [41] investigated a structural frame represented in Fig. 28, using a linear elastic material and neglecting axial deformations. Several other researchers numerically studied this structure also assuming linear elastic material [40, 42, 43] . de Souza [2 4] and Cichoń [44] analyzed this example for both elastic and inelastic materials. In the present study, a bilinear elasto -plastic material model with kinematic hardening presented in de Souza [24] is adopted for the inelastic case (see Fig. 28 for the material properties). Weak axis bending is used for both the horizontal and vertical members. Fig. 28. Lee’s frame The frame is first analyzed using the new displacement -based element with membrane locking for both elastic and inelastic cases. Fig. 29 and Fig. 30 compare the load displacement curves for the elastic case and inelastic case, respectively. It can be seen that 20 elements (10 elements per member) with membrane locking still produce over stiff results, while 40 elements can give results with sufficient accuracy. If only 10 elements with membrane locking are used, the results will be severely inaccurate (too stiff), so they are not shown in the figures. This example shows that for the displacement -based element co nsidering high order strain terms, the membrane locking problem must be addressed, otherwise more elements are needed to obtain accurate results. Fig. 29. Equilibrium path for Lee’s frame (element with membrane locking, elastic) Fig. 30. Equilibrium path for Lee’s frame (element with membrane locking, inelastic) The frame is then analyzed using the displacement -based element with membrane locking remedied. Fig. 31 and Fig. 32 illustrate the results from the present study and those from de Souza [24]. It can be concluded that the displacement -based element with membrane locking remedied can reduce the required elements amount to only 5, compared with the displacement -based element with membrane locking. For “DB5” in Fig. 31 and Fig. 32, three elements are used for the horizontal member and two elements are used for the vertical member. For “DB4”, three elements are used for the horizontal member and one element is used for the vertical member. Fig. 31. Equilibrium path for Lee’s frame (element with membrane locking remedied, elastic) Fig. 32. Equilibrium path for Lee’s frame (element with membrane locking remedied, inelastic) Conclusions Geometrically and materially nonlinear displacement -based element for structural members with asymmetric cross sections have been developed and vali dated in this work to analyze structures like transmission towers. The geometric nonlinear phenomenon is simulated by using the corotational transformation, and by adopting the total Lagrangian formulation in the basic system through employing Green-Lagran ge strains. Consequently, the newly developed element is able to take into account the axial -flexural -torsional interaction behavior of structural members in large deformation. In order to model members with asymmetric sections, the element DOFs are defined with respect to both centroid and shear center in the basic system, and then transformed to the shear center in advance of the corotational transformation. The material nonlinearity is considered by discretizing the member cross section into fibers and applying uniaxial constitutive law and numerical integration to calculate cross-section forces. Examples show that nonlinear torsion, elastic and inelastic flexural -torsional buckling, and lateral-torsional buckling can be modeled accurately with only a small number of elements. The displacement -based element yields overestimated stiffness for both elastic and inelastic inextensional bending modes unless the membrane locking phenomenon is remedied. With the approach adopted in this work to alleviate membrane locking, the improved version of the displacement -based element presents excellent performance for a comprehensive set of structural members and frames. Acknowledgement This research is supported by the National Science Foundation under Grant No. CRISP -1638234 and Northeastern University. This support is gratefully acknowledged. Appendix A. Linearization of the Governing Equation This appendix shows details of the linearization process of the governing equation to derive the tangent stiffness matrix and in ternal forces, which are used in Section 3.4. Using a Taylor series expansion to Eq. (42) and neglecting the high order terms, the following equation is obtained 𝒈𝒈𝑖𝑖+1≈𝒈𝒈𝑖𝑖+𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0𝒈𝒈�𝑫𝑫𝑏𝑏𝑖𝑖+𝑑𝑑∆𝑫𝑫𝑏𝑏�+𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0𝒈𝒈�𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖+𝑑𝑑∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡�=𝟎𝟎 (A.1) where α and β are scalers used to scale vectors ∆𝑫𝑫𝑏𝑏 and ∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡. Variables with a superscript i or (i+1) represent their values in the ith or (i+1)th iteration. The terms in the above equation are obtained as shown below 𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0𝒈𝒈�𝑫𝑫𝑏𝑏𝑖𝑖+𝑑𝑑∆𝑫𝑫𝑏𝑏� =�𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0(𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇)𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑺𝑺 𝑙𝑙0𝑑𝑑𝜕𝜕+�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0(𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇)𝑺𝑺 𝑙𝑙0𝑑𝑑𝜕𝜕 +�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0(𝑺𝑺) 𝑙𝑙0𝑑𝑑𝜕𝜕 =�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑮𝑮𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕∆𝑫𝑫𝑏𝑏+�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑲𝑲𝑠𝑠𝑵𝑵𝛿𝛿𝛿𝛿1 𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕∆𝑫𝑫𝑏𝑏 (A.2) where the following relations are used 𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇�=0 (A.3) 𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0�𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇�𝑺𝑺=𝑑𝑑𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇 𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏)𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏) 𝑑𝑑𝑑𝑑𝑺𝑺=𝑮𝑮𝑑𝑑𝒗𝒗 𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏)∆𝑫𝑫𝑏𝑏=𝑮𝑮𝑵𝑵𝛿𝛿𝛿𝛿2∆𝑫𝑫𝑏𝑏 (A.4) 𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0(𝑺𝑺)=𝑑𝑑𝑺𝑺 𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏)𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏) 𝑑𝑑𝑑𝑑 =�𝒀𝒀𝑇𝑇𝑑𝑑𝝈𝝈 𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏)𝑑𝑑𝑑𝑑∆𝑫𝑫𝑏𝑏 𝐴𝐴0 =�𝒀𝒀𝑇𝑇𝑬𝑬𝑑𝑑𝝐𝝐 𝑑𝑑(𝑑𝑑∆𝑫𝑫𝑏𝑏)𝑑𝑑𝑑𝑑∆𝑫𝑫𝑏𝑏= 𝐴𝐴0�𝒀𝒀𝑇𝑇𝑬𝑬𝒀𝒀𝑵𝑵𝛿𝛿𝛿𝛿1𝑵𝑵𝛿𝛿𝛿𝛿2𝑑𝑑𝑑𝑑∆𝑫𝑫𝑏𝑏 𝐴𝐴0 =�𝒀𝒀𝑇𝑇𝑬𝑬𝒀𝒀𝑑𝑑𝑑𝑑𝑵𝑵𝛿𝛿𝛿𝛿1𝑵𝑵𝛿𝛿𝛿𝛿2∆𝑫𝑫𝑏𝑏 𝐴𝐴0=𝑲𝑲𝑠𝑠𝑵𝑵𝛿𝛿𝛿𝛿1𝑵𝑵𝛿𝛿𝛿𝛿2∆𝑫𝑫𝑏𝑏 (A.5) in which 𝑮𝑮= ⎣⎢⎢⎢⎢⎢⎡0 0 0 0𝑁𝑁 0 0 0𝑁𝑁0 0 0 0 0 0 0 0 00 𝑁𝑁𝑧𝑧𝑠𝑠 −𝑁𝑁𝑦𝑦𝑠𝑠 0 0 0 0 0 0 0 0 00 0𝑀𝑀𝐼𝐼 0 0𝑀𝑀𝐼𝐼 𝑀𝑀𝐼𝐼𝑀𝑀𝐼𝐼 0 0 0 0 0𝑁𝑁𝑧𝑧𝑠𝑠−𝑁𝑁𝑦𝑦𝑠𝑠 0 0 0 𝑊𝑊 ⎦⎥⎥⎥⎥⎥⎤ (A.6) 𝑲𝑲𝑠𝑠=�𝒀𝒀𝑇𝑇𝑬𝑬𝒀𝒀𝑑𝑑𝑑𝑑 𝐴𝐴0=� ⎣⎢⎢⎢⎢⎡𝐸𝐸−𝑦𝑦𝐸𝐸𝑧𝑧𝐸𝐸𝑝𝑝2𝐸𝐸 0 −𝑦𝑦𝐸𝐸𝑦𝑦2𝐸𝐸−𝑦𝑦𝑧𝑧𝐸𝐸−𝑝𝑝2𝑦𝑦𝐸𝐸 0 𝑧𝑧𝐸𝐸−𝑦𝑦𝑧𝑧𝐸𝐸𝑧𝑧2𝐸𝐸𝑝𝑝2𝑧𝑧𝐸𝐸 0 𝑝𝑝2𝐸𝐸−𝑝𝑝2𝑦𝑦𝐸𝐸𝑝𝑝2𝑧𝑧𝐸𝐸𝑝𝑝4𝐸𝐸 0 0 0 0 0 4𝑛𝑛2𝐺𝐺⎦⎥⎥⎥⎥⎤ 𝑑𝑑𝑑𝑑 𝐴𝐴0 (A.7) in which 𝑝𝑝2=(𝑦𝑦−𝑦𝑦𝑠𝑠)2+(𝑧𝑧−𝑧𝑧𝑠𝑠)2 (A.8) Here, 𝑲𝑲𝑠𝑠 is the section stiffness matrix. The matrix E sets the constitutive relation between the stresses and strains. Here 𝑬𝑬=�𝐸𝐸 0 0𝐺𝐺� (A.9) where E is the tangent Young’s modulus and G is the tangent shear modulus. In addition, 𝑑𝑑 𝑑𝑑𝑑𝑑|𝛼𝛼,𝛽𝛽=0𝒈𝒈�𝑸𝑸𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖+𝑑𝑑∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡�=𝑑𝑑𝒈𝒈 𝑑𝑑(𝑑𝑑∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡)𝑑𝑑(𝑑𝑑∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡) 𝑑𝑑𝑑𝑑=−∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡 (A.10) Therefore, the following relation is obtained 𝒈𝒈𝑖𝑖+�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑮𝑮𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕∆𝑫𝑫𝑏𝑏+�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑲𝑲𝑠𝑠𝑵𝑵𝛿𝛿𝛿𝛿1 𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕∆𝑫𝑫𝑏𝑏−∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡=𝟎𝟎 (A.11) which can be expressed as the following equation with the help of Eq. (42) ��𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑮𝑮𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕+�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇𝑲𝑲𝑠𝑠𝑵𝑵𝛿𝛿𝛿𝛿1 𝑵𝑵𝛿𝛿𝛿𝛿2 𝑙𝑙0𝑑𝑑𝜕𝜕�∆𝑫𝑫𝑏𝑏 =∆𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡+𝑷𝑷𝑒𝑒𝑒𝑒𝑡𝑡𝑖𝑖−�𝑵𝑵𝛿𝛿𝛿𝛿2𝑇𝑇𝑵𝑵𝛿𝛿𝛿𝛿1𝑇𝑇 𝑙𝑙0𝑺𝑺𝑖𝑖𝑑𝑑𝜕𝜕 (A.12) This iteration equation can be expressed in short as shown in Eq. (43 ). Appendix B. The Transformation for Rigid Offsets This appendix shows details of the transformatio n used to consider eccentricities as discussed in Section 3.7. From the equilibriums on rigid offsets 1 -I and 2 -J (see Fig. 6), the following transformation for the end forces of the element reference ends P and the end forces at the shear centers 𝑷𝑷� is obtained 𝑷𝑷=𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝑇𝑇𝑷𝑷� (B.1) where the transformation matrix is given by 𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝑇𝑇= ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1 0 0 000 0 0 0 000 0 1 0 000 0 0 0 000 0 0 1 000 0 0 0 000 0 0 0 100 0 0 0 000 0−𝑑𝑑𝑍𝑍𝐼𝐼𝑑𝑑𝑌𝑌𝐼𝐼 010 0 0 0 000 𝑑𝑑𝑍𝑍𝐼𝐼 0−𝑑𝑑𝑋𝑋𝐼𝐼 001 0 0 0 000 −𝑑𝑑𝑌𝑌𝐼𝐼𝑑𝑑𝑋𝑋𝐼𝐼 0 000 1 0 0 000 0 0 0 000 0 1 0 000 0 0 0 000 0 0 1 000 0 0 0 000 0−𝑑𝑑𝑍𝑍𝐽𝐽𝑑𝑑𝑌𝑌𝐽𝐽 100 0 0 0 000𝑑𝑑𝑍𝑍𝐽𝐽 0−𝑑𝑑𝑋𝑋𝐽𝐽 010 0 0 0 000−𝑑𝑑𝑌𝑌𝐽𝐽𝑑𝑑𝑋𝑋𝐽𝐽 0 001⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ (B.2) According to the principle of virtual work, the contragredient transformation of the end displacements is obtained as this tangential relationship 𝛿𝛿𝑫𝑫�=𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝛿𝛿𝑫𝑫 (B.3) where D denotes the end displacements of the element reference ends. Small r otations are assumed in the derivations in this section. With 𝛿𝛿𝑷𝑷�=𝑲𝑲�𝛿𝛿𝑫𝑫�, as developed in Section 2.2 , we have 𝛿𝛿𝑷𝑷=𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝑇𝑇𝛿𝛿𝑷𝑷�=𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝑇𝑇𝑲𝑲�𝛿𝛿𝑫𝑫�=𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝑇𝑇𝑲𝑲�𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝛿𝛿𝑫𝑫=𝑲𝑲𝛿𝛿𝑫𝑫 (B.4) where 𝑲𝑲=𝑻𝑻𝑜𝑜𝑜𝑜𝑜𝑜𝑇𝑇𝑲𝑲�𝑇𝑇𝑜𝑜𝑜𝑜𝑜𝑜 is the tangent stiffness matrix in the global system considering rigid offsets.