A Corotational Formulation Based on Hamilton’s Principle for Geometrically Nonlinear Thin and Thick Planar Beams and Frames

Elkaranshawy, Hesham A.; Elerian, Ahmed A. H.; Hussien, Walied I.

doi:https://doi.org/10.1155/2018/2670462

Mathematical Problems in Engineering

On this page

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

Research Article | Open Access

Volume 2018 | Article ID 2670462 | https://doi.org/10.1155/2018/2670462

A Corotational Formulation Based on Hamilton’s Principle for Geometrically Nonlinear Thin and Thick Planar Beams and Frames

Hesham A. Elkaranshawy,¹Ahmed A. H. Elerian,¹and Walied I. Hussien¹

Academic Editor: Nerio Tullini

Received19 May 2018

Accepted11 Jul 2018

Published14 Aug 2018

Abstract

A corotational finite element formulation for two-dimensional beam elements with geometrically nonlinear behavior is presented. The formulation separates the rigid body motion from the pure deformation which is always small relative to the corotational element frame. The stiffness matrices and the mass matrices are evaluated using both Euler-Bernoulli and Timoshenko beam models to reveal the shear effect in thin and thick beams and frames. The nonlinear equilibrium equations are developed using Hamilton’s principle and are defined in the global coordinate system. A MATLAB code is developed for the numerical solution. In static analysis, the code employed an iterative method based on the full Newton-Raphson method without incremental loading, while, in dynamic analysis, the Newmark direct integration implicit method is also utilized. Several examples of flexible beams and frames with large displacements are presented. Not only is the method simple and time-saving, but it is also highly effective and highly accurate.

1. Introduction

Over the past few decades the geometrically nonlinear finite element analyses of flexible beams and frames with large displacements drew the attention of many authors [1–22]. Remarkably, it is still an open research subject. No wonder, new applications need to be treated like the analysis of flexible hoses in the offshore industry [23], highly flexible deployable structures which are used in architecture and civil engineering [24], functionally graded structures [25], and flexible aircraft [26, 27]. Additionally, in applications like large wind turbine blades, the structural geometric nonlinearity may affect the aerodynamic loading [28, 29]. Consequently, the module for structural geometric nonlinearity has to be coupled with plant-scale CFD simulations. Though commercial finite element programs like ANSYS and ABAQUS can correctly handle geometrically nonlinear problems, they have high computational cost that makes them intractable in such situations. Hence, there is a need for accurate and computationally inexpensive geometrically nonlinear models. Geometrically nonlinear analysis not only is used in immense applications in structural and mechanical engineering but also would apply in several fields. For instance, one can model a plant’s stem which grows vertically as a cantilever beam, so it can be studied under the affecting wind loads. Consequently, the appropriate size of plant stems of the wind related crops can be determined which is very vital for controlling plant stem failure and optimization process [30].

The key point in geometrically nonlinear analysis boils down to the description of the kinematics. In literature, there are three kinematics description formulations for a flexible beam element. They are the total Lagrangian formulation, the updated Lagrangian formulation, and the corotational formulation. The total Lagrangian formulation is the oldest kinematics description formulation and it is the most widely used one in finite element commercial programs such as ANSYS, ABAQUS, and MARC. In this formulation, the system equations are defined in terms of a fixed global coordinate system, which is not changed through analysis. Large displacements and strains are found in this formulation which needs special handling. Consequently, the system equations could be very complicated. This formulation is used by many authors like those of [1–3, 22]. In the updated Lagrangian formulation [4–6], the system equations are defined in terms of a coordinate system which is updated with the last accepted solution. This coordinate system is kept fixed during the solution cycles. Consequently, the system equations are much simpler than the corresponding equations in the total Lagrangian formulation. However, if the displacement from the current configuration to the last equilibrium configuration is large, a basic assumption is violated and consequently the results cannot be trusted. To avoid this problem, the corotational kinematical formulation was introduced. This formulation has been developed, improved, and applied by many researchers; see, for example, [7–17, 20, 23, 25, 31–34].

The corotational formulation is based on small strain beam theory. It has two coordinate systems: the fixed global coordinate system and the element corotational local coordinate system. The corotational local frame rotates and translates with each element but does not deform with it. Therefore, this formulation easily separates the rigid body motion from the pure deformation which is always small relative to the corotational element frame.

