A Unified Approach to the Timoshenko Geometric Stiffness Matrix Considering Higher-Order Terms in the Strain Tensor

Abstract

Nonlinear analyses using an updated Lagrangian formulation considering the Euler-Bernoulli 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 higher-order 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 higher-order terms in the strain tensor, as well as the Euler-Bernoulli 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 higher-order terms in the strain tensor results in larger displacements as expected.

Keywords:
Geometric Matrix; Shape Functions; Timoshenko Beam Theory; Higher-Order 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 strain-displacement relations and the interpolating (shape) functions.

The most commonly used bending solution for frame elements is the Euler-Bernoulli 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 beam-column with a moderate slenderness ratio, the Timoshenko beam theory (TBT) provides better results. Several high-order bending theories can be found in the literature, such as in Levinson (1981), Bickford (1982), Heyliger and Reddy (1988), Petrolito (1995), Reddy et al. (1997), Reddy (1997), Tessler and Gherlone (2007) and Meghare and Jadhao (2015). This work considers the Timoshenko beam theory due to its influence on beams with small slenderness, especially when subjected to axial loads (beam-columns).

Regarding the kinematic description, a geometric nonlinear analysis can be performed through the full Lagrangian, updated Lagrangian or a co-rotational approach. Basically, these descriptions differ with respect to the reference configuration considered (Felippa 2017). When consistently developed, the total Lagrangian and the updated formulation produce the same results (McGuire et al. 2000). These formulations are well developed in the literature, such as in Bathe (1996) and Bathe and Bolourchi (1979), in which the nonlinear analysis of three-dimensional structures is studied using both Lagrangian formulations. Nanakorn and Vu (2006) employ a total Lagrangian formulation to analyse a plane element subject to large displacements.

Other studies, such as Santana and Silveira (2014) and Silva et al. (2016), have built a tangent matrix for Timoshenko beam elements using the co-rotational formulation. However, in Santana and Silveira (2014), cubic (Hermitian) shape functions are used, while Silva et al. (2016) use an element proposed by Tang et al. (2015) with a consistent axial displacement field interpolation. Considering a consistent field approach and a co-rotational formulation, Kien (2012) developed a Timoshenko beam plane element for large displacements employing the co-rotational method. Silva and Silva (2010) use an element that includes the theories of Euler-Bernoulli, Timoshenko and third order linear frame analysis. Oliveira and Silva (2017) employ a unified Bernoulli-Timoshenko element.

The strain-displacement relations have an important role in a geometric nonlinear analysis. Usually, the geometric matrix is obtained considering small displacement gradients, neglecting some higher-order terms in the Green strain tensor, as done in McGuire et al. (2000).

In this work, a geometric stiffness matrix is developed considering higher-order 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 Euler-Bernoulli beam theory are well known (Bathe and Bolourchi 1979, Conci 1988, Yang and Leu 1994, Yang and Kuo 1994 and Bathe 1996); 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 beam-column 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 1979, Aguiar et al. 2014).

In some cases, only the Timoshenko elastic matrix is developed; otherwise, the geometric stiffness matrix does not consider higher-order terms in the strain tensor (Davis et al. 1972, Friedman and Kosmatka 1992, Yunhua 1998).

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 high-order finite element with an increase in the order of the polynomial of the basis functions is necessary. Therefore, some authors develop high-order expansion basis functions for structured elements aiming for the accuracy and efficiency of the method (Zheng and Dong 2011, Rodrigues et al. 2016).

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 Euler-Bernoulli 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. (1994) and Pilkey et al. (1995) 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 two-dimensional Timoshenko beam-column elements, considering higher-order 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 2013, Martha and Burgos 2014, 2015). Then, using an updated Lagrangian kinematic description and considering a higher-order 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 higher-order terms in strain tensors to predict the critical loads of framed structures, especially for small slenderness beam-columns 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):

$u\left(x,y\right)=u_0\left(x\right)-\theta \left(x\right)yv\left(x,y\right)=v_0\left(x\right)$ (1)

Figure 1
Beam displacement field

2.2 Differential equilibrium relationships in beam-columns

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.

Figure 2
Equilibrium of an undeformed beam element

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.

$\sum F_y=0\to \frac{dQ\left(x\right)}{dx}=q\left(x\right)\sum M_0=0\to \frac{dM\left(x\right)}{dx}=Q\left(x\right)$ (2)

$M\left(x\right)=EI\frac{d\theta }{dx}\to EI\frac{d^2\theta }{d{x}^2}-Q\left(x\right)=0$ (3)

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 cross-section.

2.3 Euler-Bernoulli beam theory (EBBT)

The Euler-Bernoulli 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).

$EI\frac{d^{4}v\left(x\right)}{d{x}^4}-\frac{dQ\left(x\right)}{dx}=0\to EI\frac{d^{4}v\left(x\right)}{d{x}^4}=q\left(x\right)$ (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 cross-sectional rotation and the transverse displacement are considered to be independent variables.

Figure 3
Shear deformation in the Timoshenko beam theory

According to the TBT, the total cross-section 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):

$\frac{dv}{dx}=\theta +\gamma$ (5)

The shear force acting on the section is given by the following equation:

$Q\left(x\right)=-\chi GA\gamma \left(x\right)$ (6)

where $G$ is the material shear modulus, $A$ is the cross-sectional area, and $\chi$ is the factor that defines the effective area for cross-sectional shear. Substituting equations (5) and (6) in the differential equation (3) of the infinitesimal beam element equilibrium, the following expression can be found:

$EI\frac{d^2\theta }{d{x}^2}+\chi GA\left[\frac{dv\left(x\right)}{dx}-\theta \left(x\right)\right]$ (7)

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), or the one developed in Onu (2008), 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:

$v\left(x\right)=v_b\left(x\right)+v_s\left(x\right)\Rightarrow \frac{dv\left(x\right)}{dx}=\frac{d{v}_b\left(x\right)}{dx}+\frac{d{v}_s\left(x\right)}{dx}=\theta \left(x\right)+\gamma \left(x\right)$ (8)

Using equation (4) and finding the homogeneous solution for Euler-Bernoulli 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:

$EI\frac{d^4{v}_b\left(x\right)}{d{x}^4}=0\to v_b\left(x\right)=c_0+c_{1}x+c_2\frac{x^2}{2}+c_3\frac{x^3}{6}$ (9)

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:

$EI\frac{d^2\theta }{d{x}^2}+\chi GA\left[\frac{dv\left(x\right)}{dx}-\theta \left(x\right)\right]=0\Rightarrow EI\frac{d^3{v}_b\left(x\right)}{d{x}^3}+\chi GA\gamma =0\Rightarrow \gamma =-\frac{EI{c}_3}{\chi GA}=v_s\left(x\right)$ (10)

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 (1997), in this equation.

$v\left(x\right)=c_0+c_{1}x+c_2\frac{x^2}{2}+c_3\left(\frac{x^3}{6}-\Omega L^{2}x\right)\Omega =\frac{EI}{\chi GA}\frac{1}{L^2}$ (11)

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. 2000).

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 (2018), Burgos and Martha (2013) and Onu (2008).

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 2018); therefore, for a nonlinear analysis, bar discretization is necessary. Previous work such as Burgos and Martha (2013) and Onu (2008) 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).

Figure 4
Isolated bar deformed configuration

$\begin{array}{c}{v}_0\left(x\right)=N_2^v\left(x\right)d_2^{\text{'}}+N_3^v\left(x\right)d_3^{\text{'}}+N_5^v\left(x\right)d_5^{\text{'}}+N_6^v\left(x\right)d_6^{\text{'}}\\ \theta \left(x\right)=N_2^{\theta }\left(x\right)d_2^{\text{'}}+N_3^{\theta }\left(x\right)d_3^{\text{'}}+N_5^{\theta }\left(x\right)d_5^{\text{'}}+N_6^{\theta }\left(x\right)d_6^{\text{'}}\end{array}\to \left\{\begin{array}{c}{v}_0\left(x\right)\\ \theta \left(x\right)\end{array}\right\}=\left[N\right]\left\{d\text{'}\right\}$ (12)

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:

$\begin{array}{c}{\nu }_0^h\left(x\right)=c_0+c_{1}x+c_2\frac{x^2}{2}+c_3\left(\frac{x^3}{6}-\Omega L^{2}x\right)\\ {\theta }_0^h\left(x\right)=c_1+c_{2}x+c_3\frac{x^2}{2}\end{array}\to \left\{\begin{array}{c}{v}_0\left(x\right)\\ \theta \left(x\right)\end{array}\right\}=\left[\begin{array}{cc}\begin{array}{ccc}1& x& \frac{x^2}{2}\\ 0& 1& x\end{array}& \begin{array}{c}\frac{x^3}{6}-\Omega L^{2}x\\ \frac{x^2}{2}\end{array}\end{array}\right]\left\{\begin{array}{c}{c}_0\\ c_1\\ \begin{array}{c}{c}_2\\ c_3\end{array}\end{array}\right\}=\left[X\right]\left\{C\right\}$ (13)

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):

$\{d\text{'}\}=\left\{\begin{array}{l}d{\text{'}}_1\\ d{\text{'}}_2\\ d{\text{'}}_3\\ d{\text{'}}_4\end{array}\right\}=\left\{\begin{array}{l}{v}_0(0)\\ \theta (0)\\ v_0(L)\\ \theta (L)\end{array}\right\}\to \{d\text{'}\}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 1& L& \frac{L^2}{2}& \left(\frac{1}{6}-\Omega \right)L^3\\ 0& 1& L& \frac{L^2}{2}\end{array}\right]\left\{\begin{array}{l}{c}_0\\ c_1\\ c_2\\ c_3\end{array}\right\}=[H]\left\{C\right\}$ (14)

Shape functions can be obtained using equations (12), (13) and (14), resulting in:

