Complete tangent stiffness matrix considering higher-order terms in the strain tensor and large rotations for a Euler Bernoulli - Timoshenko space beam-column element

This paper presents a unified method developed by Rodrigues et al. [1] to obtain a complete tangent stiffness matrix for spatial geometric nonlinear analysis using minimal discretization. The formulation presents four distinct important aspects to a complete analysis: interpolation (shape) functions, higher-order terms in the strain tensor and in the finite rotations, an updated Lagrangian kinematic description, and shear deformation effect (Timoshenko beam theory). Thus, the tangent stiffness matrix is calculated from the differential equation solution of deformed infinitesimal element equilibrium, considering the axial load and the shear deformation in this relation. This solution provides interpolation functions that are used in an updated Lagrangian formulation to construct the spatial tangent stiffness matrix considering higher-order terms in the strain tensor and in the finite rotations. The method provides an efficient formulation to perform geometric nonlinear analyses and predict the critical buckling load for spatial structures with moderate slenderness and with the interaction between axial and torsion effects, considering just one element in each member or a reduced discretization.• Complete expressions for a geometric nonlinear analyses considering one element per member• Spatial analyses considering higher-order terms in the strain tensor and large rotations• Shear deformation influence included


Method details
In the finite elements method (FEM), the analytical behavior of a solid can be approximated by a discrete behavior. The discrete solution is usually obtained by nodal displacements, while the continuous problem solution can be found by interpolating the nodal displacements using interpolating (shape) functions [2] . However, in general, the discrete solution using FEM introduces simplifications to the analytical mathematical idealization for the structure behavior, since the interpolation functions that define the deformed configuration of a structure are not compatible with the mathematical idealization for the response of the continuous medium [3] . An exception is for a linear analysis with constant cross-section elements, in which the cubic interpolation functions (the so-called Hermitian functions) represent the mathematical idealization of the structure behavior. However, in geometric nonlinear analysis, these cubic interpolation functions do not represent the analytical response.
Therefore, in geometric nonlinear analysis, to converge to the analytical behavior, a refined discretization or the use of a high-order finite element with an increase in the order of the polynomial of the basis functions is necessary. However, when the shape functions are obtained from the differential equation solution of the problem, the frame member discretization could be reduced [4][5][6][7] . In fact, in this case, the continuous behavior of the element is represented by nodal parameters without the consideration of any other approximation.
The differential equation that describes a linear problem is obtained from the equilibrium of an infinitesimal element in an undeformed configuration, leading to the Hermitian interpolation functions. However, for geometric nonlinear or second-order analysis, the equilibrium of an infinitesimal element in the deformed configuration should be considered, and the differential equation considers the nonlinearity caused by the axial load in the deformed infinitesimal element [8] .
Furthermore, to implement the unified method proposed by Rodrigues et al. [1] , other important effects must be considered. The influence of shear deformation in beam-columns with a moderate slenderness ratio or with a small shear to bending ratio is well known. Indeed, employing the Timoshenko beam theory (TBT) in these cases provides better results than considering the Euler Bernoulli beam theory (EBBT) [9][10] . Moreover, all expressions for the Timoshenko beam theory can be presented in a similar way to the corresponding expressions for Euler Bernoulli beam theory, differing only by auxiliary parameters that multiply the terms of the expressions. Thus the interpolation functions and the tangent stiffness matrix of a nonlinear analysis can be developed to consider either the Timoshenko or the Euler-Bernoulli beam theory, by making just a small parameter change.
In addition, to reach the expressions presented in this paper, an updated Lagrangian kinematic description considering spatial elements is used. Besides, the stiffness matrix is adjusted to consider finite rotation effects. Furthermore, the formulation employs the Green strain tensor considering all higher-order terms, which leads to accurate results and reduces the buckling load of spatial frame structures, showing that not considering these effects can overestimate the resistance of a structure.
The presented method combines three types of nonlinearities, namely, higher-order terms in the Green strain tensor, higher-order terms in the finite rotations, and the nonlinearity induced by the equilibrium of the infinitesimal element in the deformed configuration considering the axial effect. One main difference of the present formulation is that the interpolation functions for displacements and rotations in an element are directly obtained from the solution of differential equations that consider the equilibrium in the deformed configuration. As a consequence, the resulting interpolation functions depend on the axial force in the element. In addition, the present formulation considers the shear deformation influence, and expressions are developed in 3D. These considerations extend the classical tangent stiffness matrix for geometric nonlinear analysis [ 2 , 11-14 ], as explained in the sequel.
The manuscript provides explicit expressions for the proposed interpolation functions and tangent stiffness matrix. Since this formulation involves trigonometric and hyperbolic functions, for the sake of computational efficiency the tangent stiffness matrix can be conveniently rewritten in polynomial functions using a Taylor series expansion. The manuscript shows expressions for expansions using two, three, and four terms. The paper in which the formulation is deduced [1] demonstrates that the tangent stiffness matrix with four terms provides accurate results when compared to the original matrix based on trigonometric and hyperbolic functions. MATLAB and C open-source codes with expressions presented in this paper are available in the GitLab [15] repository and in the MathWorks file exchange [16] .

