# Geometrically Nonlinear Analysis for Elastic Beam Using Point Interpolation Meshless Method

Authors: Cheng He; Xinhai Wu; Tao Wang; Huan He
Journal: Shock and Vibration
 (2019)

---

## Abstract

The intrinsic beam theory, as one of the exact beam formulas, is quite suitable to describe large deformation of the flexible curved beam and has been widely used in many engineering applications. Owing to the advantages of the intrinsic beam theory, the resulted equations are expressed in first-order partial differential form with second-order nonlinear terms. In order to solve the intrinsic beam equations in a relative simple way, in this paper, the point interpolation meshless method was employed to obtain the discretization equations of motion. Different from those equations by using the finite element method, only the differential of the shape functions are needed to form the final discrete equations. Thus, the present method does not need integration process for all elements during each time step. The proposed method has been demonstrated by a numerical example, and results show that this method is highly efficient in treating this type of problem with good accuracy.


---

## Body

## 1. Introduction

Many structures can be modeled as beams in structural mechanics, such as bridge trusses, helicopter rotor blades, DNA/protein molecules, carbon nanotubes, and many others. How to precisely describe the dynamic response of beams is of great interest to many researchers [1–3].As an extension of the classical Kirchhoff–Love rod model, the equations deduced by Reissner in 1972 is regarded as one of the earliest works on geometrically exact beam and have been extended to two/three-dimensional static and dynamic large deformation problems of beam by scholars [4]. Then, the intrinsic equation was proposed for static large displacement space beam by introducing the generalized strain of the reference line [5]. For nonlinear numerical analysis, the second-order nonlinearity is the most ideal case, and the intrinsic beam theory meets this requirement. However, the basic variables in the intrinsic equations were all converted into displacement form when the solution procedure was conducted, and that makes it lack many advantages of the intrinsic formulation. Therefore, Hodges [6] presented a nonlinear intrinsic formulation for the dynamics of initially curved and twisted moving beams. Based on the intrinsic motion and constitutive equations, Hodges expressed the equations in a compact matrix form by only six generalized strain variables. Different to early works, the presented intrinsic formulation does not require reformulation while solving the equations. Thus, the intrinsic beam theory proposed by Hodges was widely used in many industries [7–11]. To evaluate different geometrically nonlinear beam theories, Bauchau et al. [12] compared the accuracy of several commonly used beam models by some standard problems. Pai [13] used a truly geometrically exact displacement-based beam theory to compare with and reveal problems of other geometrically nonlinear beam theories in the literature. In general, the main challenges and problems in geometrically exact beam modeling are (1) the problems of singularity due to the use of three independent rotation variables to describe the cross-sectional rotations, (2) the shear locking problem due to reducing the beam theory’s order by the bending-shear rotation variable, and (3) the accuracy problems caused by numerical discretization.Meanwhile, corresponding numerical methods were also proposed to solve the static and dynamic problems of the developed geometrically exact beam equations [14, 15]. Borri and Mantegazza [16] firstly presented equations of motion in the intrinsic form for dynamics of beams. The equations developed in their work are suitable to describe dynamic behavior of the initially curved and twisted beam, such as the helicopter rotor blades. Cardona and Geradin [17] developed another finite element method for the Reissner beam element, where the rotation vector appears as a dependent variable, and the singular problem of complete Lagrange formula could be avoided when the rotation angle is close to 2Pi and its multiple. In addition, Karlson and Leamy [18] developed a shooting method for computing solutions to nonlinear intrinsic beam equations. As no iteration is needed in the computing procedure, the presented approach may find application in model-based control for practical three-dimensional problems, such as the control of manipulators utilized in endoscopic surgeries and the control of spacecraft with robotic arms and long cables. Khaneh Masjedi and Ovesy [19] and Khaneh Masjedi and Maheri [20] presented a Chebyshev collocation method for the static and free vibration analyses of the geometrically exact beams with fully intrinsic formulation. Le et al. [21] presented a corotational beam element for the dynamic analysis of 3D flexible frames. Wang et al. [22] developed a novel nonlinear aeroelastic model for large wind turbine blades by combining blade element momentum theory and mixed-form formulation of geometrically exact beam theory. Considering that the beam-shaped structures have been widely used in micro- and nano-electromechanical systems, Karparvarfard et al. [23] developed a new geometrically nonlinear elastic size-dependent Euler–Bernoulli beam formulation in the framework of the second strain theory and the corresponding end conditions were also given. Kurka et al. [24] proposed a simplified numerical model of a long flexible beam with variable cross section. Jeong and Yoo [25] proposed a nonlinear modeling method to conduct the static and dynamic analyses of a flexible beam using the in-extensible beam assumption. Lenci et al. [26] compared the free vibrations of Timoshenko beams with mechanical or geometric curvature definitions. The results of different curvature definitions are different for extensible beams. Ahmed and Rhali [27] discussed the nonlinear transverse vibrations of Bernoulli–Euler beams with general boundary conditions based on the extended Hamilton’s principle.In the literatures demonstrated above, there are two main methods for solutions to dynamic problem of large deformation beam: nonlinear finite element methods and finite difference method. However, both methods require a preprocessing step where the domain is discretized, and the accuracy of the results depend on the points size, which means the more the points, the lower the computational efficiency. In this paper, the discretization formulations of the intrinsic beam model are presented by using the point interpolation meshless method. Low-order interpolation functions are used to obtain the shape functions. Since the governing equations of the intrinsic beam are in the first-order form, it is easy to get the discrete form of the intrinsic continua beam directly by substituting the shape functions into the governing partial differential equations. Then, only differential of the shape functions will appear in the final discrete first-order equations. Furthermore, no integration means no need of the Gauss quadrature, which will speed up the whole solving process.


