Utilizing Tait-Bryan Angles for Large Displacement Corotational Finite Element Static Analysis of Spatial Beams

Abstract

In this work, a corotational finite element formulation is suggested for spatial beams with geometrically nonlinear behavior subjected to static loads. We returned to the three successive rotation angle procedure, mainly the Tait-Bryan angles. By carefully defining the trigonometric rules for all rotation angles, the singularity problem, that had limited the use of these angles, is avoided. Three different types of coordinate systems are used: a fixed global coordinate system that stays fixed throughout the analysis, a fixed local coordinate system that is fixed and precisely attached to each element, and a corotational local frame for each element that moves and rotates together with the element throughout the analysis. The rigid body motion can easily be separated from the overall deformation since the deformation is always tiny relative to the corotational frame. An incremental-iterative method is used for the solution based upon the Newton-Raphson method. Different examples are solved to demonstrate the practicality, correctness, and accuracy of the proposed method. The solutions converge at a relatively quick rate.

Keywords:
Corotational formulation; geometrically nonlinear analysis; large displacement analysis; spatial beams; Tait Bryan angles

Graphical Abstract

1 INTRODUCTION

Enhancing geometrically nonlinear finite element formulations for spatial beam elements has recently become a crucial objective (Santana et al., 2022; Vo and Nanakorn, 2020; Magisano et al., 2020). It makes sense that many innovative applications (Leng et al., 2022; Trapper, 2022; Liu, 2014; Xiaohang et al., 2022; Liu and Bai, 2022) that have large displacements require accurate modelling. Leng et al. (2022) insist that geometric nonlinearity has a considerable impact on flexible offshore structures and devices and cannot be ignored. At the same time, it has become urgently necessary to increase the size of renewable energy devices as they experience extraordinary growth. As a result, Xiaohang et al. (2022) studied a flexible wind turbine blade that was 100 meters long. This study validates how crucial geometric nonlinearity is, in the analysis of such a large blade. Even in space applications, Liu and Bai (2022) conducted experimental and numerical analyses for a deployable cabin that can be used as space habitats. Most of these innovative applications have been treated using commercial finite element programs. Hence, there is a snowballing need for computationally inexpensive geometrically nonlinear finite element formulations to examine and improve such commercial programs. Such improved and flexible formulations can also deal with various engineering problems that cannot be solved, and diverge, in some cases, when using commercial finite element programs.

One pillar of geometrically nonlinear analysis is the kinematics description approach. Three formulations are usually applied to describe the kinematics of a spatial beam element. These formulations are the total Lagrangian, the updated Lagrangian, and the corotational formulations. The total Lagrangian formulation (Santana et al., 2022; Vo and Nanakorn, 2020; Mars et al., 2017; Crivelli, 1991; Simo and Vu-Quoc, 1986; Remseth, 1979) is the oldest formulation and is frequently used in commercial finite element softwares. This formulation uses the fixed global frame to define the terms of the system equation. This fixed global frame does not update through the analysis. Consequently, it causes relatively large displacements, rotations, and strains that require special techniques to handle. The updating Lagrangian formulation (Yang et al., 2002; Yang et al., 2007) employs a frame which is updated with the last accepted solution to define the terms of the system equation. Thus, the system equations are sampler than the corresponding equations in total Lagrangian formulations. The reference frame does not change during the solution cycles; hence, if the displacement from the current configuration relative to the last equilibrium configuration is large, a basic assumption is violated, and for that reason this formulation also experiences some intricacies. To avoid these complications, the corotational formulation presents an efficient kinematics description approach for large displacement analysis. The small strain theory is employed in this formulation (Belytschko and Hsieh, 1973; Elkaranshawy and Dokainish, 1995; Elkaranshawy et al., 2018; Gu, 2004; Oran, 1973; Crisfield, 1990; Le et al., 2014; Jonker and Meijaard, 2013; Bathe and Bolourchi, 1979; Benjamin, 1982; Nunes et al., 2003). With each element, the corotational frame rotates and translates with it, but it does not deform with it. This frame is continuously updated and can be used in a variety of ways (Gu, 2004). Because the deformational motion is always small with respect to the moving corotational frame, this formulation has the distinct advantage of making it easy to separate the rigid body motion from that motion.

Due to the complication of large spatial rotations, the three-dimensional formulation is not just a simple extension of the planar formulation. For instance, Remseth (1979) applied an approximate vectorial hypothesis to treat three-dimensional rotations, and he restricted his method to moderate rotations that do not exceed 15 degrees. Therefore, his method cannot be applied for large rotations. Oran (1973) utilized orthogonal axes that are rigidly attached and deformed with each joint. A joint orientation matrix was defined and used to describe these axes. The angles between the member axes and this set of orthogonal axes are used to compute the nodal rotations for each element. This procedure was improved by Crisfield (1990), Le et al. (2014), and Jonker and Meijaard (2013). The requirement for a special parametrization of the finite rotations is the cornerstone of this method. The joint orientation matrices that are stored and the parametrization of the finite rotations considerably increase the computational time needed for this method. A moving local frame with three successive rotations, similar to the Euler angle, was presented by Bathe and Bolourchi (1979). However, they did not provide trigonometric rules for all rotation angles. Benjamin (1982) adds to this method by stating rotation angles' cosines and sins in terms of kinematic variables. However, he did not determine a control sign for the cosine of the rotation angles that can be obtained using the two hypotenuses of right-angle triangles. To overcome the problem of cosine of an angle outside the interval [π/2, -π/2], Nunes et al. (2003) controlled the sign of cosine of rotation angles. However, they did not specify clearly how to determine the transformation procedure in the case of vertical members, which is very important in the modelling of three-dimensional structures. They did not analyze three-dimensional problems to thoroughly test their motion description method. Simo and Vu-Quoc (1986) pointed out the problem of singularity, in the case of adopting this method, which has to be handled carefully.

This study proposes a relatively accurate and simple corotational finite element formulation for statically loaded space frames. The material is assumed to be elastic and isotropic. The beam element cross section is uniform, and Bernoulli's hypothesis is assumed. The cross-sectional distortion, shear, and warping effect are not considered. The transformation procedure involves a regular update of the coordinates' vector with each equilibrium configuration during the analysis. In the transformation technique, three successive rotations of Tait-Bryan angles (Gu, 2004; Bathe and Bolourchi, 1979) are employed. The proposed approach uses two main steps to transform vectors and matrices from the fixed global frame to the local corotational frame. The first step is to transform from the fixed global frame to the fixed local frame, and the second step is to transform from the fixed local frame to the moving corotational local frame. All trigonometric rules of the spatial beam element with control signs are expressed, including special cases. These trigonometric rules are then used to define the rotation matrices and the relative displacement vector. Although the suggested formulation is simple and does not require special parametrization of finite rotations, it requires adjusting the load step and number of elements to reduce the relative chordal rotations that occur throughout the analysis. The system's equilibrium equation is derived using the virtual work principle. An incremental iterative procedure based on the full Newton-Raphson method is employed to solve this equation through a MATLAB code. Due to bypassing the joint orientation matrices and parametrizing of finite rotations, this code has a comparatively rapid convergence rate for equilibrium.

This section serves as an overview and highlights the value of researching geometric nonlinearity. Section 2 introduces the spatial beam element motion description approach, which involves coordinate systems. Section 3 illustrates the method of transformation between the used coordinate systems based on Tait-Bryan angles. In Section 4, the stiffness matrix and the strain energy are derived. The equilibrium equation is then derived in Section 5 using the principle of virtual work. The numerical algorithm is shown in Section 6. To show the efficacy and validate the precision of the suggested method, five numerical examples are solved and compared with the published results in Section 7. Conclusions are discussed at the end.

2 KINEMATICS DESCRIPTION

The beam material is assumed to be isotropic and elastic, and the cross section of the beam element is assumed to be uniform and doubly symmetric. Bernoulli's assumptions are considered, and cross-sectional distortion, shear, and warping effects are neglected. Always, the deformational and rotational displacements with respect to the corotational frame are assumed to be tiny. To guarantee that these requirements remain valid, and results are reliable, appropriate element sizes and load steps are carefully chosen. So, the small strain theory hypothesis is employed in the used corotational formulation.

After discretization of the structure into finite elements, the i^th beam element can be defined with two end nodes (n=1, 2). Every node has six degrees of freedom and is defined with respect to three frames, as shown in Figure 1. These coordinate systems are the fixed global coordinate system associated with the fixed global frame (X, Y, Z), the fixed local coordinate system associated with the fixed local frame ($\widehat{\mathrm{x}}$_i, $\widehat{\mathrm{y}}$_i,$\widehat{\mathrm{z}}$_i) and the moving local coordinate system associated with the corotational local frame ($\stackrel{-}{\mathrm{x}}$_i, $\stackrel{-}{\mathrm{y}}$_i,$\stackrel{-}{\mathrm{z}}$_i). This local corotational frame is updated and attached to each beam element. It also translates and rotates with the beam element but does not deform with it. Figure 1 also shows the three configurations used in the analysis: an initial configuration, the j^th equilibrium configuration, and a current configuration. The element’s initial length is L_o; after deformation in the current configuration, the element's length is equal to the arc length${\mathrm{}\mathrm{S}}_{\mathrm{i}}$, while L_c is the current chord length.

Figure 1
Kinematics description and coordinate systems of the i^th spatial beam element.

For the current configuration, as shown in Figure 2(a), the nodal displacement vector for the i^th beam element in the fixed global coordinate system is given by:

$D_{\text{i}}={\left[U_1{V}_1{W}_1\theta X_1\theta Y_1\theta Z_1{U}_2{V}_2{W}_2\theta X_2\theta Y_2\theta Z_2\right]}^{\text{T}}$ (1)