Matrix approach to geometric nonlinear analysis
In a second-order elastic analysis, the structure behavior can be stated according to the equation incremental nodal loads and reactions, and [ K t ] is a tangent stiffness matrix [2] . Usually, this matrix is composed by two components, [ K t ] = [ K e + K g ] , a linear elastic component, [ K e ] , and a nonlinear, [ K g ] , the geometric stiffness matrix. The matrix can be developed in different ways and members of the structure need to be divided into several elements to provide satisfactory results. Also, the nonlinearity included in the matrix depends on distinct aspects, as the interpolation functions, the kinematic description, and the strain-displacement relations. Moreover, the matrix considers shear deformations and the expressions depend on the so-called Timoshenko constant .
The conventional formulation of the tangent stiffness matrix is obtained from the virtual work principle considering an updated Lagrangian description using cubic (Hermitian) functions and disregard higher-order terms in the strain tensor and shear deformations (Euler-Bernoulli beam theory). However, in a geometric nonlinear analysis context, this discrete solution introduces simplifications to the mathematical idealization for the structure behavior, since these interpolation functions are not compatible with the mathematical idealization for the response of the continuous medium. Additionally, according to [10] consider higher-order terms in the strain tensor leads to smaller buckling loads, and, to predict the behavior of a beam-column with moderate slenderness, or with a small shear-to-bending rigidity ratio, the shear deformation has important effects, and in this case, the Timoshenko beam theory (TBT) provides better results.
Therefore, the method presented in this paper is for geometric nonlinear analysis and consists of calculating the tangent stiffness matrix of a 3D frame finite element, [ K t ] , according to the expressions presented below. A structure member (beam or column) can be modeled considering a reduced, or even a minimal, discretization, in which each member corresponds to just one finite element. The formulation is also capable to evaluate structures considering the Timoshenko beam theory.

Tangent stiffness matrix [ K]
Theoretical concepts of the matrix formulated are developed in [1] . The authors explain that an exponential solution can be used, or, as an alternative, hyperbolic expressions can be used in cases of tensile forces, positive axial loads ( P > 0), and trigonometric expressions in cases of compressive forces, negative axial loads ( P < 0). The separation between tension and compression can be an additional conditional statement that can influence the efficiency of the code. This separation is not necessary if the implementation works directly with exponential functions [ 8 , 17 ]. The method proposed in this article employs hyperbolic and trigonometric functions because the tangent stiffness matrix is simpler. The steps necessary to use the method are postulated: For a tensile force ( P > 0), the elastic component of the tangent stiffness matrix can be written in terms of hyperbolic functions, and the coefficients become: For a compressive force ( P < 0), the elastic component of the tangent stiffness matrix can be written in terms of trigonometric functions, and the coefficients become: (20) For a compressive force ( P < 0), the geometric component of the tangent stiffness matrix can be written in terms of trigonometric functions, and the coefficients become: Interaction matrix of axial load and torsion -[ K int ] : For a tensile force ( P > 0), the matrix that considers the interaction between the axial load and the torsion is given according: The coefficients are written in terms of hyperbolic functions: For a compressive force ( P < 0), the interaction between the axial load and the torsion is given according: The coefficients are written in terms of trigonometric functions: ϕ y = 2 y L 2