## 2. Intrinsic Continua Formulation for Beams

### 2.1. Geometric Relationship of the Deformed and Undeformed Beam

Consider an initially twisted and curved beam as shown in Figure1. The initial shape of the beam can be described by the undeformed configuration Ω0, and the configuration of the deformed beam is described by Ωf. Ωref is the reference configuration, which is the common configuration of the strait beam generally.Figure 1
Configuration of the initially curved and twisted beam undergoing finite deformation.LetR and r denote the centerline of Ωf and Ω0, respectively, and define χ1 as the running length coordinate of the material point along the center-line of the beam. biχ1 (for i=1,2,3) at material point χ1 are orthogonal unit basis vectors. Among these vectors, b1χ1 is tangent to the center-line of Ω0, while b2χ1 and b3χ1 are fixed in the cross section of the beam. Biχ1 at material point χ1 are another orthogonal unit basis vectors. Unlike b1χ1, because of the deformation of the beam, B1χ1 is not necessarily tangent to the centerline of the deformed configuration R.Defineκ and γ as the curvature vector and the strain vector of the deformed configuration Ωf, respectively. The curvature vector κ contains three components describing twist and bending relative to Bi; thus, it can be expressed by the following equation:(1)κi=κ⋅Bi.The strain vectorγ is shown as follows:(2)γ=γ112γ122γ13T,where γ11 is the axial strain measure of the beam and 2γ12 and 2γ13 are the transverse shear strain measures of the beam.Considering small deformations, the spatial rate of change of the centerline positionR can be expressed by [8](3)R′=1+γ11B1+2γ12B2+2γ13B3,where •′=∂•/∂χ1.And the spatial rate of change of the basis vectorsB can be expressed by(4)B′=κ×B.


### 2.2. Equations of Motion

