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

Authors: Hesham A. Elkaranshawy; Ahmed A. H. Elerian; Walied I. Hussien
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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 Section2. 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 (xi, zi). 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 ith 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 theith beam element.Theith beam element nodal displacement vector in the fixed global coordinate system, at current configuration, is given by(1)Di=U1W1θ1U2W2θ2Twhere Un (n=1, 2) andWn (n=1, 2) are the displacement components in X and Z directions, respectively, and θn (n =1, 2) are the counterclockwise rotations.In the kinematics of the dynamic analysis, as shown in Figure1, the nodal incremental displacement vector of the ith beam element in the global coordinate system is defined as(2)ΔDi=ΔU1ΔW1Δθ1ΔU2ΔW2Δθ2Twhere ΔUn(n=1, 2) and ΔWn (n=1, 2) are the incremental displacement components in X and Z directions, respectively, and Δθn (n=1, 2) are the counterclockwise incremental rotations. Thus, the displacement components in (1) can be identified as(3)Un=Unj+ΔUn(4)Wn=Wnj+ΔWn(5)θn=θnj+Δθnwhere Unj, Wnj, and θnj are the displacement components at the jth equilibrium configuration.For both the static and dynamic analysis, the nodal displacement vector of theith beam element in the element local coordinate system, at current configuration, is(6)di=u1w1θ1¯u2w2θ2¯Twhere un (n=1, 2) and wn (n=1, 2) are the displacement components in xi and zi directions, respectively, and θn¯(n=1, 2) are the counterclockwise rotations. The initial and current element chord lengths denoted as Lo 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 globalx-axis direction to the local xi-axis. αo and α are given by (7)sin⁡αo=Z2-Z1Lo(8)cos⁡αo=X2-X1Lo(9)sin⁡α=Z2+W2-Z1-W1L(10)cos⁡α=X2+U2-X1-U1LThe relation between θn,θn¯ of the ith element is(11)θn¯=θn–βwhere β is the rigid body rotation and is obtained as (12)β=α-αoThe relation between the displacement vectors Di, di and the rigid body displacement vector in the global coordinate system Dir takes the form(13)Di=Dir+Tidiwhere Ti is the ith beam element transformation matrix, expressed by(14)Ti=Ai00Ai,Ai=cos⁡α-sin⁡α0sin⁡αcos⁡α0001


### 2.2. The Displacement Interpolation

The classical Hermitian shape functions are used to relate the element deformation vectorvi to the element nodal displacement vector di as follows:(15)vi=NidiThe components of the shape functions Ni of the ith beam element are given by(16)N1=1-121+ξ,N2=141-ξ22+ξN3=L81-ξ21-ξ,N4=121+ξ(17)N5=141+ξ22-ξ,N6=L81+ξ21+ξwhere ξ is given by(18)ξ=-1+2xLThe transverse displacement w of the ith beam element takes the form(19)wx=N2w1+N3θ1¯+N5w2+N6θ2¯In the corotational formulation the frame (xi,zi) is updated each iteration. The xi-axis is always passing through the two end nodes of the ith element as shown in Figure 1. Therefore, the displacement components wn(n=1,2) are zeros. Thus (19) can be written as (20)wx=N3θ1¯+N6θ2¯One can get the approximate current arc length Si, as follows:(21)Si=L2∫-11w´2+1dξwhere w´ is the derivative of the function w(x) with respect to xi. This integral is evaluated using Gauss integration scheme with five Gauss points. The elongation of the ith beam element in xi direction is(22)e=e1=Si-LoThis elongation can be approximated as(23)e=e2=L-LoThe axial displacement of the first node u1is chosen to be zero while the axial displacement of the second node u2 in (6) is defined as(24)u2=e