where U_n (n=1, 2), V_n (n =1, 2) and W_n (n =1, 2) are the displacement translational components in X, Y and Z directions, respectively, and $\theta X_n.\mathrm{}\theta Y_n\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\theta Z_n$ (n =1, 2) are the counterclockwise rotations about X, Y and Z axes, respectively. The nodal incremental displacement vector of the i^th beam element in the global coordinate system is defined as:

Figure 2
Nodal displacements and forces of the i^th beam element: (a) displacements and forces positive signs in the global coordinate system and the corotaional local coordinate system, and (b) element displacements and internal forces with the attached corotaional local frame in planes ($\stackrel{-}{\mathrm{x}}$_i - $\stackrel{-}{\mathrm{y}}$_i) and ($\stackrel{-}{\mathrm{x}}$_i -$\mathrm{}\stackrel{-}{\mathrm{z}}$_i).

$\Delta D_{\text{i}}={\left[\Delta U_1\Delta V_1\Delta W_1\Delta \theta X_1\Delta \theta Y_1\Delta \theta Z_1\Delta U_2\Delta V_2\Delta W_2\Delta \theta X_2\Delta \theta Y_2\Delta \theta Z_2\right]}^{\text{T}}$ (2)

where $∆$U_n (n=1, 2), $∆$V_n (n=1, 2) and $∆$W_n (n=1, 2) are the incremental translational components of displacement in X, Y and Z directions, respectively, and $∆\theta X_n.\mathrm{}∆\theta Y_n\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}∆\theta Z_n$ (n =1, 2) are the counterclockwise incremental rotations about X, Y and Z axes, respectively. The nodal displacement vector ${\mathbf{D}}_{\mathrm{i}}$ can be updated by:

$D_{\text{i}}=D_{\text{i}}^{\text{j}}+\Delta D_{\text{i}}$ (3)

where ${\mathbf{D}}_{\mathrm{i}}^{\mathrm{j}}\mathrm{}$is the nodal displacement vector for the i^th beam element in the fixed global coordinate system, at the j^th equilibrium configuration. The nodal displacement vector ${\mathbf{D}}_{\mathrm{i}}$is divided into nodal translational displacement vector ${\mathbf{D}\mathbf{t}}_{\mathrm{i}}$ and nodal rotational displacement vector${\mathrm{}\mathbf{D}\mathbf{r}}_{\mathrm{i}}$, which can be written as:

$D{t}_{\text{i}}={\left[U_1{V}_1{W}_1{U}_2{V}_2{W}_2\right]}^{\text{T}}$ (4)

$D{r}_{\text{i}}={\left[\theta X_1\theta Y_1\theta Z_1\theta X_2\theta Y_2\theta Z_2\right]}^{\text{T}}$ (5)

Similarly, ${∆\mathbf{D}}_{\mathrm{i}}$ vector can be separated into the nodal incremental translational displacement vector $∆{\mathbf{D}\mathbf{t}}_{\mathrm{i}}$ and the nodal incremental rotational displacement vector${\mathrm{}∆\mathbf{D}\mathbf{r}}_{\mathrm{i}}$. Therefore, the nodal coordinates' vector of the i^th beam element in the fixed global system can be continually updated as follows:

$X_{\text{i}}=X_{\text{n}}^{\text{j}}+\Delta D{t}_{\text{i}},$ $X_{\text{i}}={\left[\begin{array}{ccc}{X}_i& Y_i& Z_i\end{array}\right]}^T$ (6)

where ${\mathbf{X}}_{\mathrm{i}}$is the vector of the nodal coordinates of the i^th beam element relative to the fixed global frame, at the current configuration, and ${\mathbf{X}}_{\mathrm{i}}^{\mathrm{j}}$ is the vector of the nodal coordinates relative to the fixed global frame, at the j^th equilibrium configuration. The nodal displacement vector of the i^th beam element in the element corotational local coordinate system, at the current configuration, is:

$d_{\text{i}}={\left[u_1{v}_1{w}_1\theta x_1\theta y_1\theta z_1{u}_2{v}_2{w}_2\theta x_2\theta y_2\theta z_2\right]}^{\text{T}}$ (7)

where u_n (n=1, 2), v_n (n=1, 2), and w_n (n=1, 2) are the displacement translational components in $\stackrel{-}{\mathrm{x}}$_i, $\stackrel{-}{\mathrm{y}}$_i, and $\stackrel{-}{\mathrm{z}}$_i directions, respectively, and $\theta x_n.\mathrm{}\theta y_n.\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\theta z_n$ (n =1, 2) are the counterclockwise deformational rotations after eliminating the rigid body rotations. As shown in Figures 1 and 2(b), the displacement component v_n (n=1, 2) and w_n (n=1, 2) are equal to zero because of the characteristics of the attached corotational frame. Also, the axial displacement of the first node is selected to be zero, while the axial displacement of the second node is$u_2$. As a result, the nodal displacement vector${\mathbf{}\mathbf{d}}_{\mathrm{i}}\mathrm{}$has only seven nonzero components that can simplified to:

$d_{\text{i}}={\left[000\theta x_1\theta y_1\theta z_1{u}_{2}00\theta x_2\theta y_2\theta z_2\right]}^{\text{T}}$ (8)

The axial displacement $u_2$ in equation (8) can be written as:

$u_2={\text{S}}_{\text{i}}-{\text{L}}_{\text{o}}={\text{L}}_{\text{c}}-{\text{L}}_{\text{o}}+{\text{b}}_{\text{i}}$ (9)

where ${\mathrm{b}}_{\mathrm{i}}$ is the element axial elongation due to the bowing effect, which can be determined in terms of rotations (Chan, 1992) as

$\text{b}=\frac{{\text{L}}_{\text{c}}}{30}\left[2\theta y_1{}^2-\theta y_1\theta y_2+2\theta y_2{}^2\left]+\frac{{\text{L}}_{\text{c}}}{30}\right[2\theta z_1{}^2-\theta z_1\theta z_2+2\theta z_2{}^2\right]$ (10)

The internal elastic force vector for the i^th beam element in the fixed global coordinate system, at the current configuration, can be written as:

$F_i^e={\left[F{X}_{1}F{Y}_{1}FZ{}_{1}M{X}_{1}M{Y}_{1}M{Z}_{1}F{X}_{2}F{Y}_{2}F{Z}_{2}M{X}_{2}M{Y}_{2}M{Z}_2\right]}^{\text{T}}$ (11)

The internal elastic force vector of the ith beam element in the element corotational local coordinate system, at the current configuration, is:

$f_i^e={\left[f{x}_{1}f{y}_{1}f{z}_{1}m{x}_{1}m{y}_{1}m{z}_{1}f{x}_{2}f{y}_{2}f{z}_{2}m{x}_{2}m{y}_{2}m{z}_2\right]}^{\text{T}}$ (12)

where the internal elastic force vectors' components in both the fixed global coordinate system and the corotational local coordinate system are shown in Figure 2(a) with their positive signs.

3 TRANSFORMATION PROCEDURE

The transformation procedure depends upon updating the coordinates with every equilibrium configuration during the analysis. Two main stages are employed here to perform the transformation from the fixed global frame to the moving corotational local frame. The first stage is the transformation from the fixed global frame to the fixed local frame, and the second stage is the transformation from the fixed local frame to the moving corotational local frame. Assuming that V_d is a 3D vector associated with the fixed global frame, the relation between the fixed global frame (X, Y, Z) and the fixed local frame ($\widehat{\mathrm{x}}$_i, $\widehat{\mathrm{y}}$_i, $\widehat{\mathrm{z}}$_i) can be expressed by:

${\widehat{v}}_{\text{d}}=r_{\text{o}}{V}_{\text{d}}$ (13)

where r_o is an orthogonal matrix (3 × 3) which can be determined from the direction cosines of the fixed local frame relative to the fixed global frame. For a three-dimensional frame element, this matrix turns into a (12 × 12) matrix as follows:

$T_{\text{o}}=\left[\begin{array}{cc}\begin{array}{cc}{r}_{\text{o}}& 0\\ 0& r_{\text{o}}\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}{r}_{\text{o}}& 0\\ 0& r_{\text{o}}\end{array}\end{array}\right]$ (14)

Similarly, the relation between the ($\widehat{\mathrm{x}}$_i, $\widehat{\mathrm{y}}$_i, $\widehat{\mathrm{z}}$_i) frame and the current corotational local frame ($\stackrel{-}{\mathrm{x}}$_i, $\stackrel{-}{\mathrm{y}}$_i,$\mathrm{}\stackrel{-}{\mathrm{z}}$_i) is:

${\overline{v}}_{\text{d}}=r_{\text{c}}{\widehat{v}}_{\text{d}}$ (15)

where r_c is also an orthogonal matrix (3 × 3) which can be obtained from the direction cosines of the corotational local frame relative to the fixed local frame. For a three-dimensional frame element, this matrix turns into a (12 × 12) matrix as follows:

$T_{\text{c}}=\left[\begin{array}{cc}\begin{array}{cc}{r}_{\text{c}}& 0\\ 0& r_{\text{c}}\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}{r}_{\text{c}}& 0\\ 0& r_{\text{c}}\end{array}\end{array}\right]$ (16)

One can write the vector ${\stackrel{-}{\mathbf{v}}}_{\mathrm{d}}$ in terms of V_d as

${\overline{v}}_{\text{d}}=r_{\text{r}}{V}_{\text{d}}$ (17)

where

$r_{\text{r}}=r_{\text{c}}{r}_{\text{o}}$ (18)

Therefore, the transformation matrix for the three-dimensional i^th beam element from the fixed global frame to the moving corotational local frame, at the current configuration, can be expressed by

$T_{\text{r}}=T_{\text{c}}{T}_{\text{o}}$ (19)

Both transformation matrices r_o and r_c are determined using Tait-Bryan angles, which describe the three successive rotations of the three dimensional beam element.

3.1 Transformation from the fixed global frame to the fixed local frame

3.1.1 General case