The governing equations of motion for curved and twisted beam in the intrinsic formulation presented by Hodges [28] is expressed as follows:(5)F′+κ×F+f=P˙+Ω×P,M′+κ×M+e1+γ×F+m=H+Ω×H+V×P,where •˙=∂•/∂t; V and Ω are the linear velocity and angular velocity of the cross section, respectively; P and H are the linear momentum and angular momentum per unit length, which is associated with V and Ω; and F and M are the internal force vector and moment vector, respectively. The unit vector e1=100T is measured relative to Bi.The above two equations can describe the dynamic response of the beam’s cross section located at any arbitrary material pointχ1 while undergoing external distributed force f and moment m vector per unit length.The constitutive law associating internal force vectorF and moment vector M with curvature vector κ and strain vector γ can be expressed by the following equation:(6)FM=Dγ^κ^=Dγ−γ0κ−κ0,where D is the constitutive matrix and γ0 and κ0 are the initial curvature and strain vectors, respectively.The intrinsic formula takesV, Ω, κ, and γ as the variables. There are another three equations implemented as follows to get a close form:(7)Ω′+κ×Ω=κ˙,(8)V′+κ×V+e1+γ×Ω=γ˙,(9)P1P2P3H1H2H3=μ000μχ3−μχ20μ0−μχ30000μμχ2000−μχ3μχ2i2+i300μχ3000i2i23−μχ2000i23i3V1V2V3Ω1Ω1Ω1,where Pi=P⋅Bi and Hi=H⋅Bi (for i=1,2,3); μ is the mass per unit length; χ2 and χ3 are offsets from the reference line of the cross-sectional mass centroid; and i2, i3, and i23 are cross-sectional mass moments and product of inertia.Consider that there is a symmetric cross-sectional beam and mass center of the cross section is located at the centerline. So, we will haveχ2=χ3=0 and i23=0. Then, the expression of equation (9) can be simplified as follows:(10)P=μV,H=IcΩ,where Ic=i1000i2000i3, i1=i2+i3.Substituting equation (10) into equation (3) will lead to(11)μV˙=F′+κ×F+f−μΩ×V.Then substituting equation (10) into equation (4) will give(12)M′+κ×M+e1+γ×F+m=IcΩ˙+Ω×IcΩ+μV×V.SinceV×V=0, equation (12) can be expressed by a simple form(13)IcΩ˙=M′+κ×M+e1+γ×F+m−Ω×IcΩ.Then, we have the following four equations:(14)μV˙=F′+κ×F+f−μΩ×V,IcΩ˙=M′+κ×M+e1+γ×F+m−Ω×IcΩ,Ω′κ×Ω=κ˙,V′+κ×V+e1+γ×Ω=γ˙.For isotropic materials, the constitutive law can be decoupled as(15)D=EA0000002GA/k0000002GA/k000000GJ000000EI2000000EI3.For simple expression, expressD by the following form:(16)D=DF00DM,DF=D1F000D2F000D3F,DM=D1M000D2M000D3M,then the constitutive law can be expressed in the form(17)F=DFγ,M=DMκ.By substituting equation (17) into the four equations, the following equations are obtained:(18)μV˙=DFγ′+κ×DFγ−μΩ×V+f,IcΩ˙=DMκ′+κ×DMκ+e1+γ×DFγ−Ω×IcΩ+m,κ˙=Ω′+κ×Ω,γ˙=V′+κ×V+e1+γ×Ω.Now we re-express equation (18) in another form as follows:(19)μV˙=cγ′+κ×cγ+f−μΩ×V,IcΩ˙=qκ′+κ×qκ+e1+γ×cγ+m−Ω×IcΩ,κ˙=Ω′+κ×Ω,γ˙=V′+κ×V+e1+γ×Ω,V˙=cμγ′+κ×cμγ−Ω×V+1μf,Ω˙=Ic−1qκ′+Ic−1κ×qκ+Ic−1e1+γ×cγ−Ic−1Ω×IcΩ+Ic−1m,(20)Aη˙=Qη′+fnon+fext,where(21)η=VΩκγ,A=μI0000Ic0000I0000I,Q=000DF00DM00I00I000,fext=fm00,fnon=κ×DFγ−μΩ×Vκ×DMκ+e1+γ×DFγ−Ω×IcΩκ×Ωκ×V+e1+γ×Ω.