$\left\{\begin{array}{c}{v}_0\left(x\right)\\ \theta \left(x\right)\end{array}\right\}=\left[X\right]{\left[H\right]}^{-1}\left\{d\text{'}\right\}\Rightarrow \left[N\right]=\left[X\right]{\left[H\right]}^{-1}$ (15)

Finally, Timoshenko beam element shape functions are given by:

$\left\{\begin{array}{c}{N}_2^v\left(x\right)=1+\frac{2{\left(\frac{x}{L}\right)}^3-3{\left(\frac{x}{L}\right)}^2-12\Omega \frac{x}{L}}{\left(1+12\Omega \right)}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{N}_2^{\theta }\left(x\right)=\frac{6}{L}\frac{{\left(\frac{x}{L}\right)}^2-\frac{x}{L}}{\left(1+12\Omega \right)}\\ \begin{array}{c}{N}_3^v\left(x\right)=\frac{x\left[{\left(\frac{x}{L}\right)}^2-2\left(1+3\Omega \right)\frac{x}{L}+1+6\Omega \right]}{\left(1+12\Omega \right)}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{N}_3^{\theta }\left(x\right)=1+\frac{3{\left(\frac{x}{L}\right)}^2-4\left(1+3\Omega \right)\frac{x}{L}}{\left(1+12\Omega \right)}\\ \begin{array}{c}{N}_5^v\left(x\right)=\frac{3{\left(\frac{x}{L}\right)}^2-2{\left(\frac{x}{L}\right)}^3+12\Omega \frac{x}{L}}{\left(1+12\Omega \right)}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{N}_5^{\theta }\left(x\right)=\frac{6}{L}\frac{{\frac{x}{L}-\left(\frac{x}{L}\right)}^2}{\left(1+12\Omega \right)}\\ N_6^v\left(x\right)=\frac{x\left[{\left(\frac{x}{L}\right)}^2-\left(1+6\Omega \right)\frac{x}{L}-6\Omega \right]}{\left(1+12\Omega \right)}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{N}_6^{\theta }\left(x\right)=\frac{3{\left(\frac{x}{L}\right)}^2-2\left(1-6\Omega \right)\frac{x}{L}}{\left(1+12\Omega \right)}\end{array}\end{array}\end{array}\right.$ (16)

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:

$u_0\left(x\right)=N_1\left(x\right)d_1^{\text{'}}+N_4\left(x\right)d_4^{\text{'}},N_1\left(x\right)=1-\frac{x}{L}{N}_4\left(x\right)=\frac{x}{L}$ (17)

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 co-rotational 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. (2000). The final equations are used to find the Timoshenko stiffness matrix considering higher-order terms in the strain tensor.

In this formulation, the equilibrium equations of the unknown configuration $t+\text{Δ}t$ must be written using the known variables of a reference configuration $t$. For a configuration $t+\text{Δ}t$, the virtual work of the internal forces must be equal to the virtual work of the external forces, as shown in equation (18):

${\int }_V^{}{S}_{ij}^{\left(t+∆t\right)}\delta {\epsilon }_{ij}^{\left(t+∆t\right)}dV=R^{\left(t+∆t\right)}$ (18)

where $S_{ij}^{\left(t+\Delta t\right)}$ corresponds to the second Piola-Kirchoff stress tensor, ${\epsilon }_{ij}^{\left(t+∆t\right)}$ the Green-Lagrange strain tensor, and $R^{\left(t+∆t\right)}$ the virtual work due to external loading, which could include body, surface, inertia and damping forces (Aguiar et al. 2014).

The linearized incremental equation requires small displacement increments; thus, the second Piola-Kirchoff stress tensor and the Green-Lagrange strain tensor can be written according to:

$S_{ij}^{\left(t+∆t\right)}={\tau }_{ij}^t+∆{\tau }_{ij}\mathrm{}{\epsilon }_{ij}^{\left(t+∆t\right)}={\epsilon }_{ij}^t+∆{\epsilon }_{ij}$ (19)

where ${\tau }_{ij}^t$ corresponds to the Cauchy stress tensor, $∆{\tau }_{ij}$ is the stress increment and $∆{\epsilon }_{ij}$ is the deformation increment, which can be calculated from the Green-Lagrange 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).

${\int }_V^{}{S}_{ij}^{\left(t+∆t\right)}\delta {\epsilon }_{ij}^{\left(t+∆t\right)}dV={\int }_V^{}\left({\tau }_{ij}^t+∆{\tau }_{ij}\right)\delta \left(∆{\epsilon }_{ij}\right)dV={\int }_V^{}∆{\tau }_{ij}\delta ∆{\epsilon }_{ij}dV+{\int }_V^{}{\tau }_{ij}^t\delta ∆{\epsilon }_{ij}dV$ (20)

The Green-Lagrange strain tensor in the xy plane is represented by:

${\epsilon }_{xx}=\frac{\partial u}{\partial x}+\frac{1}{2}\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial x}\right)}^2\right]{\gamma }_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial v}{\partial y}$ (21)

The strain tensor shown in relation (21) has a linear ($∆e_{ij}$) and a nonlinear part ($∆{\eta }_{ij}$); thus, the decomposition $∆{\epsilon }_{ij}=∆e_{ij}+∆{\eta }_{ij}$ can be adopted, as in equation (22):

$∆e_{xx}=\frac{\partial \mathrm{u}}{\partial \mathrm{x}}∆{\eta }_{xx}=\frac{1}{2}\left[{\left(\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{v}}{\partial \mathrm{x}}\right)}^2\right]∆e_{xy}=\frac{\partial \mathrm{u}}{\partial \mathrm{y}}+\frac{\partial v}{\partial \mathrm{x}}∆{\mathrm{\eta }}_{xy}=\frac{\partial \mathrm{u}}{\partial \mathrm{x}}\frac{\partial u}{\partial \mathrm{y}}+\frac{\partial \mathrm{v}}{\partial \mathrm{x}}\frac{\partial \mathrm{v}}{\partial \mathrm{y}}$ (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:

$∆{\tau }_{ij}=C_{ijkl}∆{\epsilon }_{kl}=C_{ijkl}∆e_{kl}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}∆{\epsilon }_{ij}=∆e_{ij}$ (23)

Finally, based on equations (20) and (23), the virtual work equation becomes:

$\begin{array}{l}\underset{V}{\overset{}{{\displaystyle \int }}}{C}_{ijkl}\Delta e_{kl}\delta \Delta e_{ij}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{ij}^t\delta \left(\Delta e_{ij}+\Delta {\eta }_{ij}\right)dV=R^{\left(t+\Delta t\right)}\\ \underset{V}{\overset{}{{\displaystyle \int }}}{C}_{ijkl}\Delta e_{kl}\delta \Delta e_{ij}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{ij}^t\delta \Delta e_{ij}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{ij}\Delta {\eta }_{ij}dV=R^{\left(t+\Delta t\right)}\end{array}$ (24)

In the xy plane, the stress vector, the constitutive matrix, and the linear and non-linear strain vectors are given by:

$\tau =\left\{\begin{array}{c}{\tau }_{xx}\\ {\tau }_{xy}\end{array}\right\}C=\left[\begin{array}{cc}E& 0\\ 0& G\end{array}\right]e=\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\gamma }_{xy}\end{array}\right\}\eta =\left\{\begin{array}{c}{\eta }_{xx}\\ {\eta }_{xy}\end{array}\right\}$ (25)

Thus, the virtual work equation, expression (24), can be expanded according to equation (26), using the relations presented next.

$\begin{array}{l}\begin{array}{c}\underset{V}{\overset{}{{\displaystyle \int }}}{C}_{ijkl}\Delta e_{kl}\delta \Delta e_{ij}dV=\underset{V}{\overset{}{{\displaystyle \int }}}{\epsilon }_{xx}.E\delta {\epsilon }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\gamma }_{xy}.G\delta {\gamma }_{xy}dV\\ \underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{ij}^t\delta \Delta e_{ij}dV=\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xx}\delta {\epsilon }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xy}\delta {\gamma }_{xy}dV\\ \underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{ij}\Delta {\eta }_{ij}dV=\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xx}\delta {\eta }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xy}\delta {\eta }_{xy}dV\end{array}\\ \underset{V}{\overset{}{{\displaystyle \int }}}{\epsilon }_{xx}.E\delta {\epsilon }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\gamma }_{xy}.G\delta {\gamma }_{xy}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xx}\delta {\epsilon }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xy}\delta {\gamma }_{xy}dV+\\ +\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xx}\delta {\eta }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xy}\delta {\eta }_{xy}dV=R^{\left(t+\Delta t\right)}\end{array}$ (26)

5 STIFFNESS MATRICES

The stiffness matrix development is performed by analysing a two-dimensional structure in the xy plane. Using an updated Lagrangian formulation, the higher-order 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 Green-Lagrange strain tensor components in equation (21), can be rewritten as expressions (27) and (28), respectively:

${\epsilon }_{xx}=\frac{\partial u}{\partial x}=\frac{\partial u_0}{\partial x}-y{\frac{\partial {\theta }_Z}{\partial x}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{}^{}{\gamma }_{xy}=\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}=\frac{\partial v_0}{\partial x}-{\theta }_Z$ (27)

${\eta }_{xx}=\frac{1}{2}\left({\frac{\partial u}{\partial x}}^2+{\frac{\partial v}{\partial x}}^2\right)=\frac{1}{2}\left({\frac{\partial u}{\partial x}}^2+{\frac{\partial v}{\partial x}}^2+y^2\frac{\partial {\theta }_z{}^2}{\partial x}\right)-y\frac{\partial u}{\partial x}{\frac{\partial {\theta }_z}{\partial x}}^{}{}^{}{\eta }_{xy}=\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial v}{\partial y}=y\frac{\partial {\theta }_z}{\partial x}{\theta }_z-{\theta }_z\frac{\partial u}{\partial x}$ (28)

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.

$\begin{array}{l}\delta W_1=\underset{V}{\overset{}{{\displaystyle \int }}}{\epsilon }_{xx}.E\delta {\epsilon }_{xx}dV={\displaystyle \underset{A}{\int }\left({\displaystyle \underset{0}{\overset{L}{\int }}\left(\frac{\partial u}{\partial x}-y\frac{\partial {\theta }_z}{\partial x}\right)E\left(\delta \frac{\partial u}{\partial x}-\delta y\frac{\partial {\theta }_z}{\partial x}\right)}dx\right)dA}\\ \delta W_1=\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial u}{\partial x}\delta \frac{\partial u}{\partial x}}dx\right)E{\displaystyle \underset{A}{\int }dA}+\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial {\theta }_z}{\partial x}\delta \frac{\partial {\theta }_z}{\partial x}}dx\right)E{\displaystyle \underset{A}{\int }{y}^{2}dA}-\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial {\theta }_z}{\partial x}\delta \frac{\partial u}{\partial x}}dx\right)E{\displaystyle \underset{A}{\int }ydA}-\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial u}{\partial x}\delta \frac{\partial {\theta }_z}{\partial x}}dx\right)E{\displaystyle \underset{A}{\int }ydA}\\ \end{array}$ (29)

$\begin{array}{l}\delta W_2=\underset{V}{\overset{}{{\displaystyle \int }}}{\gamma }_{xy}.G\delta {\gamma }_{xy}dV={\displaystyle \underset{A}{\int }\left({\displaystyle \underset{0}{\overset{L}{\int }}\left(\frac{\partial v}{\partial x}-{\theta }_Z\right)G\left(\delta \frac{\partial v}{\partial x}-\delta {\theta }_Z\right)}dx\right)dA}\\ \delta W_2=\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial v}{\partial x}\delta \frac{\partial v}{\partial x}}dx\right)G{\displaystyle \underset{A}{\int }dA}+\left({\displaystyle \underset{0}{\overset{L}{\int }}{\theta }_Z\delta {\theta }_Z}dx\right)G{\displaystyle \underset{A}{\int }dA}-\left({\displaystyle \underset{0}{\overset{L}{\int }}{\theta }_Z\delta \frac{\partial v}{\partial x}}dx\right)G{\displaystyle \underset{A}{\int }dA}-\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial v}{\partial x}\delta {\theta }_Z}dx\right)G{\displaystyle \underset{A}{\int }dA}\end{array}$ (30)