The first stage is to transform from the fixed global system to the fixed local system using the three successive rotaions ${\mathrm{\beta }}_{\mathrm{o}}$, ${\mathrm{\gamma }}_{\mathrm{o}}$ and ${\mathrm{\alpha }}_{\mathrm{o}}$, as shown in Figures 3(a), 3(b), 3(c), and 4(a)as follows:

Figure 3
Three successive rotations of a three dimensional beam element: (a) Rotation angle ${\beta }_o$of corodinate axes about Y axis, from (X, Y, X) to ($X_{{\beta }_o}$, Y, ${\mathrm{}Z}_{{\beta }_o}$), (b) Rotation angle ${\gamma }_o$ of corodinate axes about $Z_{{\beta }_o}$axis, from ($X_{{\beta }_o}$, Y, $Z_{{\beta }_o}$) to ($X_{{\gamma }_o}$,$Y_{{\gamma }_o}$, $Z_{{\beta }_o}$), and (c)Rotation angle ${\alpha }_o$ of corodinate axes about $X_{{\gamma }_o}$ axis, from ($X_{{\gamma }_o}$,$Y_{{\gamma }_o}$, $Z_{{\beta }_o}$) to ($\widehat{x}$_i, $\widehat{y}$_i,$\mathrm{}\widehat{z}$_i).

$r_{{\text{β}}_{\text{o}}}=\left[\begin{array}{ccc}\cos {\text{β}}_{\text{o}}& 0& \sin {\text{β}}_{\text{o}}\\ 0& 1& 0\\ -\sin {\text{β}}_{\text{o}}& 0& \cos {\text{β}}_{\text{o}}\end{array}\right]$ (20)

where

$\cos {\text{β}}_{\text{o}}=\frac{{\text{C}}_{\text{X}}}{{\text{C}}_{\text{XZ}}},$ $\sin {\text{β}}_{\text{o}}=\frac{{\text{C}}_{\text{Z}}}{{\text{C}}_{\text{XZ}}},$ ${\text{C}}_{\text{X}}=\frac{{\text{X}}_2-{\text{X}}_1}{{\text{L}}_{\text{o}}},$ ${\text{C}}_{\text{Y}}=\frac{{\text{Y}}_2-{\text{Y}}_1}{{\text{L}}_{\text{o}}},$ ${\text{C}}_{\text{Z}}=\frac{{\text{Z}}_2-{\text{Z}}_1}{{\text{L}}_{\text{o}}}$

