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

Elerian, Ahmed A. H.; Shebl, Saiid A.; Elkaranshawy, Hesham A.

doi:https://doi.org/10.1155/2023/6675678

Shock and Vibration

On this page

Abstract Introduction Conclusion Appendix Data Availability Conflicts of Interest References Copyright Related Articles

Research Article | Open Access

Volume 2023 | Article ID 6675678 | https://doi.org/10.1155/2023/6675678

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

Ahmed A. H. Elerian,¹Saiid A. Shebl,¹and Hesham A. Elkaranshawy¹

Academic Editor: Salvatore Caddemi

Received05 Apr 2023

Revised20 Jul 2023

Accepted02 Aug 2023

Published02 Sept 2023

Abstract

The aim of this study is to compose a corotational finite 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 filled the gaps. Three-coordinate systems are used, i.e., the fixed global coordinate system, the fixed 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. The deformation is always small relative to the corotational frame. The successive rotations between different coordinate systems are expressed using Tait–Bryan angles. Lagrange’s equation is used to derive the equation of motion, and the stiffness 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 fixed local frame approaches , and if any is greater than , 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 specifies the location of elements in all situations. Different examples of beams and frames are analysed. While the method is not intricate, it is timesaving, is highly effective, provides more stable and robust analysis, and gives sufficiently accurate results. Compared to the parametrization of the finite rotations technique, the method has a significant reduction in the convergence rate because it avoids the storage of joint orientation matrices.

1. Introduction

In recent years, the development of lightweight high-strength materials has attracted many industries [1–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. That 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 significant effect of the geometric nonlinearity on the flexible offshore structures and devices that cannot be ignored. Trapper [5] studied a pipe lay on a seafloor that experienced large deformations. Accordingly, he used a geometrically nonlinear model to calculate the maximum internal forces locations. Liu [6] investigated a geometrically nonlinear finite 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 flexible wind turbine blade. This 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. Therefore, Liu and Bai [8] considered the geometric nonlinearity during their experimental and numerical investigation of a deployable composite cabin for space habitats. The important question that can be addressed here is why shall we develop geometrically nonlinear models in the existing of various finite element commercial softwares? The answer is not complicated; there is still a need to investigate more versatile, efficient, and time-saving models that can successfully handle different engineering problems which may be exclusive or intractable in some cases using these commercial programs.

According to Crisfield [9], Turner et al. [10] were the first to study geometrically nonlinear finite 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 finite element analysis of structures was covered in some notable books by Oden [12], Bathe [13], and Crisfield [9].

In general, there are two different forms, in continuum mechanics, to describe the motion of a body, namely, Eulerian and Lagrangian formulations. The Eulerian formulations are usually used in fluid mechanics, while Lagrangian formulations are used in most other engineering fields. 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–17], the system equation terms are defined in terms of the fixed global frame that does not change through the analysis. This 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 defined relative to a frame that is updated with the last accepted solution. This 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 configuration to the last equilibrium configuration is large, a basic assumption is violated and accordingly these formulations experience also some complexities. To avoid such complexities, corotational formulations [20–31] provide a simple kinematics description method for large displacement analysis. These formulations are based on the theory of small strain. The 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 filter the rigid-body motion from the deformational motion, which is always small relative to the corotational local frame. This frame, which is continuously updated, can be defined by many different methods [23]. Remseth [17] used an approximate vectorial assumption to deal with three-dimensional rotations. Therefore, this method is not applicable for large rotations. Thus, 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 specifically, 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. The element nodal rotations are determined from the angles between this set of orthogonal axes and the member axes. This procedure is improved by Crisfield [25], Le et al. [26], Jonker and Meijaard [27], and Hsiao et al. [28]. Though 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 finite rotations. It also increases the computational time because of storing the joint orientation matrices and parametrization of the finite rotations. Bathe and Bolourchi [29] defined 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 . Nevertheless, they did not specify clearly how to determine the transformation matrices in the case of vertical members, when the rotation angle is equal or − 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] identified the problem of singularity in case of using this method. That is why many authors used parametrization of finite rotations, such as the authors in [16, 25–27].