## 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 (xi, zi). 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 ith 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 theith beam element.Theith beam element nodal displacement vector in the fixed global coordinate system, at current configuration, is given by(1)Di=U1W1θ1U2W2θ2Twhere Un (n=1, 2) andWn (n=1, 2) are the displacement components in X and Z directions, respectively, and θn (n =1, 2) are the counterclockwise rotations.In the kinematics of the dynamic analysis, as shown in Figure1, the nodal incremental displacement vector of the ith beam element in the global coordinate system is defined as(2)ΔDi=ΔU1ΔW1Δθ1ΔU2ΔW2Δθ2Twhere ΔUn(n=1, 2) and ΔWn (n=1, 2) are the incremental displacement components in X and Z directions, respectively, and Δθn (n=1, 2) are the counterclockwise incremental rotations. Thus, the displacement components in (1) can be identified as(3)Un=Unj+ΔUn(4)Wn=Wnj+ΔWn(5)θn=θnj+Δθnwhere Unj, Wnj, and θnj are the displacement components at the jth equilibrium configuration.For both the static and dynamic analysis, the nodal displacement vector of theith beam element in the element local coordinate system, at current configuration, is(6)di=u1w1θ1¯u2w2θ2¯Twhere un (n=1, 2) and wn (n=1, 2) are the displacement components in xi and zi directions, respectively, and θn¯(n=1, 2) are the counterclockwise rotations. The initial and current element chord lengths denoted as Lo 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 globalx-axis direction to the local xi-axis. αo and α are given by (7)sin⁡αo=Z2-Z1Lo(8)cos⁡αo=X2-X1Lo(9)sin⁡α=Z2+W2-Z1-W1L(10)cos⁡α=X2+U2-X1-U1LThe relation between θn,θn¯ of the ith element is(11)θn¯=θn–βwhere β is the rigid body rotation and is obtained as (12)β=α-αoThe relation between the displacement vectors Di, di and the rigid body displacement vector in the global coordinate system Dir takes the form(13)Di=Dir+Tidiwhere Ti is the ith beam element transformation matrix, expressed by(14)Ti=Ai00Ai,Ai=cos⁡α-sin⁡α0sin⁡αcos⁡α0001


## 2.2. The Displacement Interpolation

The classical Hermitian shape functions are used to relate the element deformation vectorvi to the element nodal displacement vector di as follows:(15)vi=NidiThe components of the shape functions Ni of the ith beam element are given by(16)N1=1-121+ξ,N2=141-ξ22+ξN3=L81-ξ21-ξ,N4=121+ξ(17)N5=141+ξ22-ξ,N6=L81+ξ21+ξwhere ξ is given by(18)ξ=-1+2xLThe transverse displacement w of the ith beam element takes the form(19)wx=N2w1+N3θ1¯+N5w2+N6θ2¯In the corotational formulation the frame (xi,zi) is updated each iteration. The xi-axis is always passing through the two end nodes of the ith element as shown in Figure 1. Therefore, the displacement components wn(n=1,2) are zeros. Thus (19) can be written as (20)wx=N3θ1¯+N6θ2¯One can get the approximate current arc length Si, as follows:(21)Si=L2∫-11w´2+1dξwhere w´ is the derivative of the function w(x) with respect to xi. This integral is evaluated using Gauss integration scheme with five Gauss points. The elongation of the ith beam element in xi direction is(22)e=e1=Si-LoThis elongation can be approximated as(23)e=e2=L-LoThe axial displacement of the first node u1is chosen to be zero while the axial displacement of the second node u2 in (6) is defined as(24)u2=e


## 3. The Strain Energy and Stiffness Matrices

