# Static, Vibration, and Transient Dynamic Analyses by Beam Element with Adaptive Displacement Interpolation Functions

Authors: Shuang Li; Jinjun Hu; Changhai Zhai; Lili Xie
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

An approach to analyzing structures by using beam elements is developed with adaptive displacement interpolation functions. First, the element stiffness matrix and equivalent nodal loads are derived on the basis of the equilibrium between nodal forces and section forces rather than the compatibility between nodal deformations and section deformations, which avoids discretization errors caused by the limitation of conventional polynomial interpolation functions. Then, six adaptive element displacement interpolation functions are derived and extended to include several cases, such as beams with variable cross-section, variable material properties, and many different steps in cross-section and/or material properties. To make the element usable in dynamic analyses, consistent mass matrix (CMM) and diagonally lumped mass matrix (LMM) are constructed using the presented adaptive displacement interpolation functions. All these features have made the element elegant, which is tested with a number of simple static, vibration, and dynamic examples to show its accuracy.


---

## Body

## 1. Introduction

Beams with variable cross-section and/or material properties are frequently used in aeronautical engineering (e.g., rotor shafts and functionally graded beams), mechanical engineering (e.g., robot arms and crane booms), and civil engineering (e.g., beams, columns, and steel composite floor slabs in the single direction loading case). This paper tries to look for a simple but general approach to tackling the static and dynamic responses of such beams. A brief review of prior researches follows.The beam with variable cross-section is often modeled by a large number of small uniform elements, replacing the continuous changes with a step law. This scheme is accurate for a stepped beam but approximated for a beam with continuously changed cross-section. Although in this way, it is always possible to reduce errors as much as desired and obtain acceptable results by refining meshes, the modeling and computational efforts can become excessive. Sapountzakis and coworkers [1, 2] developed boundary element methods for static torsion and torsional vibration analyses of bars with variable cross-section. Banerjee et al. [3] used the dynamic stiffness method to investigate the free bending vibration of rotating beams with linearly changed cross-section. Gimena et al. [4] investigated the static responses of curved beam with variable cross-section, in which the stiffness matrix and the equivalent nodal loads of the curved beam element were presented. Carrera and coworkers [5, 6] derived the Carrera Unified Formulation, and under that framework, they presented a method to analyze beams with arbitrary cross-sectional geometries. Firouz-Abadi et al. [7] presented a Wentzel, Kramers, Brillouin approximation-based analytical solution to free transverse vibration of a class of varied cross-section beams. Ece et al. [8] presented solutions of the vibration of beams with exponentially varying cross-section width for three different types of boundary conditions associated with simply supported, clamped, and free ends. Shin et al. [9] applied the generalized differential quadrature method and differential transformation method to vibration analysis of circular arches with variable cross-section, stating that these two methods showed fast convergence and accuracy. Exact displacement interpolation functions to beams with linearly changed cross-section were solved and then used to derive the accurate stiffness matrix [10, 11]. The use of exact displacement interpolation functions to solve varied cross-section beam problems is a straightforward way; however, they [10, 11] only focused on the beam with linearly and continuously changed cross-section.A special case is the stepped beam, a beam with abrupt changes of cross-section and/or material properties. Several works about stepped beams had been published. Naguleswaran [12] presented an analytical approach to calculating the frequencies of beams on elastic end supports and with up to three step changes in cross-section. Maurini et al. [13] presented an enhancement of assumed mode method by introducing special jump functions to catch the curvature discontinuities of the mode shapes. Kisa and Gurel [14] presented a technique to solve the free vibration problems of stepped beam with circular cross-section and an existing crack. In [15, 16], solutions to the free vibration problem of stepped beams were presented by using the properties of Green’s function.Jaworski and Dowell [17] conducted an experiment of free vibration analysis of a stepped cantilevered beam and compared the experiment results with the classical Rayleigh-Ritz method, component modal analysis, and commercial finite element software ANSYS, the local boundary conditions and nonbeam effects were discussed. Lu et al. [18] used the composite element method to analyze free and forced vibrations of stepped beams and compared with experimental results. Mao and Pietrzko [19] used the Adomian decomposition method to investigate the free vibrations of a two-stepped beam, considering different boundary conditions, step locations, and step ratios. Most recently, Zheng and Ji [20] presented an equivalent representation of a stepped beam with a uniform beam to simplify the calculation of static deformations and frequencies. The methods to analyze the stepped beams and the beams with continuously changed cross-section are generally not unified.In addition to cross-sectional geometry, other parameters such as material properties (e.g., modulus, mass density, etc.) are also changeable, which belong to the researches on beams made of functionally graded materials. For functionally graded beams, continuous material gradient variation may be orientated in the cross-section and/or in the axial direction. For the former, there have been many researches devoted to static and vibration analyses (e.g., [21–25]). For the axially functionally graded beams, few solutions are yet found for arbitrary gradient changes. Alsorbagy et al. [26] studied the free vibration characteristics of a functionally graded beam with material graduation in the axial direction by finite element method, in this method fixed displacement interpolation functions were used, leading to only an approximate representation of the effect of material graduation. To the authors’ knowledge, the research done by Huang and Li [27] may be the only available publication on axially functionally graded beams coupled with variable cross-section. In their work, high-order differential equations were solved by converting them to Fredhplm integral equations. The axially functionally graded beam is somewhat similar to material nonlinear beams in that the effect of variation of material properties can be represented by integration along an element. However, the conventional displacement interpolation functions usually used in material nonlinear beam elements are not accurate for axially functionally graded beam.The solution to structural static and dynamic problems (e.g., eigenvalue problems, time integration methods, etc.) can be more accurately computed if certain vectors and matrixes (e.g., equivalent nodal loads and mass matrices) are well established. If there are internal loads applying on an element, the equivalent nodal load is usually obtained by integrating the product of displacement interpolation functions and the applying loads along the element. Hence, the process sensitively relies on the accuracy of the displacement interpolation functions. Two schemes for construction of mass matrix are well known and widely used, leading to consistent and diagonally lumped forms. The consistent ones can be used in frequency and common dynamic analyses, and the lumped ones are usually used in explicit finite element analyses especially for impact-like problems. Archer and Whalen [28] presented a new diagonalization technique for mass matrix of beam element which contains both translational and rotational degrees of freedom (DOFs). Fried and Leong [29] treated the mass matrix as a linear matrix function, M(γ)=M1+γM2, of parameter γ such that M(γ=0) is the consistent mass matrix (CMM) and M(γ=1) is the diagonally lumped mass matrix (LMM). And a linear combination of mass matrixes obtained by different methods was also investigated in [28, 30] to achieve better frequency analysis results. Felippa [31] presented a historical notes and comments on mass matrixes. It seems that there were fewer research works on CMM for dynamic response analyses and on LMM for both frequency and dynamic response analyses, which had been performed on varied cross-section beam and functionally graded beam.So far, numerical methods such as finite element, finite difference and differential quadrature, Rayleigh-Ritz method, and analytical methods which are based on fourth-order differential equations have been used to solve problems of beams with cross-sectional geometric or material variations. Among these methods, finite element may be a better choice, because it is easy to treat boundary conditions and implement the formulas to an existed finite element framework resulting to a potential to solve large and complex problems. In the present study, a simple and unified approach is derived using Euler-Bernoulli beam theory to analyze beams with variable cross-section, variable material properties, and many different steps in cross-section and/or material properties. The paper is organized as follows. The defects of the fixed displacement interpolation functions are presented in Section2. The element stiffness matrix and equivalent nodal loads are studied in Section 3. The adaptive element displacement interpolation functions are derived in Section 4. The new element consistent and diagonally lumped mass matrixes are presented in Section 5. Numerical examples are presented to assess the performance of the presented element in Section 6. Some concluding remarks complete the study.