Two theories have widely been employed for studying beam deflections. They are Euler-Bernoulli beam theory and Timoshenko beam theory. In Euler-Bernoulli beam theory the shear effect is neglected, while in Timoshenko beam theory the shear effect is considered. Generally, according to their thickness, beams and frames can be classified into two categories: thin and thick. In thin beams and frames, the ratio of length to thickness is large. Consequently, the shear effect is slight in this case. Inversely, the shear effect analysis is very significant in thick beams and frames, which have small length to thickness ratio. The majority of researchers neglected the shear effect in the geometrically nonlinear analysis of beam elements, so they evaluated the stiffness and mass matrices using Euler-Bernoulli beam model. To model the shear effect, Timoshenko beam model was used in [1, 34–36] for the evaluation of the stiffness and mass matrices. However, no sufficient examples of beams with small length to thickness ratio have been solved.

In the present work, a fairly accurate and computationally inexpensive corotational formulation is developed. The materials are considered to be isotropic and elastic, and the element cross section is uniform. The stiffness matrices, the mass matrices, and the internal elastic force vector due to deformation are firstly computed in the local corotational coordinate system. Then, they are transformed to the fixed global coordinate system and assembled to construct the structure overall stiffness matrix, mass matrix, and overall internal elastic force vector. On the other hand, the acceleration vectors, the velocity vectors, the displacement corrector vectors, the incremental displacement vectors, and the nodal coordinates are defined in terms of the fixed global coordinate system. Consequently, the nonlinear equilibrium equations are derived using Hamilton’s principle in the fixed global coordinate system. The internal elastic force vector is evaluated using the arc length integral formula. This integral is computed using Gauss integration scheme with five Gauss points. The stiffness matrices are estimated using both Euler-Bernoulli beam model and Timoshenko beam. In static analysis, an iterative method based on the full Newton-Raphson method without incremental loading is used to obtain the solution of the nonlinear equilibrium equation, while, in dynamic analysis, an incremental-iterative method is used based on the full Newton-Raphson method and the Newmark direct integration implicit method. MATLAB codes are developed for these purposes. The beam element kinematics and the displacement interpolation are defined in Section 2. The strain energy and the stiffness matrix are formulated in Section 3. The kinetic energy and the mass matrix are expressed in Section 4. Hamilton’s principle is used to derive the nonlinear differential equation of the equilibrium of the system in Section 5. Section 6 covers the numerical algorithms. The effectiveness and the accuracy of the present method are demonstrated through working out six problems with large displacements in Section 7. The conclusions are presented in Section 8.

2. Kinematics of a Beam Element Using Corotational Finite Element Formulation

2.1. The Coordinate Systems

Two coordinate systems are used for this formulation: a global coordinate system associated with the fixed global frame (X, Z) and a local coordinate system associated with the corotational local frame (, ). This corotational frame is updated and attached to each beam element as shown in Figure 1. It translates and rotates with the beam element but does not deform with it; i.e., the corotational frame has a rigid body displacement with respect to the global frame. In the dynamic analysis, there are three configurations employed: an initial configuration, an equilibrium configuration, and a current configuration as shown in Figure 1, while, in static analysis, only two configurations are used: an initial configuration and a current configuration. After discretization of the structure into finite elements, the beam element can be defined with two end nodes (n=1, 2). Every node has three degrees of freedom and is defined in both coordinate systems.

The beam element nodal displacement vector in the fixed global coordinate system, at current configuration, is given bywhere (n=1, 2) and W_n (n=1, 2) are the displacement components in X and Z directions, respectively, and (n =1, 2) are the counterclockwise rotations.

In the kinematics of the dynamic analysis, as shown in Figure 1, the nodal incremental displacement vector of the beam element in the global coordinate system is defined aswhere (n=1, 2) and (n=1, 2) are the incremental displacement components in X and Z directions, respectively, and (n=1, 2) are the counterclockwise incremental rotations. Thus, the displacement components in (1) can be identified aswhere , , and are the displacement components at the equilibrium configuration.

For both the static and dynamic analysis, the nodal displacement vector of the beam element in the element local coordinate system, at current configuration, iswhere (n=1, 2) and (n=1, 2) are the displacement components in and directions, respectively, and (n=1, 2) are the counterclockwise rotations. The initial and current element chord lengths denoted as and L are shown in Figure 1. The initial and current locations of the corotational frame are specified by the two angles , that are measured counterclockwise from the fixed global x-axis direction to the local -axis. and are given by The relation between of the element iswhere β is the rigid body rotation and is obtained as The relation between the displacement vectors , and the rigid body displacement vector in the global coordinate system takes the formwhere is the beam element transformation matrix, expressed by

2.2. The Displacement Interpolation

The classical Hermitian shape functions are used to relate the element deformation vector to the element nodal displacement vector as follows:The components of the shape functions of the beam element are given bywhere ξ is given byThe transverse displacement of the beam element takes the formIn the corotational formulation the frame is updated each iteration. The -axis is always passing through the two end nodes of the element as shown in Figure 1. Therefore, the displacement components are zeros. Thus (19) can be written as One can get the approximate current arc length , as follows:where is the derivative of the function (x) with respect to . This integral is evaluated using Gauss integration scheme with five Gauss points. The elongation of the beam element in direction isThis elongation can be approximated asThe axial displacement of the first node is chosen to be zero while the axial displacement of the second node in (6) is defined as