To express the stiffness 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 simplified 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 effect 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 finite element formulation has been developed. The stiffness and mass matrices are evaluated using the Euler–Bernoulli beam model. The transformation procedure is based on the Tait–Bryan angles successive rotations [23, 29]. This procedure is employed in two main stages to transform vectors and matrices from the fixed global frame to the moving corotational 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. The transformation procedure also depends upon updating the coordinates with every equilibrium configuration 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. These rules are used to calculate the rotations matrices. Meanwhile, this method is simple and does not need parametrization for finite rotations, but it requires regulating the time step and number of elements in order to decrease the relative chordal rotations during the analysis. The equation of motion is formulated using Lagrange’s equation. An incremental iterative procedure is used to solve the equation of motion. This procedure is based on the full Newton–Raphson method and the Newmark direct-time integration method. The MATLAB code is developed for this purpose. This code involves a relatively rapid convergence rate for equilibrium because it avoids storing joint orientation matrices and parametrizing of finite rotations, which are often associated with parametrized formulations. This section represents a general introduction and reveals the importance of studying geometric nonlinearity. The spatial beam element motion description method, which involves the coordinate systems and the displacement interpolation, is presented in Section 2. The transformation procedure between different coordinate systems based on Tait–Bryan angles can be seen in Section 3. The strain energy and stiffness matrix are derived in Section 4. The kinetic energy and the mass matrix are presented in Section 5. Then, Lagrange’s equation is used to derive the equation of motion in Section 6. Section 7 represents the solution strategy. To expose the efficiency and validate the accuracy of the proposed model, five numerical examples are solved and compared with the published results in Section 8. Finally, the conclusions are presented in Section 9.

2. Kinematics Description of the 3D Beam Element

2.1. Basic Assumptions

The following assumptions are employed to formulate the spatial beam element:(i)The material is isotropic, elastic, and homogeneous(ii)The 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 effect 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.

2.2. Coordinate Systems