In the beam centroidal axis, ${\displaystyle \underset{A}{\int }{y}^{2}dA=I_z,}{\displaystyle \underset{A}{\int }ydA=0},$ and the expressions (29) and (30) are reduced to equations (31) and (32):

$\delta W_1=\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial u}{\partial x}\delta \frac{\partial u}{\partial x}}dx\right)EA+\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial {\theta }_Z}{\partial x}\delta \frac{\partial {\theta }_Z}{\partial x}}dx\right)E{I}_z$ (31)

$\delta W_2=\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial v}{\partial x}\delta \frac{\partial v}{\partial x}}dx\right)GA+\left({\displaystyle \underset{0}{\overset{L}{\int }}{\theta }_Z\delta {\theta }_Z}dx\right)GA-\left({\displaystyle \underset{0}{\overset{L}{\int }}{\theta }_Z\delta \frac{\partial v}{\partial x}}dx\right)GA-\left({\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial v}{\partial x}\delta {\theta }_Z}dx\right)GA$ (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 $\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:

$u_0\left(x\right)=\left\{N_u\left(x\right)\right\}\left\{u\right\}{v}_0\left(x\right)=\left\{N_v\left(x\right)\right\}\left\{v\right\}{\theta }_z\left(x\right)=\left\{N_{\theta z}\left(x\right)\right\}\left\{v\right\}$ (33)

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{'}$.

$\delta W_1={\left\{\delta u\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}EA\left\{N_u\text{'}\right\}}{\left\{N_u\text{'}\right\}}^{T}dx\left\{u\right\}+{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}E{I}_z\left\{N_{\theta z}\text{'}\right\}}{\left\{N_{\theta z}\text{'}\right\}}^{T}dx\left\{v\right\}$ (34)

$\begin{array}{l}\delta W_2={\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}GA\left\{N_v\text{'}\right\}}{\left\{N_v\text{'}\right\}}^{T}dx\left\{v\right\}+{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}GA\left\{N_{\theta z}\right\}}{\left\{N_{\theta z}\right\}}^{T}dx\left\{v\right\}\\ {{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{}^{}-{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}GA\left\{N_{\theta z}\right\}}{\left\{N_v\text{'}\right\}}^{T}dx\left\{v\right\}-{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}GA\left\{N_v\text{'}\right\}}{\left\{N_{\theta z}\right\}}^{T}dx\left\{v\right\}\end{array}$ (35)

Solving for the integrals in equations (34) and (35), the elastic stiffness matrix can be calculated, resulting in the matrix (36):

$K_{e,xy}=\left[\begin{array}{cccccc}\frac{EA}{L}& 0& 0& -\frac{EA}{L}& 0& 0\\ 0& \frac{12E{I}_Z}{L^3\left(12{\Omega }_Y+1\right)}& \frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}& 0& -\frac{12E{I}_Z}{L^3\left(12{\Omega }_Y+1\right)}& \frac{6E{I}_Z}{L^2\left(12\Omega +1\right)}\\ 0& \frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}& \frac{4E{I}_Z\left(3{\Omega }_Y+1\right)}{L\left(12{\Omega }_Y+1\right)}& 0& -\frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}& \frac{2E{I}_Z\left(1-6{\Omega }_Y\right)}{L\left(12{\Omega }_Y+1\right)}\\ -\frac{EA}{L}& 0& 0& \frac{EA}{L}& 0& 0\\ 0& -\frac{12E{I}_Z}{L^3\left(12{\Omega }_Y+1\right)}& -\frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}& 0& \frac{12E{I}_Z}{L^3\left(12{\Omega }_Y+1\right)}& -\frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}\\ 0& \frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}& \frac{2E{I}_Z\left(1-6{\Omega }_Y\right)}{L\left(12{\Omega }_Y+1\right)}& 0& -\frac{6E{I}_Z}{L^2\left(12{\Omega }_Y+1\right)}& \frac{4E{I}_Z\left(3{\Omega }_Y+1\right)}{L\left(12{\Omega }_Y+1\right)}\end{array}\right]$ (36)

5.3 Geometric stiffness matrix

From the beam displacement field, equation (1), and the Green-Lagrange strain tensor, equation (21), the nonlinear part of this tensor can be written according to:

${\eta }_{xx}=\frac{1}{2}\left({\frac{\partial u}{\partial x}}^2+{\frac{\partial v}{\partial x}}^2+y^2\frac{\partial {\theta }_z{}^2}{\partial x}\right)-y\frac{\partial u}{\partial x}\frac{\partial {\theta }_z}{\partial x}{\eta }_{xy}=y\frac{\partial {\theta }_z}{\partial x}{\theta }_z-{\theta }_z\frac{\partial u}{\partial x}$ (37)

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. (2000):

${\eta }_{xx}=\frac{1}{2}{\left({\frac{\partial u}{\partial x}}^2+{\frac{\partial v}{\partial x}}^2\right)}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{}^{}{\eta }_{xy}=0$ (38)