## 2. Defects of Fixed Element Interpolation Functions

Considering a planar Euler-Bernoulli beam with internal loads, its governing differential equations are given by(2.1)-∂∂x[E(x)A(x)∂u(x)∂x]=qx(x),∂2∂x2[E(x)I(x)∂2v(x)∂x2]=qy(x),
where x is the element reference axis connecting two element nodes; qx(x) and qy(x) are the axial and transversal internal element loads, respectively; E(x) is the material modulus; A(x) is the cross-sectional area; I(x) is the moment of inertia of cross-section; u(x) and v(x) are the axial and transversal displacements of the element, respectively.The most commonly used displacement interpolation functions for beam elements are as follows, which can be found in many textbooks [32, 33]:(2.2)N1(x)=1-xL,N2(x)=1L3(L3-3Lx2+2x3),N3(x)=1L2(L2x-2Lx2+x3),N4(x)=xL,N5(x)=1L3(3Lx2-2x3),N6(x)=1L2(x3-Lx2),
where L is the element length. These functions are linear Lagrangian polynomials and three-order Hermitian polynomials for axial and transversal displacements, respectively. The functions predefine a law that governs the variation pattern of displacements in an element. Elements adopt (2.2) all belong to the catalog of so-called three-order beam element.The functions expressed by (2.2) are fixed in shape and lead to the cross-sectional geometry and material independent displacement interpolation functions. It is well known that (2.2) are derived by (2.1) under the conditions: (1)  E(x) is constant; (2)  A(x) is constant; (3)  I(x) is constant; and (4)  qx(x) and qy(x) equal to zero. If any of these conditions is not met, (2.2) will not give accurate representations of the relationships between nodal and internal displacements. The discretization error, characterized by discrepancies between the displacements of the element and the real member, is introduced due to utilization of the inaccurate displacement patterns to approximate the real ones. Higher-order beam elements [34, 35] may be employed to improve the accuracy; however, in most cases, the discretization error also exists, computational effort increases, and the formula become more complex. Figure 1 is an example to illustrate the defects of the fixed displacement interpolation functions, in which a cantilever two-stepped beam with material modulus of E1 and E2 is subject to nodal loads. One element is used to model the beam, and the element stiffness matrix is obtained by integrating along the element to consider the effect of two stepped material properties. Considering the element state determination procedure in which the nodal displacements are obtained after a structural analysis step, we want to solve the element internal deformation field (i.e., derivatives of displacement field). The results, calculated by the known nodal displacements and (2.2), show large disparities between calculated strain and curvature (εc,κc) and real ones (εr,κr). Equation (2.2), actually, can only rationally represent the deformation with constant axial strain and linear curvature. The displacement interpolation functions are important for beam elements to derive the element equivalent nodal loads and consistent mass matrix.Comparisons of axial strain and curvature distributions. (a) Axial strain under nodal forceN, (b) Curvature under nodal force M.
(a)(b)


## 3. Element Stiffness Matrix and Equivalent Nodal Loads

### 3.1. Equilibrium Relationship in an Element

The adopted coordinate frame is presented in Figure2. It is a generalized coordinate system so that only three deformation DOFs (i.e., axial elongation and rotation deformations at the two elemental nodes) and their corresponding nodal forces are considered for the planar beam element. ox(I→J) is the axis connected geometric centroids of sections at nodes I and J. Within the element equilibrium framework, equilibrium is stated in the form(3.1)N(x)=Q1,M(x)=(xL-1)Q2+xLQ3,
where  Q=[𝒬1𝒬2𝒬3]⊤ are the nodal forces, in which the superscript ⊤ denotes the transpose; N(x) and M(x) are the section axial force and moment with respect to coordinate x.Figure 2
Coordinate frame and element forces.Because the equilibrium relationship presented by ((3.1)) is unchanged when the cross-sectional geometry and material change, the variation pattern of forces in the element is accurate. The relationship in ((3.1)) is employed to replace the displacement interpolation functions as discretization functions, hence the discretization error is avoided. For example, reconsidering the problem in Figure 1, Figure 3 shows that the calculated axial force and moment (Nc,Mc) by known nodal forces and (3.1) are same as real ones (Nr,Mr). The discretization approach directly predefines element equilibrium relationship within the generalized coordinate system rather than compatibility relationship by (2.2). Rewriting the section forces in a vector form as S(x)=[N(x)M(x)]⊤, we can get the equilibrium matrix(3.2)e(x)=[1000xL-1xL],
which connects the element nodal forces and section forces as S(x)=e(x)Q.Comparisons of axial force and moment distributions. (a) Axial forces under nodal forceN, (b) Moment under nodal force M.
(a)(b)


### 3.2. Element Stiffness Matrix