After discretization of the structure into finite elements, the i^th beam element can be defined with two end nodes (n = 1 and 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 at time = 0 (_i, _i, _i), and the moving local coordinate system associated with the corotational local frame (_i, _i, _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 dynamics: an initial configuration at , the j^th equilibrium configuration at , and a current configuration at . The element’s initial length is L_o, and after deformation in the current configuration, the element length is equal to the arc length , while L_c is the current chord length.

For the current configuration, as shown in Figure 2, the nodal displacement vector for the i^th beam element in the fixed 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 (n = 1, 2) are the counterclockwise rotations about X, Y, andZ axes, respectively.

The nodal incremental displacement vector of the i^th beam element in the global coordinate system is defined as follows: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 (n = 1, 2) are the counterclockwise incremental rotations about X, Y, and Z axes, respectively. The nodal displacement vector can be updated by the following equation:where 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 is divided into the nodal translational displacement vector and the nodal rotational displacement vector , which can be written as follows:

Similarly, the vector can be divided into the nodal incremental translational displacement vector and the nodal incremental rotational displacement vector as follows:

Thus, the nodal coordinate vector of the i^th beam element in the fixed global system can be updated continuously as follows:where is the vector of the nodal coordinates of the i^th beam element relative to the fixed global frame, at the current configuration, and 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 as follows:where , , and are the displacement translational components in _i, _i, and _i directions, respectively, and (n = 1 and 2) are the counterclockwise deformational rotations after eliminating the rigid-body rotations.

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

The internal elastic force vector of the i^th beam element in the element corotational local coordinate system, at the current configuration, is as follows:where the internal elastic force vector components in both the fixed global coordinate system and the corotational local coordinate system are shown in Figure 2 with their positive signs.

2.3. Displacement Interpolation

The classical Hermitian shape functions are used to relate the element axial elongation , lateral deflections of the centroid axis shown in Figure 3, and rotation about the centroid axis to the element nodal displacement vector as follows:

The components of the shape functions of the i^th beam element are given by the following equation:where ξ is given bywhere is the element current chord length.

The 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 and are equal to zero. Furthermore, the axial displacement of the first node is chosen to be zero while the axial displacement of the second node is . Consequently, the nodal displacement vector has only seven nonzero components which will simplify the analysis as can be seen in the next section. This vector can be written as follows:

Thus, equations (12) and (13) can be simplified to

The arc length of the i^th beam element S_i can be expressed bywhere and are the first derivatives of the functions and with respect to . This integral is evaluated using the Gaussian integration scheme. The axial elongation of the i^th beam element can be defined as follows:

This equation can also be written in terms of chord lengths as follows:where is the element initial length and is the element axial elongation due to the bowing effect, which can be determined in terms of rotations [35] as follows:

Hence, the axial displacement in equation (17) is given by

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 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 (_i, _i, _i) can be expressed by the following equation: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:

Similarly, the relation between the fixed local frame (_i, _i, _i) and the current corotational local frame (_i, _i, _i) iswhere 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:

One can write the vector in terms of V_d as follows:where

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

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

The first stage is to transform from the fixed global system to the fixed local system using the three successive rotations , , , as shown in Figures 5–8, as follows:(1)Rotation of the (X, Y, Z) coordinate axes about the Y axis: This rotation places the X and Z axis along and , respectively, while the Y axis remains the same, as shown in Figure 5. Using the direction cosines of frame (, , ) with respect to the global frame (X, Y, Z), shown in Figures 5 and 8, the rotation matrix can be determined as follows: where(2)Rotation of the (, , ) coordinate axes about the axis: This rotation places the and axis along and , respectively, while the axis remains the same, as shown in Figure 6. One can use the direction cosines of the frame (, , ) with respect to the frame (, , ), and the rotation matrix can be obtained as follows: where(3)Rotation of the (, , ) coordinate axes about axis: This rotation places the and axis along _i and _i, respectively, while the axis remains the same. It has to be noticed that _i coincides with the axis , as shown in Figure 7.

By assuming a point P lies in the (_i, _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 fixed global system coordinate (X, Y, Z) can be written as follows:

Consequently, the coordinates of point P with respect to the (, , ) frame, as shown in Figure 9, can be obtained as follows:which leads to

Now, using the direction cosines of frame (, , ), which coincides with the (_i, _i, _i) frame, with respect to previous frame (, , ), the rotation matrix can be determined as follows:where

Hence, the rotation matrix 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 is insignificant, in the case of a circular cross section element. Thus, the rotation matrix can be calculated as follows:

There is another special case where the initial position of the element is vertical, in the Y-axis direction. In order to get , there are only two successive rotations not three as in the general case. The first rotation is either 90° or 270°, as shown in Figure 10, depending on whether C_Y is +1 or −1. The other rotation is , which is shown in Figure 11, and can be determined using a reference the point P that lies in the (_i, _i) plane. In this case, equations (39) and (40) are modified as follows:

Thus, the matrix in equation (43) can be written as follows:

Substituting equation (42) or equation (45) into equation (25), the matrix can be determined.

3.2. Transformation from the Fixed Local Frame to the Moving Corotational Local Frame

The second stage is to transform from the fixed local system (_i, _i, _i) to the moving corotational local system (_i, _i, _i), at the current configuration, using the three successive rotations , , and , as shown in Figure 12. These rotations are similar to the rotations , , and in Figures 5–7; however, the calculation method depends upon the relative displacements.(1)Using the direction cosines of the frame (, , ) with respect to fixed local frame (_i, _i, _i), which is shown in Figure 12, the rotation matrix can be determined as follows: where Of particular interest is to notice that in [29], a difficulty arises when angle 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 12, , , and are the i^th beam element relative translational displacements with respect to the fixed local frame. These relative displacements of the i^th beam element can be determined from the relative displacement with respect to the fixed global system calculated from the 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 frame (, , ) with respect to frame (, , ), and the rotation matrix can be obtained as follows: When employing the technique outlined in [29], a complication rises when the angle , 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. These two trigonometric relationships can be expressed as follows: where and SN is equal +1 if and −1 if . The 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 vanishes, and the rotation is either 90° or 270°, depending on the element position. Traditionally, this case could create singularity and it had been the source of many difficulties [31], and the authors of this research work suggested preventing the rotation 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 is either 90° or 270°. The matrix can be rewritten as follows: where can be specified using the current vector of the nodal coordinates in equation (7) as follows: Hence, the value of is either +1 and the rotation is 90° or −1 and the rotation is 270°.(3)When the third rotation angle is computed by assuming a reference point P as has been done for , it is very hard to assume a point in the (_i, _i) plane every iteration with different consequence positions and directions. Accordingly, the model experienced some difficulties in converging using this method. Thus, an incremental procedure [29] is employed here to compute this angle as follows: where and are the incremental twist angles about the _i axis. These incremental angles can be determined using the following procedure. First, the incremental rotational displacement vector, with respect to the fixed global frame in equation (6), is transformed to the corresponding vector in the fixed local frame using the procedure indicated in equation (24) as follows:

Then, the incremental rotational vector relative to the fixed local frame is transformed to the corotational local frame as follows:

Therefore, in equation (60) can be determined using the following relation:

Consequently, the angle , at the current configuration, can be determined as follows:where is the twist rotation about the _i axis, at the j^th equilibrium configuration. This equation is used as a predictor between two successive time steps. Also, it is used as a corrector between two successive iterations.

Thus, the rotation matrix can be determined as follows:

The 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 in equation (27) has been specified. Then, the transformation matrix in equation (30) has been determined.

4. Strain Energy and the Stiffness Matrix

Considering isotropic elastic materials, the constitutive relation between the stress vector and the strain vector of the i^th beam element can be defined by the following equation:where is the symmetric matrix of the elastic coefficients. The strain vector is given by the following equation:where is the differential operator matrix and is the deformation vector, which can be defined using the shape functions in equation (15) as follows:where 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:

The strain energy for the i^th beam element is given by the following equation:where is the volume. Substituting equations (70) and (71) into equation (72), the strain energy can be expressed as follows:

The element local displacement vector is independent of the volume . Subsequently, equation (73) can be simplified to the following form:where k_i is the symmetric element stiffness matrix in the corotational local coordinate system, which can be defined as follows:and B is defined as follows:

Equation (74) can be written as follows: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 follows: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 A.2 and A.3.

5. The Kinetic Energy and the Mass Matrix

The kinetic energy for the i^th beam element is given by the following equation:where is the differentiation with respect to time t, is the density of the i^th beam element, and is the velocity of a general point in the i^th beam element with respect to the global coordinate system. The velocity of a general point in the i^th beam element with respect to the corotational local coordinate system can be expressed as follows: can be written in the following form:where is the nodal velocity vector of the i^th beam element with respect to the local coordinate system. The vector can be expressed as follows:where 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:

The velocity vector and the transformation matrix are independent of the volume . Subsequently, equation (83) can be simplified to the following form:where m_i is the symmetric element mass matrix in the local corotational coordinate system defined as follows:

Equation (84) can be written in the following form:where is the symmetric element mass matrix in the fixed global coordinate system defined as follows:

For the beam element used in this research work, the element mass matrix in the local coordinate system can expressed as follows:where is the ith beam element mass matrix for the translational inertia and is the mass matrix of the i^th beam element for the rotational inertia. The 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 three-coordinate 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. Therefore, the usual vector addition rules can be applied to these rotational displacements. Hence, the obtained stiffness and mass matrices are symmetric.

6. 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 is the kinetic energy, is the strain energy, is the work done by the external forces, is the nodal velocity vector in the global coordinate system, and is the nodal displacement vector in the global coordinate system. The work done by the external forces on the i^th beam element is given by the following equation:where is the vector of the external applied forces on the i^th beam element in the fixed global coordinate system.

Using equation (86), one can write the first 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 writewhere is the internal elastic force vector in the fixed global and 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 is the vector of the external applied forces of the entire structure and is the internal elastic force vector of the entire structure. Both and are defined in the fixed global coordinate system. is the mass matrix in the global coordinate system of the entire structure, and is the acceleration vector in the global coordinate system of the entire structure. Rayleigh damping is used in equation (98) to readwhere is the damping matrix in the global coordinate system of the entire structure and is the velocity vector in the global coordinate system of the entire structure. The damping matrix is defined in terms of the stiffness matrix and the mass matrix as follows:where are the damping coefficient which can be determined from the vibration modes of the system. is the assembled stiffness matrix in the global coordinate system of the entire structure.

7. Numerical Algorithms

The equation of motion is solved using an incremental iterative procedure. This 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. The iteration equilibrium convergence criterion is given by the following equation:where is the reference out of balance force, which is assumed to be the out of balance force in the first iteration, and 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 the 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 nonsymmetric matrices compared with the storage of only the triangular part in the symmetric matrices. Therefore, we are confident that the restriction we made in subsection (2.1), 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 solution, it is assumed that the displacement, velocity, and acceleration vectors are null. The values of the vectors , , and , at a known equilibrium configuration and time , are , , and , respectively. Likewise, the values of , , and , at time , are , , and , where is the time step. Using the geometric data, one can search for vertical member, which is handled differently, as stated in Section 3. The numerical solution procedure is described in the following steps at the beginning of each time step:(i)Compute , , and for each element using equations (78), (88) and (96), respectively.(ii)In the first iteration, calculate the matrices and using equations (25) and (41). In all other iterations, calculate and using equations (27) and (66).(iii)Determine the transformation matrix for each member using equation (30).(iv)Obtain , , and according to equations (77), (87), and (95).(v)Get , , and for the entire structure, by assembling the stiffness matrices, the mass matrices, and the internal elastic force vectors for all elements of the analysed structure.(vi)Calculate the out of balance force from equation (101).(vii)If the convergence condition in equation (102) 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 is calculated as follows: where is the effective matrix that can be defined as follows: where and are the Newmark parameters [38].(b)Update the incremental displacement vector as follows: where is the incremental displacement vector in the previous iteration, which is considered to be zero in the first iteration.(c)Extract vector for each element from vector . Consequently, one can update the vectors and using equations (3) and (7).(d)Using and , the relative displacements are calculated, as in equations (49)–(51).(e)From the coordinate vector , the model can check for each element position to apply either the regular rotations matrices in equations (42), (57), and (65) or the vertical member rotations matrices in equations (45), (58), and (65). Then, and are updated for each element.(f)Using the relative displacements, the rigid-body rotations can be obtained as follows:(g)Eliminate the rigid-body rotations from the rotational components and in the vector , for each element, as follows:(h)Transform the pure rotations and , which are relative to the fixed global frame, to determine the corresponding rotational components . and for vector in the current corotational local frame which are always relative to the element chord, using the procedure detailed in equations (24) and (26).(i)The axial displacement in equation (17) is obtained from equation (23).(j)Using equations (17) and (96), and can be determined, respectively.(k)Using the Newmark direct integration method [38], and can be calculated as follows:(viii)Going to the start of step iv again.(ix)The displacement vector at time can be updated as follows:(x)Start the next time step.

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 given in equation (7). 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 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 in equation (59). In addition, when the relative displacement turns out to be zero, the angle vanishes. Likewise, if the relative displacement turns out to be zero, the angle vanishes; however, the element in general is not vertical. Hence, the program deals with these special cases separately.

This function also controls the angle when the rotation is outside the interval [, −. The sign SN in equation (53) specifies the cosine of the angle to meet the corresponding element position during the motion because the terms and are always positive which cannot reflect the real sign of cosine of . At the same time, the sine of is already specified with the sign of . This function conserves the code to converge efficiently. This numerical algorithm steps are organized in Figure 13.

8. Numerical Examples

8.1. Simple Supported Beam Subjected to a Concentrated Step Load

The first example is a simply supported beam whose geometric properties are L = 12 m, A = , and and material properties are E = 210 GPa, G = 80.775 GPa, and ρ = 7850 , as shown in Figure 14. The 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 s and the error tolerance is chosen to be . Liu [6] solved this problem using a corotational formulation based on Euler’s theorem. He compared his results with the theoretical results [6]. The present results are compared with the results in [6], as shown in Figure 16. This comparison reveals that the present results remarkably agreed with the theoretical results in [6].

8.2. Clamped-Clamped Beam Subjected to a Concentrated Vertical Load at the Midspan

In this part, a clamped-clamped beam is analysed. The geometric data of the beam are shown in Figure 17. The beams modulus of elasticity equals 30,000 ksi, density is 0.098 lb/in³, and Poisson’s ratio is considered to be zero. This beam is subjected to a dynamic concentrated step load. This load is assumed to be applied suddenly to the midspan, as can be shown in Figure 18. This problem is solved by Mondkar and Powell [39] using various time steps. A secant stiffness concept is utilized to solve this problem by Chan [40]. The present results are obtained using ten element and error tolerance . The present results are compared with the results of [39] and the linear analysis results, using the same time step . The dynamic load is applied over a period of 5000 s (0.005 s). As can be shown in Figure 19, the present results are significant in agreement with the results in [39]. Figure 19 also clarifies the considerable differences 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 and the time duration is 10000 s (0.01 s). These results are compared with the so-called exact solution in [39], which is obtained by using a shorter time step of . Figure 20 reveals that the present results are closer to the so-called exact solution than the other results.

8.3. 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. The geometric data of the beam are L = 120 in (0.508 m), A = 21.9 in² (0.014 m²), and I_Y = I_Z = 100 in⁴ (4.16 × 10⁻⁵ m⁴). The material density ρ is 4.567 × 10⁻³ lb.s²/in⁴ (4.8808 × 10⁴ ), the modulus of elasticity E = 30 × 10⁶ psi (207 GPa), and Poisson’s ratio = 0.3. In this example, a viscous damping coefficient equal to 10 lb/in/s (1.7513 × 10³ N/m/s) is applied for each translational degree of freedom. This problem is classified 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 different formulations, the consistent corotational formulation and the updated Lagrangian formulation. They used a different time step in each formulation. However, they did not specify clearly the time step used in each formulation. The present results are obtained using eight elements, with the error tolerance and time step . These results are compared with the results in [6, 41]. This comparison clearly shows that the results are well consistent, as shown in Figure 23.

8.4. 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. The 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. The circular frequency of the force is = 50 rad/s, Young’s modulus E is 210 GPa, material’s density is ρ = 7,850 kg/, and Poisson’s ratio is = 0.3. The 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 s, and the error tolerance is . Le et al. [26] solved this problem using two different beam elements (cubic and linear) based on spatial spin variables and spatial rotational vector as a parametrization method for finite rotations. They obtained a reference solution with 20 beam elements. They also compared their results with the results of Simo and Vu-Quoc [42]. The 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.

8.5. A Right Angle Cantilever Beam Subjected to a Sinusoidal-Concentrated End Load

A right-angle cantilever beam is solved, as shown in Figure 29. This problem is classified as a three-dimensional large deformation problem [6]. The geometric properties of the beam are A = and and material properties are E = psi (207 GPa) and mass density ρ = (126.4 ). The beam is subjected to a concentrated vertical sinusoidal dynamic force F_Y (t) = . 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 s, and the error tolerance is . The present results are compared with the results in [6], as can be seen in Figure 30. This comparison shows that the present results are in agreement with the results in [6].

9. Conclusion

A modified corotational finite element formulation for geometrically nonlinear dynamic analysis of spatial beams and frames has been presented in this paper. The 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 stiffness and mass matrices are symmetric which significantly 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 fixed global frame to the moving corotational frame. This procedure depends essentially on Tait–Bryan angles, which are computed successively. Trigonometric rules for all rotation angles with their different cases have been defined in terms of kinematics parameters. Accordingly, the proposed method deals with the problems of the vertical members and the rotation angles greater than π/2. This contribution has been used for handling some of the reported cases which are classified 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. This code includes detection functions, which successfully control the sign of rotations during the analysis. Therefore, no convergence problems have been detected throughout the study.(v)Each iteration, the element coordinate, stiffness, and mass matrices are updated regularly. Though this updating requires time in each iteration, it significantly decreases the overall time of the analysis.(vi)The results have been compared with the published results to show the effectiveness and accuracy of the proposed formulation and the numerical algorithm. Though these problems have been solved using an adequate number of elements, accurate results have been obtained compared with the published results.(vii)The proposed formulation provides a rapid convergence rate and does not need special parametrization for finite rotations. However, it requires the adjustment of the time step and the element size in order to adjust the relative chordal rotations during each time step.

Appendix

A.1. The Shape Function Matrix

where is the first derivatives with respect to .

A.2. The Axial and Bending Stiffness Matrix

where is the modulus of elasticity, G is the modulus of rigidity, is the cross-sectional area, and are the moment of inertia about _i and _i axes, and is the polar moment of inertia.

A.3. The Geometric Stiffness Matrix

where is the axial force of the second node in _i direction, which is presented in equation (11).

A.4. The Mass Matrix for the Translational Inertia

A.5. The Mass Matrix for the Rotational Inertia

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

A. K. Naik, M. Nazeer, D. Prasad, T. Laha, and S. Roy, “Development of functionally graded ZrB2–B4C composites for lightweight ultrahigh-temperature aerospace applications,” Ceramics International, vol. 48, no. 22, pp. 33332–33339, 2022.
View at: Publisher Site | Google Scholar
Z. Zhang, L. Zhang, B. Song, Y. Yao, and Y. Shi, “Bamboo-inspired, simulation-guided design and 3D printing of light-weight and high-strength mechanical metamaterials,” Applied Materials Today, vol. 26, Article ID 101268, 2022.
View at: Publisher Site | Google Scholar
M. Shadlou, W. Sun, and T. Lou, “A novel FRP composite with high-strength, high extensibility in tension: 1-D constitutive relation,” Composite Structures, vol. 261, Article ID 113332, 2021.
View at: Publisher Site | Google Scholar
J. Leng, Q. Wang, and Y. Li, “A geometrically nonlinear analysis method for offshore renewable energy systems—examples of offshore wind and wave devices,” Ocean Engineering, vol. 250, Article ID 110930, 2022.
View at: Publisher Site | Google Scholar
P. A. Trapper, “A numerical model for geometrically nonlinear analysis of a pipe-lay on a rough seafloor,” Ocean Engineering, vol. 252, Article ID 111146, 2022.
View at: Publisher Site | Google Scholar
Y. Liu, ““The development of the co-rotational finite element for the prediction of the longitudinal load factor for a transmission line system,”,” The University of Manitoba, Winnipeg, Manitoba, 2014, PhD Thesis.
View at: Google Scholar
Q. Xiaohang, G. Zhiteng, W. Tongguang, W. Long, and K. Shitang, “Nonlinear aeroelastic response analysis of 100-meter-scale flexible wind turbine blades,” Acta Aerodynamica Sinica, vol. 40, pp. 220–230, 2022.
View at: Google Scholar
T. W. Liu and J. B. Bai, “Folding behaviour of a deployable composite cabin for space habitats-part 1: experimental and numerical investigation,” Composite Structures, vol. 302, Article ID 116244, 2022.
View at: Publisher Site | Google Scholar
M. A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, Jon Wiley & Sons, West Sussex, UK, 1st edition, 1991.
M. J. Turner, R. W. Clough, H. C. Martin, and L. J. Topp, “Large deflections of structures subjected to heating and external loads,” Journal of the Aerospace Sciences, vol. 27, pp. 99–106, 1960.
View at: Google Scholar
J. H. Argyris, Recent Advances in Matrix Methods of Structural Analysis, vol. 4, Pergamon Press, New York, NY, USA, 1964.
J. T. Oden, Finite Elements of Nonlinear Continua, McGraw-Hill, New York, NY, USA, 1972.
K. J. Bathe, Finite Element Procedures, Prentice Hall, New Jersey, NJ, USA, 1st edition, 1996.
M. V. B. Santana, C. Sansour, M. Hjiaj, and H. Somja, “An equilibrium‐based formulation with nonlinear configuration dependent interpolation for geometrically exact 3D beams,” International Journal for Numerical Methods in Engineering, vol. 123, no. 2, pp. 444–464, 2021.
View at: Publisher Site | Google Scholar
D. Vo, P. Nanakorn, and T. Q. Bui, “A total Lagrangian Timoshenko beam formulation for geometrically nonlinear isogeometric analysis of spatial beam structures,” Acta Mechanica, vol. 231, no. 9, pp. 3673–3701, 2020.
View at: Publisher Site | Google Scholar
J. C. Simo and L. Vu-Quoc, “A three-dimensional finite-strain rod model. Part II: computational aspects,” Computer Methods in Applied Mechanics and Engineering, vol. 58, no. 1, pp. 79–116, 1986.
View at: Publisher Site | Google Scholar
S. N. Remseth, “Nonlinear static and dynamic analysis of framed structures,” Computers and Structures, vol. 10, no. 6, pp. 879–897, 1979.
View at: Publisher Site | Google Scholar
Y. B. Yang, S. Kuo, and Y. S. Wu, “Incrementally small-deformation theory for nonlinear analysis of structural frames,” Engineering Structures, vol. 24, no. 6, pp. 783–798, 2002.
View at: Publisher Site | Google Scholar
Y. B. Yang, S. P. Lin, and C. S. Chen, “Rigid body concept for geometric nonlinear analysis of 3D frames, plates and shells based on the updated Lagrangian formulation,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 7, pp. 1178–1192, 2007.
View at: Publisher Site | Google Scholar
T. Belytschko and B. Hsieh, “Non-linear transient finite element analysis with convected co‐ordinates,” International Journal for Numerical Methods in Engineering, vol. 7, no. 3, pp. 255–271, 1973.
View at: Publisher Site | Google Scholar
H. A. Elkaranshawy and M. A. Dokainish, “Corotational finite element analysis of planar flexible multibody systems,” Computers and Structures, vol. 54, no. 5, pp. 881–890, 1995.
View at: Publisher Site | Google Scholar
H. A. Elkaranshawy, A. A. H. Elerian, and W. I. Hussien, “A corotational formulation based on Hamilton’s principle for geometrically nonlinear thin and thick planar beams and frames,” Mathematical Problems in Engineering, vol. 2018, Article ID 2670462, 22 pages, 2018.
View at: Publisher Site | Google Scholar
J. Gu, ““Large displacement elastic analysis of space frames allowing for flexural-torsional buckling of beams,”,” Hong Kong Polytechnic University, Hong Kong, 2004, PhD Thesis.
View at: Google Scholar
C. Oran, “Tangent stiffness in space frames,” Journal of the Structural Division, vol. 99, no. 6, pp. 987–1001, 1973.
View at: Publisher Site | Google Scholar
M. A. Crisfield, “A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements,” Computer Methods in Applied Mechanics and Engineering, vol. 81, no. 2, pp. 131–150, 1990.
View at: Publisher Site | Google Scholar
T. N. Le, J. M. Battini, and M. Hjiaj, “A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures,” Computer Methods in Applied Mechanics and Engineering, vol. 269, pp. 538–565, 2014.
View at: Publisher Site | Google Scholar
J. B. Jonker and J. P. Meijaard, “A geometrically non-linear formulation of a three-dimensional beam element for solving large deflection multibody system problems,” International Journal of Non-Linear Mechanics, vol. 53, pp. 63–74, 2013.
View at: Publisher Site | Google Scholar
K. M. Hsiao, J. Y. Lin, and W. Y. Lin, “A consistent co-rotational finite element formulation for geometrically nonlinear dynamic analysis of 3-D beams,” Computer Methods in Applied Mechanics and Engineering, vol. 169, no. 1-2, pp. 1–18, 1999.
View at: Publisher Site | Google Scholar
K. J. Bathe and S. Bolourchi, “Large displacement analysis of three‐dimensional beam structures,” International Journal for Numerical Methods in Engineering, vol. 14, no. 7, pp. 961–986, 1979.
View at: Publisher Site | Google Scholar
A. C. Benjamin, “Análise Não-Linear Geométrica De Pórticos Tridimensionais Pelo Método Dos Elementos Finitos,” Rio de Janeiro, Brasil, 1982, M.Sc. thesis.
View at: Google Scholar
C. C. Nunes, H. L. Soriano, and F. filho, “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, vol. 38, no. 8, pp. 1195–1204, 2003.
View at: Publisher Site | Google Scholar
S. Caddemi and A. Morassi, “Multi-cracked Euler–Bernoulli beams: mathematical modeling and exact solutions,” International Journal of Solids and Structures, vol. 50, no. 6, pp. 944–956, 2013.
View at: Publisher Site | Google Scholar
S. Caddemi, I. Calio, and M. Marletta, “The non-linear dynamic response of the Euler–Bernoulli beam with an arbitrary number of switching cracks,” International Journal of Non-linear Mechanics, vol. 45, no. 7, pp. 714–726, 2010.
View at: Publisher Site | Google Scholar
A. Beheshti, “Large deformation analysis of strain-gradient elastic beams,” Computers and Structures, vol. 177, pp. 162–175, 2016.
View at: Publisher Site | Google Scholar
S. L. Chan, “Large deflection kinematic formulations for three-dimensional framed structures,” Computer Methods in Applied Mechanics and Engineering, vol. 95, no. 1, pp. 17–36, 1992.
View at: Publisher Site | Google Scholar
H. Goldstein, C. Poole, and J. Safko, Classical Mechanics, Addison-Wesley, Boston, MA, USA, 1980.
A. A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, Cambridge, UK, 2020.
N. M. Newmark, “A method of computation for structural dynamics,” Journal of the Engineering Mechanics Division, vol. 85, no. 3, pp. 67–94, 1959.
View at: Publisher Site | Google Scholar
D. P. Mondkar and G. H. Powell, “Finite element analysis of non‐linear static and dynamic response,” International Journal for Numerical Methods in Engineering, vol. 11, no. 3, pp. 499–520, 1977.
View at: Publisher Site | Google Scholar
S. L. Chan, “Large deflection dynamic analysis of space frames,” Computers and Structures, vol. 58, no. 2, pp. 381–387, 1996.
View at: Publisher Site | Google Scholar
K. Behdinan, M. C. Stylianou, and B. Tabarrok, “Co-rotational dynamic analysis of flexible beams,” Computer Methods in Applied Mechanics and Engineering, vol. 154, no. 3-4, pp. 151–161, 1998.
View at: Publisher Site | Google Scholar
J. C. Simo and L. Vu-Quoc, “On the dynamics in space of rods undergoing large motions—a geometrically exact approach,” Computer Methods in Applied Mechanics and Engineering, vol. 66, no. 2, pp. 125–161, 1988.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Ahmed A. H. Elerian et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

243

Downloads

185

Citations