Corotational Finite Element Dynamic Analysis of Space Frames with Geometrically Nonlinear Behavior Based on Tait–Bryan Angles

Te aim of this study is to compose a corotational fnite element formulation for space frames with geometrically nonlinear behavior under dynamic loads. Using a moving frame through three successive rotations similar to Euler angles is one of the oldest techniques; however, there are still some gaps that require attention, mainly due to singularity. Hence, alternative techniques had been developed, sometimes elusive and computationally expensive. In this paper, we went back to the old technique and flled the gaps. Tree-coordinate systems are used, i.e., the fxed global coordinate system, the fxed local coordinate system that is attached individually to every element, and the corotational local frame for each element that moves and rotates with the element. Te deformation is always small relative to the corotational frame. Te successive rotations between diferent coordinate systems are expressed using Tait–Bryan angles. Lagrange’s equation is used to derive the equation of motion, and the stifness and mass matrices are obtained using the Euler–Bernoulli beam model. A MATLAB code is developed based on the Newton–Raphson method and the Newmark direct integration implicit method. In traditional techniques, singularity is attained when any rotation angle in the fxed local frame approaches π /2, and if any is greater than π /2, the techniques could fail to specify the location of the element. In this paper, each case is treated with a proper procedure, and special handling of trigonometric formulations prevents singularity and correctly specifes the location of elements in all situations. Diferent examples of beams and frames are analysed. While the method is not intricate, it is timesaving, is highly efective, provides more stable and robust analysis, and gives sufciently accurate results. Compared to the parametrization of the fnite rotations technique, the method has a signifcant reduction in the convergence rate because it avoids the storage of joint orientation matrices.