The section deformations of the beam element are(3.3)d(x)=[ε(x)κ(x)]⊤=[∂u(x)∂x∂2v(x)∂x2]⊤,
where ε(x) is the axial strain with respect to coordinate x; κ(x) is the section curvature with respect to coordinate x.The weak integration form of (3.3) along the beam element is(3.4)∫0L[δN(x)(∂u(x)∂x-ε(x))+δM(x)(∂2v(x)∂x2-κ(x))]dx=0,
where δ stands the virtual variation.Integration of two derivatives in (3.4) by parts, yielding(3.5)∫0LδN(x)∂u(x)∂xdx=δN(x)u(x)|0L-∫0L∂δN(x)∂xu(x)dx,(3.6)∫0LδM(x)∂2v(x)∂x2dx=δM(x)∂v(x)∂x|0L-∫0L∂δM(x)∂x∂v(x)∂xdx.
Furthermore, we repeat to integrate the second term at the right side of (3.6), yielding(3.7)∫0L∂δM(x)∂x∂v(x)∂xdx=∂δM(x)∂xv(x)|0L-∫0L∂2δM(x)∂x2v(x)dx.As shown in Figure2, there are three nodal forces in the beam generalized coordinate system. And corresponding deformation DOFs are the axial elongation and rotation deformations at the two elemental nodes. Therefore, the pattern of nodal forces and the deformations of a beam in the generalized coordinate system are the same as in a simply supported beam system. Thus, we have boundary conditions(3.8)u(0)=v(0)=v(L)=0,u(L)=D1,∂v(0)∂x=D2,∂v(L)∂x=D3,δN(L)=δQ1,δM(0)=-δQ2,δM(L)=δQ3,
where D=[D1D2D3]⊤ are generalized nodal displacements (i.e., axial elongation and rotations at nodes I and J); δ stands the virtual variation.Substituting ((3.5), (3.6)), (3.7), and (3.8) into (3.4) leads to(3.9)∫0L[δN(x)ε(x)+δM(x)κ(x)]dx+∫0L[∂δN(x)∂xu(x)-∂2δM(x)∂x2v(x)]dx=δQ1D1+δQ2D2+δQ3D3.According to the principle of virtual force,u(x) and v(x) are impossibly equal to zero with respect to all x, and one gets(3.10)∂δN(x)∂x=0,∂2δM(x)∂x2=0.Integration of (3.10) and application of the boundary conditions in (3.8), rewritten in vector form, the relationship between section force increments and generalized nodal force increments can be expressed as(3.11)δS(x)=e(x)δQ.
Substituting (3.10) and (3.11) into (3.9) and considering the arbitrary of δQ, we have(3.12)D=∫0Le(x)⊤d(x)dx.The element flexibility matrix in the generalized coordinate system can be obtained by the common finite element formula that differentiatesD with respect to Q. The stiffness matrix is then represented in the inversion form of the flexibility matrix F(3.13)K=F-1=(∂D∂Q)-1.The explicit form of element stiffness matrix in the generalized coordinate system is expressed as follows:(3.14)K=[1a1000a2a2a3-a42-a4a2a3-a420-a4a2a3-a42a2a2a3-a42],
where(3.15)a1=∫0L1E(x)A(x)dx,a2=∫0L1E(x)I(x)(xL)2dx,a3=∫0L1E(x)I(x)(1-xL)2dx,a4=-∫0L1E(x)I(x)xL(1-xL)dx.Equation (3.15) can be directly integrated to get analytical solutions or be calculated by numerical integration approaches such as Gauss quadrature formula. For the cases of beams with many different steps in cross-section and/or material properties, the integration needs to be processed in each piecewise section and then be summed up. Note that (3.15) considers variable parameters of cross-sectional geometry and material so that K is accurate for these cases if the error of numerical integration is neglected. K will be transformed to a 6 × 6 matrix in the local element coordinate system by transform matrixes before used in the finite element calculated process.


### 3.3. Internal Element Loads and Equivalent Nodal Loads

The equilibrium relationshipS(x)=e(x)Q, which is related to (3.2), is corresponding to the case without internal element loads and can be rewritten in the following form when there are internal element loads:(3.16)S(x)=e(x)Q+e̅(x)q,
where q is the internal element loads; e̅(x)q is the additional section forces caused by q. The e̅(x) and q for several common cases are given in Table 1.Table 1
Theē(x) matrixes and q for commonly used load patterns.
Load patterne̅(x)qe̅(x)=[100-(L-a)xL],x<aq=[FxFy]⊤e̅(x)=[000-a(L-x)L],x>ae̅(x)=[000xL],x<aq=[0M]⊤e̅(x)=[000-L-xL],x>ae̅(x)=[b-a00-(b-a)(2L-a-b)x2L],x<ae̅(x)=[b-x00a2(L-x)+x(b2-2Lb+Lx)2L],a<x<bq=[qxqy]⊤e̅(x)=[000-(b2-a2)(L-x)2L],x>be̅(x)=[000(b-a)xL],x<ae̅(x)=[000a(L-x)-(L-b)xL],a<x<bq=[0m]⊤e̅(x)=[000-(b-a)(L-x)L],x>bThe second term in (3.16) leads to additional section deformations with f(x)e̅(x)q. Then, the additional nodal deformations can be obtained by principle of virtual force as(3.17)D¯=∫0Le(x)⊤f(x)e̅(x)qdx,
where f(x) is the section flexibility matrix. Multiplying K by  D¯, we have the element equivalent nodal loads as follows:(3.18)Q¯=K∫0Le(x)⊤f(x)e̅(x)qdx.


## 3.1. Equilibrium Relationship in an Element

The adopted coordinate frame is presented in Figure2. It is a generalized coordinate system so that only three deformation DOFs (i.e., axial elongation and rotation deformations at the two elemental nodes) and their corresponding nodal forces are considered for the planar beam element. ox(I→J) is the axis connected geometric centroids of sections at nodes I and J. Within the element equilibrium framework, equilibrium is stated in the form(3.1)N(x)=Q1,M(x)=(xL-1)Q2+xLQ3,
where  Q=[𝒬1𝒬2𝒬3]⊤ are the nodal forces, in which the superscript ⊤ denotes the transpose; N(x) and M(x) are the section axial force and moment with respect to coordinate x.Figure 2
Coordinate frame and element forces.Because the equilibrium relationship presented by ((3.1)) is unchanged when the cross-sectional geometry and material change, the variation pattern of forces in the element is accurate. The relationship in ((3.1)) is employed to replace the displacement interpolation functions as discretization functions, hence the discretization error is avoided. For example, reconsidering the problem in Figure 1, Figure 3 shows that the calculated axial force and moment (Nc,Mc) by known nodal forces and (3.1) are same as real ones (Nr,Mr). The discretization approach directly predefines element equilibrium relationship within the generalized coordinate system rather than compatibility relationship by (2.2). Rewriting the section forces in a vector form as S(x)=[N(x)M(x)]⊤, we can get the equilibrium matrix(3.2)e(x)=[1000xL-1xL],
which connects the element nodal forces and section forces as S(x)=e(x)Q.Comparisons of axial force and moment distributions. (a) Axial forces under nodal forceN, (b) Moment under nodal force M.
(a)(b)


## 3.2. Element Stiffness Matrix