Taylor Series Expansion
The tangent stiffness matrix developed using complete expressions with hyperbolic and trigonometric functions (TBT_Large_complete) provides accurate results for nonlinear analysis and for buckling load prediction to spatial structures. However, the hyperbolic and trigonometric functions might be computationally inefficient. In addition, in some cases, mainly for initials loads steps in a nonlinear code and for reduced axial loads, the algorithm can present numerical instability. To improve computation efficiency and avoid numerical problems, a Taylor series expansion is proposed, which generates polynomial functions. Furthermore, the user can select the number of terms in the Taylor expansion according to: TBT_Large_2tr -Timoshenko beam theory, 2 terms in tangent stiffness matrix (1 elastic + 1 geometric). The geometric term corresponds to factor that multiplies the coefficient P . Conventional formulation (Rodrigues et al. [ 10 ]) TBT_Large_3tr -Timoshenko beam theory, 3 terms in tangent stiffness matrix (1 elastic + 2 geometric). The geometric term corresponds to factor that multiplies the coefficients P and P ².
TBT_Large_4tr -Timoshenko beam theory, 4 terms in tangent stiffness matrix (1 elastic + 3 geometric). The geometric term corresponds to factor that multiplies the coefficients P, P ² and P ³.
The coefficients of the complete tangent stiffness matrix of a 3D frame element, including the elastic components and the geometric components with all terms of the Taylor expansion, are shown in eq. 70 to 93 :       30240 0 0 E 3 I y 2 Iz 12 y + 1   The examples clearly evidence the efficiency of the proposed method to predict the critical buckling load for the spatial frame with no discretization and with torsion interaction. For the spatial  frame with torsion and the asymmetric spatial frame, results considering the complete formulation or the Taylor series expansion with 4 or 3 terms, and using just one element in each member provides differences with the response of the discretized structure using conventional formulation in order of 0,5%. In order to check the performance of the stiffness matrices presented in this paper, a convergence study was implemented. Fig. 4 presents the convergence rate found considering different options for the tangent stiffness matrix. The study was carried out by comparing the top displacement of a beamcolumn to the one obtained by a highly discretized solution using conventional elements. The beamcolumn was loaded according to Fig. 4 , with P = 0.6 Pcr (60% of the critical buckling load).
The presented method can be used to perform a nonlinear geometric analysis of spatial structures using minimal or a reduced discretization, whereas, while the usual tangent stiffness matrix would require a more refined mesh to predict the structure behavior. The method can be implemented in any structural analysis software. Matrix coefficients use hyperbolic and trigonometric expressions. Alternatively, polynomial expressions for geometric stiffness coefficients using a Taylor series expansion with 4 or 3 terms are proposed, and all necessary equations are available online in open source.

Analytical interpolation functions
The presented tangent stiffness was calculated by using interpolation functions for displacements and rotations in a frame element obtained directly from the solution of differential equations that consider the equilibrium in the deformed configuration. As a consequence, these interpolation functions depend on the axial force in the element, and also on the parameters of the Timoshenko beam theory. Moreover, the interpolation functions deduced in the deformed configuration can be used to improve visualization in a graphics program and the computation of the local P-delta effect using just one element per member. For a tensile force ( P > 0), the interpolation functions for displacements and rotations are written in terms of hyperbolic functions, according to: Rodrigues [17] also presents analytical interpolation functions for displacements and rotations considering exponential functions along with expressions for articulated left member end, and for articulated right member end. In cases for which the constant is very small, the Timoshenko beam theory behaves as the Euler-Bernoulli beam theory and the interpolation functions lead to the same results. Expressions for the interpolation functions considering just the Euler-Bernoulli beam theory are also featured in [17] and implemented in [15] and [16] .
Analytical interpolation functions plots are shown in Figs. 5-10 . It can be noticed that when the axial load is null ( P = 0 ), conventional cubic interpolation functions are identical to the ones developed in this work. As the axial load increases, analytical interpolations functions present a distinct behavior, as expected. This influence becomes more evident for elements with > 0 , for which the Timoshenko beam theory is used.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.