## 2.1. Geometric Relationship of the Deformed and Undeformed Beam

Consider an initially twisted and curved beam as shown in Figure1. The initial shape of the beam can be described by the undeformed configuration Ω0, and the configuration of the deformed beam is described by Ωf. Ωref is the reference configuration, which is the common configuration of the strait beam generally.Figure 1
Configuration of the initially curved and twisted beam undergoing finite deformation.LetR and r denote the centerline of Ωf and Ω0, respectively, and define χ1 as the running length coordinate of the material point along the center-line of the beam. biχ1 (for i=1,2,3) at material point χ1 are orthogonal unit basis vectors. Among these vectors, b1χ1 is tangent to the center-line of Ω0, while b2χ1 and b3χ1 are fixed in the cross section of the beam. Biχ1 at material point χ1 are another orthogonal unit basis vectors. Unlike b1χ1, because of the deformation of the beam, B1χ1 is not necessarily tangent to the centerline of the deformed configuration R.Defineκ and γ as the curvature vector and the strain vector of the deformed configuration Ωf, respectively. The curvature vector κ contains three components describing twist and bending relative to Bi; thus, it can be expressed by the following equation:(1)κi=κ⋅Bi.The strain vectorγ is shown as follows:(2)γ=γ112γ122γ13T,where γ11 is the axial strain measure of the beam and 2γ12 and 2γ13 are the transverse shear strain measures of the beam.Considering small deformations, the spatial rate of change of the centerline positionR can be expressed by [8](3)R′=1+γ11B1+2γ12B2+2γ13B3,where •′=∂•/∂χ1.And the spatial rate of change of the basis vectorsB can be expressed by(4)B′=κ×B.


## 2.2. Equations of Motion