The section deformations of the beam element are(3.3)d(x)=[ε(x)κ(x)]⊤=[∂u(x)∂x∂2v(x)∂x2]⊤,
where ε(x) is the axial strain with respect to coordinate x; κ(x) is the section curvature with respect to coordinate x.The weak integration form of (3.3) along the beam element is(3.4)∫0L[δN(x)(∂u(x)∂x-ε(x))+δM(x)(∂2v(x)∂x2-κ(x))]dx=0,
where δ stands the virtual variation.Integration of two derivatives in (3.4) by parts, yielding(3.5)∫0LδN(x)∂u(x)∂xdx=δN(x)u(x)|0L-∫0L∂δN(x)∂xu(x)dx,(3.6)∫0LδM(x)∂2v(x)∂x2dx=δM(x)∂v(x)∂x|0L-∫0L∂δM(x)∂x∂v(x)∂xdx.
Furthermore, we repeat to integrate the second term at the right side of (3.6), yielding(3.7)∫0L∂δM(x)∂x∂v(x)∂xdx=∂δM(x)∂xv(x)|0L-∫0L∂2δM(x)∂x2v(x)dx.As shown in Figure2, there are three nodal forces in the beam generalized coordinate system. And corresponding deformation DOFs are the axial elongation and rotation deformations at the two elemental nodes. Therefore, the pattern of nodal forces and the deformations of a beam in the generalized coordinate system are the same as in a simply supported beam system. Thus, we have boundary conditions(3.8)u(0)=v(0)=v(L)=0,u(L)=D1,∂v(0)∂x=D2,∂v(L)∂x=D3,δN(L)=δQ1,δM(0)=-δQ2,δM(L)=δQ3,
where D=[D1D2D3]⊤ are generalized nodal displacements (i.e., axial elongation and rotations at nodes I and J); δ stands the virtual variation.Substituting ((3.5), (3.6)), (3.7), and (3.8) into (3.4) leads to(3.9)∫0L[δN(x)ε(x)+δM(x)κ(x)]dx+∫0L[∂δN(x)∂xu(x)-∂2δM(x)∂x2v(x)]dx=δQ1D1+δQ2D2+δQ3D3.According to the principle of virtual force,u(x) and v(x) are impossibly equal to zero with respect to all x, and one gets(3.10)∂δN(x)∂x=0,∂2δM(x)∂x2=0.Integration of (3.10) and application of the boundary conditions in (3.8), rewritten in vector form, the relationship between section force increments and generalized nodal force increments can be expressed as(3.11)δS(x)=e(x)δQ.
Substituting (3.10) and (3.11) into (3.9) and considering the arbitrary of δQ, we have(3.12)D=∫0Le(x)⊤d(x)dx.The element flexibility matrix in the generalized coordinate system can be obtained by the common finite element formula that differentiatesD with respect to Q. The stiffness matrix is then represented in the inversion form of the flexibility matrix F(3.13)K=F-1=(∂D∂Q)-1.The explicit form of element stiffness matrix in the generalized coordinate system is expressed as follows:(3.14)K=[1a1000a2a2a3-a42-a4a2a3-a420-a4a2a3-a42a2a2a3-a42],
where(3.15)a1=∫0L1E(x)A(x)dx,a2=∫0L1E(x)I(x)(xL)2dx,a3=∫0L1E(x)I(x)(1-xL)2dx,a4=-∫0L1E(x)I(x)xL(1-xL)dx.Equation (3.15) can be directly integrated to get analytical solutions or be calculated by numerical integration approaches such as Gauss quadrature formula. For the cases of beams with many different steps in cross-section and/or material properties, the integration needs to be processed in each piecewise section and then be summed up. Note that (3.15) considers variable parameters of cross-sectional geometry and material so that K is accurate for these cases if the error of numerical integration is neglected. K will be transformed to a 6 × 6 matrix in the local element coordinate system by transform matrixes before used in the finite element calculated process.


## 3.3. Internal Element Loads and Equivalent Nodal Loads

The equilibrium relationshipS(x)=e(x)Q, which is related to (3.2), is corresponding to the case without internal element loads and can be rewritten in the following form when there are internal element loads:(3.16)S(x)=e(x)Q+e̅(x)q,
where q is the internal element loads; e̅(x)q is the additional section forces caused by q. The e̅(x) and q for several common cases are given in Table 1.Table 1
Theē(x) matrixes and q for commonly used load patterns.
Load patterne̅(x)qe̅(x)=[100-(L-a)xL],x<aq=[FxFy]⊤e̅(x)=[000-a(L-x)L],x>ae̅(x)=[000xL],x<aq=[0M]⊤e̅(x)=[000-L-xL],x>ae̅(x)=[b-a00-(b-a)(2L-a-b)x2L],x<ae̅(x)=[b-x00a2(L-x)+x(b2-2Lb+Lx)2L],a<x<bq=[qxqy]⊤e̅(x)=[000-(b2-a2)(L-x)2L],x>be̅(x)=[000(b-a)xL],x<ae̅(x)=[000a(L-x)-(L-b)xL],a<x<bq=[0m]⊤e̅(x)=[000-(b-a)(L-x)L],x>bThe second term in (3.16) leads to additional section deformations with f(x)e̅(x)q. Then, the additional nodal deformations can be obtained by principle of virtual force as(3.17)D¯=∫0Le(x)⊤f(x)e̅(x)qdx,
where f(x) is the section flexibility matrix. Multiplying K by  D¯, we have the element equivalent nodal loads as follows:(3.18)Q¯=K∫0Le(x)⊤f(x)e̅(x)qdx.


## 4. Adaptive Element Displacement Interpolation Functions