${\text{C}}_{\text{XZ}}={\left({\text{C}}_{\text{X}}{}^2+{\text{C}}_{\text{Z}}{}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},$ ${\text{L}}_{\text{o}}={\left({\left({\text{X}}_2-{\text{X}}_1\right)}^2+{\left({\text{Y}}_2-{\text{Y}}_1\right)}^2+{\left({\text{Z}}_2-{\text{Z}}_1\right)}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ (21)

$r_{{\text{γ}}_{\text{o}}}=\left[\begin{array}{ccc}\cos {\text{γ}}_{\text{o}}& \sin {\text{γ}}_{\text{o}}& 0\\ -\sin {\text{γ}}_{\text{o}}& \cos {\text{γ}}_{\text{o}}& 0\\ 0& 0& 1\end{array}\right]$ (22)

where

$\cos {\text{γ}}_{\text{o}}={\text{C}}_{\text{XZ}},$ $\sin {\text{γ}}_{\text{o}}={\text{C}}_{\text{Y}}$ (23)

$r_{{\text{α}}_{\text{o}}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\text{α}}_{\text{o}}& \sin {\text{α}}_{\text{o}}\\ 0& -\sin {\text{α}}_{\text{o}}& \cos {\text{α}}_{\text{o}}\end{array}\right]$ (24)

where

$\cos {\text{α}}_{\text{o}}=\frac{{\text{Y}}_{{\text{Pγ}}_{\text{o}}}}{{\left({\left({\text{Y}}_{{\text{Pγ}}_{\text{o}}}\right)}^2+{\left({\text{Z}}_{{\text{Pγ}}_{\text{o}}}\right)}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}},$$\sin {\text{α}}_{\text{o}}=\frac{{\text{Z}}_{{\text{Pγ}}_{\text{o}}}}{{\left({\left({\text{Y}}_{{\text{Pγ}}_{\text{o}}}\right)}^2+{\left({\text{Z}}_{{\text{Pγ}}_{\text{o}}}\right)}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$.(25)

where ${\mathrm{X}}_{\mathrm{P}{\mathrm{\gamma }}_{\mathrm{o}}}$, ${\mathrm{Y}}_{\mathrm{P}{\mathrm{\gamma }}_{\mathrm{o}}}$, and ${\mathrm{Z}}_{\mathrm{P}{\mathrm{\gamma }}_{\mathrm{o}}}$are the coordinates of an assumed point P relative to ${\mathrm{\gamma }}_{\mathrm{o}}$ frame as shown in Figure 4(b), which can be obtained from the following equation:

Figure 4
Transformation from the fixed global frame to the fixed local frame: (a) three successive rotations (${\mathrm{\beta }}_{\mathrm{o}}$, ${\mathrm{\gamma }}_{\mathrm{o}}$, ${\mathrm{\alpha }}_{\mathrm{o}}$) of a three dimensional beam element, and (b) corodinates of a reference point P relative to frame (${\mathrm{X}}_{{\mathrm{\gamma }}_{\mathrm{o}}}$,${\mathrm{Y}}_{{\mathrm{\gamma }}_{\mathrm{o}}}$, ${\mathrm{Z}}_{{\mathrm{\gamma }}_{\mathrm{o}}}$).

${\left[\begin{array}{ccc}{\text{X}}_{{\text{Pγ}}_{\text{o}}}& {\text{Y}}_{{\text{Pγ}}_{\text{o}}}& {\text{Z}}_{{\text{Pγ}}_{\text{o}}}\end{array}\right]}^{\text{T}}=r_{{\text{γ}}_{\text{o}}}{r}_{{\text{β}}_{\text{o}}}{\left[\begin{array}{ccc}{\text{X}}_{\text{P}1}& {\text{Y}}_{\text{P}1}& {\text{Z}}_{\text{P}1}\end{array}\right]}^{\text{T}}$ (26)

Hence, the rotation matrix ${\mathbf{r}}_{\mathrm{o}}$ can be obtained as:

$r_{\text{o}}=r_{{\text{α}}_{\text{o}}}{r}_{{\text{γ}}_{\text{o}}}{r}_{{\text{β}}_{\text{o}}}$ (27)

By substitution, this matrix takes the form:

$$r_{\text{o}}=\left[\begin{array}{ccc}{C}_X& C_Y& C_Z\\ \frac{-C_X{C}_Y\cos {\alpha }_o-C_Z\sin {\alpha }_o}{C_X{}_Z}& C_X{}_Z\cos {\alpha }_o& \frac{-C_Y{C}_Z\cos {\alpha }_o+C_X\sin {\alpha }_o}{C_X{}_Z}\\ \frac{C_X{C}_Y\sin {\alpha }_o-C_Z\cos {\alpha }_o}{C_X{}_Z}& -C_X{}_Z\sin {\alpha }_o& \frac{C_Y{C}_Z\sin {\alpha }_o+C_X\cos {\alpha }_o}{C_X{}_Z}\end{array}\right]$$ (28)

3.1.2 Special cases

It should be noted that the rotation angle ${\mathrm{\alpha }}_{\mathrm{o}}$ is insignificant, in the case of a circular cross-section element. Thus, the rotation matrix ${\mathbf{r}}_{\mathrm{o}}$ can be calculated as:

$r_{\text{o}}=r_{{\text{γ}}_{\text{o}}}{r}_{{\text{β}}_{\text{o}}}$ (29)

There is another special case where the initial position of the element is vertical, in the Y-axis direction. In order to get${\mathbf{}\mathbf{r}}_{\mathrm{o}}$, there are only two successive rotations not three as in the general case. The first rotation ${\mathrm{\gamma }}_{\mathrm{o}}$is either 90^o or 270^o, as shown in Figure 5, depending on whether C_Y is +1 or -1. The other rotation${\mathrm{}\mathrm{\alpha }}_{\mathrm{o}}$, which is shown in Figure 6, can be determined using a reference point P that lies in the ($\widehat{\mathrm{x}}$_i, $\widehat{\mathrm{y}}$_i) plane. In this case, Eq. (25) is modified as follows:

Figure 5
The rotation angle ${\mathrm{\gamma }}_{\mathrm{o}}$ for vertical member case.

Figure 6
The rotation angle ${\mathrm{\alpha }}_{\mathrm{o}}$ with the reference point for vertical member case.

$\cos {\text{α}}_{\text{o}}=-\frac{{\text{X}}_{\text{P}1}}{\sqrt{{\left({\text{X}}_{\text{P}1}\right)}^2+{\left({\text{Z}}_{\text{P}1}\right)}^2}}\left({\text{C}}_{\text{Y}}\right),\sin {\text{α}}_{\text{o}}=\frac{{\text{Z}}_{\text{P}1}}{\sqrt{{\left({\text{X}}_{\text{P}1}\right)}^2+{\left({\text{Z}}_{\text{P}1}\right)}^2}}$ (30)

Thus, the matrix ${\mathbf{r}}_{\mathrm{o}}$ in Eq. (29), can be written as:

$r_{\text{o}}=\left[\begin{array}{ccc}0& {\text{C}}_{\text{Y}}& 0\\ -{\text{C}}_{\text{Y}}\cos {\text{α}}_{\text{o}}& 0& \sin {\text{α}}_{\text{o}}\\ {\text{C}}_{\text{Y}}\sin {\text{α}}_{\text{o}}& 0& \cos {\text{α}}_{\text{o}}\end{array}\right]$ (31)

Substituting Eq. (28) or Eq. (31) into Eq. (14), the matrix ${\mathbf{T}}_{\mathrm{o}}$can be determined.

3.2 Transformation from the fixed local frame to the moving corotational local frame

3.2.1 General case

The second stage is to transform from the fixed local system ($\widehat{\mathrm{x}}$_i, $\widehat{\mathrm{y}}$_I, $\widehat{\mathrm{z}}$_i) to the moving corotational local system ($\stackrel{-}{\mathrm{x}}$_i,$\stackrel{-}{\mathrm{}\mathrm{y}}$_i,$\stackrel{-}{\mathrm{z}}$_i), at the current configuration, using the three successive rotaions ${\mathrm{\beta }}_{\mathrm{c}}$, ${\mathrm{\gamma }}_{\mathrm{c}}$ and ${\mathrm{\alpha }}_{\mathrm{c}}$, as shown in Figure 7. The transformation procedure is similar to the rotations ${\mathrm{\beta }}_{\mathrm{o}}$, ${\mathrm{\gamma }}_{\mathrm{o}}$ and ${\mathrm{\alpha }}_{\mathrm{o}}$ in Figures 3(a), 3(b), 3(c), and 4(a), but the trigonomatric rules for rotation angles are based on the relative displcacements between the two end points of each element, as follows:

$r_{{\text{β}}_{\text{c}}}=\left[\begin{array}{ccc}\cos {\text{β}}_{\text{c}}& 0& \sin {\text{β}}_{\text{c}}\\ 0& 1& 0\\ -\sin {\text{β}}_{\text{c}}& 0& \cos {\text{β}}_{\text{c}}\end{array}\right]$ (32)

where

Figure 7
Transformation from fixed local frame to moving corotational local frame using three successive rotations (${\mathrm{\beta }}_{\mathrm{c}}$, ${\mathrm{\gamma }}_{\mathrm{c}}$, ${\mathrm{\alpha }}_{\mathrm{c}}$) of a three dimensional beam element.

$\cos {\text{β}}_{\text{c}}=\frac{{\text{L}}_{\text{o}}+{\widehat{U}}_i^r}{\overline{P_1{P}_2^{\text{'}}}},$ $\sin {\text{β}}_{\text{c}}=\frac{{\widehat{W}}_i^r}{\overline{P_1{P}_2^{\text{'}}}},$ $\overline{P_1{P}_2^{\text{'}}}={\left({\left({\text{L}}_{\text{o}}+{\widehat{U}}_i^r\right)}^2+{\left({\widehat{W}}_i^r\right)}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ (33)

It is worth noting that Bathe and Bolourchi (1979) provided only an expression for the cosine of the angle. Hence, when the angle${\mathrm{}\mathrm{\beta }}_{\mathrm{c}}>90°$a problem appears. Therefore, in this work an expression for the sine is provided. The location of the element can be precisely specified by both trigonometric relations. As shown in Figure 7,${\widehat{U}}_i^r$,${\widehat{V}}_i^r$ and${\widehat{W}}_i^r$ are the i^th beam element relative translational displacements with respect to the fixed local frame. These relative displacements can be obtained from the corresponding relative displacement with respect to the fixed global system that can be determined from the vector ${\mathbf{D}\mathbf{t}}_{\mathrm{i}}$in Eq. (4), and the matrix r_o in Eq. (28), or Eq. (31) for the vertical member, as follows:

$\left[\begin{array}{ccc}{\widehat{U}}_i^r& {\widehat{V}}_i^r& {\widehat{W}}_i^r\end{array}\right]=\left[\begin{array}{cc}-r_{\text{o}}& r_{\text{o}}\end{array}\right]D{t}_{\text{i}}$ (34)

Consequently, the rotation matrix ${\mathbf{r}}_{{\mathrm{\gamma }}_{\mathrm{c}\mathrm{}}}$can be determined as:

$r_{{\text{γ}}_{\text{c}}}=\left[\begin{array}{ccc}\cos {\text{γ}}_{\text{c}}& \sin {\text{γ}}_{\text{c}}& 0\\ -\sin {\text{γ}}_{\text{c}}& \cos {\text{γ}}_{\text{c}}& 0\\ 0& 0& 1\end{array}\right]$ (35)

where

$\cos {\text{γ}}_{\text{c}}=\left(SN\right)\frac{\overline{P_1{P}_2^{\text{'}}}}{{\text{L}}_{\text{c}}},$ $\sin {\text{γ}}_{\text{c}}=\frac{{\widehat{V}}_i^r}{{\text{L}}_{\text{c}}},$ ${\text{L}}_{\text{c}}={\left({\left(\overline{P_1{P}_2^{\text{'}}}\right)}^2+{\left({\widehat{V}}_i^r\right)}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ (36)

and SN is equal +1 if $\left({\mathrm{L}}_{\mathrm{o}}+{\widehat{U}}_i^r\right)\ge 0$ and -1 if $\left({\mathrm{L}}_{\mathrm{o}}+{\widehat{U}}_i^r\right)<0$. Both trigonometric rules are defined properly with controlled signs, enabling the precise determination of the element's position. Hence, the model can deal with angle${\mathrm{}\mathrm{\gamma }}_{\mathrm{c}}>90°$ contrary to the technique of Bathe and Bolourchi (1979).

The rotation matrix for relative translational displacements can be written as

$r_{{\text{d}}_{\text{c}}}=r_{{\text{γ}}_{\text{c}}}{r}_{{\text{β}}_{\text{c}}}$ (37)

$r_{{\text{d}}_{\text{c}}}=\left[-\begin{array}{ccc}\cos {\text{β}}_{\text{c}}\cos {\text{γ}}_{\text{c}}& \sin {\text{γ}}_{\text{c}}& \sin {\text{β}}_{\text{c}}\cos {\text{γ}}_{\text{c}}\\ \cos {\text{β}}_{\text{c}}\sin {\text{γ}}_{\text{c}}& \cos {\text{γ}}_{\text{c}}& -\sin {\text{β}}_{\text{c}}\sin {\text{γ}}_{\text{c}}\\ -\sin {\text{β}}_{\text{c}}& 0& \cos {\text{β}}_{\text{c}}\end{array}\right]$ (38)

Then, the rotation matrix ${\mathbf{r}}_{{\mathrm{\alpha }}_{\mathrm{c}}}$ can be obtained as:

$r_{{\text{α}}_{\text{c}}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\text{α}}_{\text{c}}& \sin {\text{α}}_{\text{c}}\\ 0& -\sin {\text{α}}_{\text{c}}& \cos {\text{α}}_{\text{c}}\end{array}\right]$ (39)

The angle${\mathrm{\alpha }}_{\mathrm{c}}$ is computed using the incremental procedure described in Bathe and Bolourchi (1979). The rotation matrix r_c in Eq. (16) can be expressed as:

$r_{\text{c}}=r_{{\text{α}}_{\text{c}}}{r}_{{\text{d}}_{\text{c}}}$ (40)

3.2.2 Special cases

In case of a vertical member that is parallel to the Y-axis, the rotation ${\mathrm{\beta }}_{\mathrm{c}}$ vanishes and the rotation${\mathrm{\gamma }}_{\mathrm{c}}$is either 90^o or 270^o, depending on the member position. Traditionally, this case could create singularity and it had been the source of many difficulties (Nunes et al., 2003), and the authors of this research work suggested to avoid the rotation ${\mathrm{\gamma }}_{\mathrm{c}}$ to be either 90^o or 270^o. In this work, this problem is solved by letting the code search for the alignment of the element, which means specifying if the rotation${\mathrm{\gamma }}_{\mathrm{c}}$ is either 90^o or 270^o. Thus, the matrix ${\mathbf{r}}_{{\mathrm{d}}_{\mathrm{c}\mathrm{}}}$can be rewritten as:

$r_{{\text{d}}_{\text{c}}}=\left[\begin{array}{ccc}0& {\text{C}}_{\text{Y}}{}^{\text{'}}& 0\\ -{\text{C}}_{\text{Y}}{}^{\text{'}}& 0& 0\\ 0& 0& 1\end{array}\right]$ (41)

where ${{\mathrm{C}}_{\mathrm{Y}\mathrm{}}}^{\text{'}}$ can be specified using the current vector of the nodal coordinates in Eq. (8) as follows:

${\text{C}}_{\text{Y}}{}^{\text{'}}=\frac{{\text{Y}}_2-{\text{Y}}_1}{{\text{L}}_{\text{c}}}$ (42)

Thus, the value of ${{\mathrm{C}}_{\mathrm{Y}\mathrm{}}}^{\text{'}}$ is either +1 and the rotation ${\mathrm{\gamma }}_{\mathrm{c}}$ is 90^o or -1 and the rotation ${\mathrm{\gamma }}_{\mathrm{c}}$ is 270^o.

Substituting Eqs. (38) or (41) and (40) into Eq. (36), the matrix r_c can be determined. Hence, the transformation matrix ${\mathbf{T}}_{\mathrm{c}}$ in Eq. (16) has been specified. Then, the transformation matrix ${\mathbf{T}}_{\mathrm{r}}$ in Eq. (19) has been determined.

4 STRAIN ENERGY AND STIFFNESS MATRIX

Considering isotropic elastic materials, the constitutive relation between the stress vector ${\mathbf{\sigma }}_{\mathrm{i}}$ and the strain vector ${\mathbf{ϵ}}_{\mathrm{i}}$ of the i^th beam element can be expressed by:

${\sigma }_{\text{i}}=E_{\text{i}}{ϵ}_{\text{i}}$ (43)

where ${\mathbf{E}}_{\mathrm{i}}$ is the symmetric matrix of the elastic coefficients. The strain vector is given by:

${ϵ}_{\text{i}}=D_f{v}_q$ (44)

where ${\mathbf{D}}_f$is the differential operator matrix and ${\mathit{v}}_{\mathit{q}}$is the deformation vector, which can be defined as follows:

$v_q=N_{\text{i}}{d}_{\text{i}}$ (45)

where ${\mathbf{N}}_{\mathrm{i}}$ is the shape functions matrix, which is given in ^{Appendix 1}. Substituting Eq. (45) into Eq. (44), the strain vector can be expressed as follows:

${ϵ}_{\text{i}}=D_f{N}_{\text{i}}{d}_{\text{i}}$ (46)

Combining Eqs. (43) and (46), the stress vector can be rewritten as follows:

${\sigma }_i=E_i{D}_f{N}_{\text{i}}{d}_{\text{i}}$ (47)

The strain energy for the i^th beam element ${}_{\mathrm{i}}$ is given by:

${\Pi }_{\text{i}}=\frac{1}{2}{{\displaystyle \int }}_{{\text{V}}_{\text{i}}}^{}{\sigma }_{\text{i}}^{\text{T}}{ϵ}_{\text{i}}d{\text{V}}_{\text{i}}$ (48)

where ${\mathrm{V}}_{\mathrm{i}}$ is the volume. Substituting Eqs. (46) and (47) into Eq. (48), the strain energy can be expressed as

${\Pi }_{\text{i}}=\frac{1}{2}{{\displaystyle \int }}_{{\text{V}}_{\text{i}}}^{}{\left(E_{\text{i}}{D}_f{N}_{\text{i}}{d}_{\text{i}}\right)}^{\text{T}}{D}_f{N}_{\text{i}}{d}_{i}d{\text{V}}_{\text{i}}$ (49)

The element local displacement vector ${\mathbf{d}}_{\mathrm{i}}$is independent of the volume${\mathrm{V}}_{\mathrm{i}}$. Subsequently, Eq. (49) can be simplified to the following form:

${\Pi }_{\text{i}}=\frac{1}{2}{d}_{\text{i}}^{\text{T}}{k}_{\text{i}}{d}_{\text{i}}$ (50)

where k_i is the symmetric element stiffness matrix in the corotational local coordinate system, which can be defined as:

$k_{\text{i}}={{\displaystyle \int }}_{{\text{V}}_{\text{i}}}^{}{B}^T{E}_{\text{i}}Bd{\text{V}}_{\text{i}}$ (51)

and B is defined as:

$B=D_f{N}_{\text{i}}$ (52)

Equation (50) can be written as:

${\Pi }_{\text{i}}=\frac{1}{2}{d}_{\text{i}}^{\text{T}}{k}_{\text{i}}{d}_{\text{i}}=\frac{1}{2}{\left(T_{\text{r}}{D}_{\text{i}}\right)}^{\text{T}}{k}_{\text{i}}\left(T_{\text{r}}{D}_{\text{i}}\right)=D_{\text{i}}^{\text{T}}{K}_{\text{i}}{D}_{\text{i}}$ (53)

$K_{\text{i}}=T_{\text{r}}{}^{\text{T}}{k}_{\text{i}}{T}_{\text{r}}$ (54)

where K_i is the symmetric element stiffness matrix in the fixed global coordinate system.

For the beam element used in this research work, the element stiffness matrix in the local coordinate system k_i can be expressed as

$k_{\text{i}}=k_1+k_2$ (55)

where k₁ is the axial and bending stiffness matrix and k₂ is the geometric stiffness matrix for the i^th beam element. Stiffness matrices are attached in ^{Appendix 2} and ^{Appendix 3}.

It is important to note that, in general, the rotation is non-commutative. However, the kinematic description employed in this paper depends on the continuous updating of the corotational frame to separate the large rigid-body motion from the deformational motion. We limit our analysis to small strain theory as explained in section 2, therefore the local rotational displacements relative to the element chord in each load step are always small. Hence, the vectorial addition can be applied to these relatively small rotational displacements, which offers simplicity and convenience. Also, the resulted stiffness matrices are symmetric.

5 THE EQUILIBRIUM EQUATION

The virtual work principle (Goldstein et al., 1980; Shabana, 2020) can be classified as the basis for variational principles of mechanics. This principle takes the form:

$\delta {\text{W}}_{\text{T}}={\displaystyle \sum }_{i=1}^n\left(\delta Q_i^{int}-\delta Q_i^{ext}\right)=0$ ^. (56)

where ${\mathrm{W}}_{\mathrm{T}}$ is the total virtual work, $Q_i^{int}$is internal virtual work for the i^th beam element and $Q_i^{ext}$ is the external virtual work for the i^th beam element. The internal virtual work $Q_i^{int}$ can be written as

$\delta Q_i^{int}={{\displaystyle \int }}_{{\text{V}}_{\text{i}}}^{}\delta {ϵ}_{\text{i}}^{\text{T}}{\sigma }_{\text{i}}d{\text{V}}_{\text{i}}$ (57)

Substituting Eqs. (46) and (47) in Eq. (57) one obtains:

$\delta Q_i^{int}={{\displaystyle \int }}_{{\text{V}}_{\text{i}}}^{}\delta {\left(D_f{N}_{\text{i}}{d}_i\right)}^{\text{T}}{E}_{\text{i}}{D}_f{N}_{\text{i}}{d}_{\text{i}}d{\text{V}}_{\text{i}}$ ^. (58)

Using Eqs. (51), (52) and (53) one can write Eq. (58) as:

$\delta Q_i^{int}=\delta D_{\text{i}}^{\text{T}}{K}_{\text{i}}{D}_{\text{i}}$ (59)

The external virtual work for the i^th beam element in static analysis can be written as

$\delta Q_i^{ext}=\delta D_{\text{i}}^{\text{T}}{F}_{\text{i}}^{\text{P}}$ (60)

Substituting Eqs. (59) and (60) in Eq. (56) one obtains:

$\delta {\text{W}}_{\text{T}}={\displaystyle \sum }_{i=1}^n\left(\delta D_{\text{i}}^{\text{T}}{K}_{\text{i}}{D}_{\text{i}}-\delta D_{\text{i}}^{\text{T}}{F}_{\text{i}}^{\text{P}}\right)=0$ (61)

Thus,

$K_{\text{i}}{D}_{\text{i}}-F_{\text{i}}^{\text{P}}=0$ (62)

One can write

$F_{\text{i}}^{\text{e}}=K_{\text{i}}{D}_{\text{i}}=T_{\text{r}}{}^{\text{T}}{f}_{\text{i}}^{\text{e}}$ (63)

where ${\mathbf{F}}_{\mathrm{i}}^{\mathrm{e}}\mathrm{}$is the internal elastic force vector in the fixed global and ${\mathbf{f}}_{\mathrm{i}}^{\mathrm{e}}$ is the internal elastic force vector in the local coordinate systems, which is given by

$f_{\text{i}}^{\text{e}}=k_{\text{i}}{d}_{\text{i}}$ (64)

Accordingly, the nonlinear equilibrium equation of the entire structure is

$F^{\text{e}}-F^{\text{P}}=0$ (65)

where ${\mathbf{F}}^{\mathrm{P}}$ is the vector of the external applied forces of the entire structure and ${\mathbf{F}}^{\mathrm{e}}$ is the internal elastic force vector of the entire structure. Both ${\mathbf{F}}^{\mathrm{P}}$ and ${\mathbf{F}}^{\mathrm{e}}$ are defined in the fixed global coordinate system.

6 NUMERICAL ALGORITHM

The equilibrium equation is solved using an incremental iterative procedure. This procedure is based on the full Newton-Raphson method (Bathe, 1996). Eq. (65) can be rewritten as

$\psi =F^{\text{e}}-{\lambda }_c{F}^{\text{P}}=0$ (66)

where $\mathbf{\psi }$ is the out of balance force and ${}_c$ is the loading factor, which can be determined as follows:

${\lambda }_c=\frac{\tau }{\omega }$ (67)

where $\tau$ is the current number of load increment and $\omega$is the total number of load increments. The load increment ${∆\mathbf{F}}^{\mathrm{P}}$ can be written as:

$\Delta F^{\text{P}}=\frac{F^{\text{P}}}{\omega }$ (68)

The iteration equilibrium convergence criterion is given by:

$\parallel \psi \parallel \le e_r\parallel {\psi }_f\parallel$ (69)

where ${\mathbf{\psi }}_{\mathbf{f}}$is the reference out of balance force, which is assumed to be the out of balance force in the first iteration, and$\mathrm{}{e}_r$is the error tolerance.

It is worth mentioning that it is well known that element matrices are used only in the iterative process for the incremental solution, and they do not have to be exact. They are required to allow the solution to converge and satisfy the specified tolerance during solution iterations. That is why many authors have used tangent, secant, or even initial stiffness matrices in their nonlinear FE formulations. Using exact matrices typically costs more time because of the storage of non-symmetric matrices compared with the storage of only the triangular part in the symmetric matrices. Therefore, we are confident that the restriction we made in section (2), which produced symmetric stiffness matrices, allows the whole solution to go faster without affecting the overall accuracy of the results. Through the solution of various numerical examples in the next section, we have assessed the effectiveness and the accuracy of this method.

At the beginning of the iteration procedure, null displacement vectors are assumed, and then the following steps are performed at the beginning of each load increment:

• I. Computing${\mathbf{k}}_{\mathrm{i}}$ and ${\mathbf{f}}_{\mathrm{i}}^{\mathrm{e}}$ for each element using Eqs. (51), (55) and (64), respectively.

• II. In the first iteration, calculating the matrices${\mathbf{}\mathbf{r}}_{\mathrm{o}}$ and${\mathbf{}\mathbf{T}}_{\mathrm{o}}$ using Eqs. (28) and (14). In all other iterations, calculating${\mathbf{}\mathbf{r}}_{\mathrm{c}}$ and ${\mathbf{T}}_{\mathrm{c}}$ using Eqs. (40) and (16).

• III. Determining$$the transformation matrix ${\mathbf{T}}_{\mathrm{r}}$ for each member using Eq. (19).

• IV. Obtaining${\mathbf{K}}_{\mathrm{i}}$ and ${\mathbf{F}}_{\mathrm{i}}^{\mathrm{e}}$ according to Eqs. (54) and (63).

• V. Getting$\mathbf{K}$ and ${\mathbf{F}}^{\mathrm{e}}$ for the entire structure, by assembling the stiffness matrices and the internal elastic force vectors for all elements of the analyzed structure.

• VI. Calculating the out of balance force $\mathbf{\psi }$ from Eq. (66).

• VII. If the convergence condition in Eq. (69) is satisfied, stop the iteration and go to step IX. Otherwise, start the following iteration:

  a. Using the Newton-Raphson method, a displacement corrector vector $\mathbf{R}$ is calculated as:

  $R=-K^{-1}\psi$(70)

  b. Updating the incremental displacement vector as follows:

  $\Delta D=\Delta D_0+R$(71)

where ${∆\mathbf{D}}_0$ is the incremental displacement vector in the previous iteration, which is considered to be zero in the first iteration.

c. Extracting vector $∆{\mathbf{D}}_{\mathrm{i}}$for each element from vector$∆\mathbf{D}$. Consequently, one can update the vectors ${\mathbf{D}}_{\mathrm{i}}$ and ${\mathbf{X}}_{\mathrm{i}}$using Eqs. (3) and (6).

d. Using ${\mathbf{D}}_{\mathrm{i}}$ and${\mathbf{}\mathbf{r}}_{\mathrm{o}}$, the relative displacements are calculated, as in Eq. (34).

e. From the coordinate vector${\mathbf{X}}_{\mathrm{i}}$, the model can check for each element position to apply either the regular rotations matrices in Eqs. (28), (38) and (39), or the vertical member rotations matrices in Eqs. (31), (41) and (39). Then, ${\mathbf{T}}_{\mathrm{c}}$ and ${\mathbf{T}}_{\mathrm{r}}$ are updated for each element.

f. Using the relative displacements, the rigid body rotations can be obtained as:

${\mu }_Y={\text{tan}}^{-1}\left(\frac{{\widehat{W}}_i^r}{\sqrt{{\left({\text{L}}_{\text{o}}+{\widehat{U}}_i^r\right)}^2}}\right)$ (72)

${\mu }_Z={\text{tan}}^{-1}\left(\frac{{\widehat{V}}_i^r}{\sqrt{{\left({\text{L}}_{\text{o}}+{\widehat{U}}_i^r\right)}^2+{\left({\widehat{W}}_i^r\right)}^2}}\right)$ ^. (73)

g. Eliminating the rigid body rotations from the rotational components ${\theta Y}_n$and ${\theta Z}_n$ in the vector ${\mathbf{D}}_{\mathrm{i}\mathrm{}}$, for each element, as follows:

$$\stackrel{`}{\theta Y_n}=\theta Y_n+{\mu }_Y$$ (74)

$\stackrel{`}{\theta Z_n}=\theta Z_n-{\mu }_Z$ (75)

h. Transforming the pure rotations ${\theta X}_n.$ $\stackrel{`}{{\theta Y}_n}$ and $\stackrel{`}{{\theta Z}_n}$, which are relative to the fixed global frame, to determine the corresponding rotational components ${\theta x}_{n.}$ ${\theta y}_n$and ${\theta z}_n$ for vector${\mathbf{d}}_{\mathrm{i}\mathrm{}}$in the current corotational local frame, which are always relative to the element chord, using the procedure detailed in Eqs. (13) and (15).

i. The axial displacement $u_2$ in Eq. (8) is obtained from Eq. (9).

j. Using Eqs. (8) and (64), ${\mathbf{d}}_{\mathrm{i}}$ and ${\mathbf{f}}_{\mathrm{i}}^{\mathrm{e}}$ can be determined, respectively.

• VIII. Going to the start of step IV again.

• IX. The displacement vector at current configuration ${\mathbf{D}}_{\mathrm{N}+1}$can be updated using the displacement vector of the previous configuration ${\mathbf{D}}_{\mathrm{N}}$as:

  $D_{\text{N}+1}=D_{\text{N}}+\Delta D$(76)

• X. Starting the next load step.

The flowchart for the numerical solution is shown in Figure 8. It should be noted that this code includes a detection function in the MATLAB code, which accurately determines the position of each element using the nodal coordinate vector. Thus, it can easily select the appropriate trigonometric rules and the sign of rotations. The detection function specifies the position of the element using the nodal coordinate vector ${\mathbf{X}}_{\mathrm{i}}$ given in Eq. (6). When the relative difference between the coordinate in X and Z direction of the element end nodes approaches zero together, the code detects that the angle ${\mathrm{\beta }}_{\mathrm{c}}$ turns to be zero and the element is in the direction of Y axis. In this case, the code searches for the alignment of the element based on the sign of${{\mathrm{C}}_{\mathrm{Y}\mathrm{}}}^{\text{'}}$in Eq. (42). Additionally, when the relative displacement ${\widehat{V}}_i^r$ turns out to be zero, the angle ${\mathrm{\gamma }}_{\mathrm{c}}$ vanishes. Likewise, if the relative displacement ${\widehat{W}}_i^r$ turns out to be zero, the angle ${\mathrm{\beta }}_{\mathrm{c}}$ vanishes, however, the element in general is not vertical. Hence, the program deals with these special cases separately.

Figure 8
The flowchart of the numerical solution.

This function also controls the angle ${\mathrm{\gamma }}_{\mathrm{c}}$ when the rotation is outside the interval [$\mathrm{\pi }/2$, -$\mathrm{\pi }/2]$. The sign SN in Eq. (36) specifies the cosine of the angle ${\mathrm{\gamma }}_{\mathrm{c}}$ to meet the corresponding element position during the motion because the terms $\stackrel{-}{P_1{P}_2^{\mathrm{\text{'}}}}$ and ${\mathrm{L}}_{\mathrm{c}}$ are always positive which cannot reflect the real sign of cosine of ${\mathrm{\gamma }}_{\mathrm{c}}$. At the same time, the sine of ${\mathrm{\gamma }}_{\mathrm{c}}$ is already specified with the sign of ${\mathrm{}\widehat{V}}_{i\mathrm{}}^r$. This function conserves the code to converge efficiently.

The geometric stiffness matrix is essential in improving convergence. Updating of element matrices in every iteration is also crucial for convergence acceleration. Hence, the proposed formulation experiences a rapid convergence rate. Providing that the load step and the element size are adjusted to satisfy the requirements of the small strain theory.

7 NUMERICAL EXAMPLES

7.1 Clamped- clamped beam subjected to a concentrated vertical load at the mid-span

A clamped-clamped beam is solved, in this part. The beam geometric data are shown in Figure 9. The beam's material properties are given in Table 1. This beam is subjected to a concentrated vertical load at mid-span. This problem is analysed by Mondkar and Powell (1977) using different load stepping procedures. They used 16 elements to solve this problem. The present results are obtained using only ten elements and one loading increment. The used error tolerance is${\mathrm{}10}^{-2}$. The present results are compared with the results of Mondkar and Powell (1977), which used seven loading increments, and the linear analysis results. The present results are significantly in agreement with the results in Mondkar and Powell (1977), as can be shown in Figure 10. Figure 10 also shows that the linear displacement is considerably larger than the corresponding nonlinear displacement.

Figure 9
Geometrical data of the clamped-clamped beam.

Table 1
The beam's material properties (Mondkar and Powell, 1977).

Modulus of elasticity E(ksi)	Density ρ(lb/in3)	Poisson's ratio 𝑣
30,000	0.098	0

Figure 10
Load-displacement curves of the clamped-clamped beam.

7.2 Cantilever beam subjected to a vertical end load at the free end

In this example, a cantilever beam subjected to a concentrated end load is analysed, as shown in Figure 11. The geometric data and material properties of the beam are shown in Table 2. The applied concentrated vertical force P_Y = $-600$ KN. Liu (2014) and Crivelli (1991) solved this problem with different number of elements. The present results are obtained using one loading increment and the error tolerance${\mathrm{e}}_{\mathrm{r}}$ $={10}^{-2}$. The free end vertical displacement and rotation are compared with the results of Liu (2014) and Crivelli (1991), as shown in Table 3. It can be noticed from Table 3 that our results converge to the more accurate values with only 2 elements compared to 8 elements for the results of both Liu (2014) and Crivelli (1991).

Figure 11
Load-displacement curves of the clamped-clamped beam.

Table 2
The geometrical data and material properties of the cantilever beam.

Length L(m)	Cross-sectional area A(m2)	IY (m4)	IZ (m4)	Modulus of elasticity E (GPa)	Modulus of rigidity G(GPa)
5.0	4.8 x 10-3	4.45 x 10-5	4.45 x 10-5	207	80.775

Table 3
Comparison of results for the cantilever beam subjected to vertical end load.

Number of elements	Displacement (m)			Rotation (rad)
	Liu (2014)	Crivelli (1991)	Present	Liu (2014)	Crivelli (1991)	Present
1	2.261	1.833	2.232	0.701	0.747	0.698
2	2.167	2.078	2.155	0.674	0.688	0.673
4	2.157	2.139	2.154	0.673	0.676	0.673
8	2.154	2.154	2.153	0.673	0.673	0.672
16	2.153	2.157	2.152	0.672	0.672	0.672

7.3 Cantilever beam subjected to a concentrated end moment at the free end

The cantilever beam shown in Figure 12 is solved in this example. Table 4 provides the material properties and geometric data of the beam. This beam is subjected to a concentrated end moment M_Z$=10\mathrm{}000\mathrm{}\mathrm{N}.\mathrm{m}\mathrm{}$at the free end. This problem is classified as a large displacement and large rotation problem (Liu, 2014). Wang and Rattanawangcharoen (2008) used 20 elements and 10 loading increments to solve this problem. Liu (2014) solved this problem with one loading increment and ten beam elements. The present results are obtained using one loading increment, ten beam elements and the error tolerance$e_r$ $={10}^{-2}$. The present results are compared with the results in (Liu, 2014; Wang and Rattanawangcharoen, 2008) in Table 5. The comparison shows that the present results agree sufficiently with the referred results.

Figure 12
Load-displacement curves of the cantilever beam.

Table 4
The geometrical data and material properties of the cantilever beam.

Length L(m)	Cross-sectional area A(m2)	IY (m4)	IZ (m4)	Modulus of elasticity E (GPa)	Poisson's ratio 𝑣
1.0	1 x 10-4	5 x 10-8	5 x 10-8	210	0.3

Table 5
Comparison of results for the cantilever beam subjected to an end moment.

Displacement (cm)	Liu (2014)	Wang and Rattanawangcharoen (2008)	Present
UTip	-15.82	-14.70	-14.45
VTip	-45.99	-47.19	-45.13

7.4 Cantilever beam subjected to a concentrated inclined load at the free end

In this numerical example, a cantilever beam subjected to an inclined end load is analyzed, as shown in Figure 13(a). Dowell and Traybar (1975) conducted an experimental study on this cantilever to invistigave its geometrically nonlinear behaviour. Rosen et al. (1987a,b) proposed some mathematical models to solve this problem and used the geometric data and material properties provided in Table 6. Noting that they used different values for EI_Y in both papers, the value of EI_Y in Table 6 is that given in Rosen et al. (1987b), which matches the experiment of Dowell and Traybar (1975). Ten beam elements and four loading increments are used in the present analysis. The error tolerance ${\mathrm{e}}_{\mathrm{r}}$ is chosen to be${10}^{-3}$. The results for the tip displacements V_Tip, W_Tip vs the inclination angle $\mathrm{\gamma }$_v for three different cases of end load P_TiP are shown in Figures 13(b), 13(c). The tip twist angle $\mathrm{\theta }$_Tip vs the inclination angle $\mathrm{\gamma }$_v, for three different cases of end load P_TiP is shown in Figure 13(d). The proposed results are compared with corresponding numerical results (model D) in Rosen et al. (1987a), and with the experimental results in Dowell and Traybar (1975). This comparison shows that the present results are in good agreement with the previously published results, especially with the experimental results.

Figure 13
Three-dimensional cantilever beam subjected to an inclined end load: (a) geometrical data, (b) V_Tip displcament curves, (c)W_Tip displcament curves, and (d) twist angle $\mathrm{\theta }$_Tip curves.

Table 6
The geometrical data and material properties of the beam.

Length L(m)	Cross-sectional dimensions		EIY (N.m2)	EIZ (N.m2)	GJ (N.m2)
	b(m)	h(m)
0.508	0.3175 x 10-2	0.127 x 10-1	2.4783	36.2695	2.9623

7.5 Circular cantilever beam subjected to a vertical out of Plane end load

In this example, a 45^o circular cantilever beam is considered as shown in Figure 14(a). The material and geometric data of the bend are given in Table 7. This example is a three-dimensional large displacement and rotation problem (Liu, 2014; Chan, 1992). Bathe and Bolourchi (1979) solved this problem using a moving frame through three successive rotations similar to Euler angles. They used 16 three-dimensional solid elements and eight straight beam elements. They also employed 60 load steps in the analysis of this bend. Chan (1992) analysed this problem using various formulations for comparison. However, he did not determine the number of elements used in this example. In the present study, four load increments and eight straight beam elements are used. The chosen value for the error tolerance ${\mathrm{e}}_{\mathrm{r}}$ is${10}^{-2}$. The present results are compared with the linear analysis results, the results in Bathe and Bolourchi (1979), and the results of the joint orientation approach in Chan (1992), as can be seen in Figures 14(b), 14(c), and 14(d). The present results are highly consistent with the compared results, especially those of Bathe and Bolourchi (1979) with a smaller number of load steps and elements.

Figure 14
45^o circular cantilever beam: (a) Geometrical data, (b) dimensionless load-displcament U/R curves, (c) dimensionless load-displcament V/R curves, and (d) dimensionless load-displcament W/R curves.

Table 7
The geometrical data and material properties of the circular beam (Bathe and Bolourchi, 1979).

Radius R(in)	Cross-sectional dimensions		Modulus of elasticity E (psi)	Poisson's ratio 𝑣
	b(in)	h(in)
100	1.0	1.0	107	0

8 CONCLUSIONS

This paper has presented an adapted corotational finite element formulation for spatial beams and frames with geometrically nonlinear behavior, subjected to static loads. Elastic and isotropic properties of the materials have been postulated. The distortion of cross-section, shear, and warping effects have all been disregarded. Through the concept of the used corotational frame, which continually updates and moves with each element, this approach significantly decreases the active degrees of freedom in the local displacement. The corotational frame additionally distinguishes between the deformational displacement and the rigid body displacement. The small strain hypothesis has been used since the deformational displacement is always small relative to the corotational frame. Also, the load steps and the number of elements both have been adjusted to keep the rotational displacement relative to the element chord relatively small. As a result, the vectorial addition can be applied to the rotational displacements and the stiffness matrices are symmetric. Consequently, storage of only the triangular part has been counted, which significantly reduces the computational time and accelerates the rate of convergence. The equilibrium equation has been derived using the virtual work principle. Two steps are suggested for the transformation of matrices and vectors from the fixed global frame to the corotational moving frame depending on Tait-Bryan angles, which have been computed in successive manner. The trigonometric formulas for all rotation angles including the special cases of the beam element have been investigated. Thus, the proposed method has addressed the issues with the rotation angles beyond the range [π/2, -π/2] and the vertical members. This contribution has been used to deal with some of these special cases, which have been classified as a singularity problem by many researchers in the field.

An incremental-iterative method based on the Newton-Raphson method has been used to solve the equilibrium equation. A MATLAB code has been constructed to accomplish this. This code provides a detecting function designed to control the rotational angle sign and locate the position of each member during the analysis. Each iteration, the stiffness matrices and the element coordinates' vector are continuously updated. Even though this update takes time in each iteration, it considerably reduces the total time of analysis. The proposed numerical algorithm's effectiveness and precision have been demonstrated by comparing the obtained results with published results of analytical formulations and experiments. Five numerical examples of large-displacement frames and beams have been solved and analyzed. Although in the analysis of these problems a reasonable number of elements were used, accurate results have been obtained in comparison with the published results.

The kinematic representation, the transformation method introduced, the symmetric stiffness matrices, and the mathematical derivation of the equilibrium equations constitute a distinctive blend of corotational formulation, setting it apart from other formulations found in the existing literature. A relatively rapid convergence rate has been observed in the proposed method since it does not depend on the well-known joint orientation method that requires a special parametrization of the finite rotations which considerably increase the computational time. The proposed method is not only simple and timesaving, but also it is highly effective and accurate. It seems that the method presented in this research has the potential to become a valuable tool for the analysis and resolution of numerous engineering applications due to its simplicity, accuracy, and effectiveness.

APPENDIX 1 - THE SHAPE FUNCTIONS MATRIX

${\mathbf{N}}_{\mathrm{i}}=\mathrm{}\left[\begin{array}{c}\begin{array}{cccccc}{\mathrm{}N}_1& 0& 0& 0& 0& 0\\ 0& {\mathrm{}N}_2& 0& 0& 0& {\mathrm{}N}_3\\ 0& 0& {\mathrm{}N}_2\mathrm{}& 0& -N_3& 0\\ 0& 0& 0& {\mathrm{}N}_1\mathrm{}& 0& 0\\ 0& 0& {{\mathrm{}N}_2}^{\mathrm{\text{'}}}& 0& {{-N}_3}^{\mathrm{\text{'}}}& 0\\ 0& {{\mathrm{}N}_2}^{\mathrm{\text{'}}}& 0& 0& 0& {{\mathrm{}N}_3}^{\mathrm{\text{'}}}\end{array}\begin{array}{cccccc}{\mathrm{}N}_4& 0& 0& 0& 0& 0\\ 0& {\mathrm{}N}_5& 0& 0& 0& {\mathrm{}N}_6\\ 0& 0& {\mathrm{}N}_5\mathrm{}& 0& -N_6& 0\\ 0& 0\mathrm{}& 0& {\mathrm{}N}_4& 0& 0\\ 0& 0& {{\mathrm{}N}_5}^{\mathrm{\text{'}}}& 0& {{-N}_6}^{\mathrm{\text{'}}}& 0\\ 0& {{\mathrm{}N}_5}^{\mathrm{\text{'}}}& 0& 0& 0& {{\mathrm{}N}_6}^{\mathrm{\text{'}}}\mathrm{}\end{array}\\ \\ \end{array}\right]$

where and ${\left(\right)}^{\text{'}}$ is the first derivatives with respect to ${\stackrel{-}{\mathrm{x}}}_{\mathrm{i}}$.

APPENDIX 2 - THE AXIAL AND BENDING STIFFNESS MATRIX

${\mathbf{k}}_1=\left[\begin{array}{c}\begin{array}{cccccc}\frac{\mathrm{E}{a}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0& 0& 0& 0\\ 0& \frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& 0& 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}\\ 0& 0& \frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0\\ 0& 0& 0& \frac{\mathrm{G}{\mathrm{J}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0\\ 0& 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& \frac{4\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0\\ 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& 0& 0& \frac{4\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}\end{array}\begin{array}{cccccc}-\frac{\mathrm{E}{a}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0& 0& 0& 0\\ 0& -\frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& 0& 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}\\ 0& 0& -\frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0\\ 0& 0& 0& -\frac{\mathrm{G}{\mathrm{J}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0\\ 0& 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& \frac{2\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0\\ 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& 0& 0& \frac{2\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}\end{array}\\ \begin{array}{cccccc}-\frac{\mathrm{E}{a}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0& 0& 0& 0\\ 0& -\frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& 0& 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}\\ 0& 0& -\frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0\\ 0& 0& 0& -\frac{\mathrm{G}{\mathrm{J}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0\\ 0& 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& \frac{2\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0\\ 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& 0& 0& \frac{2\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}\end{array}\begin{array}{cccccc}\frac{\mathrm{E}{a}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0& 0& 0& 0\\ 0& \frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& 0& 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}\\ 0& 0& \frac{12\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^3}& 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0\\ 0& 0& 0& \frac{\mathrm{G}{\mathrm{J}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0& 0\\ 0& 0& \frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& \frac{4\mathrm{E}{{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}& 0\\ 0& -\frac{6\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}^2}& 0& 0& 0& \frac{4\mathrm{E}{{\mathrm{I}}_{\mathrm{z}}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{o}}}\end{array}\\ \end{array}\right]$

where $\mathrm{E}$ is the modulus of elasticity, G is the modulus of rigidity, $a_i$is the cross-sectional area, ${{\mathrm{I}}_{\mathrm{y}}}_{\mathrm{i}}$and ${{\mathrm{I}}_{\mathrm{Z}}}_{\mathrm{i}}$are the moment of inertia about $\widehat{\mathrm{y}}$_i and $\widehat{\mathrm{z}}$_i axes, and ${\mathrm{J}}_{\mathrm{i}}$ is the polar moment of inertia.

APPENDIX 3 - THE GEOMETRIC STIFFNESS MATRIX

${\mathbf{k}}_2=\frac{{fx}_2}{L_{\mathrm{o}}}\mathrm{}\left[\begin{array}{c}\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& \frac{6}{5}& 0& 0& 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}\mathrm{}\\ 0& 0& \frac{6}{5}\mathrm{}& 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0\\ 0& 0& 0& 0\mathrm{}& 0& 0\\ 0& 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& -\frac{2\mathrm{}}{15\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2& 0\\ 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& 0& 0& \frac{2\mathrm{}}{15\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2\end{array}\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& -\frac{6}{5}& 0& 0& 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}\\ 0& 0& -\frac{6}{5}\mathrm{}& 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0\\ 0& 0\mathrm{}& 0& 0& 0& 0\\ 0& 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& -\frac{1\mathrm{}}{30\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2& 0\\ 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& 0& 0& -\frac{1\mathrm{}}{30\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2\end{array}\\ \begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& -\frac{6}{5}& 0& 0& 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}\mathrm{}\\ 0& 0& -\frac{6}{5}\mathrm{}& 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0\\ 0& 0& 0& 0\mathrm{}& 0& 0\\ 0& 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& -\frac{1}{30\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2& 0\\ 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& 0& 0& -\frac{1\mathrm{}}{30\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2\end{array}\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& \frac{6}{5}& 0& 0& 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}\\ 0& 0& \frac{6}{5}\mathrm{}& 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0\\ 0& 0\mathrm{}& 0& 0& 0& 0\\ 0& 0& \frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& \frac{2\mathrm{}}{15\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2& 0\\ 0& -\frac{1}{10\mathrm{}}{\mathrm{L}}_{\mathrm{o}}& 0& 0& 0& \frac{2\mathrm{}}{15\mathrm{}}{\mathrm{L}}_{\mathrm{o}}^2\end{array}\\ \end{array}\right]$

where ${\mathrm{f}\mathrm{x}}_{2\mathrm{}}$is the axial force of the second node in $\stackrel{-}{\mathrm{x}}$_i direction, which is presented in Eq. (12).

References

• Bathe, K. (1996). Finite Element Procedures. New Jersey, Prentice Hall.
• Bathe, K., Bolourchi, S., (1979). Large displacement analysis of three‐dimensional beam structures. International Journal for Numerical Methods in Engineering 14(7): 961-986.
• Belytschko, T., Hsieh, B., (1973). Non-linear transient finite element analysis with connected co‐ordinates.” International Journal for Numerical Methods in Engineering 7(3): 255-271.
• Benjamin, A. (1982). Análise não-linear geométrica de pórticos tridimensionais pelo método dos elementos finitos, M.Sc. Thesis (in Portuguese), Federal University of Rio de Janeiro, Rio de Janeiro, Brasil.
• Chan, S. (1992). Large deflection kinematic formulations for three-dimensional framed structures. Computer Methods in Applied Mechanics 95(1): 17-36.
• Crisfield, M. (1990). A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements. Computer Methods in Applied Mechanics and Engineering 81(2): 131-150.
• Crivelli, L. (1991). A Total-Lagrangian Beam Element for Analysis of Nonlinear Space Structures. Ph.D. Thesis, University of Colorado Boulder, USA.
• Dowell, E., Traybar, J., (1975). An experimental study of the nonlinear stiffness of a rotor blade undergoing flap, lag and twist deformations, U. S. Army Air Mobility Research Development Laboratory (AMES), Technical Report No. 1257, Princeton University, USA.
• Elkaranshawy, H., Dokainish, M., (1995). Corotational finite element analysis of planar flexible multibody systems. Computers & Structures 54(5): 881-890.
• Elkaranshawy, H., Elerian, A., Hussien, W., (2018). A corotational formulation based on Hamilton’s principle for geometrically nonlinear thin and thick planar beams and frames. Mathematical Problems in Engineering 2018.
• Goldstein, H., Poole, C., Safko, J., (1980). Classical mechanics, Addison-Wesley (Boston).
• Gu, J. (2004). Large displacement elastic analysis of space frames allowing for flexural-torsional buckling of beams, Ph. D. Thesis, Hong Kong Polytechnic University, Hong Kong.
• Jonker, J., Meijaard, J., (2013). A geometrically non-linear formulation of a three-dimensional beam element for solving large deflection multibody system problems. International Journal of Non-linear Mechanics 53: 63-74.
• Le, T., Battini, J., Hjiaj, M., (2014). A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures. Computer Methods in Applied Mechanics Engineering 269: 538-565.
• Leng, J., Wang, Q., Li, Y., (2022). A geometrically nonlinear analysis method for offshore renewable energy systems—Examples of offshore wind and wave devices. Ocean Engineering 250: 110930.
• Liu, T., Bai, J., (2022). Folding behaviour of a deployable composite cabin for space habitats-part 1: Experimental and numerical investigation. Composite Structures 302: 116244.
• Liu, Y. (2014). The development of the co-rotational finite element for the prediction of the longitudinal load factor for a transmission line system, Ph.D. Thesis, University of Manitoba, Canada.
• Magisano, D., Leonetti, L., Madeo, A., Garcea, G., (2020). A large rotation finite element analysis of 3D beams by incremental rotation vector and exact strain measure with all the desirable features. Computer Methods in Applied Mechanics 361: 112811.
• Mars, J., Koubaa, S., Wali, M., Dammak, F., (2017). Numerical analysis of geometrically non-linear behavior of functionally graded shells. Latin American Journal of Solids Structures 14: 1952-1978.
• Mondkar, D., Powell, G., (1977). Finite element analysis of non‐linear static and dynamic response. International Journal for Numerical Methods in Engineering 11(3): 499-520.
• Nunes, C., Soriano, H., Venancio Filho, F., (2003). Geometric non-linear analysis of space frame with rotation greater than 90°, with Euler angles and quasi-fixed local axes system. International Journal of Non-linear Mechanics 38(8): 1195-1204.
• Oran, C. (1973). Tangent stiffness in space frames. Journal of the Structural Division 99(6): 987-1001.
• Remseth, S. (1979). Nonlinear static and dynamic analysis of framed structures. Computers & Structures 10(6): 879-897.
• Rosen, A., Loewy, R., Mathew, M., (1987a). Nonlinear analysis of pretwisted rods using “principal curvature transformation,” Part II: Numerical results”. AIAA journal 25(4): 598-604.
• Rosen, A., Loewy, R., Mathew, M., (1987b). Nonlinear dynamics of slender rods. AIAA journal 25(4): 611-619.
• Santana, M., Sansour, C., Hjiaj, M., Somja, H., (2022). An equilibrium‐based formulation with nonlinear configuration dependent interpolation for geometrically exact 3D beams. International Journal for Numerical Methods in Engineering 123(2): 444-464.
• Shabana, A. (2020). Dynamics of multibody systems, Cambridge University press.
• Simo, J., Vu-Quoc, L., (1986). A three-dimensional finite-strain rod model. Part II: Computational aspects. Computer Methods in Applied Mechanics 58(1): 79-116.
• Trapper, P. (2022). A numerical model for geometrically nonlinear analysis of a pipe-lay on a rough seafloor. Ocean Engineering 252: 111146.
• Vo, D., Nanakorn, P., (2020). “A total Lagrangian Timoshenko beam formulation for geometrically nonlinear isogeometric analysis of planar curved beams.” Acta Mechanica 231: 2827-2847.
• Wang, Z., Rattanawangcharoen, N., (2008). The analysis of a transmission guyed tower line system subjected to downburst wind. Technical Report, University of Manitoba, Canada.
• Xiaohang, Q., Zhiteng, G., Tongguang, W., Long, W., Shitang, K., (2022). Nonlinear aeroelastic response analysis of 100-meter-scale flexible wind turbine blades. Acta Aerodynamica Sinica 40(4): 220-230.
• Yang, Y., Kuo, S., Wu, Y., (2002). Incrementally small-deformation theory for nonlinear analysis of structural frames. Engineering Structures 24(6): 783-798.
• Yang, Y., Lin, S., Chen, C., (2007). Rigid body concept for geometric nonlinear analysis of 3D frames, plates and shells based on the updated Lagrangian formulation. Computer methods in applied mechanics engineering 196(7): 1178-1192.

Edited by

Editor: Marco L. Bittencourt

Publication Dates

• Publication in this collection
  12 Feb 2024
• Date of issue
  2024

History

• Received
  17 Aug 2023
• Reviewed
  13 Nov 2023
• Accepted
  21 Dec 2023
• Published
  12 Jan 2024