The governing equations of motion for curved and twisted beam in the intrinsic formulation presented by Hodges [28] is expressed as follows:(5)F′+κ×F+f=P˙+Ω×P,M′+κ×M+e1+γ×F+m=H+Ω×H+V×P,where •˙=∂•/∂t; V and Ω are the linear velocity and angular velocity of the cross section, respectively; P and H are the linear momentum and angular momentum per unit length, which is associated with V and Ω; and F and M are the internal force vector and moment vector, respectively. The unit vector e1=100T is measured relative to Bi.The above two equations can describe the dynamic response of the beam’s cross section located at any arbitrary material pointχ1 while undergoing external distributed force f and moment m vector per unit length.The constitutive law associating internal force vectorF and moment vector M with curvature vector κ and strain vector γ can be expressed by the following equation:(6)FM=Dγ^κ^=Dγ−γ0κ−κ0,where D is the constitutive matrix and γ0 and κ0 are the initial curvature and strain vectors, respectively.The intrinsic formula takesV, Ω, κ, and γ as the variables. There are another three equations implemented as follows to get a close form:(7)Ω′+κ×Ω=κ˙,(8)V′+κ×V+e1+γ×Ω=γ˙,(9)P1P2P3H1H2H3=μ000μχ3−μχ20μ0−μχ30000μμχ2000−μχ3μχ2i2+i300μχ3000i2i23−μχ2000i23i3V1V2V3Ω1Ω1Ω1,where Pi=P⋅Bi and Hi=H⋅Bi (for i=1,2,3); μ is the mass per unit length; χ2 and χ3 are offsets from the reference line of the cross-sectional mass centroid; and i2, i3, and i23 are cross-sectional mass moments and product of inertia.Consider that there is a symmetric cross-sectional beam and mass center of the cross section is located at the centerline. So, we will haveχ2=χ3=0 and i23=0. Then, the expression of equation (9) can be simplified as follows:(10)P=μV,H=IcΩ,where Ic=i1000i2000i3, i1=i2+i3.Substituting equation (10) into equation (3) will lead to(11)μV˙=F′+κ×F+f−μΩ×V.Then substituting equation (10) into equation (4) will give(12)M′+κ×M+e1+γ×F+m=IcΩ˙+Ω×IcΩ+μV×V.SinceV×V=0, equation (12) can be expressed by a simple form(13)IcΩ˙=M′+κ×M+e1+γ×F+m−Ω×IcΩ.Then, we have the following four equations:(14)μV˙=F′+κ×F+f−μΩ×V,IcΩ˙=M′+κ×M+e1+γ×F+m−Ω×IcΩ,Ω′κ×Ω=κ˙,V′+κ×V+e1+γ×Ω=γ˙.For isotropic materials, the constitutive law can be decoupled as(15)D=EA0000002GA/k0000002GA/k000000GJ000000EI2000000EI3.For simple expression, expressD by the following form:(16)D=DF00DM,DF=D1F000D2F000D3F,DM=D1M000D2M000D3M,then the constitutive law can be expressed in the form(17)F=DFγ,M=DMκ.By substituting equation (17) into the four equations, the following equations are obtained:(18)μV˙=DFγ′+κ×DFγ−μΩ×V+f,IcΩ˙=DMκ′+κ×DMκ+e1+γ×DFγ−Ω×IcΩ+m,κ˙=Ω′+κ×Ω,γ˙=V′+κ×V+e1+γ×Ω.Now we re-express equation (18) in another form as follows:(19)μV˙=cγ′+κ×cγ+f−μΩ×V,IcΩ˙=qκ′+κ×qκ+e1+γ×cγ+m−Ω×IcΩ,κ˙=Ω′+κ×Ω,γ˙=V′+κ×V+e1+γ×Ω,V˙=cμγ′+κ×cμγ−Ω×V+1μf,Ω˙=Ic−1qκ′+Ic−1κ×qκ+Ic−1e1+γ×cγ−Ic−1Ω×IcΩ+Ic−1m,(20)Aη˙=Qη′+fnon+fext,where(21)η=VΩκγ,A=μI0000Ic0000I0000I,Q=000DF00DM00I00I000,fext=fm00,fnon=κ×DFγ−μΩ×Vκ×DMκ+e1+γ×DFγ−Ω×IcΩκ×Ωκ×V+e1+γ×Ω.


## 3. Discretization for the Intrinsic Beam Formula