Introduction
In recent years, the development of lightweight highstrength materials has attracted many industries [1][2][3]. No wonder, industry always faces new challenges and needs to reduce the cost of designs. Such designs inevitably experience large displacements but with small strain. Tat is why geometrically nonlinear analysis plays a crucial role in such designs, which cannot be analysed using the traditional linear analysis. Remarkably, it remains a fertile ground for research due to the continuous demand for accurate, adaptable, and computationally inexpensive geometrically nonlinear formulations for treating innovative applications.
To elucidate this point, Leng et al. [4] pointed to the signifcant efect of the geometric nonlinearity on the fexible ofshore structures and devices that cannot be ignored. Trapper [5] studied a pipe lay on a seafoor that experienced large deformations. Accordingly, he used a geometrically nonlinear model to calculate the maximum internal forces locations. Liu [6] investigated a geometrically nonlinear fnite element formulation to determine the dynamic response of a guyed transmission line system that contains large displacements from the vibration of cables. With the unprecedented growth in renewable energy, increasing the scale of renewable energy devices becomes an urgent need. As a result, Xiaohang et al. [7] analysed a 100-meter fexible wind turbine blade. Tis investigation reveals that geometric nonlinearity plays an indispensable role in the computation of the dynamic response of such giant blades. Due to the increasing interest in deep space exploration, there is a continuous need for space applications and habitats. It is worthy to say that such applications of light-weight materials contain large displacements. Terefore, Liu and Bai [8] considered the geometric nonlinearity during their experimental and numerical investigation of a deployable composite cabin for space habitats. Te important question that can be addressed here is why shall we develop geometrically nonlinear models in the existing of various fnite element commercial softwares? Te answer is not complicated; there is still a need to investigate more versatile, efcient, and time-saving models that can successfully handle diferent engineering problems which may be exclusive or intractable in some cases using these commercial programs.
According to Crisfeld [9], Turner et al. [10] were the frst to study geometrically nonlinear fnite element analysis in the sixtieths of the last century. In fact, Argyris [11] made considerable premature contributions in genuine geometrically nonlinear analysis procedures. Essentially, the geometrically nonlinear fnite element analysis of structures was covered in some notable books by Oden [12], Bathe [13], and Crisfeld [9].
In general, there are two diferent forms, in continuum mechanics, to describe the motion of a body, namely, Eulerian and Lagrangian formulations. Te Eulerian formulations are usually used in fuid mechanics, while Lagrangian formulations are used in most other engineering felds. With regard to geometrically nonlinear analysis, the Lagrangian formulations are commonly used in the form of total Lagrangian, updated Lagrangian, and corotational formulations. In total Lagrangian formulations [14][15][16][17], the system equation terms are defned in terms of the fxed global frame that does not change through the analysis. Tis generates relatively large strains, displacements, and rotations that need special procedures to handle.
While in updating Lagrangian formulations [18,19], the terms of the equation are defned relative to a frame that is updated with the last accepted solution. Tis reference frame does not change during the solution cycles. As a result, the system equations are much simpler than the corresponding equations in total Lagrangian formulations. However, if the displacement from the current confguration to the last equilibrium confguration is large, a basic assumption is violated and accordingly these formulations experience also some complexities. To avoid such complexities, corotational formulations [20][21][22][23][24][25][26][27][28][29][30][31] provide a simple kinematics description method for large displacement analysis. Tese formulations are based on the theory of small strain. Te corotational local frame translates and moves with each element but does not deform with it. Consequently, one obvious advantage of this formulation is that it can easily flter the rigid-body motion from the deformational motion, which is always small relative to the corotational local frame. Tis frame, which is continuously updated, can be defned by many diferent methods [23]. Remseth [17] used an approximate vectorial assumption to deal with three-dimensional rotations. Terefore, this method is not applicable for large rotations. Tus, he limited his approach to rotations to the range of 12-15 degrees. However, the formulation of the three-dimensional beam element is not just a simple extension of the two-dimensional formulation, mainly because of the complexity of the three-dimensional large rotations. More specifcally, the three-dimensional large rotations are noncommutative and nonadditive. To handle this problem, Oran [24] used a joint orientation matrix to describe a set of orthogonal axes that are rigidly attached and deformed with the joints of a structure. Te element nodal rotations are determined from the angles between this set of orthogonal axes and the member axes. Tis procedure is improved by Crisfeld [25], Le et al. [26], Jonker and Meijaard [27], and Hsiao et al. [28]. Tough this method does not put restriction on the size of the time step and can use smaller number of elements, a key stone for the method is the need for a special parametrization of the fnite rotations. It also increases the computational time because of storing the joint orientation matrices and parametrization of the fnite rotations. Bathe and Bolourchi [29] defned a moving local frame through three successive rotations similar to Euler angles, sometimes called Tait-Bryan angles, Bunge, or other conventions. However, they did not formulate trigonometric rules for all rotation angles. Benjamin [30] contributed to the solution of this problem by obtaining cosines and sins of all rotation's angles. However, he did not specify a control sign for the cosine of rotation angles which can be obtained using two hypotenuses of right-angle triangles. Nunes et al. [31] controlled the sign of the cosine of rotation angles to overcome the problem of the cosine of an angle greater than π/2. Nevertheless, they did not specify clearly how to determine the transformation matrices in the case of vertical members, when the rotation angle is equal π/2 or −π/2 exactly, which is very important in the modeling of three-dimensional structures. Furthermore, they did not solve any three-dimensional problem by examining their motion description method. Due to such problems, Simo and Vu-Quoc [16] identifed the problem of singularity in case of using this method. Tat is why many authors used parametrization of fnite rotations, such as the authors in [16,[25][26][27].
To express the stifness and mass matrices, the Euler-Bernoulli beam theory was extensively used in the formulations by many researchers (see, for example, [21,22,27,28,32,33]) since it simplifed the analysis and in the same time obtains adequate results. In the Euler-Bernoulli beam theory, the material is considered isotropic and elastic, and the cross section of the elements is uniform. A normal plane section on the centroid axis before deformation remains plane after deformation and normal to the axis is employed. Warping and cross-sectional distortion are not considered. To investigate the efect of shear formulation, the Timoshenko beam model has been used by other researchers (see, for example, [15,22,34]).
In this paper, a relatively accurate and simple corotational formulation for three-dimensional fnite element formulation has been developed. Te stifness and mass matrices are evaluated using the Euler-Bernoulli beam model. Te transformation procedure is based on the Tait-Bryan angles successive rotations [23,29]. Tis procedure is employed in two main stages to transform vectors and matrices from the fxed global frame to the moving corotational frame. Te frst stage is the transformation from the fxed global frame to the fxed local frame, and the second stage is the transformation from the fxed local frame to the moving corotational local frame. Te transformation procedure also depends upon updating the coordinates with every equilibrium confguration during the analysis. In order to formulate a consistent model, the trigonometric rules for special cases of spatial beam elements are considered with the control sign. Tese rules are used to calculate the rotations matrices. Meanwhile, this method is simple and does not need parametrization for fnite rotations, but it requires regulating the time step and number of elements in order to decrease the relative chordal rotations during the analysis. Te equation of motion is formulated using Lagrange's equation. An incremental iterative procedure is used to solve the equation of motion. Tis procedure is based on the full Newton-Raphson method and the Newmark direct-time integration method. Te MATLAB code is developed for this purpose. Tis code involves a relatively rapid convergence rate for equilibrium because it avoids storing joint orientation matrices and parametrizing of fnite rotations, which are often associated with parametrized formulations. Tis section represents a general introduction and reveals the importance of studying geometric nonlinearity. Te spatial beam element motion description method, which involves the coordinate systems and the displacement interpolation, is presented in Section 2. Te transformation procedure between diferent coordinate systems based on Tait-Bryan angles can be seen in Section 3. Te strain energy and stifness matrix are derived in Section 4. Te kinetic energy and the mass matrix are presented in Section 5. Ten, Lagrange's equation is used to derive the equation of motion in Section 6. Section 7 represents the solution strategy. To expose the efciency and validate the accuracy of the proposed model, fve numerical examples are solved and compared with the published results in Section 8. Finally, the conclusions are presented in Section 9.