Considering isotropic elastic materials, the relation between the stress vectorσi and the strain vector ϵi for the ith beam element is given by(25)σi=Eiϵiwhere Ei is the symmetric matrix of the elastic coefficients. One can write the strain vector as follows:(26)ϵi=Dfviwhere Df is the differential operator matrix. Substituting (15) into (26), the strain vector can be expressed as follows:(27)ϵi=DfNidiCombining (25) and (27), the stress vector can be rewritten as(28)σi=EiDfNidiThe strain energy for the ith beam element Πi is given by(29)Πi=12∫ViσiTϵidViwhere Vi is the volume. Substituting (27) and (28) into (29), the strain energy can be expressed by(30)Πi=12∫ViEiDfNidiTDfNididViThe element local displacement vector di is independent of the volume Vi. Subsequently, (30) can be simplified to the following form:(31)Πi=12diTkidiwhere ki is the symmetric element stiffness matrix in the local coordinate system, defined as follows:(32)ki=∫ViBTEiBdViandB is defined as(33)B=DfNiSubstituting (13) into (31), one can write (31) in the following form:(34)Πi=12TiTDi-DirTkiTiTDi-Dir=12D-iTKiD-iwhere D-i=Di-Dir is the nodal pure displacement vector of the ith beam element and Ki is the symmetric element stiffness matrix in the fixed global coordinate system, defined as follows:(35)Ki=TikiTiTFor the beam element used in this research work, the element stiffness matrix in the local coordinate system ki can be expressed as(36)ki=k1+k2where k1 is the axial and bending stiffness matrix and k2 is the geometric stiffness matrix for the ith 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 theith beam element KEi is given by(37)KEi=12∫ViρR˙iTR˙idViwhere ˙ is the differentiation with respect to time t, ρ is the density in the ith beam element, and R˙i is the velocity of a general point in the ith beam element with respect to the global coordinate system. For Ai in (14), one can write(38)AiAiT=IwhereI is the identify matrix. Consequently, one can write (37) as(39)KEi=12∫ViρAiTR˙iTAiTR˙idViThe velocity of a general point in the ith beam element with respect to the local coordinate system r˙i can be expressed as(40)r˙i=AiTR˙ir˙i can be written in the form(41)r˙i=Nid˙iwhere d˙i is the nodal velocity vector of the ith beam element with respect to the local coordinate system. Furthermore, d˙i can be expressed as(42)d˙i=TiTD˙iwhere D˙i is the nodal velocity vector of the ith beam element in the global coordinate system. Combining (39) − (42), the kinetic energy can be written as follows:(43)KEi=12∫ViρNiTiTD˙iTNiTiTD˙idViThe velocity vector D˙i and transformation matrix Ti are independent of the volume Vi. Subsequently, (43) can be simplified to the following form:(44)KEi=12D˙iTTimiTiTD˙iwhere mi is the symmetric element mass matrix in the local corotational coordinate system, defined as follows:(45)mi=12∫ViρNiTNidViEquation (44) can be written in the form(46)KEi=12D˙iTMiD˙iwhere Mi is the symmetric element mass matrix in the fixed global coordinate system, defined as follows: (47)Mi=TimiTiTFor the beam element used in this research work, the element mass matrix in the local coordinate system mi can be expressed as(48)mi=m1+m2where m1 is the ith beam element mass matrix of the translational inertia and m2 is the mass matrix of ith 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(49)δ∫t1t2Lgdt=0where Lg is the Lagrangian functional. Lg is defined as(50)Lg=KEi-PEi+WDiwhere KEi is the kinetic energy, PEi is the potential energy, and WDi is the work done by the external forces. In the structural applications, the potential energy PEi is the elastic strain energy Πi which is given in (34). WDi is defined by(51)WDi=DiTFiPwhere FiP is the vector of the external applied forces on the ith 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:(52)∫t1t2δKEi-Πi+WDidt=0Using (46), kinetic energy term in (52) can be written as(53)δKEi=δ12D˙iTMiD˙i=δD˙iTMiD˙iUsing (34), the strain energy term in (52) can be expressed as(54)δΠi=δ12D-iTKiD-i=δD-iTKiD-iUsing (51), one can write the work done term in (52) as(55)δWDi=δDiTFiP=δDiTFiPSubstituting (53), (54), and (55) in (52) gives(56)∫t1t2δD˙iTMiD˙i-δD-iTKiD-i+δDiTFiPdt=0The differentiation with respect to time and the variation are interchangeable. Therefore, the term δD˙iT in (56) can be written in the following form:(57)δD˙iT=δdDiTdt=ddtδDiTConsequently, substituting (57) in the kinetic energy term in (56) and integrating the same term by parts, putting in mind that Mi is constant during the time interval from t1 to t2, one obtains(58)∫t1t2δD˙iTMiD˙idt=∫t1t2ddtδDiTMiD˙idt(59)∫t1t2δD˙iTMiD˙idt=δDiTMiD˙it1t2-∫t1t2δDiTMiD¨idtwhere D¨i is the nodal acceleration vector of the ith beam element in the global coordinate system. The variation δDi at t1 and t2 is zero. Thus, the first term on the right hand side in (60) vanishes, Equation (60) can be written as(60)∫t1t2δD˙iTMiD˙idt=-∫t1t2δDiTMiD¨idtSubstituting (60) in (56), one obtains(61)∫t1t2-δDiTMiD¨i-δD-iTKiD-i+δDiTFiPdt=0For small time step t2-t1, one can approximate δD-iT to δDiT; hence (61) can be expressed as(62)∫t1t2δDiT-MiD¨i-KiD-i+FiPdt=0The only way for the integration in (62) to be zero is the vanishing of the integrand itself. Therefore, (62) can be rewritten as(63)δDiT-MiD¨i-KiD-i+FiP=0As mentioned before, the displacement variation is arbitrary; then(64)KiD-i+MiD¨i-FiP=0According to (13), (34), (35), and (64), one can write(65)Fie=KiD-i=Tifiewhere Fie is the internal elastic force vector in the fixed global and fie is the internal elastic force vector in the local coordinate systems, which is given by(66)fie=kidiSubstituting (65) in (64) gives(67)Fie+MiD¨i-FiP=0Equation (67) is the equation of motion of the ith 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(68)Fe+MD¨-FP=0where FP is the vector of the external applied forces of the entire structure and Fe is the internal elastic force vector of the entire structure. Both FP and Fe 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(69)Fe+CD˙+MD¨-FP=0where C is the damping matrix in the global coordinate system of the entire structure and D˙ is the velocity vector in the global coordinate system of the entire structure. The damping matrix C is defined in terms of the stiffness matrix and the mass matrix as(70)C=μM+ηKwhere μ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 ith beam element can be expressed as(71)Fie-FiP=0Therefore, the nonlinear equilibrium equation of the overall structure can be written as follows:(72)Fe-FP=0