We use a series of basis functionsφ1χ1φ2χ1⋯φnχ1 to interpolating Vi, Ωi, κi, and γ¯i, i=1,2,3, where γ¯1=γ1, γ¯2=2γ2, γ¯3=2γ3, and γ=γ¯1γ¯2γ¯3T. For simplicity of the expressions, we use x to replace χ1. Then, the expressions of the interpolation will be(22)Vix=∑j=1najiφjx=φ1xφ2x⋯φnxa1ia2i⋮ani=φTxai,i=1,2,3,(23)Ωix=∑j=1nbjiφjx=φ1xφ2x⋯φnxb1ib2i⋮bni=φTxbi,i=1,2,3,(24)κix=∑j=1ncjiφjx=φ1xφ2x⋯φnxc1ic2i⋮cni=φTxci,i=1,2,3,(25)γ¯ix=∑j=1ndjiφjx=φ1xφ2x⋯φnxd1id2i⋮dni=φTxdi,i=1,2,3.To each interior point located atxr, which is also the rth influence domain, we have(26)Vixr=φTxrai,Ωixr=φTxrbi,κixr=φTxrci,γ¯ixr=φTxrdi,i=1,2,3.As toVi, if each influence domain contains 3 grid points and n=3, for example, then, the number of the interior point will be r=2,3,…,n−1. The following expressions will be satisfied as(27)Vixr−1=φ1xr−1φ2xr−1φ3xr−1ai,Vixr=φ1xrφ2xrφ3xrai,Vixr+1=φ1xr+1φ2xr+1φ3xr+1ai,i=1,2,3.Writing equation (27) in a matrix form, the following equation is obtained:(28)Vixr−1VixrVixr+1=φ1xr−1φ2xr−1φ3xr−1φ1xrφ2xrφ3xrφ1xr+1φ2xr+1φ3xr+1ai,i=1,2,3.Since the expression ofφTx is known to us, the unknown ai can be determined by(29)ai=φ1xr−1φ2xr−1φ3xr−1φ1xrφ2xrφ3xrφ1xr+1φ2xr+1φ3xr+1−1Vixr−1VixrVixr+1,i=1,2,3.Substituting equation (29) into equation (22) will lead to(30)Vix=φ1xφ2xφ3xφ1xr−1φ2xr−1φ3xr−1φ1xrφ2xrφ3xrφ1xr+1φ2xr+1φ3xr+1−1Vixr−1VixrVixr+1,i=1,2,3.S(x) is denoted by the following equation:(31)Sx=φ1xφ2xφ3xφ1xr−1φ2xr−1φ3xr−1φ1xrφ2xrφ3xrφ1xr+1φ2xr+1φ3xr+1−1=S1xS2xS3x.Then,V within the rth influence domain can be expressed by(32)Vix=SxVixr−1VixrVixr+1,i=1,2,3.Following the same process, we will have(33)Ωix=SxΩixr−1ΩixrΩixr+1,i=1,2,3,κix=Sxκixr−1κixrκixr+1,i=1,2,3,γ¯ix=Sxγ¯ixr−1γ¯ixrγ¯ixr+1,i=1,2,3.Align all of the variables in the following form:(34)δr=ηr−1TηrTηr+1TT,where δr denotes generalized grid point displacement for the rth influence domain, ηr−1, ηr, and ηr+1 are the generalized grid point displacements of the r−1th, rth, and r+1th grid point, respectively, and(35)ηk=V1kV2kV3kΩ1kΩ2kΩ3kκ1kκ2kκ3kγ¯1kγ¯2kγ¯3kT,k=r−1,r,r+1,where the superscript k denotes the kth grid point. Or, equation (35) can be expressed by(36)ηk=VkTΩkTκkTγkTT,k=r−1,r,r+1,where VkT=V1kV2kV3k, ΩkT=Ω1kΩ2kΩ3k, κkT=κ1kκ2kκ3k, and γkT=γ¯1kγ¯2kγ¯3k.Thenηx within xr−1 to xr+1 can be expressed by(37)ηx=N1xN2xN3xηr−1ηrηr+1=Nxδr,where the shape function Nx is a 12×36 matrix. N1x, N2x, and N3x are all 12×12 diagonal matrices, and the corresponding expressions are(38)Nix=SixI,i=1,2,3,where the matrix I is a 12×12 identity matrix.Substituting equation (37) into equation (20) will result in(39)ANxδ˙r=QN′xδr+fnon+fext.As to the interpolation meshless method, equation (39) must be satisfied at each grid point xr; thus, the following equations are obtained:(40)ANx2δ˙2=QN′x2δ2+fnon2+fext2,ANx3δ˙3=QN′x3δ3+fnon3+fext3,⋮ANxn−1δ˙n−1=QN′xn−1δn−1+fnonn−1+fextn−1,where(41)fnonr=κr×DFγr−μΩr×Vrκr×DMκr+e1+γr×DFγr−Ωr×IcΩrκr×Ωrκr×Vr+e1+γr×Ωr,r=2,3,…,n−1.Equation (40) is a series of first-order ordinary differential equations which describe motion of each point except the tips.If we depart the shape functionNx as follows:(42)Nx=NVxNΩxNκxNγx,then we will have(43)V=NVxδr,Ω=NΩxδr,κ=Nκxδr,γ=Nγxδr,where NVx, NΩx, Nκx, and Nγx∈ℝ3×36. Denote NVr=NVxr, NΩr=NΩxr, Nκr=Nκxr, and Nγr=Nγxr.SubstitutingVr=NVrδr,Ωr=NΩrδr,κr=Nκrδr,andγr=Nγrδr,r=2,3,…,n−1, into equation (41), the following is obtained:(44)fnonr=NVrδr×DFNγrδr−μNΩrδr×NVrδrNκrδr×DMNκrδr+e1+Nγrδr×DFNγrδr−NΩrδr×IcNΩrδrNκrδr×NΩrδrNκrδr×NVrδr+e1+Nγrδr×NΩrδr,r=2,3,…,n−1.To each tip point, we consider that each influence domain contains 2 grid points andn=2. Following sequence similar to the previous steps, for the first point and the last point, the generalized deformation can be interpolated by the following expressions:(45)Vix=φ1xφ2xφ1x1φ2x1φ1x2φ2x2−1Vix1Vix2,i=1,2,3atx1,Vix=φ1xφ2xφ1xn−1φ2xn−1φ1xnφ2xn−1Vixn−1Vixn,i=1,2,3atxn.Denoting(46)S1x=φ1xφ2xφ1x1φ2x1φ1x2φ2x2−1=S1,1xS1,2x,Snx=φ1xφ2xφ1xn−1φ2xn−1φ1xnφ2xn−1=Sn,1xSn,2x.All of the variables at x1 can be expressed by(47)Vix=S1xVix1Vix2,i=1,2,3,Ωix=S1xΩix1Ωix2,i=1,2,3,κix=S1xκix1κix2,i=1,2,3,γ¯ix=S1xγ¯ix1γ¯ix2,i=1,2,3.And all of the variables atxn can be expressed by(48)Vix=SnxVixn−1Vixn,i=1,2,3,Ωix=SnxΩixn−1Ωixn,i=1,2,3,κix=Snxκixn−1κixn,i=1,2,3,γ¯ix=Snxγ¯ixn−1γ¯ixn,i=1,2,3.Then, the generalized response vector at each tip can be expressed by(49)ηx=N11xN12xη1η2=N1xδ1,atx1,ηx=Nn1xNn2xηn−1ηn=Nnxδn,atxn,where N11x, N12x, Nn1x, and Nn2x are all 12×12 diagonal matrices; thus, N1x and Nnx∈ℝ12×24. The expressions of N1ix and Nnix are listed in the following:(50)N1ix=S1,ixI,Nnix=Sn,ixI,i=1,2.Substituting equation (49) into equation (20) and replacing x as x1 and xn, respectively, will give(51)AN1x1δ˙1=QN1′x1δ1+fnon1+fext1,ANnxnδ˙n=QNn′xnδn+fnonn+fextn,where fnonr has the same form as shown in equation (41).By the same steps, we departN1x and Nnx as follows:(52)Nrx=NVrxNΩrxNκrxNγrx,r=1,n,where NVrx, NΩrx, Nκrx, and Nγrx∈ℝ3×24, r=1,n. Denote NVr=NVxr, NΩr=NΩxr, Nκr=Nκxr, and Nγr=Nγxr,r=1,n.SubstitutingVr=NVrδr,Ωr=NΩrδr,κr=Nκrδr,andγr=Nγrδr, r=1,n, into equation (41) will give the same form as shown in equation (42).Denote(53)δ=δ1Tδ2T⋯δnTT.To equations (40) and (51), we expand δr, r=1,3,…,n, to δ and put the coefficient matrices on the appropriate position and then assembly each equation. Finally, we will obtain a series of first-order nonlinear differential equations as follows:(54)A¯δ˙=Q¯δ+f¯non+f¯ext,where(55)f¯non=fnon1fnon2⋮fnonn,f¯ext=fext1fext2⋮fextn.The expression offnonr has been illustrated as equation (44), and fextr is the external force which is known at the beginning of the simulation.From above, it is obvious that equation (54) is a series of first-order ordinary differential equations of time, which can be solved by using ordinary time-domain integral algorithm.After calculating all the intrinsic parameters, the radius vector of the beam at any cross section can be expressed as(56)R′=1+γ1B1+γ2B2+γ3B3,(57)Bi′=κ×Bi,i=1,2,3.Deformation of any point on the beam is obtained from equations (56) and (57).