Keep in mind that boundary conditions of the beam element in the generalized coordinate system, which expressed by (3.8), are the same as those of a simply supported beam. The internal elemental axis displacement u(x) and transverse displacement v(x) in the generalized coordinate system can be obtained by dummy unit force method [36]. A beam with an axial load F̅x and transversal load F̅y is shown in Figure 4(a). When using the dummy unit force method, F̅x and F̅y are set with one to calculate the section forces. For giving a more clear description to the element deformation pattern, simply supported beam like supports are also plotted in Figure 4 and later in Figure 5.Beam virtual and real deformation states. (a) With unit virtual loads, (b) With real nodal forces.
(a)(b)Figure 5
Transform between generalized and local coordinate system.Considering the equilibrium relationship, we can get the section axial forceN̅(x0,x) and moment M̅(x0,x) with respect to coordinate x0 in Figure 4(a):(4.1)N̅(x0,x)=H(x-x0),(4.2)M̅(x0,x)=(-1+xL)x0H(x-x0)+(-1+x0L)xH(-x+x0),
where H(·) is the Heaviside step function. Because H(0)=1, (4.2) has a discontinuity at x0=x, equals to two times as its real value. So, the definition of M̅(x0,x)|x0=x=M̅(x,x)/2 is employed to get a rational value at the point x0=x. This definition is only to get a continuous expression of  M̅(x0,x) and has no influences on the following integrations.Consider the beam with nodal forcesQ=[𝒬1𝒬2𝒬3] shown in Figure 4(b). According to the dummy unit force method, the element axial and transversal displacements with respect to coordinate x are calculated by the following formulas:(4.3)u(x)=∫0LN̅(x0,x)N(x0)E(x0)A(x0)dx0,v(x)=∫0LM̅(x0,x)M(x0)E(x0)I(x0)dx0,
where N(x0) and M(x0) are the section axial force and moment with respect to coordinate x0 under real nodal loads in Figure 4(b).ReplaceQ by K·D and consider ((3.1)) and (4.3), yielding(4.4)u(x)=[∫0LN̅(x0,x)e1(x0)E(x0)A(x0)dx0K]D,v(x)=[∫0LM̅(x0,x)e2(x0)E(x0)I(x0)dx0K]D,
where e1(x0) and e2(x0) are the first and second rows of (3.2), respectively.The element axial and transversal displacements can also be defined with respect to local element coordinate system. Assuming that the displacement and rotation between the generalized and local coordinate systems are small, according to Figure5, the geometric transform relationship can be obtained as
(4.5)ue(x)=[∫0LN̅(x0,x)e1(x0)E(x0)A(x0)dx0K⋅T]De+D1e,ve(x)=[∫0LM̅(x0,x)e2(x0)E(x0)I(x0)dx0K⋅T]De+D5e-D2eLx+D2e,
where De=[D1eD2eD3eD4eD5eD6e]⊤ is the nodal displacements in the local element coordinate system, in which Die are the axial, transversal, and rotation displacements at nodes I and J; T·De=D, where T is a transform matrix as follows. And T is also the matrix that used to transform (3.14) and (3.18) into the local element coordinate system(4.6)T=[-10010001L10-1L001L00-1L1].Equation (4.5) can be rewritten in the following form:(4.7)ue(x)=N1(x)D1e+N4(x)D4e,ve(x)=N2(x)D2e+N3(x)D3e+N5(x)D5e+N6(x)D6e,
where(4.8)N1(x)=1-K(1,1)γ1(x),N2(x)=1-xL+[K(2,2)+K(2,3)]γ2(x)+[K(3,2)+K(3,3)]γ3(x),N3(x)=K(2,2)γ2(x)+K(3,2)γ3(x),N4(x)=1-N1(x),N5(x)=1-N2(x),N6(x)=K(2,3)γ2(x)+K(3,3)γ3(x),(4.9)γ1(x)=∫0LH(x-x0)E(x0)A(x0)dx0,γ2(x)=∫0L(1-x0/L)[x0/L(1-x/L)H(x-x0)+x/L(1-x0/L)H(-x+x0)]E(x0)I(x0)dx0,γ3(x)=∫0L-x0/L[x0/L(1-x/L)H(x-x0)+x/L(1-x0/L)H(-x+x0)]E(x0)I(x0)dx0.It can be easily seen that the six functions in (4.8), which present the relationships between element nodal and internal displacements, have functions similar to the element displacement interpolation functions shown in (2.2). However, (4.8) has some superior performance. Compared with the fixed functions, (4.8) are adaptive element displacement interpolation functions. The functions are accurate for beams with variable cross-section, variable material properties, and many different steps in cross-section and/or material properties.


## 5. Consistent and Lumped Mass Matrixes

Here, we use the presented adaptive displacement interpolation functions to derive a CMM that can cope with the beams with variable cross-sectional geometry and material properties. The element mass matrix in local element coordinate system is given by the common finite element formula(5.1)M=∫0Lρ(x)A(x)N(x)N(x)⊤dx,
where ρ(x) is the material mass density; N(x) is the interpolation function matrix as follows, in which N1(x) to N6(x) are expressed by (4.8):(5.2)N(x)=[N1(x)00N4(x)000N2(x)N3(x)0N5(x)N6(x)]⊤.It is worth noting the equivalent nodal loads in local element coordinate system can also be derived by the adaptive displacement interpolation functions. The expression given by (5.3) is the same as the result presented by (3.18) when they are transformed to the local element coordinate system(5.3)Q̅=∫0LN(x)qdx+∑mN(xm)Fm+∑ndN(xn)dxnMn,
where q is the element axial and transversal distribution loads; Fm is the element axial and transversal concentrated loads with respect to coordinate xm; Mn is the concentrated moment with respect to coordinate xn. The q, Fm and Mn are expressed as(5.4)q=[qxqy]⊤,Fm=[FmxFmy]⊤,Mn=[0Mn]⊤.The explicit finite element method has been extensively developed for dynamic analyses to meet the special demands of engineering applications. In this case, the LMM is preferable. The LMM can be created from the CMM. One suggestion, called row sum approach, is to calculate the LMM from diagonalization by adding the off-diagonal entries in each row to diagonal entry and setting the rotational mass terms with zero [28]. An alternative, called diagonal scaling approach, is to ignore the off-diagonal terms of consistent mass matrix and scale the diagonal entries to satisfy the conservation requirement [33]. The third approach uses a numerical integration to obtain a diagonal matrix in which the integration points coincide with the element nodes [33]. These three approaches make sense when all DOFs of the problem have the same physical interpretation; that is, the element has only translational DOFs. Assuming that the total rotational inertia of the diagonal rotational mass terms is same the as the consistent one, [28, 34] developed rotationally consistent lumped mass matrixes for beam element in which the rotational mass terms are -MbeamL2/24. In [28], a value of -11MbeamL2/156, which is the average of -MbeamL2/24 and the rotational mass term obtained from diagonal scaling approach, had also been used to get more accurate vibration frequencies. The drawback of mass matrixes in [28, 34] is the introduction of negative mass entries to rotational DOFs and negative definite mass matrix. In practice, shifting and condensation approaches have to be usually used to overcome the numerical convergence problem. In this study, a simple approach is presented which uses the row sum approach, and the rotational mass terms are calculated by scale factors determined for the conservation requirement of translational mass. The presented LMM  is
(5.5)M=[m1,1+m1,4000000m2,2+m2,5000000α1m3,3000000m4,1+m4,4000000m5,2+m5,5000000α2m6,6],
where mi,j is the elements in (5.1) at the position i, j. α1 and α2 are scale factors which can be calculated by the diagonal scaling approach according to mass conservation [28]. Here, we modify them as α1=α2=1. Some verifications had shown that the modification can be a little more accurate than directly using the factors calculated by the diagonal scaling approach.Equation (5.5) is rational for a uniform beam element and gives good results for frequency and dynamic analyses. However, for nonuniform beam the two translational terms at one node are not equal, leading to a nondiagonal mass matrix when transforming it to global coordinate system. Hence, a scheme to average the two translational terms at one node is used to form a new LMM as follows: (5.6)M=[m̅1,1000000m̅2,2000000m3,3000000m̅4,4000000m̅5,5000000m6,6],
where
(5.7)m̅1,1=m̅2,2=m1,1+m1,4+m2,2+m2,52,m̅4,4=m̅5,5=m4,1+m4,4+m5,2+m5,52.The presented LMM can also consider the variation of cross-sectional geometry and material properties, because it is created from the presented LMM. Furthermore, the matrix has better convergence performance than commonly used ones which set the rotational mass terms with zero or negative values.


