Abstract
Nonlinear analyses using an updated Lagrangian formulation considering the EulerBernoulli beam theory have been developed with consistency in the literature, with different geometric matrices depending on the nonlinear displacement parts considered in the strain tensor. When performing this type of analysis using the Timoshenko beam theory, in general, the stiffness and the geometric matrices present additional degrees of freedom. This work presents a unified approach for the development of a geometric matrix employing the Timoshenko beam theory and considering higherorder terms in the strain tensor. This matrix is obtained using shape functions calculated directly from the solution of the differential equation of the problem. The matrix is implemented in the Ftool software, and its results are compared against several matrices found in the literature, with or without higherorder terms in the strain tensor, as well as the EulerBernoulli or Timoshenko beam theories. Examples show that the use of the Timoshenko beam theory has a strong influence, especially when the structure has small slenderness (short members). For high axial load values, the consideration of higherorder terms in the strain tensor results in larger displacements as expected.
Keywords:
Geometric Matrix; Shape Functions; Timoshenko Beam Theory; HigherOrder Terms in Strain Tensor
1 INTRODUCTION
A structural geometric nonlinear analysis using the finite element method (FEM) depends on the consideration of four aspects: the bending theory, the kinematic description, the straindisplacement relations and the interpolating (shape) functions.
The most commonly used bending solution for frame elements is the EulerBernoulli beam theory (EBBT), which is the one implemented in most structural analysis software, with a large number of applications.
However, in some cases, this theory cannot predict the correct behaviour of the structure. To study a laminated beam or a beamcolumn with a moderate slenderness ratio, the Timoshenko beam theory (TBT) provides better results. Several highorder bending theories can be found in the literature, such as in Levinson (1981Levinson, M. (1981). A new rectangular beam theory, Journal of Sound and Vibration. 74: 8187.), Bickford (1982Bickford, W.B. (1982) A consistent higher order beam theory, Developments in Theoretical and Applied Mechanics 11: 137.), Heyliger and Reddy (1988Heyliger, P.R., Reddy, J.N. (1988). A higher order beam finite element for bending and vibration problems, Journal of Sound and Vibration 126: 309326.), Petrolito (1995Petrolito, J. (1995). Stiffness analysis of beams using a higherorder theory, Computers and Structures 55: 3339.), Reddy et al. (1997Reddy, J.N. (1997). On lockingfree shear deformable beam finite elements, Computer Methods in Applied Mechanics and Engineering 149: 113132.), Reddy (1997Reddy, J.N., Wang, C.M., Lee, K.H. (1997). Relationships between bending solutions of classical and shear deformation beam theories, International Journal Solids Structures 34(26): 33733384.), Tessler and Gherlone (2007Tessler, A., Gherlone, M.C.M. (2007). Refinement of Timoshenko beam theory for composite and sandwich beams using zigzag kinematics. DC: NASA, Washington. (NASA TP2007215086).) and Meghare and Jadhao (2015Meghare, T.K., Jadhao, P.D. (2015). A simple higher order theory for bending analysis of steel beams, International Journal of Civil Engineering 2: 3138.). This work considers the Timoshenko beam theory due to its influence on beams with small slenderness, especially when subjected to axial loads (beamcolumns).
Regarding the kinematic description, a geometric nonlinear analysis can be performed through the full Lagrangian, updated Lagrangian or a corotational approach. Basically, these descriptions differ with respect to the reference configuration considered (Felippa 2017Felippa, C.A., (2017), Nonlinear finite element methods / NFEM, Lecture notes for the course nonlinear finite element methods, Center for Aerospace Structures, University of Colorado, Boulder/USA.). When consistently developed, the total Lagrangian and the updated formulation produce the same results (McGuire et al. 2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.). These formulations are well developed in the literature, such as in Bathe (1996Bathe, K.J. (1996). Finite Element Procedures, PrenticeHall, Englewood Cliffs, NJ, USA.) and Bathe and Bolourchi (1979), in which the nonlinear analysis of threedimensional structures is studied using both Lagrangian formulations. Nanakorn and Vu (2006Nanakorn, P., Vu, L.N. (2006). A 2D fieldconsistent beam element for large displacement analysis using the total Lagrangian formulation, Finite Elements in Analysis and Design 42: 12401247.) employ a total Lagrangian formulation to analyse a plane element subject to large displacements.
Other studies, such as Santana and Silveira (2014Santana, M.V.B., Silveira, R.A.M. (2014). Sistema computacional gráfico interativo para problemas de instabilidade em pórticos planos, in: XXXV Iberian Latin American Congress on Computational Methods in Engineering, Fortaleza, Brazil.) and Silva et al. (2016Silva, J.L., Lemes, I.J.M., Silveira, R.A.M., Silva, A.R.D. (2016). Influência da teoria de viga na análise geometricamente não linear de estruturas reticuladas, in: XXXVII Iberian Latin American Congress on Computational Methods in Engineering, Brasília, Brazil.), have built a tangent matrix for Timoshenko beam elements using the corotational formulation. However, in Santana and Silveira (2014Santana, M.V.B., Silveira, R.A.M. (2014). Sistema computacional gráfico interativo para problemas de instabilidade em pórticos planos, in: XXXV Iberian Latin American Congress on Computational Methods in Engineering, Fortaleza, Brazil.), cubic (Hermitian) shape functions are used, while Silva et al. (2016Silva, J.L., Lemes, I.J.M., Silveira, R.A.M., Silva, A.R.D. (2016). Influência da teoria de viga na análise geometricamente não linear de estruturas reticuladas, in: XXXVII Iberian Latin American Congress on Computational Methods in Engineering, Brasília, Brazil.) use an element proposed by Tang et al. (2015Tang, Y.Q., Zhou, Z.H., Chan, S.L. (2015). Nonlinear BeamColumn Element Under Consistent Deformation, International Journal of Structural Stability and Dynamics 15(5).) with a consistent axial displacement field interpolation. Considering a consistent field approach and a corotational formulation, Kien (2012Kien, N.D. (2012). A Timoshenko beam element for large displacement analysis of planar beams and frames, International Journal of Structural Stability and Dynamics 12(6).) developed a Timoshenko beam plane element for large displacements employing the corotational method. Silva and Silva (2010Silva, S.S., Silva, W.T.M. (2010). Estudo de pórticos planos utilizando um elemento finito de viga unificado em um programa de análise linear, in: Mecánica Computacional Vol XXIX, Buenos Aires, Argentina. [in Portuguese].) use an element that includes the theories of EulerBernoulli, Timoshenko and third order linear frame analysis. Oliveira and Silva (2017Oliveira, G.C., Silva, W.T.M. (2017). Análise nãolinear de arcos utilizando o elemento de viga unificado BernoulliTimoshenko e a formulação corotacional, Civil Engineering Electronic Journal. 13(2): 116.) employ a unified BernoulliTimoshenko element.
The straindisplacement relations have an important role in a geometric nonlinear analysis. Usually, the geometric matrix is obtained considering small displacement gradients, neglecting some higherorder terms in the Green strain tensor, as done in McGuire et al. (2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.).
In this work, a geometric stiffness matrix is developed considering higherorder terms in the strain tensor based on an updated Lagrangian formulation for the equilibrium equations, in which the complete Green strain tensor is employed. The steps to build this matrix considering the EulerBernoulli beam theory are well known (Bathe and Bolourchi 1979Bathe, K.J., Bolourchi, S. (1979). Large displacement analysis of threedimensional beam structures, International Journal for Numerical Methods in Engineering 14: 961986., Conci 1988Conci, A. (1988) Análise de estruturas reticuladas de aço com consideração de empenamento e nãolinearidades geométrica e material, Ph.D. Dissertation, Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, PUCRio, Rio de Janeiro, Brazil., Yang and Leu 1994Yang, Y. B., Leu, L.J. (1994). Nonlinear stiffnesses in analysis of planar frames, Computer Methods in Applied Mechanics and Engineering 117: 233247., Yang and Kuo 1994 Yang, Y. B., Kuo, S.R. (1994). Theory & analysis of nonlinear framed structures, Prentice Hall, Simon & Schuster (Asia) Pte ltd, Singapore.and Bathe 1996Bathe, K.J. (1996). Finite Element Procedures, PrenticeHall, Englewood Cliffs, NJ, USA.); however, this work employs these steps using a Timoshenko beam element, i.e., it also considers shear deformations.
One way to formulate the elastic and geometric stiffness matrices of the Timoshenko beamcolumn elements is achieved by using the usual beam nodal displacements and leaving the shear distortions as independent variables. This formulation results in additional stiffness terms leading to an element stiffness matrix of order 14, and static condensation is used to reduce the matrix order to 12 (Bathe and Bolourchi 1979Bathe, K.J., Bolourchi, S. (1979). Large displacement analysis of threedimensional beam structures, International Journal for Numerical Methods in Engineering 14: 961986., Aguiar et al. 2014Aguiar, L.L., Almeida, C.A., Paulino, G.H. (2014). A threedimensional multilayered pipe beam element: nonlinear analysis, Computers and Structures 138: 142161.).
In some cases, only the Timoshenko elastic matrix is developed; otherwise, the geometric stiffness matrix does not consider higherorder terms in the strain tensor (Davis et al. 1972Davis, R., Henshell, R.D., Warburton G.B. (1972). A Timoshenko beam element, Journal of Sound and Vibration 22 (4): 475487., Friedman and Kosmatka 1992Friedman, Z., Kosmatka, J.B. (1992). An improved twonode Timoshenko beam finite element, Computers and Structures 47 (3): 473481., Yunhua 1998Yunhua, L. (1998). Explanation and elimination of shear locking and membrane locking with field consistence approach, Computer Methods in Applied Mechanics and Engineering 162: 249269.).
In the FEM, the continuous (analytical) behaviour of a solid is approximated by a discrete behaviour. Usually, the discrete solution is obtained by nodal displacements, while the continuous solution can be found by means of interpolating functions, which are polynomials in general. To converge to the analytical behaviour, a refined discretization or the use of a highorder finite element with an increase in the order of the polynomial of the basis functions is necessary. Therefore, some authors develop highorder expansion basis functions for structured elements aiming for the accuracy and efficiency of the method (Zheng and Dong 2011Zheng, X., Dong, S. (2011). An eigenbased highorder expansion basis for structured spectral elements, Journal of Computational Physics 230: 85738602., Rodrigues et al. 2016Rodrigues, C. F., Suzuki, J. L., Bittencourt, M. L. (2016). Construction of minimum energy highorder Helmholtz bases for structured elements, Journal of Computational Physics 306: 269290.).
However, when shape functions are obtained from the solution of the differential equations of an infinitesimal element, without the consideration of any other approximation except for those already covered in the analytical idealization of the element behaviour, discretization is not necessary for a linear analysis. When an undeformed configuration is considered for this infinitesimal element, the axial force influence on transversal behaviour is not considered, leading to cubic functions. When shear deformations are not considered, the EulerBernoulli beam theory leads to Hermitian functions.
In the formulation presented in this work, no additional terms are necessary to develop the geometric stiffness matrix since shape functions are obtained directly from the solution of the differential equations of the problem, considering the Timoshenko beam theory.
Schramm et al. (1994Schramm, U., Kitis, L., Kang, W., Pilkey W.D. (1994). On the shear deformation coefficient in beam theory, Finite Elements in Analysis and Design 16: 141162.) and Pilkey et al. (1995Pilkey, W.D., Kang, W., Schramm U. (1995). New structural matrices for a beam element with shear deformation, Finite Elements in Analysis and Design 19: 2544.) also developed a stiffness matrix from the differential equations considering shear effects; however, shape functions were not developed.
Thus, the main contribution of this research is to provide a unified approach to building the geometric stiffness matrices of twodimensional Timoshenko beamcolumn elements, considering higherorder terms in the Green strain tensor. The shape functions are obtained directly from the solution of the differential equation of an infinitesimal Timoshenko beam element (Burgos and Martha 2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil., Martha and Burgos 2014Martha, L.F., Burgos, R.B. (2014). Diferenças na consideração da distorção no modelo de Timoshenko de uma viga submetida a carregamento axial, in: South American Workshops on Structural Engineering, Montevideo, Uruguay. [in Portuguese]., 2015Martha, L.F., Burgos, R.B. (2015). Possíveis inconsistências na consideração da distorção por cisalhamento numa viga submetida a carregamento axial, in: Brazilian Concrete Congress, Bonito, Brazil. [in Portuguese].). Then, using an updated Lagrangian kinematic description and considering a higherorder Green strain tensor, the elastic and geometric stiffness matrices are obtained.
The results clearly show the influence of the beam theory used and the importance of considering higherorder terms in strain tensors to predict the critical loads of framed structures, especially for small slenderness beamcolumns and high axial loads.
Research on the extension of this formulation is being developed to fully exploit its potential and reduce the influence of bar discretization. Future work will consider shape functions obtained from the equilibrium of a deformed infinitesimal element, that includes the influence of axial force and 3D elements.
2 BEAM BEHAVIOR IDEALIZATION
This section presents and solves the differential equations that define the analytical behaviour of an undeformed infinitesimal beam element considering the Timoshenko beam theory, obtained from the equilibrium conditions, compatibility relations and constitutive material laws.
2.1 Displacement field
As shown in Figure 1, the displacement field of a beam is defined according to relation (1):
2.2 Differential equilibrium relationships in beamcolumns
Equilibrium conditions must be satisfied for the entire structure, and the analysis includes the equilibrium conditions of an infinitesimal beam element. Figure 2 shows an infinitesimal beam element subjected to a distributed transversal load q and a distributed longitudinal load p.
From the equilibrium of an infinitesimal beam element, equation (2) can be obtained. From the approximate relation between bending moment and curvature, expression (3) can be written.
where $v\left(x\right)$ is the infinitesimal element transversal displacement, $q\left(x\right)$ is the transverse distributed load, $Q\left(x\right)$ is the transverse shear force and $M\left(x\right)$ is the bending moment in the crosssection.
2.3 EulerBernoulli beam theory (EBBT)
The EulerBernoulli beam theory considers the rotation as the derivative of the transverse displacement ($\theta =dv/dx$). Thus, using equations (2) and (3), the differential equation that governs the element behaviour can be written according to relation (4).
2.4 Timoshenko beam theory (TBT)
In the Timoshenko beam theory, shear distortion is considered as an additional rotation of the cross section, according to Figure 3. Therefore, the crosssectional rotation and the transverse displacement are considered to be independent variables.
According to the TBT, the total crosssection rotational $\left(dv/dx\right)$ is composed of bending rotation $\left(\theta \right)$ increased by the shear distortion $\left(\gamma \right)$, according to the expression (5):
The shear force acting on the section is given by the following equation:
where $G$ is the material shear modulus, $A$ is the crosssectional area, and $\chi $ is the factor that defines the effective area for crosssectional shear. Substituting equations (5) and (6) in the differential equation (3) of the infinitesimal beam element equilibrium, the following expression can be found:
2.5 Differential equation solution
There are different forms of the solution for the differential equation (7) of a Timoshenko infinitesimal beam element, such as the auxiliary function, presented in Shirima and Giger (1992)Shrima, L.M., Giger, M.W. (1992). Timoshenko beam element resting on twoparameter elastic foundation, Journal of Engineering Mechanics 118(2): 280295., or the one developed in Onu (2008Onu, G. (2008). Finite elements on generalized elastic foundation in Timoshenko beam theory, Journal of Engineering Mechanics 134(9): 763776.), which will be used in this work. The transverse displacement of the structure will have a bending and a shear contribution, ${v}_{b}\left(x\right)$ and ${v}_{s}\left(x\right)$, respectively, according to:
Using equation (4) and finding the homogeneous solution for EulerBernoulli beams, expression (9) for the displacement ${v}_{b}\left(x\right)$, i.e., the bending contribution, can be found. The homogeneous solution corresponds to an unloaded element situation depending only on the boundary conditions and obtained as:
Using the result obtained for ${v}_{b}\left(x\right)$ from equation (9), in the differential relation (7), the expression for the shear distortion becomes:
Finally, based on equations (9) and (10), the equation of the total transverse displacement of the beam, equation (8), is described by the following equation (11), and it is usual to adopt a dimensionless factor Ω, introduced by Reddy (1997Reddy, J.N., Wang, C.M., Lee, K.H. (1997). Relationships between bending solutions of classical and shear deformation beam theories, International Journal Solids Structures 34(26): 33733384.), in this equation.
3 TIMOSHENKO SHAPE FUNCTIONS
In a finite element analysis, the analytical behaviour of a solid can be approximated by a discrete behaviour. The discrete solution is usually obtained by nodal displacements, while the continuous problem solution can be found by interpolating the nodal displacements using shape functions (McGuire et al. 2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.).
When the shape functions are obtained from the differential equation solution of the problem, the bar discretization is unnecessary because continuous behaviour of the element can be represented by nodal parameters without the consideration of any other approximation, except for those already contained in the analytical idealization of the element behaviour. In this section, these shape functions for a Timoshenko beam element are calculated as in Martha (2018Martha, L.F. (2018). Análise matricial de estruturas com orientação a objetos, Elsevier, Rio de Janeiro, Brazil.), Burgos and Martha (2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.) and Onu (2008Onu, G. (2008). Finite elements on generalized elastic foundation in Timoshenko beam theory, Journal of Engineering Mechanics 134(9): 763776.).
However, in this work, the differential equation used to obtain the shape functions is the one for an undeformed element, Figure 2, without the influence of the axial load in transversal behaviour (Martha 2018Martha, L.F. (2018). Análise matricial de estruturas com orientação a objetos, Elsevier, Rio de Janeiro, Brazil.); therefore, for a nonlinear analysis, bar discretization is necessary. Previous work such as Burgos and Martha (2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.) and Onu (2008Onu, G. (2008). Finite elements on generalized elastic foundation in Timoshenko beam theory, Journal of Engineering Mechanics 134(9): 763776.) consider this influence, but only for small displacement gradients.
3.1 Timoshenko shape functions development
Shape functions shown in equation (12) interpolate nodal displacements (Figure 4).
As previously observed, for an infinitesimal beam element, the homogeneous solution of the problem gives the transversal displacement (${\nu}_{0}^{h}$) from equation (11), and the rotation of the cross section (${\theta}_{0}^{h}$) can be obtained by substituting this transversal displacement in equation (8). Thus, both displacement solutions can be written in matrix form, according to the following equation:
The boundary conditions are obtained by evaluating the homogeneous solution of these displacements (${\nu}_{0}^{h}$ and ${\theta}_{0}^{h}$) at the extreme nodes of the bar, according to equation (14):
Shape functions can be obtained using equations (12), (13) and (14), resulting in:
Finally, Timoshenko beam element shape functions are given by:
The axial displacement of the element is given according to equation (17), and in general, the interpolation functions are linear without the interaction between the axial and transversal formulations:
4 UPDATED LAGRANGIAN FORMULATION
To predict the geometric nonlinear behaviour of structures using the finite element method, three different descriptions can be used: total Lagrangian, updated Lagrangian and, more recently, the corotational approach. These formulations differ essentially with respect to the reference configuration. Both Lagrangian formulations, when used consistently, lead to the same results; in this work, the stiffness matrices are calculated considering the updated Lagrangian formulation, and the steps shown have been presented in McGuire et al. (2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.). The final equations are used to find the Timoshenko stiffness matrix considering higherorder terms in the strain tensor.
In this formulation, the equilibrium equations of the unknown configuration $t+\text{\Delta}t$ must be written using the known variables of a reference configuration $t$. For a configuration $t+\text{\Delta}t$, the virtual work of the internal forces must be equal to the virtual work of the external forces, as shown in equation (18):
where ${S}_{ij}^{\left(t+\Delta t\right)}$ corresponds to the second PiolaKirchoff stress tensor, ${\epsilon}_{ij}^{\left(t+\u2206t\right)}$ the GreenLagrange strain tensor, and ${R}^{\left(t+\u2206t\right)}$ the virtual work due to external loading, which could include body, surface, inertia and damping forces (Aguiar et al. 2014Aguiar, L.L., Almeida, C.A., Paulino, G.H. (2014). A threedimensional multilayered pipe beam element: nonlinear analysis, Computers and Structures 138: 142161.).
The linearized incremental equation requires small displacement increments; thus, the second PiolaKirchoff stress tensor and the GreenLagrange strain tensor can be written according to:
where ${\tau}_{ij}^{t}$ corresponds to the Cauchy stress tensor, $\u2206{\tau}_{ij}$ is the stress increment and $\u2206{\epsilon}_{ij}$ is the deformation increment, which can be calculated from the GreenLagrange strain tensor.
Knowing that in the reference configuration, the element is subjected only to rigid body motion, $\delta {\epsilon}_{ij}^{t}$ = 0, the left side of expression (18) can be written as equation (20).
The GreenLagrange strain tensor in the xy plane is represented by:
The strain tensor shown in relation (21) has a linear ($\u2206{e}_{ij}$) and a nonlinear part ($\u2206{\eta}_{ij}$); thus, the decomposition $\u2206{\epsilon}_{ij}=\u2206{e}_{ij}+\u2206{\eta}_{ij}$ can be adopted, as in equation (22):
The stress increment is obtained from the material constitutive relation. Considering also a linear approximation for stress and strain increment leads to the expression:
Finally, based on equations (20) and (23), the virtual work equation becomes:
In the xy plane, the stress vector, the constitutive matrix, and the linear and nonlinear strain vectors are given by:
Thus, the virtual work equation, expression (24), can be expanded according to equation (26), using the relations presented next.
5 STIFFNESS MATRICES
The stiffness matrix development is performed by analysing a twodimensional structure in the xy plane. Using an updated Lagrangian formulation, the higherorder terms in the Green strain tensor and the shape functions are developed, and the elastic and geometric stiffness matrices are calculated.
5.1 Timoshenko plane element
According to the beam displacement field equation (1), the linear and nonlinear parts of the GreenLagrange strain tensor components in equation (21), can be rewritten as expressions (27) and (28), respectively:
5.2 Elastic stiffness matrix
The linear part of equation (26) provides the elastic stiffness matrix. For a better understanding, this part is divided into two: $\delta {W}_{1}$, the first integral, and $\delta {W}_{2}$, the second integral of equation (26). Employing relations (27), $\delta {W}_{1}$ and $\delta {W}_{2}$ can be written as equations (29) and (30), respectively.
In the beam centroidal axis, $\underset{A}{\int}{y}^{2}dA={I}_{z},\phantom{\rule{0.25em}{0ex}}}{\displaystyle \underset{A}{\int}ydA=0},$ and the expressions (29) and (30) are reduced to equations (31) and (32):
According to expressions (12), (17) and Figure 4, the rotation, transverse and axial displacement in any section of an element can be calculated by interpolating the nodal displacements of the element, $\left\{u\right\}=\left\{\begin{array}{cc}{d}_{1}^{\text{'}}& {d}_{4}^{\text{'}}\end{array}\right\}$ and $\phantom{\rule{0.25em}{0ex}}\left\{v\right\}=\left\{\begin{array}{cc}\begin{array}{cc}{d}_{2}^{\text{'}}& {d}_{3}^{\text{'}}\end{array}& \begin{array}{cc}{d}_{5}^{\text{'}}& {d}_{6}^{\text{'}}\end{array}\end{array}\right\}$, with the shape functions developed in equations (16) and (17), i.e., ${N}_{v}=\left\{\begin{array}{cc}\begin{array}{cc}{N}_{2}^{v}\left(x\right)& {N}_{3}^{v}\left(x\right)\end{array}& \begin{array}{cc}{N}_{5}^{v}\left(x\right)& {N}_{6}^{v}\left(x\right)\end{array}\end{array}\right\}$, ${N}_{\theta}=\left\{\begin{array}{cc}\begin{array}{cc}{N}_{2}^{\theta}\left(x\right)& {N}_{3}^{\theta}\left(x\right)\end{array}& \begin{array}{cc}{N}_{5}^{\theta}\left(x\right)& {N}_{6}^{\theta}\left(x\right)\end{array}\end{array}\right\}$ and ${N}_{u}=\left\{\begin{array}{cc}{N}_{1}\left(x\right)& {N}_{4}\left(x\right)\end{array}\right\}$. In this manner, the element displacements can be represented by the following relations:
Finally, substituting relation (33) in the displacements of equations (31) and (32), expressions (34) and (35) can be found considering the notation $\partial /\partial x=\left(\right)\text{'}$.
Solving for the integrals in equations (34) and (35), the elastic stiffness matrix can be calculated, resulting in the matrix (36):
5.3 Geometric stiffness matrix
From the beam displacement field, equation (1), and the GreenLagrange strain tensor, equation (21), the nonlinear part of this tensor can be written according to:
It is important to observe that in expression (37), it is usual to neglect the flexural part of the strain in ${\eta}_{xx}$. Usually, ${\eta}_{xy}$ is also neglected, and equation (38) is obtained as in McGuire et al. (2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.):
In this work, however, the complete strain tensor is used (equation (37)), as in Yang and Kuo (1994Yang, Y. B., Kuo, S.R. (1994). Theory & analysis of nonlinear framed structures, Prentice Hall, Simon & Schuster (Asia) Pte ltd, Singapore.), considering the Timoshenko and not the EulerBernoulli beam theory. Therefore, the geometric stiffness matrix is calculated with the nonlinear part of the virtual work principle, equation (26), and can be written as expression (39).
Applying relations, $\underset{A}{\int}{t}_{xx}dA}=P{;}^{}{\displaystyle \underset{A}{\int}{t}_{xx}ydA}={M}_{z}{;}^{}{\displaystyle \underset{A}{\int}{t}_{xy}dA}={Q}_{y},$expression (39) becomes equation (40). With equation (33), the displacements can be written using shape functions, resulting in equations (41) and (42).
According to Figure 5, which shows a planar framework, and considering a constant shear force, the bending moment and the shear force equations of a planar frame element can be calculated as presented in (43).
Frame element  (adapted from Pereira 2002Pereira, A. (2002). Projeto ótimo de pórticos planos com restrição à flambagem, Msc. Dissertation, Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, PUCRio, Rio de Janeiro, Brazil.)
Substituting the shape functions, presented in equations (16) and (17), solving the integrals of the problem leads to the geometric stiffness matrix, considering the Timoshenko beam theory and higherorder terms in strain tensor, shown by matrix (44) below:
In this work, this matrix is called TBT_large. Without the consideration of shear strain, i.e.$,\phantom{\rule{0.25em}{0ex}}GA\chi \to \infty $, the problem is solved by EulerBernoulli beam theory. Considering higherorder terms in the strain tensor in this work this matrix is called EBBT_large and is the same as developed in Yang and Leu (1994Yang, Y. B., Kuo, S.R. (1994). Theory & analysis of nonlinear framed structures, Prentice Hall, Simon & Schuster (Asia) Pte ltd, Singapore.).
In this paper, matrices considering small displacement gradients were also used for comparison. These matrices use the EulerBernoulli beam theory and are developed in many references, such as McGuire et al. (2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.). The elements using these matrices were named in this work as EBBT_small. The stiffness matrix considering Timoshenko beam theory and small displacement gradients can be obtained by neglecting higherorder terms in the Green strain tensor. In this work, the examples that were modelled with these matrices were named TBT_small. This matrix is also presented in Burgos and Martha (2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.) using a Taylor series expansion.
6 NUMERICAL ANALYSIS
This section aims to evaluate the different solutions for the nonlinear analysis of structures exposed in this work, using simple numerical examples, by considering the EulerBernoulli and Timoshenko beam theories. All stiffness matrices were implemented in the Ftool computer software (Martha 1999Martha, L.F. (1999). Ftool: A structural analysis educational interactive tool, Proceedings of Workshop in Multimedia Computer Techniques in Engineering Education, Institute for Structural Analysis, Technical University of Graz, Austria, 5165.), and their results were compared with results developed in the literature and with the Mastan2 v3.5 software developed by McGuire et al. (2000McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.), by considering both beam theories, i.e., EBBT_Mastan (EulerBernoulli beam theory) and TBT_Mastan (Timoshenko beam theory).
6.1 Critical loads of columns
The first example presented studies of the buckling load of columns with different boundary conditions. To verify the influence of the Timoshenko beam theory in the nonlinear analysis and the influence of considering higherorder terms in the strain tensor, different slenderness ratios were studied. The columns presented in Figure 6 have length L=1 m, Young’s modulus E=107 kN/m2, a section form factor $\chi $=1 and a null Poisson’s ratio $v$.
Analysed Columns  (adapted from Silva et al. 2016Silva, J.L., Lemes, I.J.M., Silveira, R.A.M., Silva, A.R.D. (2016). Influência da teoria de viga na análise geometricamente não linear de estruturas reticuladas, in: XXXVII Iberian Latin American Congress on Computational Methods in Engineering, Brasília, Brazil.)
The columns were modelled with 5 elements, and the critical load is normalised by the flexural rigidity of the cross section (EI). The equilibrium paths of the structure were studied, as shown in Figure 7 (clamped column), considering different slenderness ratios, using different geometric matrices and comparing the results with the software Mastan2 and the analytical Euler’s critical load. Neglecting shear deformation, the critical load is given by ${P}_{cr}={\pi}^{2}EI/{\left(2L\right)}^{2}$.
Analysing the equilibrium paths, it can be observed that for high values of slenderness, the beam behaviour is similar for both the EulerBernoulli and Timoshenko theories. No relevant differences are observed with or without the consideration of higherorder terms in the strain tensor.
However, as the slenderness decreases, $\lambda <4.0$, the influence of the beam theory considered becomes evident. Moreover, a clear difference is noticed when considering higherorder terms in the strain tensor. Thus, for high loads and small slenderness, it is important to consider Timoshenko beam theory and higherorder terms in strain tensors.
In this case, the element developed in this work (TBT_Large) reduces the critical buckling load, which means that not considering higherorder terms and beam theory influence can be against safety in some specific cases (for small slenderness).
For a simply supported column, the buckling load is given by ${P}_{cr}={\pi}^{2}EI/{\left(L\right)}^{2}$, and the equilibrium path using different geometric matrices is shown in Figure 8.
The same conclusions can be observed from these equilibrium paths: for high values of slenderness ratios, the behaviour considering EulerBernoulli and Timoshenko theories is similar. Additionally, relevant differences are observed regarding the consideration of higherorder terms in the strain tensor. However, as the slenderness decreases, $\lambda <4.0$, the influence of the beam theory and the consideration of higherorder terms in strain tensor (TBT_Large) considerably reduces the critical buckling load of the column.
For a fixed and simply supported column, the buckling load is given by ${P}_{cr}={\pi}^{2}EI/{\left(0.7L\right)}^{2}$, and the equilibrium path using different geometric matrices is shown in Figure 9.
Even for high values of slenderness ratios, considering Timoshenko beam theory results in a substantial reduction of the critical buckling load. The use of a geometric matrix that considers higherorder terms in the strain tensor (TBT_Large) can also significantly change the results, providing a reduction of the critical buckling load, and again it can be concluded that not considering higherorder terms and beam theory influence can be against safety in the prediction of the stability of columns, even for high values of slenderness ratios.
6.2 Continuous beam column critical load
The second example examines the influence of beam theory and the consideration of higherorder terms in the strain tensor in a continuous beam column subjected to axial load at its vertex. This problem is presented in Timoshenko and Gere (1963Timoshenko, S.P., Gere, J.M (1963). Theory of elastic stability, International student edition, second edition, McGrawHill, Singapore.) and is shown in Figure 10.
Continuous beam column  (adapted from Timoshenko and Gere 1963Timoshenko, S.P., Gere, J.M (1963). Theory of elastic stability, International student edition, second edition, McGrawHill, Singapore.)
The geometry of the structure, material and section proprieties are the same as in the first example. To verify the influence of the Timoshenko beam theory in the nonlinear analysis and the influence of the consideration of higherorder terms in strain tensor, different slenderness ratios were considered, modifying the section height (h).
Each span was discretized using 5 elements, and the critical load was normalised by the flexural rigidity of the cross section (EI). According to Timoshenko and Gere (1963Timoshenko, S.P., Gere, J.M (1963). Theory of elastic stability, International student edition, second edition, McGrawHill, Singapore.), when the span lengths are equal, the critical buckling load is given by ${P}_{cr}={\pi}^{2}EI/{\left(L\right)}^{2}$. The equilibrium paths of the structure were studied and are shown in Figure 11 and Figure 12.
As in example 1, for high values of slenderness ratios, the beam behaviour considering the EulerBernoulli and Timoshenko theories is similar, and the equilibrium paths overlap. No relevant differences are observed with or without the consideration of higherorder terms in the strain tensor.
As the slenderness ratio is reduced ($\lambda <4.0$), the beam theory influence becomes more evident, and a clear difference is noticed when considering the higherorder terms in the strain tensor (TBT_Large). Not considering these influences can be against safety in the prediction of the stability of these structures since a considerable reduction in the critical buckling load can be seen in Figure 12.
6.3 Roorda’s frame analysis
The third example examines the influence of beam theory and the consideration of higherorder terms in the strain tensor in Roorda’s frame, subjected to a vertical load at its vertex. The studied frame is shown in Figure 13. The structure has a section form factor $\chi $=5/6, Poisson’s ratio $v$=0.3, Young’s modulus E=107 kN/m2 and length L=1 m.
Neglecting shear deformation, the critical load has an analytical value given by:
However, analytically, the critical load considering shear deformation can be obtained by solving the following equation (Burgos and Martha 2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.):
Burgos and Martha (2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.) studied the influence of Ω on the critical load; however, this expression gives better results for smaller values of Ω since nonlinearity leads to changes in the compressive force applied in the column. This effect is greater in beamcolumns with small slenderness ratios.
The results obtained with the different formulations for stiffness and geometric matrices are shown in Figure 14. The bars are discretized with 5 finite elements each.
For Ω = 0, the equilibrium paths for all formulations are similar, the equilibrium paths overlap, and the result for critical load is close to the analytical one. Meanwhile, for other values of Ω, the TBT_Large element presents results closer to the analytical solution proposed in Burgos and Martha (2013Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.). Moreover, the influence of the beam theory in the problem solution becomes evident with the reduction in the critical buckling load that can be seen in Figure 14. This figure also shows that employing the EulerBernoulli beam theory in structures with high Ω values leads to an incorrect prediction of their behaviour.
From these analyses, as presented in Figure 14, it can also be concluded that considering higherorder terms in the strain tensor provides a small reduction in the critical load buckling compared with neglecting these effects, for the same bar discretization. Therefore, combining accurate beam theory and a more complete stiffness matrix leads to lower critical loads in frames.
6.4 Unbraced portal frame analysis
The last example examines the beam theory influence and the consideration of higherorder terms in the strain tensor in an unbraced portal frame, Figure 15. The structure has the same properties as example 6.3, section form factor $\chi $=5/6, Poisson’s ratio $v$=0.3, Young’s modulus E=107 kN/m2 and length L=1 m.
The frame is loaded by vertical loads $\mu P$ and by a small lateral disturbing load $H$. Bazant and Cedolin (2010Bazant, Z., Cedolin, L. (2010). Stability of Structures, World Scientific Publishing Co. Pte. Ltd., Singapore.) present a solution for this problem, considering the EulerBernoulli beam theory, based on Euler’s critical load according to:
The results obtained for horizontal displacement of the upper left node using the different formulations for the stiffness matrix are exposed in Figure 16. The bars were discretized with 5 finite elements each.
Fixed Frame  (adapted from Bazant and Cedolin 2010Bazant, Z., Cedolin, L. (2010). Stability of Structures, World Scientific Publishing Co. Pte. Ltd., Singapore.)
As expected, for large slenderness ratios, the results using different beam theories are identical, and the consideration of higherorder terms in the strain tensor is also irrelevant.
However, based on Figure 16, the beam theory considered in the analysis has a considerable influence on the results for small slenderness ratios with a reduction in the critical buckling load. Moreover, when considering higherorder terms in the strain tensor, this reduction in the frame buckling load for both beam theories is amplified.
7 CONCLUSIONS
This research presented a unified approach to developing the tangent stiffness matrix for a Timoshenko beam element in the geometric nonlinear analysis of structures, using an updated Lagrangian description and considering higherorder terms in the strain tensor. These matrices are obtained using shape functions calculated from the solution of the equilibrium differential equation of the problem, with no additional approximations necessary. The matrices developed can consider either the Timoshenko or the EulerBernoulli beam theory, by making a small parameter change.
Estimates for the buckling loads of columns and frames using the proposed nonlinear geometric matrix were found using the Timoshenko or EulerBernoulli beam theory. The reported results clearly illustrate the relevance of considering the Timoshenko beam theory. The influence is more evident in elements with small slenderness ratios. It is also shown that this difference becomes more evident close to the buckling load, that is, when high axial loads are applied.
The results also show that considering higherorder terms in the strain tensor leads to smaller buckling loads than using a nonlinear geometric matrix that does not take these effects into account. Therefore, for cases with high axial loads and small slenderness ratios, these considerations lead to larger displacements in columns and in framed structures. Thus, structural design neglecting these effects can be against safety, due to underestimated results.
The developed geometric stiffness matrix can accurately describe the behaviour for both beam theories. The matrices consider higherorder terms in the strain tensor, leading to accurate results with no additional bar discretization than what is usually used in FE analysis. Considering this development, only this matrix needs to be implemented since for slender beams, it will naturally converge into the usual EulerBernoulli geometric stiffness matrix.
It is also worth mentioning that, in this paper, shape functions do not consider axial load influence; however, research on the extension of the basis function is being developed to fully exploit the potential of the proposed formulation. Future work will combine element nonlinearity, using higherorder terms in the Green strain tensor, with the infinitesimal element nonlinearity considering axial load influence. To consider this behaviour, shape functions obtained directly from the equilibrium of a deformed infinitesimal element will be used. Additionally, 3D geometric matrices are being developed to verify the contribution of torsion and axial force to the analyses.
Acknowledgement
This work has been partially supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and FAPERJ.
References
 Aguiar, L.L., Almeida, C.A., Paulino, G.H. (2014). A threedimensional multilayered pipe beam element: nonlinear analysis, Computers and Structures 138: 142161.
 Bathe, K.J. (1996). Finite Element Procedures, PrenticeHall, Englewood Cliffs, NJ, USA.
 Bathe, K.J., Bolourchi, S. (1979). Large displacement analysis of threedimensional beam structures, International Journal for Numerical Methods in Engineering 14: 961986.
 Bazant, Z., Cedolin, L. (2010). Stability of Structures, World Scientific Publishing Co. Pte. Ltd., Singapore.
 Bickford, W.B. (1982) A consistent higher order beam theory, Developments in Theoretical and Applied Mechanics 11: 137.
 Burgos, R.B., Martha, L.F. (2013). Exact shape functions and tangent stiffness matrix for the buckling of beamcolumns considering shear deformation. in: XXXIV Iberian Latin American Congress on Computational Methods in Engineering, Pirenópolis, Brazil.
 Conci, A. (1988) Análise de estruturas reticuladas de aço com consideração de empenamento e nãolinearidades geométrica e material, Ph.D. Dissertation, Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, PUCRio, Rio de Janeiro, Brazil.
 Davis, R., Henshell, R.D., Warburton G.B. (1972). A Timoshenko beam element, Journal of Sound and Vibration 22 (4): 475487.
 Felippa, C.A., (2017), Nonlinear finite element methods / NFEM, Lecture notes for the course nonlinear finite element methods, Center for Aerospace Structures, University of Colorado, Boulder/USA.
 Friedman, Z., Kosmatka, J.B. (1992). An improved twonode Timoshenko beam finite element, Computers and Structures 47 (3): 473481.
 Heyliger, P.R., Reddy, J.N. (1988). A higher order beam finite element for bending and vibration problems, Journal of Sound and Vibration 126: 309326.
 Kien, N.D. (2012). A Timoshenko beam element for large displacement analysis of planar beams and frames, International Journal of Structural Stability and Dynamics 12(6).
 Levinson, M. (1981). A new rectangular beam theory, Journal of Sound and Vibration. 74: 8187.
 McGuire, W., Gallagher, R.H., Ziemian, R.D. (2000). Matrix structural analysis, John Wiley & Sons Inc, New York, NY, USA.
 Martha, L.F. (2018). Análise matricial de estruturas com orientação a objetos, Elsevier, Rio de Janeiro, Brazil.
 Martha, L.F. (1999). Ftool: A structural analysis educational interactive tool, Proceedings of Workshop in Multimedia Computer Techniques in Engineering Education, Institute for Structural Analysis, Technical University of Graz, Austria, 5165.
 Martha, L.F., Burgos, R.B. (2014). Diferenças na consideração da distorção no modelo de Timoshenko de uma viga submetida a carregamento axial, in: South American Workshops on Structural Engineering, Montevideo, Uruguay. [in Portuguese].
 Martha, L.F., Burgos, R.B. (2015). Possíveis inconsistências na consideração da distorção por cisalhamento numa viga submetida a carregamento axial, in: Brazilian Concrete Congress, Bonito, Brazil. [in Portuguese].
 Meghare, T.K., Jadhao, P.D. (2015). A simple higher order theory for bending analysis of steel beams, International Journal of Civil Engineering 2: 3138.
 Nanakorn, P., Vu, L.N. (2006). A 2D fieldconsistent beam element for large displacement analysis using the total Lagrangian formulation, Finite Elements in Analysis and Design 42: 12401247.
 Oliveira, G.C., Silva, W.T.M. (2017). Análise nãolinear de arcos utilizando o elemento de viga unificado BernoulliTimoshenko e a formulação corotacional, Civil Engineering Electronic Journal. 13(2): 116.
 Onu, G. (2008). Finite elements on generalized elastic foundation in Timoshenko beam theory, Journal of Engineering Mechanics 134(9): 763776.
 Pereira, A. (2002). Projeto ótimo de pórticos planos com restrição à flambagem, Msc. Dissertation, Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, PUCRio, Rio de Janeiro, Brazil.
 Petrolito, J. (1995). Stiffness analysis of beams using a higherorder theory, Computers and Structures 55: 3339.
 Pilkey, W.D., Kang, W., Schramm U. (1995). New structural matrices for a beam element with shear deformation, Finite Elements in Analysis and Design 19: 2544.
 Reddy, J.N. (1997). On lockingfree shear deformable beam finite elements, Computer Methods in Applied Mechanics and Engineering 149: 113132.
 Reddy, J.N., Wang, C.M., Lee, K.H. (1997). Relationships between bending solutions of classical and shear deformation beam theories, International Journal Solids Structures 34(26): 33733384.
 Rodrigues, C. F., Suzuki, J. L., Bittencourt, M. L. (2016). Construction of minimum energy highorder Helmholtz bases for structured elements, Journal of Computational Physics 306: 269290.
 Santana, M.V.B., Silveira, R.A.M. (2014). Sistema computacional gráfico interativo para problemas de instabilidade em pórticos planos, in: XXXV Iberian Latin American Congress on Computational Methods in Engineering, Fortaleza, Brazil.
 Schramm, U., Kitis, L., Kang, W., Pilkey W.D. (1994). On the shear deformation coefficient in beam theory, Finite Elements in Analysis and Design 16: 141162.
 Shrima, L.M., Giger, M.W. (1992). Timoshenko beam element resting on twoparameter elastic foundation, Journal of Engineering Mechanics 118(2): 280295.
 Silva, J.L., Lemes, I.J.M., Silveira, R.A.M., Silva, A.R.D. (2016). Influência da teoria de viga na análise geometricamente não linear de estruturas reticuladas, in: XXXVII Iberian Latin American Congress on Computational Methods in Engineering, Brasília, Brazil.
 Silva, S.S., Silva, W.T.M. (2010). Estudo de pórticos planos utilizando um elemento finito de viga unificado em um programa de análise linear, in: Mecánica Computacional Vol XXIX, Buenos Aires, Argentina. [in Portuguese].
 Tang, Y.Q., Zhou, Z.H., Chan, S.L. (2015). Nonlinear BeamColumn Element Under Consistent Deformation, International Journal of Structural Stability and Dynamics 15(5).
 Tessler, A., Gherlone, M.C.M. (2007). Refinement of Timoshenko beam theory for composite and sandwich beams using zigzag kinematics. DC: NASA, Washington. (NASA TP2007215086).
 Timoshenko, S.P., Gere, J.M (1963). Theory of elastic stability, International student edition, second edition, McGrawHill, Singapore.
 Yang, Y. B., Leu, L.J. (1994). Nonlinear stiffnesses in analysis of planar frames, Computer Methods in Applied Mechanics and Engineering 117: 233247.
 Yang, Y. B., Kuo, S.R. (1994). Theory & analysis of nonlinear framed structures, Prentice Hall, Simon & Schuster (Asia) Pte ltd, Singapore.
 Yunhua, L. (1998). Explanation and elimination of shear locking and membrane locking with field consistence approach, Computer Methods in Applied Mechanics and Engineering 162: 249269.
 Zheng, X., Dong, S. (2011). An eigenbased highorder expansion basis for structured spectral elements, Journal of Computational Physics 230: 85738602.

Available online: March 14, 2019.
Publication Dates

Publication in this collection
25 Apr 2019 
Date of issue
2019
History

Received
11 Sept 2018 
Reviewed
21 Feb 2019 
Accepted
13 Mar 2019