3. The Strain Energy and Stiffness Matrices

Considering isotropic elastic materials, the relation between the stress vector and the strain vector for the beam element is given bywhere is the symmetric matrix of the elastic coefficients. One can write the strain vector as follows:where is the differential operator matrix. Substituting (15) into (26), the strain vector can be expressed as follows:Combining (25) and (27), the stress vector can be rewritten asThe strain energy for the beam element is given bywhere is the volume. Substituting (27) and (28) into (29), the strain energy can be expressed byThe element local displacement vector is independent of the volume . Subsequently, (30) can be simplified to the following form:where is the symmetric element stiffness matrix in the local coordinate system, defined as follows:and B is defined asSubstituting (13) into (31), one can write (31) in the following form:where is the nodal pure displacement vector of the beam element and is the symmetric element stiffness matrix in the fixed global coordinate system, defined as follows:For the beam element used in this research work, the element stiffness matrix in the local coordinate system can be expressed aswhere is the axial and bending stiffness matrix and is the geometric stiffness matrix for the beam element. For Euler-Bernoulli beam model and Timoshenko beam model, these stiffness matrices can be found, for example, in [16, 17, 34] and [34–36], respectively. For completeness, they are stated in Appendix A.

4. The Kinetic Energy and Mass Matrices

The kinetic energy for the beam element is given bywhere is the differentiation with respect to time t, is the density in the beam element, and is the velocity of a general point in the beam element with respect to the global coordinate system. For in (14), one can writewhere I is the identify matrix. Consequently, one can write (37) asThe velocity of a general point in the beam element with respect to the local coordinate system can be expressed as can be written in the formwhere is the nodal velocity vector of the beam element with respect to the local coordinate system. Furthermore, can be expressed aswhere is the nodal velocity vector of the beam element in the global coordinate system. Combining (39) − (42), the kinetic energy can be written as follows:The velocity vector and transformation matrix are independent of the volume . Subsequently, (43) can be simplified to the following form:where is the symmetric element mass matrix in the local corotational coordinate system, defined as follows:Equation (44) can be written in the formwhere 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 be expressed aswhere is the beam element mass matrix of the translational inertia and is the mass matrix of beam element for the rotational inertia. For Euler-Bernoulli beam model and Timoshenko beam model, these mass matrices can be attained, for instance, from [16, 17, 34] and [30, 34, 35], respectively. For completeness, they are given in Appendix B.

5. The Nonlinear Equilibrium Equations

Hamilton’s principle is used to derive the equilibrium equations of the system. Mathematically, Hamilton’s principle takes the formwhere is the Lagrangian functional. is defined aswhere is the kinetic energy, is the potential energy, and is the work done by the external forces. In the structural applications, the potential energy is the elastic strain energy which is given in (34). is defined bywhere is the vector of the external applied forces on the beam element in the fixed global coordinate system. The variation and the integration operators are interchangeable. Therefore, one can write (49) in the following form:Using (46), kinetic energy term in (52) can be written asUsing (34), the strain energy term in (52) can be expressed asUsing (51), one can write the work done term in (52) asSubstituting (53), (54), and (55) in (52) givesThe differentiation with respect to time and the variation are interchangeable. Therefore, the term in (56) can be written in the following form:Consequently, substituting (57) in the kinetic energy term in (56) and integrating the same term by parts, putting in mind that is constant during the time interval from to , one obtainswhere is the nodal acceleration vector of the beam element in the global coordinate system. The variation at and is zero. Thus, the first term on the right hand side in (60) vanishes, Equation (60) can be written asSubstituting (60) in (56), one obtainsFor small time step , one can approximate to ; hence (61) can be expressed asThe only way for the integration in (62) to be zero is the vanishing of the integrand itself. Therefore, (62) can be rewritten asAs mentioned before, the displacement variation is arbitrary; thenAccording to (13), (34), (35), and (64), 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 bySubstituting (65) in (64) givesEquation (67) is the equation of motion of the 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 formwhere 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. One can introduce Rayleigh damping in (68) 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 aswhere are the damping coefficients which can be obtained from vibration modes of the system. is the stiffness matrix in the global coordinate system of the entire structure.

For static analysis, the second term in (68) is identically zero. Thus, the nonlinear equilibrium equation of the beam element can be expressed asTherefore, the nonlinear equilibrium equation of the overall structure can be written as follows:

6. Numerical Algorithms

6.1. Static Analysis

The full Newton-Raphson iterative method without incremental loading [37] is employed to solve the nonlinear equilibrium equation of the overall structure. Equation (72) can be written aswhere is the out of balance force. The iteration equilibrium convergence criterion is determined bywhere is the Euclidean norm, is the reference out of balance force, and is the error tolerance. is considered to be the out of balance force obtained in the first iteration. The error tolerance is chosen according to the problem. A MATLAB code is developed for the numerical solution. At the beginning of the iteration procedure, null displacement vectors are assumed, and then the following pursues:(1)Compute and for each element using (36) and (66), respectively.(2)Determine and according to (35) and (65), respectively.(3)Assemble the stiffness matrices and the internal elastic force vectors to obtain and for the entire structure.(4)Obtain the out of balance force using (73).(5)Iterations stop, if the convergence criterion satisfies the condition in (74).(6)Otherwise, a displacement corrector is calculated using full Newton-Raphson method as(7) is added to the displacement vector of the entire structure, . Consequently, in (1) can be extracted from . Therefore, the current configuration can be determined using (9) − (12).(8)Consequently, in (6) can be specified. It has to be noticed that, for solution stability purposes, in (6) is computed from (24), (22), and (23).(9)Start new iteration.

The flowchart for the strategy is given in Appendix C.

6.2. Dynamic Analysis

An incremental-iterative method is employed to solve the nonlinear dynamic equilibrium equations. This method is based on the full Newton-Raphson method and the Newmark direct integration implicit method [37, 38]. Equation (69) can be written asAt the beginning of the solution, displacement, velocity, and acceleration vectors are assumed to be nulls. , , and are the values of , , and at a known equilibrium configuration, at time . Likewise, , , and are the values of , , and at time , where is the time step. A MATLAB code is developed for the numerical solution. At the beginning of each time step, the procedure is outlined as follows: (1)Compute , , and for each element using (36), (48), and (66), respectively.(2)Determine , , and according to (35), (47), and (65), respectively.(3)Assemble the elements stiffness and mass matrices and the internal elastic force vectors to obtain , , and for the entire structure.(4)Obtain the out of balance force using (76).(5)If the convergence criterion satisfies the condition in (74), stop the iteration and go to step . Otherwise, start the iteration:(a)a displacement corrector is calculated using full Newton-Raphson method as where is the effective matrix, and it can be defined as where and are the Newmark parameters.(b)Calculate the incremental displacement vector as where is the incremental displacement vector of the previous iteration. values are considered to be zeros in the first iteration.(c) can be extracted for each element from . Then, in (1) and consequently the current configuration can be specified using (3)−(5) and (9)−(12).(d)Using (6), can be determined. It has to be noticed that in (6) is equated to in (22) and (24).(e)Using the Newmark direct integration method, and can be expressed as It has to be noticed that the Newmark implicit method is unconditionally stable if γ ≥ 0.5 and ≥ (2γ + 1)² /16. In the present work, Newmark parameters are chosen to be 0.25 and 0.5.(6)Go to the beginning of step again.(7)The displacement vector of the overall structure at time can be determined as(8)Start the next time step.

The flowchart for the strategy is given in Appendix C.

7. Numerical Examples

In this section several problems are analyzed and compared with the published results to demonstrate the accuracy and the efficiency of the proposed formulation and algorithm. Both Euler-Bernoulli beam model and Timoshenko beam model are used to evaluate the stiffness and mass matrices and the results of both models are compared. If they are distinct a separate curve for each model appears in the graphs. Otherwise, only one curve appears to represent the present study.

7.1. Static Analysis

7.1.1. Cantilever Beam Subjected to a Concentrated Vertical End Load

The first numerical example is a cantilever beam subjected to a vertical concentrated point load P at the free end as shown in Figure 2. Five beam elements are used in the analysis, , and the error tolerance is chosen to be . Young’s modulus E is considered to be 200 GPa and Poisson’s ratio is = 0.3. Two cases are investigated; in the first case, the beam geometric data are L = 10 m, b = 0.5 m, and h = 0.25 m. Therefore, the length to thickness ratio is L/h =40. The results of Timoshenko and Euler-Bernoulli beam elements are identical in this case. The present results are compared with the results of the total Lagrangian formulation [2] and the elliptic integral numerical solution [18] in Figure 3 and Table 1. U/L and W/L are the normalized displacements and PL²/EI is the normalized load. Considering the semianalytical results of [18] as the reference results, the comparison shows that the results of the proposed formulation is more accurate than the results of the formulation in [2]. Figure 4 and Table 2 show the normalized displacement errors. Clearly, the present results give the smallest errors comparing to the total Lagrangian formulation results in [2].