## 4. Numerical Example

Consider a rectangular cross-sectional cantilever beam with a length of 75 cm, a cross-sectional size of10cm×10cm, a line density of 0.1 kg/m, an elastic modulus of 1 MPa, and a Poisson’s ratio of 0.25. At the initial moment, two constant concentrated forces of the same magnitude were applied at the free end of the beam along the Y and Z direction, respectively, as shown in Figure 2.Figure 2
Sketch of the cantilever beam subjected to the follower force.The dynamic responses of the cantilever beam in 1 s were calculated by the proposed method, and the displacement responses of the selected five nodes are plotted as shown by the solid curve in Figure3. In this figure, x1 = 0 to x1 = 75 represent the coordinates of the selected point in the natural coordinate system of the beam.Figure 3
Comparison of the deformation at several given locations based on the ABAQUS model and the proposed model.
(a)
(b)
(c)As a comparison, the dynamic model of the cantilever beam was built by using 480 C3D8 hexahedral elements of same size in the nonlinear transient finite element analysis program ABAQUS. The displacement response at the corresponding positions calculated by ABAQUS is shown by the dotted curve in Figure3.Figure3 shows that the cantilever beam has a great deformation, and the three displacement components at the free end have reached above 40 cm. Comparing the results calculated by the above two methods in Figure 3, the dynamic response curves of the three displacement components are in good agreement with each other, which verify the accuracy of the proposed method.As shown in Figure3, the cantilever beam oscillates under the concentrated force and the displacement curve exhibits a “U” shape. Since the concentrated force magnitudes in the Y and Z directions are the same, the corresponding response curves are also symmetrical. At about 0.4 s, the beam has a rebound.Figure4 shows the time-history responses of the curvature and strain, which are the basic variables of the intrinsic beam formula. The different colored curves in Figure 4 represent dynamic responses at different locations on the beam (for example, the red solid curve represents the curvature and strain response of the free end (x = 75 cm)). Since only concentrated forces are applied at the free end, the two curvature components are zero, as shown in the red lines in Figures 4(a) and 4(b). In addition, it shows that the κ2,κ3 and γ2,γ3 of the displacement components u2,u3 are symmetrical due to the symmetry of the load. This is consistent with physical phenomenon.Figure 4
Curvature and strain response of the beam at several different locations.
(a)
(b)
(c)
(d)
(e)


## 5. Conclusion

Different from traditional nonlinear beam theory, the intrinsic beam formula only has first derivative term and second-order nonlinear term. In this paper, by using the intrinsic beam theory, the point interpolation meshless method was applied to calculate the dynamic response of large flexible beams. Conclusions are as follows:(1)
Only a set of discrete interpolate points are needed by using the presented method, which makes it easy to establish the discrete nonlinear dynamic equations for the large flexible beams.(2)
Compared with the traditional nonlinear finite element methods, the presented method avoids the integration process in each calculation step. Because of this, the computational efficiency for solving the large flexible beam equation is greatly improved. Due to the geometrically exact characteristics of the intrinsic beam theory, dynamic response of the flexible beam can be calculated very accurately with less discrete points. Numerical results show that the proposed method is highly efficient in treating this type of problem with good accuracy.(3)
The governing equation of the large flexible beam is discretized by using the point interpolation meshless method. If the compactly supported functions are chosen as the basis functions, the discrete governing equation can be expressed in the sparse matrices form, which will benefit the computational speed.



---

*Source: 9065365-2019-05-15.xml*