## 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(73)φ=Fe-FP=0where φ is the out of balance force. The iteration equilibrium convergence criterion is determined by(74)φ≤erφfwhere  is the Euclidean norm, φf is the reference out of balance force, and er is the error tolerance. φf is considered to be the out of balance force obtained in the first iteration. The error tolerance er 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)
Computeki and fie for each element using (36) and (66), respectively.(2)
DetermineKi and Fie according to (35) and (65), respectively.(3)
Assemble the stiffness matrices and the internal elastic force vectors to obtainK and Fe 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 correctorR is calculated using full Newton-Raphson method as(75)R=-K-1φ(7)
R is added to the displacement vector of the entire structure, D. Consequently, Di in (1) can be extracted from D. Therefore, the current configuration can be determined using (9) − (12).(8)
Consequently,di in (6) can be specified. It has to be noticed that, for solution stability purposes, u2 in (6) is computed from (24), (22), and (23).(9)
Start new iteration.The flowchart for the strategy is given in AppendixC.


### 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(76)φ=Fe+CD˙+MD¨-FP=0At the beginning of the solution, displacement, velocity, and acceleration vectors are assumed to be nulls. DN, D˙N, and D¨N are the values of D, D˙, and D¨ at a known equilibrium configuration, at time t=tN. Likewise, DN+1, D˙N+1, and D¨N+1 are the values of D, D˙, and D¨ at time t=tN+Δt, where Δt 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)
Computeki, mi, and fie for each element using (36), (48), and (66), respectively.(2)
DetermineKi, Mi, and Fie according to (35), (47), and (65), respectively.(3)
Assemble the elements stiffness and mass matrices and the internal elastic force vectors to obtainK, M, and Fe 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 (7). Otherwise, start the iteration:
(a)
a displacement correctorR is calculated using full Newton-Raphson method as(77)R=-K--1φwhereK- is the effective matrix, and it can be defined as(78)K-=1τΔt2M+γτΔtC+Kwhereτ and γ are the Newmark parameters.(b)
Calculate the incremental displacement vectorΔD as(79)ΔD=ΔD0+RwhereΔD0 is the incremental displacement vector of the previous iteration. ΔD0 values are considered to be zeros in the first iteration.(c)
ΔDi can be extracted for each element from ΔD. Then, Di in (1) and consequently the current configuration can be specified using (3)−(5) and (9)−(12).(d)
Using (6), di can be determined. It has to be noticed that u2 in (6) is equated to e1 in (22) and (24).(e)
Using the Newmark direct integration method,D˙N+1 and D¨N+1 can be expressed as(80)D¨N+1=ΔDτΔt2-D˙NτΔt-12τ-1D¨N(81)D˙N+1=D˙N+Δt1-γD¨N+γD¨N+1It 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 timet=tN+Δt can be determined as(82)DN+1=DN+ΔD(8)
Start the next time step.The flowchart for the strategy is given in AppendixC.


## 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(73)φ=Fe-FP=0where φ is the out of balance force. The iteration equilibrium convergence criterion is determined by(74)φ≤erφfwhere  is the Euclidean norm, φf is the reference out of balance force, and er is the error tolerance. φf is considered to be the out of balance force obtained in the first iteration. The error tolerance er 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)
Computeki and fie for each element using (36) and (66), respectively.(2)
DetermineKi and Fie according to (35) and (65), respectively.(3)
Assemble the stiffness matrices and the internal elastic force vectors to obtainK and Fe 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 correctorR is calculated using full Newton-Raphson method as(75)R=-K-1φ(7)
R is added to the displacement vector of the entire structure, D. Consequently, Di in (1) can be extracted from D. Therefore, the current configuration can be determined using (9) − (12).(8)
Consequently,di in (6) can be specified. It has to be noticed that, for solution stability purposes, u2 in (6) is computed from (24), (22), and (23).(9)
Start new iteration.The flowchart for the strategy is given in AppendixC.