In the second case, a thick beam of thickness h = 2.5 m is considered. Therefore, the length to thickness ratio is L/h = 4. The shear effect cannot be ignored in this case. Consequently, Timoshenko and Euler-Bernoulli beam elements have dissimilar results. Both results are compared with the results of the commercial finite element program ANSYS using 100 eight-noded solid elements, in Figure 5 and Table 3. Considering the ANSYS program results as the reference results, Timoshenko beam model gives more reasonable results than Euler-Bernoulli beam mode.

7.1.2. Pinned-Fixed Square Diamond Frame

A pinned-fixed diamond frame of side length L and flexural rigidity EI loaded in tension as shown in Figure 6 is considered in this example. Two elements per member are used in the analysis, , and the error tolerance is chosen to be . In the first case study, the beam geometric data are L = 14.1421 m, b = 0.5 m, and h = 0.1 m. Therefore, the length to thickness ratio is L/h = 141.4214. The results of Timoshenko and Euler-Bernoulli beam elements are identical in this case. The results of the present formulation are compared with the total Lagrangian formulation [2] and the elliptic integral numerical solutions [18] in Figure 7 and Table 4. U/L and W/L are the normalized displacements and PL²/EI is the normalized load. Good agreement can be noticed between the present results and the semianalytical results in [18]. Figure 8 and Table 5 show the normalized displacement errors.

In the second case study, a thick beam of thickness h = 2.8284 m is considered. Therefore, the length to thickness ratio is L/h = 5. The shear effect exists in this case. It can be noticed that both ends for each side of the diamond frame moves. Therefore, Timoshenko and Euler-Bernoulli beam elements have slight different results. Both results are compared with the results of the commercial finite element program ANSYS using 200 eight-noded solid elements as the reference results. Timoshenko beam model gives closed results to ANSYS results compared to Euler-Bernoulli beam model; see Figure 9 and Table 6.

7.1.3. Portal Frame Subjected to Concentrated Vertical and Horizontal Loads

The portal frame ABCD shown in Figure 10 is considered in this example. The frame consists of three equal members of 120 in. The beam cross section data are as follows: b = 0.6619 in. and h = 17.7809 in. Therefore, the length to thickness ratio is L/h = 6.7488. The present results are obtained using both Euler-Bernoulli beam model and Timoshenko beam model to demonstrate the shear effect.

This problem has been solved before in [1, 10, 14, 19]. The solution of James et al. [19] failed to converge approximately after for and after for as can be seen in Figures 11 and 12.

Beheshti [1] compared his results and the results in [19] for loading case n = 0.5. He did not give results for loads behind the load-limits in [19] as shown in Figure 12. Oran and Kassimali [10] and Hsiao and Hou [14] succeeded in breaking the load-limits of [19]. Reasonable agreement can be observed between their results (Figures 11 and 12). Nevertheless, for case, a significant divergence can be observed due to the presence of a limit point in [14]. The present results of the Euler-Bernoulli beam model are in consistency with all of these results especially with results in [14]. The iterative technique used not only produces a limit point for case but also for case. Four elements for each member are used, , and the used error tolerance is . As expected the Timoshenko model gives a higher deflection than Euler-Bernoulli model for the same load.

7.2. Dynamic Analysis

7.2.1. Cantilever Beam Subjected to a Sinusoidal Concentrated Vertical end Load

In this example, a cantilever beam is subjected to a vertical concentrated sinusoidal dynamic load at the free end as shown in Figure 13. Young’s modulus E is considered 210 GPa, the material’s density ρ = 7,850 kg/, and Poisson’s ratio = 0.3. and its frequency = 50 rad/s. Ten beam elements are used in the analysis, , time step s, and the error tolerance . Figures 14 and 15 show the present solutions using Timoshenko and Euler-Bernoulli beam elements.

Thanh-Nam et al. [20] solved this problem using various formulations and obtained a reference solution with 48 elements using the total Lagrangian formulation. The present results using Euler-Bernoulli model are compared with the results in [20], as shown in Figures 16 and 17. The comparison with the so-called reference solution in [20] shows that the results of the proposed formulation are more accurate than the results of the other formulations of [20].

7.2.2. Damped Finite Vibration of a Cantilever Beam Subjected to a Concentrated Vertical End Load

In this example, a cantilever beam was subjected to a point vertical dynamic load at the free end, as shown in Figure 18. This load is completely removed at t = 1.5, as can be seen in Figure 19; consequently the beam undergoes free vibration. The geometric and material properties data are , , , , and and is the damping coefficient. Only four beam elements are used in the analysis, , time step , and the error tolerance is chosen to be . The results of Timoshenko and Euler-Bernoulli beam elements are identical in this problem. The present results are compared with the results in [15, 22], as shown in Figure 20. This comparison shows the good agreement between the proposed formulation results and the solutions reported in [15, 22].

