Mathematical Problems in Engineering
Volume 2018 (2018), Article ID 2670462, 22 pages
https://doi.org/10.1155/2018/2670462

Research Article

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

Department of Engineering Mathematics and Physics, Alexandria University, Alexandria 21544, Egypt

Correspondence should be addressed to Hesham A. Elkaranshawy; hesham_elk@alexu.edu.eg

Received 19 May 2018; Accepted 11 July 2018; Published 14 August 2018

Academic Editor: Nerio Tullini

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.

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.

Figure 1: Corotational coordinate systems for the beam element.

The beam element nodal displacement vector in the fixed global coordinate system, at current configuration, is given by

where (n=1, 2) and Wn (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 as

where (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 as

where , , 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, is

where (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 is

where β 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 form

where 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 by

where ξ is given by

The transverse displacement of the beam element takes the form

In 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 is

This elongation can be approximated as

The 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 by

where 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 as

The strain energy for the beam element is given by

where is the volume. Substituting (27) and (28) into (29), the strain energy can be expressed by

The 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 as

Substituting (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 as

where 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 by

where 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 write

where I is the identify matrix. Consequently, one can write (37) as

The 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 form

where is the nodal velocity vector of the beam element with respect to the local coordinate system. Furthermore, can be expressed as

where 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 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 be expressed as

where 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 form

where is the Lagrangian functional. is defined as

where 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 by

where 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 as

Using (34), the strain energy term in (52) can be expressed as

Using (51), one can write the work done term in (52) as

Substituting (53), (54), and (55) in (52) gives

The 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 obtains

where 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 as

Substituting (60) in (56), one obtains

For small time step , one can approximate to ; hence (61) can be expressed as

The only way for the integration in (62) to be zero is the vanishing of the integrand itself. Therefore, (62) can be rewritten as

As mentioned before, the displacement variation is arbitrary; then

According to (13), (34), (35), and (64), one can write

where 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 (65) in (64) gives

Equation (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 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. One can introduce Rayleigh damping in (68) to read

where 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

where 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 as

Therefore, 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 as

where is the out of balance force. The iteration equilibrium convergence criterion is determined by

where 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 as

At 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].

Table 1: Comparison of results for cantilever beam subjected to end vertical load, h=0.25m.

U / L
W / L

Mattiasson [18]
Nanakorn and Vu [2]
Present
Mattiasson [18]
Nanakorn and Vu [2]
Present

1
0.05643
0.05559
0.05598
0.30172
0.29946
0.30189

2
0.16064
0.15672
0.15948
0.49346
0.48748
0.49418

3
0.25442
0.24757
0.25269
0.60325
0.59534
0.60459

4
0.32894
0.32003
0.32677
0.66996
0.66126
0.67186

5
0.38763
0.37733
0.38511
0.71379
0.70479
0.71616

6
0.43459
0.42335
0.43177
0.74457
0.73550
0.74736

7
0.47293
0.46103
0.46986
0.76737
0.75831
0.77053

8
0.50483
0.49245
0.50153
0.78498
0.77597
0.78848

9
0.53182
0.51909
0.52832
0.79906
0.79011
0.80029

10
0.55500
0.54199
0.55131
0.81061
0.80173
0.81471

Table 2: Displacement errors for cantilever beam subjected to end vertical load, h=0.25m.

U / L
W / L

Error Nanakorn and Vu [2]
Error present
Error Nanakorn and Vu [2]
Error present

1
0.00084
0.00045
0.00226
0.00017

2
0.00392
0.00116
0.00598
0.00072

3
0.00685
0.00173
0.00791
0.00134

4
0.00891
0.00217
0.00870
0.00190

5
0.01030
0.00252
0.00900
0.00237

6
0.01124
0.00282
0.00907
0.00279

7
0.01190
0.00307
0.00906
0.00316

8
0.01238
0.00330
0.00901
0.00350

9
0.01273
0.00350
0.00895
0.00123

10
0.01301
0.00369
0.00888
0.00410

Figure 2: Cantilever beam with an end vertical point load.

Figure 3: Dimensionless load-displacement curves of the cantilever beam subjected to end vertical point load, h=0.25m.

Figure 4: The displacement error curves of the cantilever beam subjected to end vertical point load, h=0.25m.

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.

Table 3: Comparison of the present results for cantilever beam subjected to end vertical load, h=2.5m.

U / L
W / L

Euler Bernoulli
Timoshenko
ANSYS
Euler Bernoulli
Timoshenko
ANSYS

1
0.0545
0.0573
0.0566
0.3024
0.3115
0.3141

2
0.1554
0.1616
0.1568
0.4971
0.5096
0.5140

3
0.2462
0.2545
0.2526
0.6111
0.6250
0.6332

4
0.3184
0.3278
0.3253
0.6824
0.6975
0.7081

5
0.3751
0.3854
0.3863
0.7309
0.7473
0.7595

6
0.4204
0.4314
0.4324
0.7665
0.7840
0.7978

7
0.4573
0.4689
0.4735
0.7941
0.8128
0.8279

8
0.4879
0.5001
0.5091
0.8165
0.8364
0.8527

9
0.5137
0.5265
0.5383
0.8353
0.8564
0.8739

10
0.5358
0.5491
0.5634
0.8517
0.8739
0.8927

Figure 5: Dimensionless load-displacement curves of the cantilever beam subjected to end vertical point load, h=2.5m.

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.

Table 4: Comparison of results for pinned-fixed square diamond frame loaded in tension, h=0.1m.

U / L
W / L

Mattiasson [18]
Nanakorn and Vu [2]
Present
Mattiasson [18]
Nanakorn and Vu [2]
Present

1
0.13960
0.15592
0.13908
0.11252
0.12276
0.11331

2
0.23184
0.24894
0.23073
0.16429
0.17228
0.16637

3
0.29447
0.30958
0.29295
0.19183
0.19750
0.19507

4
0.33940
0.35211
0.33757
0.20839
0.21238
0.21259

5
0.37322
0.38367
0.37116
0.21931
0.22211
0.22429

6
0.39966
0.40810
0.39739
0.22703
0.22895
0.23265

7
0.42097
0.42765
0.41852
0.23279
0.23405
0.23893

8
0.43855
0.44370
0.43594
0.23726
0.23801
0.24384

9
0.45335
0.45716
0.45059
0.24084
0.24120
0.24781

10
0.46601
0.46865
0.46401
0.24380
0.24383
0.24951

Table 5: Displacement errors for pinned-fixed square diamond frame loaded in tension, h=0.1m.

U / L
W / L

Error Nanakorn and Vu [2]
Error present
Error Nanakorn and Vu  [2]
Error present

1
0.01632
0.00052
0.01024
0.00079

2
0.01710
0.00111
0.00799
0.00208

3
0.01511
0.00152
0.00567
0.00324

4
0.01271
0.00183
0.00399
0.00420

5
0.01045
0.00206
0.00280
0.00498

6
0.00844
0.00227
0.00192
0.00562

7
0.00668
0.00245
0.00126
0.00614

8
0.00515
0.00261
0.00075
0.00658

9
0.00381
0.00276
0.00036
0.00697

10
0.00264
0.00200
0.00003
0.00571

Figure 6: Pinned-fixed square diamond frame loaded in tension.

Figure 7: Dimensionless load-displacement curves of the pinned-fixed square diamond frame loaded in tension, h=0.1m.

Figure 8: Displacement error curves of the pinned-fixed square diamond frame loaded in tension, h = 0.1 m.

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.

Table 6: Comparison of the present results for pinned-fixed square diamond frame loaded in tension, h = 2.8284 m.

U / L
W / L

Euler Bernoulli
Timoshenko
ANSYS
Euler Bernoulli
Timoshenko
ANSYS

1
0.1375
0.1402
0.1398
0.1156
0.1174
0.1169

2
0.2280
0.2319
0.2302
0.1715
0.1739
0.1729

3
0.2893
0.2940
0.2910
0.2032
0.2058
0.2051

4
0.3331
0.3383
0.3381
0.2239
0.2266
0.2264

5
0.3661
0.3716
0.3715
0.2388
0.2417
0.2418

6
0.3917
0.3976
0.3974
0.2503
0.2534
0.2535

7
0.4123
0.4185
0.4184
0.2598
0.2630
0.2632

8
0.4292
0.4357
0.4355
0.2680
0.2712
0.2714

9
0.4434
0.4501
0.4499
0.2752
0.2785
0.2787

10
0.4555
0.4625
0.4624
0.2817
0.2851
0.2853

Figure 9: Dimensionless load-displacement curves of the pinned-fixed square diamond frame loaded in tension, h = 2.8284 m.

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.

Figure 10: Portal frame subjected to concentrated vertical and horizontal loads.

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.

Figure 11: Load-displacement curves of portal frame for loading case n = 0.1.

Figure 12: Load-displacement curves of portal frame for loading case n = 0.5.

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.

Figure 13: Cantilever beam with an end vertical dynamic sinusoidal point load.

Figure 14: Time-horizontal displacement curve for cantilever beam with an end vertical sinusoidal point load (present results).

Figure 15: Time-vertical displacement curve for cantilever beam with an end vertical sinusoidal point load (present results).

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

Figure 16: Time-horizontal displacement curve for cantilever beam with an end vertical sinusoidal point load (results comparison).

Figure 17: Time-vertical displacement curve for cantilever beam with an end vertical sinusoidal point load (results comparison).

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

Figure 18: Cantilever beam with an end vertical concentrated dynamic load.

Figure 19: The dynamic load history.

Figure 20: Dynamic response of cantilever beam.

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.

Figure 21: Fixed-fixed beam with a vertical dynamic point load at the midspan.

(a)

(b)

Figure 22: Dynamic Load history in the two case studies; (a) first case and (b) second case.

Figure 23: Time-vertical displacement curve for fixed-fixed thin beam.

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.

Figure 24: Time-vertical displacement curve for fixed-fixed thick beam.

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 are

where is the modulus of elasticity, is the cross-sectional area, is the moment of inertia, and

is the internal elastic force of the second node in direction. Stiffness matrices for Timoshenko beam model are

where Q is defined as

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

Figure 25: Solution strategy for the geometrically nonlinear static analysis.

Figure 26: Solution strategy for the geometrically nonlinear 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.

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.

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.

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.

Y.-B. Yang and L.-J. Leu, “Force recovery procedures in nonlinear analysis,” Computers & Structures, vol. 41, no. 6, pp. 1255–1261, 1991.

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.

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.

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.

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.

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.

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.

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.

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.

K. M. Hsiao and F. Y. Hou, “Nonlinear finite element analysis of elastic frames,” Computers & Structures, vol. 26, no. 4, pp. 693–701, 1987.

H. Kuo-Mo and J. Jing-Yuh, “Nonlinear dynamic analysis of elastic frames,” Computers & Structures, vol. 33, no. 4, pp. 1057–1063, 1989.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.