## 6. Numerical Results

The objective of the numerical applications in this section is to assess the performance of the presented element. The examples illustrate superior performance of the presented element than conventional beam elements in several cases.


### 6.1. Linearly Tapered Cantilever Beam

A cantilever beam shown in Figure6 is taken, with section width of 0.5 m and the section height varies as h(x)=h0-2x/25, where h0=1 m is the section height at the fixed end. The material modulus E is 210 GPa. The mass density ρ is 7.8 × 103 kg/m3. ANSYS beam54 element [37] (belongs to three-order beam element) is used for comparison. The beam54 element can model both beam with uniform and variable cross-section (denoted as ANSYS-1 and ANSYS-2 in Table 2). Figure 7 shows the comparison of the displacement interpolation functions between one presented element (4.8) and one three-order beam element (2.2). As seen in Figure 7, the presented displacement interpolation functions give a more accurate representation to the deformations of the beam (e.g., considering rotation DOFs, the section height of the beam is larger at the left end and unit rotation results in larger beam deformation. However, three-order beam element gives the same deformations at the two ends). Figure 8 shows the comparison of CMM and LMM between one presented element and one three-order beam element. The LMM of the three-order beam element is calculated by the common method in the textbook [32, 33], in which  LMM=Mbeam/2×diag[110110]. Similar to the observation from comparison of displacement interpolation functions, the presented mass matrixes can also represent the effect of the varied cross-section. The calculated equivalent nodal loads by one presented element and one three-order beam element are [0130735870693-518] and [01000166701000-1667], respectively, (Units: kN, kN/m). It is obvious that the three-order beam element gives irrational element equivalent nodal loads, because the loads are independent of cross-sectional geometry and material which have the same values at the two nodes. Table 2 shows the static displacement, rotation, and frequency results. It can be seen that to achieve comparable accuracy, many more ANSYS elements are required than the presented beam elements. It is seen that when 5 elements are used to calculate the frequency, the presented LMM only has an error of 1.86% relative to the presented CMM. For the dynamic analysis, a 0.1 s (about 1/10 of 1st natural period) pulse loads with amplitudes of the equivalent nodal loads are applied to the free end. Figure 9 illustrates that the presented CMM and LMM are more accurate than the conventional method. Only the results of ANSYS-2 and the presented element are plotted in Figure 9 for brevity.Table 2
The tip displacement, tip rotation, and 1st natural vibration frequency.
Tip disp. (×10−3m)Tip rotation (×10−3rad)1st-freq. (Hz)Ele.This paperANSYS-1ANSYS-2This paperANSYS-1ANSYS-2PresentedCMMANSYS-1,CMMANSYS-2,CMMPresentedLMMANSYS-1,LMMANSYS-2,LMM160.92132.28×10.9317.64×10.6395.049×6.2592.885×465.4764.2011.8611.5110.2389.58810.1789.9299.0289.110563.7863.0111.5311.3210.2359.81110.19610.0459.4439.4941061.6161.4411.0711.0310.11910.2252061.0961.0510.9610.9510.19810.2334060.9660.9510.9310.9310.21810.2358060.9360.9310.22312060.9360.9210.224Note: ×, If 1–3 elements are used for the ANSYS beam54 with variable cross-section (ANSYS-2), ANSYS reports error information of losing accuracy for the sharp change of the beam cross-section.Figure 6
Cantilever beam under the distributed load.Figure 7
Comparison of the displacement interpolation functions between one presented element and one three-order beam element.Comparison of the mass matrixes between one presented element and one three-order beam element. (a) CMM of this paper, (b) CMM of three-order beam, (c) LMM of this paper, and (d) LMM of three-order beam.
(a)(b)(c)(d)Figure 9
Dynamic responses of the beam at right end (results for 1 element of ANSYS-2 are not presented for the reason stated in Table2). (a) Nodal displacement in y direction using CMM, (b) Nodal displacement in y direction using LMM, (c) Nodal rotation using CMM, and (d) Nodal rotation using LMM.


### 6.2. Stepped Beam with Multiple Sectional Geometry

In [17], a free vibration test of an aluminum stepped beam (ρ=2664kg/m3,E=60.6GPa) with dimensions of 463.55 mm × 25.4 mm × 3.175 mm was conducted. The detailed dimensions are shown in Figure 10. The experimental (EXP.) frequencies in x-y plane (in-plane) and x-z plane (out-of-plane) are given in Table 3 together with the calculated solutions, including the Rayleigh-Ritz (RR) method and component modal analysis (CMA), analytical method which is based on the Euler-Bernoulli beam model, FEM with the Euler beam model, Timoshenko beam model, 2D shell, 3D solid model, and the presented beam model. The results of the presented beam element are obtained by one element and the results of beam elements in [17] are obtained by 500 elements. The error 10.8% is due to the nonbeam effects [17] of the test beam in its x-y plane. However, from the results in Table 3, it can be observed that using only one presented element can give acceptable accurate results compared with other beam elements. Another advantage of the presented element is that it can model the stepped beam with only one element rather than ANSYS which cannot model the stepped beam until 13 elements are used.Table 3
Comparison of the frequencies obtained from different methods.
MethodJaworski and Dowell [17]Zheng and  Ji [20]This paperEXP.RRCMAFE/ANSYSANALYTICALFEEuler beamEuler beamEuler beamTimoshenko beam2D shell3D solidEulerbeamEuler beam1st-freq. (in-plane)49.3854.80(+11.0%)54.99(+11.4%)54.47(+10.3%)54.43(+10.2%)49.62(+0.5%)49.8(+0.9%)50.96 (+3.2%)*54.71 (+10.8%)1st-freq. (out-of-plane)10.6310.75 (+1.1%)10.82(+1.8%)10.75(+1.1%)10.75(+1.1)10.44(−1.8%)10.4(−1.6%)—*10.77 (+1.3%)Note: * 54.71 and 10.77 are obtained by 1 presented element. In [17], the results of beam model are obtained by 500 elements. All results are obtained using CMM.Figure 10
Stepped aluminum beam test in [17]. (Unit: mm).


### 6.3. Functionally Graded Beam with Continuous Variable Material and Density