In this work, however, the complete strain tensor is used (equation (37)), as in Yang and Kuo (1994), considering the Timoshenko and not the Euler-Bernoulli 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).

$\begin{array}{l}\delta W_{NL}=\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xx}\delta {\eta }_{xx}dV+\underset{V}{\overset{}{{\displaystyle \int }}}{\tau }_{xy}\delta {\eta }_{xy}dV\\ \delta W_{NL}={\displaystyle \underset{A}{\int }\left({\displaystyle \underset{0}{\overset{L}{\int }}{t}_{xx}\delta \left(\frac{1}{2}\left({\frac{\partial u}{\partial x}}^2+{\frac{\partial v}{\partial x}}^2+y^2{\frac{\partial {\theta }_z}{\partial x}}^2\right)-y\frac{\partial {\theta }_z}{\partial x}\frac{\partial u}{\partial x}\right)}dx\right)dA}+{\displaystyle \underset{A}{\int }\left({\displaystyle \underset{0}{\overset{L}{\int }}{t}_{xy}\delta \left(y\frac{\partial {\theta }_z}{\partial x}{\theta }_z-{\theta }_z\frac{\partial u}{\partial x}\right)}dx\right)dA}\\ \delta W_{NL}=\frac{1}{2}\left({\displaystyle \underset{0}{\overset{L}{\int }}\delta {\frac{\partial u}{\partial x}}^2}dx\right){\displaystyle \underset{A}{\int }{t}_{xx}dA}+\frac{1}{2}\left({\displaystyle \underset{0}{\overset{L}{\int }}\delta {\frac{\partial v}{\partial x}}^2}dx\right){\displaystyle \underset{A}{\int }{t}_{xx}dA}+\frac{1}{2}\left({\displaystyle \underset{0}{\overset{L}{\int }}\delta {\frac{\partial {\theta }_z}{\partial x}}^2}dx\right){\displaystyle \underset{A}{\int }{y}^2{t}_{xx}dA}-\left({\displaystyle \underset{0}{\overset{L}{\int }}\delta \frac{\partial {\theta }_z}{\partial x}\frac{\partial u}{\partial x}}dx\right){\displaystyle \underset{A}{\int }{t}_{xx}ydA}\\ {{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}+\left({\displaystyle \underset{0}{\overset{L}{\int }}\delta \frac{\partial {\theta }_z}{\partial x}{\theta }_z}dx\right){\displaystyle \underset{A}{\int }{t}_{xy}ydA}-\left({\displaystyle \underset{0}{\overset{L}{\int }}\delta {\theta }_z\frac{\partial u}{\partial x}}dx\right){\displaystyle \underset{A}{\int }{t}_{xy}dA}\end{array}$ (39)

Applying relations, ${\displaystyle \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).

$\delta W_{NL}=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}\left[P\delta \left({\frac{\partial u}{\partial x}}^2+{\frac{\partial v}{\partial x}}^2\right)+P\frac{I_z}{A}\delta \left({\frac{\partial {\theta }_z}{\partial x}}^2\right)\right]}dx+{\displaystyle \underset{0}{\overset{L}{\int }}\left[M_z\delta \left(\frac{\partial {\theta }_z}{\partial x}\frac{\partial u}{\partial x}\right)-Q_y\delta \left({\theta }_z\frac{\partial u}{\partial x}\right)\right]}dx$ (40)