## 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(76)φ=Fe+CD˙+MD¨-FP=0At the beginning of the solution, displacement, velocity, and acceleration vectors are assumed to be nulls. DN, D˙N, and D¨N are the values of D, D˙, and D¨ at a known equilibrium configuration, at time t=tN. Likewise, DN+1, D˙N+1, and D¨N+1 are the values of D, D˙, and D¨ at time t=tN+Δt, where Δt 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)
Computeki, mi, and fie for each element using (36), (48), and (66), respectively.(2)
DetermineKi, Mi, and Fie according to (35), (47), and (65), respectively.(3)
Assemble the elements stiffness and mass matrices and the internal elastic force vectors to obtainK, M, and Fe 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 (7). Otherwise, start the iteration:
(a)
a displacement correctorR is calculated using full Newton-Raphson method as(77)R=-K--1φwhereK- is the effective matrix, and it can be defined as(78)K-=1τΔt2M+γτΔtC+Kwhereτ and γ are the Newmark parameters.(b)
Calculate the incremental displacement vectorΔD as(79)ΔD=ΔD0+RwhereΔD0 is the incremental displacement vector of the previous iteration. ΔD0 values are considered to be zeros in the first iteration.(c)
ΔDi can be extracted for each element from ΔD. Then, Di in (1) and consequently the current configuration can be specified using (3)−(5) and (9)−(12).(d)
Using (6), di can be determined. It has to be noticed that u2 in (6) is equated to e1 in (22) and (24).(e)
Using the Newmark direct integration method,D˙N+1 and D¨N+1 can be expressed as(80)D¨N+1=ΔDτΔt2-D˙NτΔt-12τ-1D¨N(81)D˙N+1=D˙N+Δt1-γD¨N+γD¨N+1It 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 timet=tN+Δt can be determined as(82)DN+1=DN+ΔD(8)
Start the next time step.The flowchart for the strategy is given in AppendixC.


## 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 Figure2. Five beam elements are used in the analysis, u2=e1, and the error tolerance er 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 PL2/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.
P L 2 E I 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.81471Table 2
Displacement errors for cantilever beam subjected to end vertical load, h=0.25m.
P L 2 E I 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.00410Figure 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 Figure5 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.
P L 2 E I 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.8927Figure 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 Figure6 is considered in this example. Two elements per member are used in the analysis, u2=e1, and the error tolerance er is chosen to be 10-7. 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 PL2/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.
P L 2 E I 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.24951Table 5
Displacement errors for pinned-fixed square diamond frame loaded in tension, h=0.1m.
P L 2 E I 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.00571Figure 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 Figure9 and Table 6.Table 6
Comparison of the present results for pinned-fixed square diamond frame loaded in tension, h = 2.8284 m.
P L 2 E I 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.2853Figure 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 Figure10 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 P=4000kips for n=0.1 and after P=3000kips for n=0.5 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 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, u2=e2, and the used error tolerance is er=10-6. 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 loadP=Posin(⁡ωt) at the free end as shown in Figure 13. Young’s modulus E is considered 210 GPa, the material’s density ρ = 7,850 kg/m3, and Poisson’s ratio v = 0.3. Po=10MN and its frequency ω= 50 rad/s. Ten beam elements are used in the analysis, u2=e1, time step Δt=10-4s, and the error tolerance er=10-5. 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 loadP 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=106, EI=1,000, ρA=1, ρI=10, and v=0 and μ is the damping coefficient. Only four beam elements are used in the analysis, u2=e1, time step Δt=0.05, and the error tolerance er 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].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 Figure21 is considered in this example. Young’s modulus E is 30,000 ksi, the material’s density ρ is 2.536X10-4lbs2/in4, and Poisson’s ratio v 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 640lb as shown in Figure 22(a). Ten beam elements are used in the analysis, u2=e1, time step Δt=50μs, and the error tolerance er=10-5. 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.Figure 22
Dynamic Load history in the two case studies; (a) first case and (b) second case.
(a) (b)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 to4X105lb, as shown in Figure 22(b). Ten beam elements are used in the analysis, u2=e1, time step Δt=25μs, and the error tolerance er=10-5. 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.