An axially functionally graded cantilever beam shown in Figure11 (ρ=5.7×103×e5(1-x)kg/m3,E=200×e5(1-x)GPa) with dimensions of 1 m × 0.04 m × 0.04 m is analyzed. The ANSYS beam 54 element [36] is used for comparison. Figure 12 shows the comparison of the displacement interpolation functions between one presented element and one three-order beam element. Figure 13 shows the comparison of CMM and LMM between one presented element and one three-order beam element. Similar to the example in Section 6.1, the Figures 12 and 13 illustrate that the presented element can represent the effects of the variable material modulus and mass density. Table 4 shows the static displacement, rotation, and frequency results. It can be seen that many ANSYS elements are required to achieve comparable accuracy, while only a few of presented beam elements are needed. The presented LMM only has an error of 2.73% relative to the presented CMM to calculate the frequency when 5 elements are used. The results of frequency by the presented element in Table 4, and also in Table 2, demonstrate some known knowledge observed by previous researches; for example, the rate is same of using CMM and LMM to achieve a nearly convergence solution. An observation in the paper is that although the same element meshes are used to achieve a nearly convergence frequency solution, the CMM is more accurate when less element meshes are used, because its convergence gradient is smaller than the LMM. For the dynamic analysis, a 0.007 s (about 1/10 of 1st natural period) pulse load with an amplitude of the nodal load is applied to the free end. Figure 14 illustrates that the presented CMM and LMM are more accurate than the conventional method. It should be noted that the ANSYS software gives a poor modeling to this example. The CMM of each ANSYS beam element mesh is calculated by the software itself, in which the midpoint material modulus and mass density of the each mesh are inputted for the calculation. This means that constant material modulus and mass density are used for each element and the exponential variations of them within the element cannot be considered by the software. The LMM of each ANSYS beam element mesh is calculated by Mbeam/2×diag[110110] and then inputted into the software as nodal lumped masses, in which Mbeam is manually calculated to consider the variation of material mass density within the element. These lead to that, as shown in Table 4, ANSYS presents inaccurate static responses and even lower accuracy of the frequencies obtained using CMM (compared with LMM) when a small number of element meshes are used.Table 4
The tip displacement, tip rotation, and 1st natural vibration frequency.
Ele.Tip disp. (×10−3m)Tip rotation (×10−3rad)1st-freq. (Hz)This paperANSYSThis paperANSYSPresented CMMANSYSCMMPresentedLMMANSYSLMM116.4132.0344.9848.02153.1117.1681.2538.45222.0352.25142.9564.88120.3894.04318.8749.32142.1895.92131.68117.28417.7747.67142.04112.62136.03127.25517.2746.77142.01122.00138.13132.281016.6445.42136.622016.4845.05140.734016.4544.96141.798016.44142.0512016.43Figure 11
Functionally graded beam under the nodal load.Figure 12
Comparison of the displacement interpolation functions between one presented element and one three-order beam element.Comparison of the mass matrixes between one presented element and one three-order beam element. (a) CMM of this paper, (b) CMM of three-order beam, (c) LMM of this paper, and (d) LMM of three-order beam.
(a)(b)(c)(d)Figure 14
Dynamic responses of the beam at right end. (a) Nodal displacement iny direction using CMM, (b) Nodal displacement in y direction using LMM, (c) Nodal rotation using CMM, and (d) Nodal rotation using LMM.


## 6.1. Linearly Tapered Cantilever Beam

A cantilever beam shown in Figure6 is taken, with section width of 0.5 m and the section height varies as h(x)=h0-2x/25, where h0=1 m is the section height at the fixed end. The material modulus E is 210 GPa. The mass density ρ is 7.8 × 103 kg/m3. ANSYS beam54 element [37] (belongs to three-order beam element) is used for comparison. The beam54 element can model both beam with uniform and variable cross-section (denoted as ANSYS-1 and ANSYS-2 in Table 2). Figure 7 shows the comparison of the displacement interpolation functions between one presented element (4.8) and one three-order beam element (2.2). As seen in Figure 7, the presented displacement interpolation functions give a more accurate representation to the deformations of the beam (e.g., considering rotation DOFs, the section height of the beam is larger at the left end and unit rotation results in larger beam deformation. However, three-order beam element gives the same deformations at the two ends). Figure 8 shows the comparison of CMM and LMM between one presented element and one three-order beam element. The LMM of the three-order beam element is calculated by the common method in the textbook [32, 33], in which  LMM=Mbeam/2×diag[110110]. Similar to the observation from comparison of displacement interpolation functions, the presented mass matrixes can also represent the effect of the varied cross-section. The calculated equivalent nodal loads by one presented element and one three-order beam element are [0130735870693-518] and [01000166701000-1667], respectively, (Units: kN, kN/m). It is obvious that the three-order beam element gives irrational element equivalent nodal loads, because the loads are independent of cross-sectional geometry and material which have the same values at the two nodes. Table 2 shows the static displacement, rotation, and frequency results. It can be seen that to achieve comparable accuracy, many more ANSYS elements are required than the presented beam elements. It is seen that when 5 elements are used to calculate the frequency, the presented LMM only has an error of 1.86% relative to the presented CMM. For the dynamic analysis, a 0.1 s (about 1/10 of 1st natural period) pulse loads with amplitudes of the equivalent nodal loads are applied to the free end. Figure 9 illustrates that the presented CMM and LMM are more accurate than the conventional method. Only the results of ANSYS-2 and the presented element are plotted in Figure 9 for brevity.Table 2
The tip displacement, tip rotation, and 1st natural vibration frequency.
Tip disp. (×10−3m)Tip rotation (×10−3rad)1st-freq. (Hz)Ele.This paperANSYS-1ANSYS-2This paperANSYS-1ANSYS-2PresentedCMMANSYS-1,CMMANSYS-2,CMMPresentedLMMANSYS-1,LMMANSYS-2,LMM160.92132.28×10.9317.64×10.6395.049×6.2592.885×465.4764.2011.8611.5110.2389.58810.1789.9299.0289.110563.7863.0111.5311.3210.2359.81110.19610.0459.4439.4941061.6161.4411.0711.0310.11910.2252061.0961.0510.9610.9510.19810.2334060.9660.9510.9310.9310.21810.2358060.9360.9310.22312060.9360.9210.224Note: ×, If 1–3 elements are used for the ANSYS beam54 with variable cross-section (ANSYS-2), ANSYS reports error information of losing accuracy for the sharp change of the beam cross-section.Figure 6
Cantilever beam under the distributed load.Figure 7
Comparison of the displacement interpolation functions between one presented element and one three-order beam element.Comparison of the mass matrixes between one presented element and one three-order beam element. (a) CMM of this paper, (b) CMM of three-order beam, (c) LMM of this paper, and (d) LMM of three-order beam.
(a)(b)(c)(d)Figure 9
Dynamic responses of the beam at right end (results for 1 element of ANSYS-2 are not presented for the reason stated in Table2). (a) Nodal displacement in y direction using CMM, (b) Nodal displacement in y direction using LMM, (c) Nodal rotation using CMM, and (d) Nodal rotation using LMM.