Kinematics Description of the 3D
Beam Element 2.1. Basic Assumptions. Te following assumptions are employed to formulate the spatial beam element: (i) Te material is isotropic, elastic, and homogeneous (ii) Te cross section of the beam element is symmetric about both axes (iii) Euler-Bernoulli beam assumptions, which state that a normal plane section on the centroid axis before deformation remains plane after deformation and is normal to that axis, are employed. (iv) Warping, cross-sectional distortion, and shear efect are not taken into account Hence, the small strain theory is the basis for the corotational formulation used in the analysis. Accordingly, deformational and rotational displacements are always small with respect to the corotational frame. Appropriate element sizes and time steps are chosen to ensure that these conditions remain valid and the results are accurate.

Coordinate Systems.
After discretization of the structure into fnite elements, the i th beam element can be defned with two end nodes (n = 1 and 2). Every node has six degrees of freedom and is defned with respect to three frames as shown in Figure 1. Tese coordinate systems are the fxed global coordinate system associated with the fxed global frame (X, Y, Z), the fxed local coordinate system associated with the fxed local frame at time = 0 (x , and the moving local coordinate system associated with the corotational local frame (x i , y i , z i ). Tis 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 confgurations used in dynamics: an initial confguration at time � 0, the j th equilibrium confguration at time � t, and a current confguration at time � t + ∆t. Te element's initial length is L o , and after deformation in the current confguration, the element length is equal to the arc length S i , while L c is the current chord length.
For the current confguration, as shown in Figure 2, the nodal displacement vector for the i th beam element in the fxed global coordinate system is given by the following equation: 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 θX n . θY n and θZ n (n = 1, 2) are the counterclockwise rotations about X, Y, andZ axes, respectively. Te nodal incremental displacement vector of the i th beam element in the global coordinate system is defned as follows: Shock and Vibration 3 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 ∆θX n .∆θY n and ∆θZ n (n = 1, 2) are the counterclockwise incremental rotations about X, Y, and Z axes, respectively. Te nodal displacement vector D i can be updated by the following equation: where D j i is the nodal displacement vector for the i th beam element in the fxed global coordinate system at the j th equilibrium confguration. Te nodal displacement vector D i is divided into the nodal translational displacement

Shock and Vibration
vector Dt i and the nodal rotational displacement vector Dr i , which can be written as follows: Similarly, the ∆D i vector can be divided into the nodal incremental translational displacement vector ∆Dt i and the nodal incremental rotational displacement vector ∆Dr i as follows: Tus, the nodal coordinate vector of the i th beam element in the fxed global system can be updated continuously as follows: where X i is the vector of the nodal coordinates of the i th beam element relative to the fxed global frame, at the current confguration, and X j i is the vector of the nodal coordinates relative to the fxed global frame, at the j th equilibrium confguration.
Te nodal displacement vector of the i th beam element in the element corotational local coordinate system, at the current confguration, is as follows: where u n (n � 1, 2), v n (n � 1, 2), and w n (n � 1, 2) are the displacement translational components in x i , y i , and z i directions, respectively, and θx n . θy n . and θz n (n = 1 and 2) are the counterclockwise deformational rotations after eliminating the rigid-body rotations. Te internal elastic force vector for the i th beam element in the fxed global coordinate system, at the current confguration, can be written as follows:

Shock and Vibration 5
Te internal elastic force vector of the i th beam element in the element corotational local coordinate system, at the current confguration, is as follows: where the internal elastic force vector components in both the fxed global coordinate system and the corotational local coordinate system are shown in Figure 2 with their positive signs.

Displacement
Interpolation. Te classical Hermitian shape functions are used to relate the element axial elongation u (x i ), lateral defections of the centroid axis v(x i ) and w(x i ) shown in Figure 3, and rotation about the centroid axis θ x (x i ) to the element nodal displacement vector d i as follows: Te components of the shape functions N i of the i th beam element are given by the following equation: where ξ is given by where L c is the element current chord length. Te function ξ nodal values for the i th beam element are shown in Figure 4. Transverse displacements can be determined, at any point along the element, using shape functions and the nodal displacement values in the local coordinate system.
Due to the nature of the attached corotational local frame, as shown in Figures 1 and 3, the displacement components v n (n � 1,2) and w n (n � 1,2) are equal to zero. Furthermore, the axial displacement of the frst node is chosen to be zero while the axial displacement of the second node is u 2 . Consequently, the nodal displacement vector d i has only seven nonzero components which will simplify the analysis as can be seen in the next section. Tis vector can be written as follows: Tus, equations (12) and (13) can be simplifed to Te arc length of the i th beam element S i can be expressed by where v ′ and w ′ are the frst derivatives of the functions v(x i ) and w(x i ) with respect to x i . Tis integral is evaluated using the Gaussian integration scheme. Te axial elongation of the i th beam element can be defned as follows: Tis equation can also be written in terms of chord lengths as follows:

Shock and Vibration
where L o is the element initial length and b i is the element axial elongation due to the bowing efect, which can be determined in terms of rotations [35] as follows: Hence, the axial displacement u 2 in equation (17) is given by

Transformation Procedure
Te transformation procedure depends upon updating the coordinates with every equilibrium confguration during the analysis. Two main stages are employed here to perform transformation from the fxed global frame to the moving corotational local frame. Te frst stage is the transformation from the fxed global frame to the fxed local frame, and the second stage is the transformation from the fxed local frame to the moving corotational local frame. Assuming that V d is a 3D vector associated with the fxed global frame, the relation between the fxed global frame (X, Y, Z) and the fxed local frame ( where r o is an orthogonal matrix (3 × 3) which can be determined from the direction cosines of the fxed local frame relative to the fxed global frame. For a three-dimensional frame element, this matrix turns into a (12 × 12) matrix as follows: Similarly, the relation between the fxed local frame ( 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 fxed local frame. For a threedimensional frame element, this matrix turns into a (12 × 12) matrix as follows: One can write the vector v d in terms of V d as follows: where Terefore, the transformation matrix for the threedimensional i th beam element from the fxed global frame to the moving corotational local frame, at the current confguration, can be expressed by 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.   (1) Rotation β o of the (X, Y, Z) coordinate axes about the Y axis: Tis rotation places the X and Z axis along X β o and Z β o , respectively, while the Y axis remains the same, as shown in Figure 5.

Transformation from the Fixed Global
with respect to the global frame (X, Y, Z), shown in Figures 5 and 8, the rotation matrix r β o can be determined as follows: where , and the rotation matrix r c o can be obtained as follows: where axes about X c o axis: Tis rotation places the Y c o and Z c o axis along y ∧ i and z ∧ i , respectively, while the X c o axis remains the same. It has to be noticed that x ∧ i coincides with the axis X c o , as shown in Figure 7.
By assuming a point P lies in the (x ∧ i , y ∧ i ) plane, as shown in Figure 8, the point P coordinates relative to the starting point of the i th beam element with respect to the fxed global system coordinate (X, Y, Z) can be written as follows: Consequently, the coordinates of point P with respect to the (X c o , Y c o , Z c o ) frame, as shown in Figure 9, can be obtained as follows: which leads to Now, using the direction cosines of , the rotation matrix r α o can be determined as follows: where Hence, the rotation matrix r o can be obtained as follows: By substitution, this matrix takes the form as follows: 3.1.1. Special Cases. It should be noted that the rotation angle α o is insignifcant, in the case of a circular cross section element. Tus, the rotation matrix r o can be calculated as follows: Tere is another special case where the initial position of the element is vertical, in the Y-axis direction. In order to get r o , there are only two successive rotations not three as in the general case. Te frst rotation c o is either 90°or 270°, as shown in Figure 10, depending on whether C Y is +1 or −1. Te other rotation is α o , which is shown in Figure 11, and can be determined using a reference the point P that lies in the (x ∧ i , y ∧ i ) plane. In this case, equations (39) and (40) are modifed as follows: Tus, the matrix r o in equation (43) can be written as follows: Substituting equation (42) or equation (45) into equation (25), the matrix T o can be determined. (1) Using the direction cosines of the β c frame (

Transformation from the Fixed Local Frame to the Moving
, which is shown in Figure 12, the rotation matrix r β c can be determined as follows: where Of particular interest is to notice that in [29], a difculty arises when angle β c > 90°since this reference gave only an expression of the cosine of the angle. In this work, an expression for sine is also given. Both trigonometric relations can specify the location of the element exactly. As shown in Figure i are the i th beam element relative translational displacements with respect to the fxed local frame. Tese relative displacements of the i th beam element can be determined from the relative displacement with respect to the fxed global system calculated from the Dt i vector in equation (4) and the transformation matrix r o in equation (41) or equation (43) for the vertical member as follows: where (2) Similarly, one can use the direction cosines of the c c frame ( , and the rotation matrix r c c can be obtained as follows: When employing the technique outlined in [29], a complication rises when the angle c c > 90°, as this reference solely provides an expression for the sine of the angle. In this study, however, we present an expression for the cosine as well, enabling precise determination of the element's position. Tese two trigonometric relationships can be expressed as follows: Shock and Vibration cos c c � (SN) where Te rotation matrix for relative translational displacements can be written as follows: In case when an element becomes vertical, that means it is parallel to the Y-axis, the rotation β c vanishes, and the rotation c c is either 90°or 270°, depending on the element position. Traditionally, this case could create singularity and it had been the source of many difculties [31], and the authors of this research work suggested preventing the rotation c c to be either 90°or 270°. In this work, this problem is solved by letting the code search for the alignment of the element, that means to specify if the rotation c c is either 90°or 270°. Te matrix r d c can be rewritten as follows: where C Y ′ can be specifed using the current vector of the nodal coordinates in equation (7) as follows: Hence, the value of C Y ′ is either +1 and the rotation c c is 90°or −1 and the rotation c c is 270°.
(3) When the third rotation angle α c is computed by assuming a reference point P as has been done for α o , it is very hard to assume a point in the (x i , y i ) plane every iteration with diferent consequence positions and directions. Accordingly, the model experienced some difculties in converging using this method. Tus, an incremental procedure [29] is employed here to compute this angle as follows: where ∆θx 1 and ∆θx 2 are the incremental twist angles about the x i axis. Tese incremental angles can be determined using the following procedure. First, the incremental rotational displacement vector, with respect to the fxed global frame in equation (6), is transformed to the corresponding vector in the fxed local frame using the procedure indicated in equation (24) as follows: Ten, the incremental rotational vector relative to the fxed local frame is transformed to the corotational local frame as follows: Terefore, ∆α c in equation (60) can be determined using the following relation: Consequently, the angle α c , at the current confguration, can be determined as follows: where α j c is the twist rotation about the x i axis, at the j th equilibrium confguration. Tis equation is used as a predictor between two successive time steps. Also, it is used as a corrector between two successive iterations.
Tus, the rotation matrix r α c can be determined as follows: Te rotation matrix r c in equation (27) can be expressed as follows: Substituting equations (65) and (57) (or equation (58) in case of vertical member) into equation (66), the matrix r c can be determined. Hence, the transformation matrix T c in equation (27) has been specifed. Ten, the transformation matrix T r in equation (30) has been determined.

Strain Energy and the Stiffness Matrix
Considering isotropic elastic materials, the constitutive relation between the stress vector σ i and the strain vector ϵ i of the i th beam element can be defned by the following equation: where E i is the symmetric matrix of the elastic coefcients. Te strain vector is given by the following equation: where D f is the diferential operator matrix and v q is the deformation vector, which can be defned using the shape functions in equation (15) as follows: where N i is the shape functions matrix, which is given in Appendix A.1. Substituting equation (69) into equation (68), the strain vector can be expressed as follows: Combining equations (67) and (70), the stress vector can be rewritten as follows: Te strain energy for the i th beam element Π i is given by the following equation: where V i is the volume. Substituting equations (70) and (71) into equation (72), the strain energy can be expressed as follows: Shock and Vibration Te element local displacement vector d i is independent of the volume V i . Subsequently, equation (73) can be simplifed to the following form: where k i is the symmetric element stifness matrix in the corotational local coordinate system, which can be defned as follows: and B is defned as follows: Equation (74) can be written as follows: where K i is the symmetric element stifness matrix in the fxed global coordinate system. For the beam element used in this research work, the element stifness matrix in the local coordinate system k i can be expressed as follows: where k 1 is the axial and bending stifness matrix and k 2 is the geometric stifness matrix for the i th beam element. Stifness matrices are attached in Appendix A.2 and A.3.

The Kinetic Energy and the Mass Matrix
Te kinetic energy for the i th beam element KE i is given by the following equation: where _ () is the diferentiation with respect to time t, ρ is the density of the i th beam element, and _ R i is the velocity of a general point in the i th beam element with respect to the global coordinate system. Te velocity of a general point in the i th beam element with respect to the corotational local coordinate system _ r i can be expressed as follows: where _ r i can be written in the following form: where _ d i is the nodal velocity vector of the i th beam element with respect to the local coordinate system. Te vector _ d i can be expressed as follows: where _ D i is the nodal velocity vector of the i th beam element in the global coordinate system. Combining equations (79) and (80), the kinetic energy can be written as follows: Te velocity vector _ D i and the transformation matrix T r are independent of the volume V i . Subsequently, equation (83) can be simplifed to the following form: where m i is the symmetric element mass matrix in the local corotational coordinate system defned as follows: Equation (84) can be written in the following form: where M i is the symmetric element mass matrix in the fxed global coordinate system defned as follows: For the beam element used in this research work, the element mass matrix in the local coordinate system m i can expressed as follows: where m 1 is the ith beam element mass matrix for the translational inertia and m 2 is the mass matrix of the i th beam element for the rotational inertia. Te mass matrices are attached in Appendices A. 4

and A.5.
It is worth mentioning that the rotation is noncommutative. However, in our analysis, we used the threecoordinate system that separates the large rigid-body motion from the internal deformation. We restrict our analysis to the small strain theory, as we mention in subsection (2.1), and the internal rotational displacements within each element in each time step are always small. Terefore, the usual vector addition rules can be applied to these rotational displacements. Hence, the obtained stifness and mass matrices are symmetric.

The Equation of Motion
Lagrange's equation reported by [36,37] is used to derive the equation of motion in this section. It can be written for the i th beam element as follows: where KE i is the kinetic energy, Π i is the strain energy, WD i is the work done by the external forces, _ D i is the nodal velocity vector in the global coordinate system, and D i is the nodal displacement vector in the global coordinate system. Te work done by the external forces on the i th beam element is given by the following equation: where F P i is the vector of the external applied forces on the i th beam element in the fxed global coordinate system.
Using equation (86), one can write the frst two terms in the left-hand side of the previous relation as follows: Using equation (77), the third term in the left-hand side of equation (89) can be written as follows: Also, one can write Substituting equations (91)-(93) in equation (89), one obtains One can write where F e i is the internal elastic force vector in the fxed global and f e i is the internal elastic force vector in the local coordinate systems, which is given by Substituting equation (95) in equation (94) gives Equation (97) is the equation of motion of the i th beam element. Assembling the element force vectors, acceleration vectors, and mass matrices leads to the equation of motion of the overall structure, which takes the following form: where F P is the vector of the external applied forces of the entire structure and F e is the internal elastic force vector of the entire structure. Both F P and F e are defned in the fxed global coordinate system. M is the mass matrix in the global coordinate system of the entire structure, and € D is the acceleration vector in the global coordinate system of the entire structure. Rayleigh damping is used in equation (98) to read where C is the damping matrix in the global coordinate system of the entire structure and _ D is the velocity vector in the global coordinate system of the entire structure. Te damping matrix C is defned in terms of the stifness matrix and the mass matrix as follows: where μ and η are the damping coefcient which can be determined from the vibration modes of the system. K is the assembled stifness matrix in the global coordinate system of the entire structure.

Numerical Algorithms
Te equation of motion is solved using an incremental iterative procedure. Tis procedure is based on the full Newton-Raphson and the Newmark direct integration implicit method [13,38]. Equation (99) can be rewritten as follows: where ψ is the out of balance force. Te iteration equilibrium convergence criterion is given by the following equation: where ψ f is the reference out of balance force, which is assumed to be the out of balance force in the frst iteration, and 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. Tey are Shock and Vibration required to allow the solution to converge and satisfy the specifed tolerance during the solution iterations. Tat is why many authors have used tangent, secant, or even initial stifness matrices in their nonlinear FE formulations. Using exact matrices typically costs more time because of the storage of nonsymmetric matrices compared with the storage of only the triangular part in the symmetric matrices. Terefore, we are confdent that the restriction we made in subsection (2.1), which produced symmetric stifness matrices, allows the whole solution to go faster without afecting the overall accuracy of the results. Trough the solution of various numerical examples in the next section, we have assessed the efectiveness and the accuracy of this method.
At the beginning of the solution, it is assumed that the displacement, velocity, and acceleration vectors are null.  (25) and (41). In all other iterations, calculate r c and T c using equations (27) and (66). (iii) Determine the transformation matrix T r for each member using equation (30).  (102) is satisfed, stop the iteration and go to step ix. Otherwise, start the following iteration: (a) Using the Newton-Raphson method, a displacement corrector vector R is calculated as follows: where K is the efective matrix that can be defned as follows: where τ and ∃ are the Newmark parameters [38]. (b) Update the incremental displacement vector as follows: where ∆D 0 is the incremental displacement vector in the previous iteration, which is considered to be zero in the frst iteration. (c) Extract vector ∆D i for each element from vector ∆D. Consequently, one can update the vectors D i and X i using equations (3) and (7).
(g) Eliminate the rigid-body rotations from the rotational components θY n and θZ n in the vector D i , for each element, as follows: (h) Transform the pure rotations θX n .θY ̖ n and θZ ̖ n , which are relative to the fxed global frame, to determine the corresponding rotational components θx n . θy n and θz n for vector d i in the current corotational local frame which are always relative to the element chord, using the procedure detailed in equations (24) and (26). (i) Te axial displacement u 2 in equation (17) is obtained from equation (23). (j) Using equations (17) and (96), d i and f e i can be determined, respectively. (k) Using the Newmark direct integration method [38], _ D N+1 and € D N+1 can be calculated as follows: (viii) Going to the start of step iv again. 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. Tus, it can easily select the appropriate trigonometric rules and the sign of rotations. Te detection function specifes the position of the element using the nodal coordinate vector X i given in equation (7). When the relative diference between the coordinate in X and Z direction of the element end nodes approaches zero together, the code detects that the angle β c turns to be zero and the element is in the direction of the Y axis. In this case, the code searches for the alignment of the element based on the sign of C Y ′ in equation (59). In addition, when the relative displacement V ∧ r i turns out to be zero, the angle c c vanishes. Likewise, if the relative displacement W ∧ r i turns out to be zero, the angle β c vanishes; however, the element in general is not vertical. Hence, the program deals with these special cases separately.
Tis function also controls the angle c c when the rotation is outside the interval [π/2, −π/2]. Te sign SN in equation (53) specifes the cosine of the angle c c to meet the corresponding element position during the motion because the terms P 1 P 2 ′ and L c are always positive which cannot refect the real sign of cosine of c c . At the same time, the sine of c c is already specifed with the sign of V ∧ r i . Tis function conserves the code to converge efciently. Tis numerical algorithm steps are organized in Figure 13.

Simple Supported Beam Subjected to a Concentrated
Step Load. Te frst example is a simply supported beam whose geometric properties are L � 12 m, A � 8.06 X 10 − 3 m 2 , and I Y � I Z � 1.858 X 10 − 4 m 4 and material properties are E � 210 GPa, G � 80.775 GPa, and ρ � 7850 kg/m 3 , as shown in Figure 14. Te beam is subjected to a vertical concentrated load, at point C, which increases linearly to reach the value of 10 KN at 0.15 second before being steady, as can be seen in Figure 15. Six beam elements are used in the analysis, and time step ∆t � 5 X 10 − 4 s and the error tolerance e r is chosen to be 10 − 4 . Liu [6] solved this problem using a corotational formulation based on Euler's theorem. He compared his results with the theoretical results [6]. Te present results are compared with the results in [6], as shown in Figure 16. Tis comparison reveals that the present results remarkably agreed with the theoretical results in [6].

Clamped-Clamped Beam Subjected to a Concentrated
Vertical Load at the Midspan. In this part, a clamped-clamped beam is analysed. Te geometric data of the beam are shown in Figure 17. Te beams modulus of elasticity equals 30,000 ksi, density is 0.098 lb/in 3 , and Poisson's ratio is considered to be zero. Tis beam is subjected to a dynamic concentrated step load. Tis load is assumed to be applied suddenly to the midspan, as can be shown in Figure 18. Tis problem is solved by Mondkar and Powell [39] using various time steps. A secant stifness concept is utilized to solve this problem by Chan [40]. Te present results are obtained using ten element and error tolerance e r � 10 − 2 . Te present results are compared with the results of [39] and the linear analysis results, using the same time step ∆t � 50μs. Te dynamic load is applied over a period of 5000 μs (0.005 s). As can be shown in Figure 19, the present results are signifcant in agreement with the results in [39]. Figure 19 also clarifes the considerable diferences between linear and the nonlinear results. Another comparison is made between the present results, the results of Chan [40], and the results of Mondkar and Powell [39]. In this comparison, the time step is ∆t � 100μs and the time duration is 10000 μs (0.01 s). Tese results are compared with the so-called exact solution in [39], which is obtained by using a shorter time step of ∆t � 25 μs. Figure 20 reveals that the present results are closer to the so-called exact solution than the other results.

Damped Cantilever Beam Subjected to a Ramp-Ramp
Load at the Free End. In this example, a cantilever beam subjected to a ramp-ramp dynamic load is presented, as shown in Figures 21 and 22. Te geometric data of the beam are L � 120 in (0.508 m), A � 21.9 in 2 (0.014 m 2 ), and I Y � I Z � 100 in 4 (4.16 × 10 −5 m 4 ). Te material density ρ is 4.567 × 10 −3 lb.s 2 /in 4 (4.8808 × 10 4 kg/m 3 ), the modulus of elasticity E � 30 × 10 6 psi (207 GPa), and Poisson's ratio v � 0.3. In this example, a viscous damping coefcient equal to 10 lb/in/s (1.7513 × 10 3 N/m/s) is applied for each translational degree of freedom. Tis problem is classifed as a large rotation and large displacement problem [6]. Liu [6] solved this problem using a lumped mass matrix using eight elements. However, he did not clarify the time step used in the analysis. Behdinan et al. [41] also analysed this problem using two diferent formulations, the consistent corotational formulation and the updated Lagrangian formulation. Tey used a diferent time step in each formulation. However, they did not specify clearly the time step used in each formulation. Te present results are obtained using eight elements, with the error tolerance e r � 10 − 2 and time step ∆t � 1 x 10 − 3 . Tese results are compared with the results in [6,41]. Tis comparison clearly shows that the results are well consistent, as shown in Figure 23.

Cantilever Beam Subjected to a Sinusoidal-Concentrated End Load and a Concentrated End Moment.
A cantilever beam of length L = 10 m with uniform cross section, as can be seen in Figure 24, is solved in this section. Te beam is subjected to a vertical concentrated sinusoidal out of plane dynamic force F Z (t) and a bending moment M Z (t) at the free end, as shown in Figure 25. Te circular frequency of the force is ω = 50 rad/s, Young's modulus E is 210 GPa, material's density is ρ = 7,850 kg/m 3 , and Poisson's ratio is v = 0.3. Te dimensions of the cross section are e = 0.25 m and a = 0.3 m. Ten beam elements are used in the analysis, time step is ∆t � 5 x 10 − 4 s, and the error tolerance is e r � 10 − 2 . Le et al. [26] solved this problem using two diferent beam elements (cubic and linear) based on spatial spin variables and spatial rotational vector as a parametrization method for fnite rotations. Tey

18
Shock and Vibration Liu [6] Teoritical [6] Present study obtained a reference solution with 20 beam elements. Tey also compared their results with the results of Simo and Vu-Quoc [42]. Te results of the proposed formulation are closer to the so-called reference solution than the results of other formulations given in [26], as shown in Figures 26-28.

A Right Angle Cantilever Beam Subjected to a Sinusoidal-Concentrated End Load.
A right-angle cantilever beam is solved, as shown in Figure 29. Tis problem is classifed as a three-dimensional large deformation problem [6]. Te geometric properties of the beam are A � 21.9in 2 (0.014m 2 ) and I Y � I Z � 100in 4 (4.16 X 10 − 5 m 4 ) and material properties are E � 30 X 10 6 psi (207 GPa) and mass density ρ � 4.567 X 10 − 3 lb.s 2 /in 3 (126.4 kg.s 2 /m 3 ). Te beam is subjected to a concentrated vertical sinusoidal dynamic force F Y (t) � 500000 sin(50t)lb (2224 sin(50t) kN). Liu [6] solved this problem using eight elements with lumped mass. He compared his results with the results of ANSYS program using the same number of elements. Eight beam elements are used in the proposed analysis, time step is ∆t � 3 x 10 − 4 s, and the error tolerance is e r � 10 − 2 . Te present results are compared with the results in [6], as can be seen in Figure 30. Tis comparison shows that the present results are in agreement with the results in [6].         Shock and Vibration 23

Conclusion
A modifed corotational fnite element formulation for geometrically nonlinear dynamic analysis of spatial beams and frames has been presented in this paper. Te following points are the main outcomes of this study: (i) Owing to the nature of the used corotational frame, which continuously updates and translates with each element, the used formulation condenses the computed local displacement vector compared to other formulations. In addition, it can smoothly separate the rigid-body rotations from the deformational rotations. (ii) Since the deformation is always small with respect to the corotational frame, the small strain theory has been applied. Accordingly, the stifness and mass matrices are symmetric which signifcantly reduces the required storage memory and consequently reduces the computational time and improves the convergence rate. (iii) A two-stage procedure is proposed to transform vectors and matrices from the fxed global frame to the moving corotational frame. Tis procedure depends essentially on Tait-Bryan angles, which are computed successively. Trigonometric rules for all rotation angles with their diferent cases have been defned in terms of kinematics parameters. Accordingly, the proposed method deals with the problems of the vertical members and the rotation angles greater than π/2. Tis contribution has been used for handling some of the reported cases which are classifed as singularity problems. (iv) For the numerical solution, an incremental iterative method based on the full Newton-Raphson method and the Newmark direct integration implicit method has been used to solve the equations of motion. A MATLAB code has been written for this purpose. Tis code includes detection functions, which successfully control the sign of rotations during the analysis. Terefore, no convergence problems have been detected throughout the study.  Liu [6] ANSYS [6] Present study Figure 30: Te dynamic response of the right-angle frame. 24 Shock and Vibration the adjustment of the time step and the element size in order to adjust the relative chordal rotations during each time step.

A.1. The Shape Function Matrix
where ( ) ′ is the frst derivatives with respect to x i .

Shock and Vibration 25
where E is the modulus of elasticity, G is the modulus of rigidity, a i is the cross-sectional area, I y i and I Zi are the moment of inertia about y ∧ i and z ∧ i axes, and J i is the polar moment of inertia.  where fx 2 is the axial force of the second node in x i direction, which is presented in equation (11).

Data Availability
Te data used to support the fndings of this study are available from the corresponding author upon request.

Conflicts of Interest
Te authors declare that they have no conficts of interest.