## 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 Figure2. Five beam elements are used in the analysis, u2=e1, and the error tolerance er 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 PL2/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.
P L 2 E I 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.81471Table 2
Displacement errors for cantilever beam subjected to end vertical load, h=0.25m.
P L 2 E I 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.00410Figure 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 Figure5 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.
P L 2 E I 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.8927Figure 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 Figure6 is considered in this example. Two elements per member are used in the analysis, u2=e1, and the error tolerance er is chosen to be 10-7. 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 PL2/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.
P L 2 E I 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.24951Table 5
Displacement errors for pinned-fixed square diamond frame loaded in tension, h=0.1m.
P L 2 E I 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.00571Figure 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 Figure9 and Table 6.Table 6
Comparison of the present results for pinned-fixed square diamond frame loaded in tension, h = 2.8284 m.
P L 2 E I 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.2853Figure 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 Figure10 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 P=4000kips for n=0.1 and after P=3000kips for n=0.5 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 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, u2=e2, and the used error tolerance is er=10-6. As expected the Timoshenko model gives a higher deflection than Euler-Bernoulli model for the same load.


## 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 Figure2. Five beam elements are used in the analysis, u2=e1, and the error tolerance er 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 PL2/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.
P L 2 E I 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.81471Table 2
Displacement errors for cantilever beam subjected to end vertical load, h=0.25m.
P L 2 E I 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.00410Figure 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 Figure5 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.
P L 2 E I 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.8927Figure 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 Figure6 is considered in this example. Two elements per member are used in the analysis, u2=e1, and the error tolerance er is chosen to be 10-7. 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 PL2/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.
P L 2 E I 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.24951Table 5
Displacement errors for pinned-fixed square diamond frame loaded in tension, h=0.1m.
P L 2 E I 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.00571Figure 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 Figure9 and Table 6.Table 6
Comparison of the present results for pinned-fixed square diamond frame loaded in tension, h = 2.8284 m.
P L 2 E I 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.2853Figure 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 Figure10 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 P=4000kips for n=0.1 and after P=3000kips for n=0.5 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 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, u2=e2, and the used error tolerance is er=10-6. 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 loadP=Posin(⁡ωt) at the free end as shown in Figure 13. Young’s modulus E is considered 210 GPa, the material’s density ρ = 7,850 kg/m3, and Poisson’s ratio v = 0.3. Po=10MN and its frequency ω= 50 rad/s. Ten beam elements are used in the analysis, u2=e1, time step Δt=10-4s, and the error tolerance er=10-5. 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 loadP 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=106, EI=1,000, ρA=1, ρI=10, and v=0 and μ is the damping coefficient. Only four beam elements are used in the analysis, u2=e1, time step Δt=0.05, and the error tolerance er 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].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 Figure21 is considered in this example. Young’s modulus E is 30,000 ksi, the material’s density ρ is 2.536X10-4lbs2/in4, and Poisson’s ratio v 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 640lb as shown in Figure 22(a). Ten beam elements are used in the analysis, u2=e1, time step Δt=50μs, and the error tolerance er=10-5. 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.Figure 22
Dynamic Load history in the two case studies; (a) first case and (b) second case.
(a) (b)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 to4X105lb, as shown in Figure 22(b). Ten beam elements are used in the analysis, u2=e1, time step Δt=25μs, and the error tolerance er=10-5. 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.


## 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 loadP=Posin(⁡ωt) at the free end as shown in Figure 13. Young’s modulus E is considered 210 GPa, the material’s density ρ = 7,850 kg/m3, and Poisson’s ratio v = 0.3. Po=10MN and its frequency ω= 50 rad/s. Ten beam elements are used in the analysis, u2=e1, time step Δt=10-4s, and the error tolerance er=10-5. 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 loadP 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=106, EI=1,000, ρA=1, ρI=10, and v=0 and μ is the damping coefficient. Only four beam elements are used in the analysis, u2=e1, time step Δt=0.05, and the error tolerance er 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].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 Figure21 is considered in this example. Young’s modulus E is 30,000 ksi, the material’s density ρ is 2.536X10-4lbs2/in4, and Poisson’s ratio v 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 640lb as shown in Figure 22(a). Ten beam elements are used in the analysis, u2=e1, time step Δt=50μs, and the error tolerance er=10-5. 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.Figure 22
Dynamic Load history in the two case studies; (a) first case and (b) second case.
(a) (b)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 to4X105lb, as shown in Figure 22(b). Ten beam elements are used in the analysis, u2=e1, time step Δt=25μs, and the error tolerance er=10-5. 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.



---

*Source: 2670462-2018-08-14.xml*