## 6.2. Stepped Beam with Multiple Sectional Geometry

In [17], a free vibration test of an aluminum stepped beam (ρ=2664kg/m3,E=60.6GPa) with dimensions of 463.55 mm × 25.4 mm × 3.175 mm was conducted. The detailed dimensions are shown in Figure 10. The experimental (EXP.) frequencies in x-y plane (in-plane) and x-z plane (out-of-plane) are given in Table 3 together with the calculated solutions, including the Rayleigh-Ritz (RR) method and component modal analysis (CMA), analytical method which is based on the Euler-Bernoulli beam model, FEM with the Euler beam model, Timoshenko beam model, 2D shell, 3D solid model, and the presented beam model. The results of the presented beam element are obtained by one element and the results of beam elements in [17] are obtained by 500 elements. The error 10.8% is due to the nonbeam effects [17] of the test beam in its x-y plane. However, from the results in Table 3, it can be observed that using only one presented element can give acceptable accurate results compared with other beam elements. Another advantage of the presented element is that it can model the stepped beam with only one element rather than ANSYS which cannot model the stepped beam until 13 elements are used.Table 3
Comparison of the frequencies obtained from different methods.
MethodJaworski and Dowell [17]Zheng and  Ji [20]This paperEXP.RRCMAFE/ANSYSANALYTICALFEEuler beamEuler beamEuler beamTimoshenko beam2D shell3D solidEulerbeamEuler beam1st-freq. (in-plane)49.3854.80(+11.0%)54.99(+11.4%)54.47(+10.3%)54.43(+10.2%)49.62(+0.5%)49.8(+0.9%)50.96 (+3.2%)*54.71 (+10.8%)1st-freq. (out-of-plane)10.6310.75 (+1.1%)10.82(+1.8%)10.75(+1.1%)10.75(+1.1)10.44(−1.8%)10.4(−1.6%)—*10.77 (+1.3%)Note: * 54.71 and 10.77 are obtained by 1 presented element. In [17], the results of beam model are obtained by 500 elements. All results are obtained using CMM.Figure 10
Stepped aluminum beam test in [17]. (Unit: mm).


## 6.3. Functionally Graded Beam with Continuous Variable Material and Density

An axially functionally graded cantilever beam shown in Figure11 (ρ=5.7×103×e5(1-x)kg/m3,E=200×e5(1-x)GPa) with dimensions of 1 m × 0.04 m × 0.04 m is analyzed. The ANSYS beam 54 element [36] is used for comparison. Figure 12 shows the comparison of the displacement interpolation functions between one presented element and one three-order beam element. Figure 13 shows the comparison of CMM and LMM between one presented element and one three-order beam element. Similar to the example in Section 6.1, the Figures 12 and 13 illustrate that the presented element can represent the effects of the variable material modulus and mass density. Table 4 shows the static displacement, rotation, and frequency results. It can be seen that many ANSYS elements are required to achieve comparable accuracy, while only a few of presented beam elements are needed. The presented LMM only has an error of 2.73% relative to the presented CMM to calculate the frequency when 5 elements are used. The results of frequency by the presented element in Table 4, and also in Table 2, demonstrate some known knowledge observed by previous researches; for example, the rate is same of using CMM and LMM to achieve a nearly convergence solution. An observation in the paper is that although the same element meshes are used to achieve a nearly convergence frequency solution, the CMM is more accurate when less element meshes are used, because its convergence gradient is smaller than the LMM. For the dynamic analysis, a 0.007 s (about 1/10 of 1st natural period) pulse load with an amplitude of the nodal load is applied to the free end. Figure 14 illustrates that the presented CMM and LMM are more accurate than the conventional method. It should be noted that the ANSYS software gives a poor modeling to this example. The CMM of each ANSYS beam element mesh is calculated by the software itself, in which the midpoint material modulus and mass density of the each mesh are inputted for the calculation. This means that constant material modulus and mass density are used for each element and the exponential variations of them within the element cannot be considered by the software. The LMM of each ANSYS beam element mesh is calculated by Mbeam/2×diag[110110] and then inputted into the software as nodal lumped masses, in which Mbeam is manually calculated to consider the variation of material mass density within the element. These lead to that, as shown in Table 4, ANSYS presents inaccurate static responses and even lower accuracy of the frequencies obtained using CMM (compared with LMM) when a small number of element meshes are used.Table 4
The tip displacement, tip rotation, and 1st natural vibration frequency.
Ele.Tip disp. (×10−3m)Tip rotation (×10−3rad)1st-freq. (Hz)This paperANSYSThis paperANSYSPresented CMMANSYSCMMPresentedLMMANSYSLMM116.4132.0344.9848.02153.1117.1681.2538.45222.0352.25142.9564.88120.3894.04318.8749.32142.1895.92131.68117.28417.7747.67142.04112.62136.03127.25517.2746.77142.01122.00138.13132.281016.6445.42136.622016.4845.05140.734016.4544.96141.798016.44142.0512016.43Figure 11
Functionally graded beam under the nodal load.Figure 12
Comparison of the displacement interpolation functions between one presented element and one three-order beam element.Comparison of the mass matrixes between one presented element and one three-order beam element. (a) CMM of this paper, (b) CMM of three-order beam, (c) LMM of this paper, and (d) LMM of three-order beam.
(a)(b)(c)(d)Figure 14
Dynamic responses of the beam at right end. (a) Nodal displacement iny direction using CMM, (b) Nodal displacement in y direction using LMM, (c) Nodal rotation using CMM, and (d) Nodal rotation using LMM.


## 7. Conclusions

In the paper, an approach has been presented to solve structural static, vibration, and dynamic problems by beam element. The main feature of the approach is the use of adaptive element displacement interpolation functions, which makes the approach quite suitable to solve the responses of varied cross-section beam and functionally graded beam, such as beams with variable cross-section, variable material properties, many different steps in cross-section and/or material properties, and their coupled problems. The element stiffness matrix and equivalent nodal loads are derived. A consistent mass matrix is constructed using the presented adaptive displacement interpolation functions. Then, a diagonally lumped mass scheme which considers the rotational terms is presented by the condensation of the consistent mass matrix. Compared with previously published research works, especially on functionally graded beam, the presented formula is simpler. The accuracy and convergence of this approach has been compared satisfactorily with the existing methods and the experimental results with a number of simple static, eigenvalue, and dynamic examples. One advantage of the presented approach is that the element is of a one-element-one-member configuration, modeling with this type of element would not need to take into account the variation and discontinuity between different parts of the member. Another advantage of the approach is that it is more accurate than conventional finite beam elements in most cases.



---

*Source: 212709-2012-02-16.xml*