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

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 stiffnessmatrices and themassmatrices are evaluated using both Euler-Bernoulli andTimoshenko beammodels 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 iterativemethod based on the full Newton-Raphsonmethod 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.


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 .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.
Mathematical Problems in Engineering 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][2][3]22].In the updated Lagrangian formulation [4][5][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][35][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.

Kinematics of a Beam Element Using
Corotational Finite Element Formulation 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 i th 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 i th beam element nodal displacement vector in the fixed global coordinate system, at current configuration, is given by where   (n=1, 2) and W  (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 i th 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 j th equilibrium configuration.

Mathematical Problems in Engineering
For both the static and dynamic analysis, the nodal displacement vector of the i th 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  i and  i directions, respectively, and   (n=1, 2) are the counterclockwise rotations.The initial and current element chord lengths denoted as L o and L are shown in Figure 1.The initial and current locations of the corotational frame are specified by the two angles  o ,  that are measured counterclockwise from the fixed global x-axis direction to the local  i -axis. o and  are given by The relation between   ,   of the i th element is where  is the rigid body rotation and is obtained as The relation between the displacement vectors D i , d i and the rigid body displacement vector in the global coordinate system D r i takes the form where T i is the i th beam element transformation matrix, expressed by 2.2.The Displacement Interpolation.The classical Hermitian shape functions are used to relate the element deformation vector  i to the element nodal displacement vector d i as follows: The components of the shape functions N i of the i th beam element are given by where  is given by The transverse displacement  of the i th beam element takes the form In the corotational formulation the frame ( i ,  i ) is updated each iteration.The  i -axis is always passing through the two end nodes of the i th element as shown in Figure 1.Therefore, the displacement components   ( = 1, 2) are zeros.Thus (19) can be written as One can get the approximate current arc length S i , as follows: where ẃ is the derivative of the function (x) with respect to  i .This integral is evaluated using Gauss integration scheme with five Gauss points.The elongation of the i th beam element in   direction is This elongation can be approximated as The axial displacement of the first node  1 is chosen to be zero while the axial displacement of the second node  2 in ( 6) is defined as

The Strain Energy and Stiffness Matrices
Considering isotropic elastic materials, the relation between the stress vector  i and the strain vector  i for the i th beam element is given by where E i is the symmetric matrix of the elastic coefficients.
One can write the strain vector as follows: where D  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 i th beam element Π i is given by where V i is the volume.Substituting ( 27) and ( 28) into ( 29), the strain energy can be expressed by The element local displacement vector d i is independent of the volume V i .Subsequently, (30) can be simplified to the following form: where k i 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 D i =D i − D r i is the nodal pure displacement vector of the i th beam element and K i 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 k i can be expressed as where k 1 is the axial and bending stiffness matrix and k 2 is the geometric stiffness matrix for the i th 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][35][36], respectively.For completeness, they are stated in Appendix A.

The Kinetic Energy and Mass Matrices
The kinetic energy for the i th beam element KE i is given by where ( ̇) is the differentiation with respect to time t,  is the density in the i th beam element, and Ṙi is the velocity of a general point in the i th beam element with respect to the global coordinate system.For  i 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 i th beam element with respect to the local coordinate system ṙ i can be expressed as ṙ i can be written in the form where ḋ i is the nodal velocity vector of the i th beam element with respect to the local coordinate system.Furthermore, ḋ i can be expressed as where Ḋi is the nodal velocity vector of the i th beam element in the global coordinate system.Combining (39) − (42), the kinetic energy can be written as follows: The velocity vector Ḋi and transformation matrix T i are independent of the volume V i .Subsequently, (43) can be simplified to the following form: Mathematical Problems in Engineering where m i is the symmetric element mass matrix in the local corotational coordinate system, defined as follows: Equation ( 44) can be written in the form where M i 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 m i can be expressed as where m 1 is the i th beam element mass matrix of the translational inertia and m 2 is the mass matrix of i th 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.

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 L g is the Lagrangian functional.L g is defined as where KE i is the kinetic energy, PE i is the potential energy, and WD i is the work done by the external forces.In the structural applications, the potential energy PE i is the elastic strain energy Π i which is given in (34).WD i is defined by where    is the vector of the external applied forces on the i th 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  ḊT i 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 M i is constant during the time interval from t 1 to t 2 , one obtains where Di is the nodal acceleration vector of the i th beam element in the global coordinate system.The variation D i at t 1 and t 2 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 t 2 − t 1 , one can approximate D T i to D T i ; 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 i th beam element.Assembling the element force vectors, acceleration vectors, and mass matrices leads to the equation of motion of the overall structure which takes the form where F P is the vector of the external applied forces of the entire structure and F e is the internal elastic force vector of the entire structure.Both F P and F e are defined in the fixed global coordinate system.M is the mass matrix in the global coordinate system of the entire structure and D is the acceleration vector in the global coordinate system of the entire structure.One can introduce Rayleigh damping in (68) to read where C 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 C is defined in terms of the stiffness matrix and the mass matrix as where  and  are the damping coefficients which can be obtained from vibration modes of the system.K 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 i th beam element can be expressed as Therefore, the nonlinear equilibrium equation of the overall structure can be written as follows: F e − F P = 0 (72) 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:
(3) Assemble the stiffness matrices and the internal elastic force vectors to obtain K 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 R is calculated using full Newton-Raphson method as (7) R is added to the displacement vector of the entire structure, D. Consequently, D i in (1) can be extracted from D. Therefore, the current configuration can be determined using ( 9) − (12).( 8) Consequently, d i in (6) can be specified.It has to be noticed that, for solution stability purposes,  2 in ( 6) is computed from ( 24), (22), and ( 23). ( 9) Start new iteration.
The flowchart for the strategy is given in Appendix C.

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 (1) Compute k i , m i , and    for each element using ( 36), (48), and (66), respectively.
(3) Assemble the elements stiffness and mass matrices and the internal elastic force vectors to obtain K, M, and F e for the entire structure.
(5) If the convergence criterion satisfies the condition in (74), stop the iteration and go to step (7).Otherwise, start the iteration: (a) a displacement corrector R is calculated using full Newton-Raphson method as where K is the effective matrix, and it can be defined as where  and  are the Newmark parameters.(b) Calculate the incremental displacement vector ΔD as where ΔD 0 is the incremental displacement vector of the previous iteration.ΔD 0 values are considered to be zeros in the first iteration.(c) ΔD i can be extracted for each element from ΔD.
(e) Using the Newmark direct integration method, Ḋ+1 and DN+1 can be expressed as It has to be noticed that the Newmark implicit method is unconditionally stable if  ≥ 0.5 and  ≥ (2 + 1) 2 /16.In the present work, Newmark parameters are chosen to be  = 0.25 and  = 0.5.
(6) Go to the beginning of step (5) again.
(7) The displacement vector of the overall structure at time t = t N + Δ can be determined as (8) Start the next time step.
The flowchart for the strategy is given in Appendix C.

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.

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,  2 =  1 , and the error tolerance   is chosen to be 10 −6 .Young's modulus E is considered to be 200 GPa and Poisson's ratio is V = 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 2 /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 U/L Mattiasson [18] U/L Nanakorn and Vu   [18] in Figure 7 and Table 4. U/L and W/L are the normalized displacements and PL 2 /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.

Mathematical Problems in Engineering
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.This problem has been solved before in [1,10,14,19].The solution of James et al. [19] failed to converge approximately after P = 4000 kips for n = 0.1 and after P = 3000 kips for n = 0.5 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 n = 0.5 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 n = 0.5 case but also for n = 0.1 case.Four elements for each member are used, u 2 = e 2 , and the used error tolerance is e r = 10 −6 .
As expected the Timoshenko model gives a higher deflection than Euler-Bernoulli model for the same load.

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

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 P 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 EA = 10 6 , EI = 1, 000, A = 1, I = 10, and V = 0 and  is the damping coefficient.Only four beam elements are used in the analysis,  2 =  1 , time step Δ = 0.05, and the error tolerance   is chosen to be 10 −5 .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]  Reference solution in [20] Formulation of [20] Lumped mass matrix of [20] Present study  Reference solution in [20] Formulation of [20] Lumped mass matrix of [20] Present study  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.
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 4 X 10 5 lb, as shown in Figure 22   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.

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.

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  i direction.
Stiffness matrices for Timoshenko beam model are  The shear correction factor   is a function of the crosssectional shape, e.g.,   =5/6, for rectangular cross section and V 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

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 1 :
Figure 1: Corotational coordinate systems for the i th beam element.

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

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

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

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

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

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

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

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

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

Figure 17 :
Figure 17: Time-vertical displacement curve for cantilever beam with an end vertical sinusoidal point load (results comparison).
(b).Ten beam elements are used in the analysis,  2 =  1 , time step Δ = 25 s, and the error tolerance   = 10 −5 .The shear effect cannot be ignored

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

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

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

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

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

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

Table 3 .
Considering the ANSYS program results as the reference results, Timoshenko beam model gives more reasonable results than Euler-Bernoulli beam mode.

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

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

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