$\begin{array}{l}{•}^{}{}^{}{\displaystyle \underset{0}{\overset{L}{\int }}\left[P\delta \left(u{\text{'}}^2+v{\text{'}}^2\right)+P\frac{I_z}{A}\delta \left({\theta }_z{\text{'}}^2\right)\right]}dx=\\ ={\left\{\delta u\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}P\left\{N_u\text{'}\right\}}{\left\{N_u\text{'}\right\}}^{T}dx\left\{u\right\}+{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}P\left\{N_v\text{'}\right\}}{\left\{N_v\text{'}\right\}}^{T}dx\left\{v\right\}+{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}P\frac{I_z}{A}\left\{N_{\theta z}\text{'}\right\}}{\left\{N_{\theta z}\text{'}\right\}}^{T}dx\left\{v\right\}\end{array}$ (41)

$\begin{array}{l}{•}^{}{}^{}{\displaystyle \underset{0}{\overset{L}{\int }}\left[M_z\delta \left({\theta }_z\text{'}u\text{'}\right)-Q_y\delta \left({\theta }_{z}u\text{'}\right)\right]}dx=\\ ={\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}{M}_z\left\{N_{\theta }\text{'}\right\}}{\left\{N_u\text{'}\right\}}^{T}dx\left\{u\right\}+{\left\{\delta u\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}{M}_z\left\{N_u\text{'}\right\}}{\left\{N_{\theta }\text{'}\right\}}^{T}dx\left\{v\right\}-{\left\{\delta v\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}{Q}_y\left\{N_{\theta }\right\}}{\left\{N_u\text{'}\right\}}^{T}dx\left\{u\right\}-\\ {{}^{}}^{}{{}^{}}^{}{{}^{}}^{}-{\left\{\delta u\right\}}^T{\displaystyle \underset{0}{\overset{L}{\int }}{Q}_y\left\{N_u\text{'}\right\}}{\left\{N_{\theta }\right\}}^{T}dx\left\{v\right\}\end{array}$ (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).

Figure 5
Frame element - (adapted from Pereira 2002)

$M_Z=-M_{Z1}+{\frac{\left(M_{Z1}+M_{Z2}\right)x}{L}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{}^{}{Q}_y=-\frac{\left(M_{Z1}+M_{Z2}\right)}{L}$ (43)

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 higher-order terms in strain tensor, shown by matrix (44) below:

$\begin{array}{l}{K}_{xy}^g=[\begin{array}{ccc}\frac{P}{L}& 0& -\frac{M_{Z1}}{L}\\ 0& \frac{6P\left(120{\Omega }_Y{}^2+20{\Omega }_Y+1\right)}{5L{\left(12{\Omega }_Y+1\right)}^2}+\frac{12P{I}_Z}{A{L}^3{\left(12{\Omega }_Y+1\right)}^2}& \frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}+\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}\\ -\frac{M_{z1}}{L}& \frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}+\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}& \frac{2LP\left(90{\Omega }_Y{}^2+15{\Omega }_Y+1\right)}{15{\left(12{\Omega }_Y+1\right)}^2}+\frac{4P{I}_Z\left(36{\Omega }_Y{}^2+6{\Omega }_Y+1\right)}{AL{\left(12{\Omega }_Y+1\right)}^2}\\ -\frac{P}{L}& 0& \frac{M_{Z1}}{L}\\ 0& -\frac{6P\left(120{\Omega }_Y{}^2+20{\Omega }_Y+1\right)}{5L{\left(12{\Omega }_Y+1\right)}^2}-\frac{12P{I}_Z}{A{L}^3{\left(12{\Omega }_Y+1\right)}^2}& -\frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}-\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}\\ -\frac{M_{z2}}{L}& \frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}+\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}& -\frac{LP\left(360{\Omega }_Y{}^2+60{\Omega }_Y+1\right)}{30{\left(12{\Omega }_Y+1\right)}^2}-\frac{2P{I}_Z\left(72{\Omega }_Y{}^2+12{\Omega }_Y-1\right)}{AL{\left(12{\Omega }_Y+1\right)}^2}\end{array}\\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ {}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}{{}^{}}^{}\begin{array}{ccc}-\frac{P}{L}& 0& -\frac{M_{Z2}}{L}\\ 0& -\frac{6P\left(120{\Omega }_Y{}^2+20{\Omega }_Y+1\right)}{5L{\left(12{\Omega }_Y+1\right)}^2}-\frac{12P{I}_Z}{A{L}^3{\left(12{\Omega }_Y+1\right)}^2}& \frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}+\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}\\ \frac{M_{Z1}}{L}& -\frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}-\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}& -\frac{LP\left(360{\Omega }_Y{}^2+60{\Omega }_Y+1\right)}{30{\left(12{\Omega }_Y+1\right)}^2}-\frac{2P{I}_Z\left(72{\Omega }_Y{}^2+12{\Omega }_Y-1\right)}{AL{\left(12{\Omega }_Y+1\right)}^2}\\ \frac{P}{L}& 0& \frac{M_{z2}}{L}\\ 0& \frac{6P\left(120{\Omega }_Y{}^2+20{\Omega }_Y+1\right)}{5L{\left(12{\Omega }_Y+1\right)}^2}+\frac{12P{I}_Z}{A{L}^3{\left(12{\Omega }_Y+1\right)}^2}& -\frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}-\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}\\ \frac{M_{z2}}{L}& -\frac{P}{10{\left(12{\Omega }_Y+1\right)}^2}-\frac{6P{I}_Z}{A{L}^2{\left(12{\Omega }_Y+1\right)}^2}& \frac{2LP\left(90{\Omega }_Y{}^2+15{\Omega }_Y+1\right)}{15{\left(12{\Omega }_Y+1\right)}^2}+\frac{4P{I}_Z\left(36{\Omega }_Y{}^2+6{\Omega }_Y+1\right)}{AL{\left(12{\Omega }_Y+1\right)}^2}\end{array}]\end{array}$ (44)

In this work, this matrix is called TBT_large. Without the consideration of shear strain, i.e.$,GA\chi \to \infty$, the problem is solved by Euler-Bernoulli beam theory. Considering higher-order terms in the strain tensor in this work this matrix is called EBBT_large and is the same as developed in Yang and Leu (1994).