7.2.3. Fixed-Fixed Beam Subjected to a Concentrated Vertical Point Load at the Midspan

A fixed-fixed beam subjected to a vertical concentrated dynamic point load P at the midspan as shown in Figure 21 is considered in this example. Young’s modulus E is 30,000 ksi, the material’s density ρ is , and Poisson’s ratio is 0. Due to the symmetry of the problem, only one-half of the beam is considered and the horizontal displacement of the midpoint vanishes. To show the effect of shear, two case studies are considered. In the first case study, the beam geometric data are b = 1 in and h = 0.125 in is solved. Consequently, the length to thickness ratio for this thin beam is L/h = 160. The dynamic load P is equal to as shown in Figure 22(a). Ten beam elements are used in the analysis, , time step s, and the error tolerance . The results of Timoshenko and Euler-Bernoulli beam elements are identical in this case. A very good agreement can be noticed between the present results and the solutions reported in [15, 21], as shown in Figure 23.

(a)

(b)

In the second case study, the beam geometric data are b = 1 in and h = 2.5 in is solved. Thus, the length to thickness ratio for this thick beam is L/h = 8. The dynamic load P is equal to , as shown in Figure 22(b). Ten beam elements are used in the analysis, , time step s, and the error tolerance . The shear effect cannot be ignored in this case. Timoshenko and Euler-Bernoulli beam models have different results as can be seen in Figure 24. As expected, the Timoshenko model gives a higher deflection than Euler-Bernoulli model for the same load.

8. Conclusions

In this paper, a corotational finite element formulation for geometrically nonlinear static and dynamic analysis of planar beams and frames has been developed. The corotational frame translates and rotates with the beam element but does not deform with it. Hence, the rigid body motion is separated from the pure deformation which is always small with respect to the corotational frame. Hamilton’s principle has been utilized to derive the nonlinear equilibrium equations. In static analysis, a sequential numerical algorithm has been developed using an iterative method based on the full Newton-Raphson method without incremental loading. Furthermore, for dynamic analysis an incremental-iterative method based on the full Newton-Raphson method and the Newmark direct integration implicit method has been used to solve the nonlinear equilibrium equations. MATLAB codes have been written for both cases. Euler-Bernoulli beam model and Timoshenko beam model have been used to evaluate the stiffness matrices and the mass matrices. An efficient technique has been employed to determine the internal elastic force vector.

A wide range of numerical examples of beam and frame structures with large displacements have been solved and analyzed. The results have been compared with the published results to show the effectiveness and the accuracy of the proposed formulation and numerical algorithm. Though these problems have been solved using smaller number of elements than the published cases, accurate results have been obtained. The maximum number of iterations to reach the converged solution in this work is nine, which is reasonable. Superior to some published work which could not obtain acceptable solutions for some problems, in this work adequate solutions for all problems have been obtained.

Both thin and thick beams and frames have been investigated. For thin beams with large length to thickness ratio, both Euler-Bernoulli beam model and Timoshenko beam model are identical, while for thick beams with small length to thickness ratio meaningful differences can be perceived. Therefore, some results using Euler-Bernoulli beam model and Timoshenko beam model have been compared with the corresponding result of the commercial finite element program ANSYS using enormous number of eight-noded solid elements as the reference result. The results revealed that, for thick beams and frames, Timoshenko beam model produces more accurate results due to the inclusion of the shear effects.

It should be mentioned that using the geometric stiffness matrix in equilibrium equations is very important to improve the convergence properties, especially in dynamic analysis. More specifically the method to determine the axial force in this matrix is crucial. Also, the regular updating of stiffness and mass matrices, in each iteration, is very essential for convergence acceleration. With confidence, the proposed formulation is adequately accurate, simple, and computationally inexpensive.

Appendix

A. Stiffness Matrices

In this appendix, the stiffness matrices in (36) are specified. For Euler-Bernoulli beam model, the stiffness matrices arewhere is the modulus of elasticity, is the cross-sectional area, is the moment of inertia, andis the internal elastic force of the second node in direction. Stiffness matrices for Timoshenko beam model arewhere Q is defined asThe shear correction factor is a function of the cross-sectional shape, e.g., =5/6, for rectangular cross section and is Poisson’s ratio.

B. Mass Matrices

In this appendix, the stiffness matrices in (48) are specified. For Euler-Bernoulli beam model, the mass matrices are

For Timoshenko beam model, the mass matrices can be written as

C. Solution Strategies

The flowcharts for the strategies of the numerical algorithms in Section 6 are presented in this appendix: Figure 25 for static analysis and Figure 26 for dynamic analysis.