In this paper, matrices considering small displacement gradients were also used for comparison. These matrices use the Euler-Bernoulli beam theory and are developed in many references, such as McGuire et al. (2000). 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 higher-order 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 (2013) 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 Euler-Bernoulli and Timoshenko beam theories. All stiffness matrices were implemented in the Ftool computer software (Martha 1999), and their results were compared with results developed in the literature and with the Mastan2 v3.5 software developed by McGuire et al. (2000), by considering both beam theories, i.e., EBBT_Mastan (Euler-Bernoulli 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 non-linear analysis and the influence of considering higher-order 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$.

Figure 6
Analysed Columns - (adapted from Silva et al. 2016)

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$.

Figure 7
Equilibrium paths for a clamped column

Analysing the equilibrium paths, it can be observed that for high values of slenderness, the beam behaviour is similar for both the Euler-Bernoulli and Timoshenko theories. No relevant differences are observed with or without the consideration of higher-order 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 higher-order terms in the strain tensor. Thus, for high loads and small slenderness, it is important to consider Timoshenko beam theory and higher-order 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 higher-order 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.

Figure 8
Equilibrium paths for a simply supported column

The same conclusions can be observed from these equilibrium paths: for high values of slenderness ratios, the behaviour considering Euler-Bernoulli and Timoshenko theories is similar. Additionally, relevant differences are observed regarding the consideration of higher-order terms in the strain tensor. However, as the slenderness decreases, $\lambda <4.0$, the influence of the beam theory and the consideration of higher-order 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.

Figure 9
Equilibrium paths for a fixed and simply supported column

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 higher-order 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 higher-order 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 higher-order 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 (1963) and is shown in Figure 10.

Figure 10
Continuous beam column - (adapted from Timoshenko and Gere 1963)

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 non-linear analysis and the influence of the consideration of higher-order 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 (1963), 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.

Figure 11
Equilibrium paths for continuous beam column, for λ = 10 and λ = 6.6

Figure 12
Equilibrium paths for continuous beam column, for λ = 4.0 and λ = 2.0

As in example 1, for high values of slenderness ratios, the beam behaviour considering the Euler-Bernoulli and Timoshenko theories is similar, and the equilibrium paths overlap. No relevant differences are observed with or without the consideration of higher-order 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 higher-order 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 higher-order 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:

$P_{cr}=-1.40694\frac{{\pi }^{2}EI}{L^2}$ (45)

However, analytically, the critical load considering shear deformation can be obtained by solving the following equation (Burgos and Martha 2013):

$\left[{\mu }^2{L}^2\left(1+6\Omega \right)+3\right]\sin \left(\mu L\right)-3\mu L\cos \left(\mu L\right)=\mathrm{0,}{\mu }^2=-P_{cr}/EI\left(1+P_{cr}/G\chi A\right)$ (46)

Burgos and Martha (2013) 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 beam-columns 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.

Figure 13
Roorda’s frame

Figure 14
Roorda’s Frame equilibrium path for different values of $\Omega$

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 (2013). 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 Euler-Bernoulli 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 higher-order 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 higher-order 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 (2010) present a solution for this problem, considering the Euler-Bernoulli beam theory, based on Euler’s critical load according to:

${\mu }_{cr}P=0.744\frac{{\pi }^{2}EI}{L^2}$ (47)

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.

Figure 15
Fixed Frame - (adapted from Bazant and Cedolin 2010)

Figure 16
Frame equilibrium path

As expected, for large slenderness ratios, the results using different beam theories are identical, and the consideration of higher-order 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 higher-order 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 higher-order 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 Euler-Bernoulli 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 Euler-Bernoulli 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 higher-order 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 higher-order 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 Euler-Bernoulli 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 higher-order 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 three-dimensional multilayered pipe beam element: nonlinear analysis, Computers and Structures 138: 142-161.
• Bathe, K.J. (1996). Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, USA.
• Bathe, K.J., Bolourchi, S. (1979). Large displacement analysis of three-dimensional beam structures, International Journal for Numerical Methods in Engineering 14: 961-986.
• 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 beam-columns 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ão-linearidades geométrica e material, Ph.D. Dissertation, Department of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, PUC-Rio, Rio de Janeiro, Brazil.
• Davis, R., Henshell, R.D., Warburton G.B. (1972). A Timoshenko beam element, Journal of Sound and Vibration 22 (4): 475-487.
• Felippa, C.A., (2017), Non-linear finite element methods / NFEM, Lecture notes for the course non-linear finite element methods, Center for Aerospace Structures, University of Colorado, Boulder/USA.
• Friedman, Z., Kosmatka, J.B. (1992). An improved two-node Timoshenko beam finite element, Computers and Structures 47 (3): 473-481.
• Heyliger, P.R., Reddy, J.N. (1988). A higher order beam finite element for bending and vibration problems, Journal of Sound and Vibration 126: 309-326.
• 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: 81-87.
• 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, 51-65.
• 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: 31-38.
• Nanakorn, P., Vu, L.N. (2006). A 2D field-consistent beam element for large displacement analysis using the total Lagrangian formulation, Finite Elements in Analysis and Design 42: 1240-1247.
• Oliveira, G.C., Silva, W.T.M. (2017). Análise não-linear de arcos utilizando o elemento de viga unificado Bernoulli-Timoshenko e a formulação co-rotacional, Civil Engineering Electronic Journal. 13(2): 1-16.
• Onu, G. (2008). Finite elements on generalized elastic foundation in Timoshenko beam theory, Journal of Engineering Mechanics 134(9): 763-776.
• 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, PUC-Rio, Rio de Janeiro, Brazil.
• Petrolito, J. (1995). Stiffness analysis of beams using a higher-order theory, Computers and Structures 55: 33-39.
• 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: 25-44.
• Reddy, J.N. (1997). On locking-free shear deformable beam finite elements, Computer Methods in Applied Mechanics and Engineering 149: 113-132.
• 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): 3373-3384.
• Rodrigues, C. F., Suzuki, J. L., Bittencourt, M. L. (2016). Construction of minimum energy high-order Helmholtz bases for structured elements, Journal of Computational Physics 306: 269-290.
• 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: 141-162.
• Shrima, L.M., Giger, M.W. (1992). Timoshenko beam element resting on two-parameter elastic foundation, Journal of Engineering Mechanics 118(2): 280-295.
• 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 Beam-Column 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 TP-2007-215086).
• Timoshenko, S.P., Gere, J.M (1963). Theory of elastic stability, International student edition, second edition, McGraw-Hill, Singapore.
• Yang, Y. B., Leu, L.J. (1994). Non-linear stiffnesses in analysis of planar frames, Computer Methods in Applied Mechanics and Engineering 117: 233-247.
• 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: 249-269.
• Zheng, X., Dong, S. (2011). An eigen-based high-order expansion basis for structured spectral elements, Journal of Computational Physics 230: 8573-8602.

• 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