Data Availability

All the data used in the submitted manuscript are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

References

A. Beheshti, “Large deformation analysis of strain-gradient elastic beams,” Computers & Structures, vol. 177, pp. 162–175, 2016.
View at: Publisher Site | Google Scholar
P. Nanakorn and L. N. Vu, “A 2D field-consistent beam element for large displacement analysis using the total Lagrangian formulation,” Finite Elements in Analysis and Design, vol. 42, pp. 1240–1247, 2006.
View at: Publisher Site | Google Scholar
K. Meftah, W. Zouari, L. Sedira, and R. Ayad, “Geometric non-linear hexahedral elements with rotational DOFs,” Computational Mechanics, vol. 57, no. 1, pp. 37–53, 2016.
View at: Publisher Site | Google Scholar | MathSciNet
G. Shi and G. Z. Voyiadjis, “Geometrically nonlinear analysis of plates by assumed strain element with explicit tangent stiffness matrix,” Computers & Structures, vol. 41, no. 4, pp. 757–763, 1991.
View at: Publisher Site | Google Scholar
Y.-B. Yang and L.-J. Leu, “Force recovery procedures in nonlinear analysis,” Computers & Structures, vol. 41, no. 6, pp. 1255–1261, 1991.
View at: Publisher Site | Google Scholar
Y. B. Yang, S. R. 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
G. Wempner, “Finite elements, finite rotations and small strains of flexible shells,” International Journal of Solids and Structures, vol. 5, no. 2, pp. 117–153, 1969.
View at: Publisher Site | Google Scholar
M. Mostafa, M. V. Sivaselvan, and C. A. Felippa, “A solid-shell corotational element based on ANDES, ANS and EAS for geometrically nonlinear structural analysis,” International Journal for Numerical Methods in Engineering, vol. 95, no. 2, pp. 145–180, 2013.
View at: Publisher Site | Google Scholar | MathSciNet
X. Long, K. H. Tan, and C. K. Lee, “A 3D co-rotational beam element for steel and RC framed structures,” Structural Engineering and Mechanics, vol. 48, no. 5, pp. 587–613, 2013.
View at: Publisher Site | Google Scholar
C. Oran and A. Kassimali, “Large deformations of framed structures under static and dynamic loads,” Computers & Structures, vol. 6, no. 6, pp. 539–547, 1976.
View at: Publisher Site | Google Scholar
A. Alhasawi, P. Heng, M. Hjiaj, S. Guezouli, and J.-M. Battini, “Co-rotational planar beam element with generalized elasto-plastic hinges,” Engineering Structures, vol. 151, pp. 188–205, 2017.
View at: Publisher Site | Google Scholar
T. Wu, Z. Liu, and J. Hong, “A recursive formulation based on corotational frame for flexible planar beams with large displacement,” Journal of Central South University, vol. 25, no. 1, pp. 208–217, 2018.
View at: Publisher Site | Google Scholar
S. Faroughi, H. H. Khodaparast, and M. I. Friswell, “Non-linear dynamic analysis of tensegrity structures using a co-rotational method,” International Journal of Non-Linear Mechanics, vol. 69, pp. 55–65, 2015.
View at: Publisher Site | Google Scholar
K. M. Hsiao and F. Y. Hou, “Nonlinear finite element analysis of elastic frames,” Computers & Structures, vol. 26, no. 4, pp. 693–701, 1987.
View at: Publisher Site | Google Scholar
H. Kuo-Mo and J. Jing-Yuh, “Nonlinear dynamic analysis of elastic frames,” Computers & Structures, vol. 33, no. 4, pp. 1057–1063, 1989.
View at: Publisher Site | Google Scholar
H. A. Elkaranshawy and M. A. Dokainish, “Corotational finite element analysis of planar flexible multibody systems,” Computers & Structures, vol. 54, no. 5, pp. 881–890, 1995.
View at: Publisher Site | Google Scholar
H. A. Elkaranshawy, Dynamic Modelling of Flexible Multibody Systems Subjected to Frictional Impact, Ph.D Thesis, McMaster University, Hamilton, Canada, 1995.
K. Mattiasson, “Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals,” International Journal for Numerical Methods in Engineering, vol. 17, no. 1, pp. 145–153, 1981.
View at: Publisher Site | Google Scholar
M. E. James, T. J. Kozik, and J. E. Martinex, “Effect of curvature on nonlinear frame analyses,” Journal of the Structural Division, vol. 100, 1974.
View at: Google Scholar
T.-N. Le, J.-M. Battini, and M. Hjiaj, “Efficient formulation for dynamics of corotational 2D beams,” Computational Mechanics, vol. 48, no. 2, pp. 153–161, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
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
J. C. Simo and L. Vu-Quoc, “On the dynamics of flexible beams under large overall motions-The plane case: Part I and part II,” Journal of Applied Mechanics, vol. 53, no. 4, pp. 849–854, 1986.
View at: Publisher Site | Google Scholar
F. E. Valdivia and J. P. Ávila, “Nonlinear Dynamic Analysis of a Submerged Hose Using Co-rotational Formulation,” in Proceedings of the Conference: DINCON 2015 – Conferência Brasileira de Dinâmica, Controle e Aplicações, vol. 4, pp. 1–7, 2016.
View at: Publisher Site | Google Scholar
N. Friedman and A. Ibrahimbegovic, “Overview of Highly Flexible, Deployable Lattice Structures Used in Architecture and Civil Engineering Undergoing Large Displacements,” YBL Journal of Built Environment, vol. 1, no. 1, pp. 85–103, 2013.
View at: Publisher Site | Google Scholar
D. K. Nguyen and T. T. Tran, “A Corotational Formulation for Large Displacement Analysis of Functionally Graded Sandwich Beam and Frame Structures,” Mathematical Problems in Engineering, vol. 2016, Article ID 5698351, 12 pages, 2016.
View at: Publisher Site | Google Scholar | MathSciNet
W. Su, “Nonlinear aeroelastic analysis of aircraft with strut-braced highly flexible wings,” in Proceedings of the 58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, p. 1351, Grapevine, TX, USA, 2017.
View at: Google Scholar
C. Riso, F. Mastroddi, and C. E. Cesnik, “Coupled Flight Dynamics and Aeroelasticity of Very Flexible Aircraft Based on Commercial Finite Element Solvers,” in Proceedings of the AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, p. 1685, Kissimmee, Florida, USA, 2018.
View at: Publisher Site | Google Scholar
I. Fleming and D. J. Luscher, “A model for the structural dynamic response of the CX-100 wind turbine blade,” Wind Energy, vol. 17, no. 6, pp. 877–900, 2014.
View at: Publisher Site | Google Scholar
L. Dangora, J. Harvie, K. Wichman, and D. J. Luscher, “Experiment-based assessment of NLBeam for modeling geometrically nonlinear dynamic deformations,” in Proceedings of the Conference Proceedings of the Society for Experimental Mechanics Series, pp. 47–64, Springer, 2013.
View at: Google Scholar
P. Sivanagendra, Geometrically Nonlinear Elastic Analysis of Frames with Application to Vision-Based Force-Sensing and Mechanics of Plant Stems, Master Thesis, Indian Institute of Science Bangalore, India, 2006.
J. Rungamornrat and P. Tangnovarad, “Analysis of linearly elastic inextensible frames undergoing large displacement and rotation,” Mathematical Problems in Engineering, vol. 2011, Article ID 592958, 37 pages, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
J. Rungamornrat, S. Sihanartkatakul, and P. Kanchanakitcharoen, “Adaptive, Small-Rotation-Based, Corotational Technique for Analysis of 2D Nonlinear Elastic Frames,” Mathematical Problems in Engineering, vol. 2014, Article ID 640983, 15 pages, 2014.
View at: Publisher Site | Google Scholar
A. A. Elerian, W. I. Hussien, and H. A. Elkaranshawy, “A corotational formulation for geometrically nonlinear analysis of 2d beams and frames,” in Proceedings of the 26th CANCAM, vol. 29, Victoria, British Columbia, Canada, 2017.
View at: Google Scholar
A. A. Elerian, A Corotational Finite Element Formulation for Geometrically Nonlinear Static and Dynamic Analysis of Planar Beams and Frames, M.Sc. thesis, Alexandria university, Alexandria, Egypt, 2017.
D. K. Nguyen, “A Timoshenko beam element for large displacement analysis of planar beams and frames,” International Journal of Structural Stability and Dynamics, vol. 12, no. 6, Article ID 1250048, 2012.
View at: Publisher Site | Google Scholar | MathSciNet
B. S. Gan, T. T. Huong, and N. D. Kien, “Post-buckling Behaviour of Axially FGM Planar Beams and Frames,” Procedia Engineering, vol. 171, pp. 147–158, 2017.
View at: Google Scholar
K. J. Bathe, Finite Element Procedures, Prentice Hall Inc, New Jersey, USA, 2nd edition, 2014.
N. M. Newmark, “A method of computation for structural dynamics,” Journal of the Engineering Mechanical Division, vol. 85, pp. 67–94, 1959.
View at: Google Scholar

Copyright

Copyright © 2018 Hesham A. Elkaranshawy 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

3906

Downloads

1677

Citations