A node-based smoothed point interpolation method for dynamic analysis of rotating flexible beams

C.F.Du·D.G.Zhang·L.Li·G.R.Liu

1 Introduction

Discretization of the deformation field of a flexible body in rigid-flexible coupled systems requires careful consideration to accurately model the dynamical behavior of the systems.Currently,the discrete methods including the assumed mode method(AMM)[1–3]and the finite element method(FEM)[4–7]are widely used.The AMM is originated from the method of using natural modes of vibration.The advantage is that a good approximate result can be achieved by using few modes,and hence the computation efficiency is often very high.However,how to find the natural modes can be difficult when the flexible body has an irregular shape.The FEM on the other hand is versatile for complex geometry.It is a robust and thoroughly developed discrete method that is widely used to solve many types of linear and non-linear engineering problems.Most practical engineering problems related to solids and structures are currently solved by using FEM packages that are commercially available.However,the FEM has also inherent shortcomings.It uses the Galerkin weak form and the essential formulations are confined with elements,which leads to discontinuous of derivatives for the displacement functions arose the interface of the elements,so that the stress solutions is low accuracy[8].Therefore,researchers have begun to look for new theory and discrete method for effective modeling the rigid-flexible coupled systems dynamics,such as the absolute nodal coordinate formulation[9,10]and the B-spline interpolation method[11,12].

Recently,mesh-free methods have been proposed and applied to solve complex problems in engineering and sciences.Mesh-free methods use shape functions created using nodes,including the smoothed particle hydrodynamic method(SPH)[13–15],the element-free Galerkin method(EFG)[16],there producing kernel particle methods(RKPM)[17],the meshless local Petrov–Galerkin method(MLPG)[18],the point interpolation method(PIM)[19,20],and the radial point interpolation method(RPIM)[21,22].The nodebased smoothed PIM(NS-PIM or LC-PIM termed originally[23])has been recently developed and found very stable,is free from volumetric locking,and capable ofproducing upper bound solutions for force driving static problems[24–26].The discretized system equations in an NS-PIM are established using weakened weak(W2)forms based on a set of background cells,and the numerical operations are not confined within the cells,but across the neighboring cells.The W2 form introduces the so-called softening effects to the discretized system equation,and hence can effectively overcome the disadvantage of the FEM[27].This work aims to extend the mesh-free NS-PIM to solve problems of rigid flexible coupled systems dynamics.

The behavior of beams based on the Euler–Bernoulli theory are governed by a typical 4th order differential equation with respect to coordinates.In the early work of solving 4th order boundary value problems,both displacement and slope(derivative of the displacement)are treated as independent variables.In this paper,due to the use of W2 concept,we solve the 4th order boundary value problems considering only displacement as the unknown variable,and its approximation is performed using the simple linear point interpolation method[28].In other words,we solve 4th order differential equations using only 1st order approximation.In our examination for static problems,the beams with three types of boundary conditions are analyzed,and the results are compared with the analytic solution,which shows the effectiveness of this method.The NS-PIM is then further extended for solving a rigid-flexible coupled system dynamics problem,considering a rotating flexible cantilever beam.In this case,we consider not only the transverse deformations of the rotating flexible cantilever beam,but also the longitudinal deformations.The rigid-flexible coupled dynamic equations of the system are derived via employing Lagrange’s equations of the second kind.Simulation results of the NS-PIM are compared with those obtained using FEM and AMM.It is found that compared with FEM,the NS-PIM makes the system soften and has anti-ill solving ability under the same calculation conditions.

2 Static analysis

2.1 Brief on the weak form for the Euler–Bernoulli beam

Considering the small deformation situation,for a giving bending stiffnessEI,the strong form of governing equation of an Euler–Bernoulli beam is expressed by the well-known fourth-order differential equation[27]

wherevis transverse deflection,andb yis the distributed load in the transverse direction along the beam.Depending on the problem,the boundary conditions can be given using a complementary pair of equations listed as follows

whereMandQdenote the bending moment and shear force acting at the cross-section of the beam,respectively.ΓvandΓθare essential boundaries where deflection and slope are specified;ΓMandΓQare natural boundaries where bending moment and shear force are specified,respectively.The complementary condition requires that when the deflection(or slope)is specified at an end of the beam,the shear force(or moment)cannot be specified.This is to avoid the possibility of the Cauchy boundary value problem not being well posed[29].For example,Eq.(2)may pair with Eq.(3)(for fully clamped),Eq.(2)with Eq.(4)(simply support),Eq.(3)with Eq.(5)(sliding),and Eq.(4)with Eq.(5)(free).

The standard Galerkin weak form used in the FEM formulation can be rewritten as

in which we may note that the assumed deflection functionvmust be of at least 2nd order consistency.It is weaker compared to the strong form give in Eq.(1)which requires the assumed deflection function is at least 4th order consistency.The “weaker”here means that the assumed deflection function has lower consistency.

2.2 Formulation of NS-PIM for Euler–Bernoulli beam

2.2.1 Smoothed Galerkin weak form:a typical W2

In a W2 formulation,we consider a beam discretized withNnnodes,as shown in Fig.1.Every two neighboring nodes form a cell known as background cell.For each node,a smoothing domain is formed that consists of the nearer halves of the two

Fig.1 Discretization of the problem domain with nodes and the no debased smoothing domains

Using the gradient smoothing operations[30],the smoothed first derivative ofvcan be expressed as

in which we may note:(1)the deflection functionvcan now only be of the first order consistency(a weakened weak form);(2)the domain integration on the left hand side is now converted to a simple summation(for better computational efficiency).

2.2.2 NS-PIM equations for Euler–Bernoulli beams

Because of the W2 formulation,in our NS-PIM,each node has only one degree of freedom(DOF)which is the deflection.We now use a linear interpolation of the deflectionvin each cell(cellΩcnin Fig.1)to approximate the deflection,which can be expressed as

where the linear shape function can be simply given as

wherelcn=xn+1?xnis the length of the cellΩcn,v n,andv n+1denote the nodal values of the deflection at nodesnandn+1.Substituting Eq.(11)into Eq.(9),the smoothed curvature for the smoothing domainΩsnwith interior nodes can be given as

wheren=2,3,...,Nn?1.

For simply supported boundary conditions,the zero deflection is set at the corresponding boundary point,just like in the FEM.For clamped boundary conditions,however,the following treatment is needed to impose the zero-slope conditions.The reason is the formulation of NS-PIM does not have rotation as a DOF.For the clamped boundary on the left end(n=1),we used the smoothed gradient that can be given by

Similarly,for the clamped boundary on the right end(n=Nn),the smoothed gradient is given by

It is clear that the zero-slope condition is imposed using the relation of two nodes on the boundary,which can be conveniently done without using two DOFs at a node.

Next,substituting Eq.(13)into Eq.(10),a set of discretized algebraic system equations can be obtained in the following matrix form

whereΦ(x)is the global matrix of all these nodal shape functions.The global smoothed stiffness matrix is assembled in the form of

whereK nwhich is the stiffness matrix associated withΩsnis computed using

The summation in Eq.(18)means an assembly(node marched summary).The assembling process is similar as that in the FEM,but must be based on the smoothing domains.We may note that in the FEM for beam,there are two DOFs(deflection and rotation)at each node,which is double of that in NS-PIM when the same number of nodes is used.

2.3 Numerical examination of NS-PIM for static problems

The NS-PIM is coded and examination has been performed to verify the formulation for static problems.Figure 2 shows three types of boundary conditions with distributed load that are considered in this examination.The parameters are taken asL=1 m,E I=1 N ·m2,andq=1 N/m.Figure3a–c shows the results of deflection under three types of boundary conditions using NS-PIM,FEM,AMM,together with the exact solution.The NS-PIM uses 22 DOFs,the FEM also uses 22 DOFs,and the AMM uses the first three modes.It is shown that these numerical results agree very well,and the differences are not distinguishable.These examples indicate that the NS-PIM has high accuracy even if using linear interpolation.

The results of deflection are computed and listed in Tables1–3.It can be seen that the NS-PIM provides an upper bound solution for the deflection.This demonstrates that the NS-PIM makes the system soften.

Figure 4 shows a cantilever beam subjected to pure bending moment at free end.The bending momentMe=1 N·m.Figure 5 shows the results of deflection.The NS-PIM uses 22 DOFs.It is shown that the numerical results agree very well with the exact solution,and the differences are not distinguishable.It also indicates that the pure bending boundary condition is able to be addressed by NS-PIM even if there is no rotation as the degree of freedom.

3 NS-PIM of rigid-flexible coupled system dynamics

3.1 NS-PIM for rotating flexible beams

For rotating flexible beams,we assume that for the bending deformations,the rotating flexible beam obeys the Euler–Bernoulli beam theory.This means that the cross-section of the beam is perpendicular to the central axis of the beam and remains in a plane after deformation.

Fig.2 Three types of boundary conditions of beams with distributed load:cantilever beam,simply–simply supported beam,and clamped-simply supported beam,respectively

Consider a flexible beam attached to a rigid hub of radiusa,and rotating about theZaxis,following the right-hand rule convention in the inertial coordinate systemOXYZ,as shown in Fig.6.A planar floating coordinate system,which rotates with the beam is denoted byO?x yz.The rotary inertia of the hub isJh.Physical parameters of the beam are as follows:lengthL,cross-section areaS,Young’s modulusE,area moment of inertiaI,and mass densityρ.The position vector of arbitrary pointPon the flexible beam in the inertial frameOXYZafter deformation can be given by

Fig.3 The deflection under different boundary conditions:a cantilever beam;b simply–simply supported beam;c clamped-simply supported beam

Table 1 The deflection of cantilever beam

Table 2 The deflection of simply–simply supported beam

Table 3 The deflection of clamped-simply supported beam

In Eq.(21),Θis the direction cosine matrix of floating coordinate system relative to the inertial coordinate system.The

Fig.4 A cantilever beam subjected to pure bending moment at free end

Fig.5 The deflection with pure bending moment at free end

Fig.6 The deformation of rotating flexible beam

deformation vectoruin the floating coordinate systemO?x yzcan be expressed by

wherew1is the axial stretch deformation andw2is the transverse deformations of the flexible beam.w cis the coupling term of the deformation which is caused by the transverse deformation.It can be given by

The velocity vector of arbitrary pointPin the inertial coordinate system can be obtained through the first derivative of Eq.(20)versus time.It can be given by

Then the kinetic energy of the system is

The elastic potential energy of the system is

3.2 Dynamic equations of the rotating flexible beam system

The flexible beam is divided intoNccells shown in Fig.1.As discussed in Sect.2.3,using Eqs.(11)and(12),the axial stretch deformationw1and the transverse deformationsw2inith cell can be expressed

whereAi(t)andB i(t)are the axial stretch deformation vector and the transverse deformations vector of nodesiandi+1 with the change of time,respectively.The generalized coordinates of the beam are

So the nodal coordinate vector can be written as

whereR iis a Boolean matrix given by

Substituting Eq.(31)into Eq.(22),the deformation matrix of arbitrary pointPon the flexible beam can be expressed by

whereH(x)is the coupling shape functions,it can be given by

tor.Substituting Eqs.(25)and(26)into Lagrange’s equations of the second kind

whereJohandJobare constant;S,M,C,andDare constant matrixes.The underlined term is additional coupling term which caused by transverse bending deformation and it makes the longitudinal deformation shorten.This term appears in the generalized mass matrix and the generalized force matrix.

With the discussion of Sect.2.3,using the smoothing operations,the stiffness matrixK2can be expressed by

Using Eqs.(18)and(19),the global smoothed stiffness matrix can be obtained.

3.3 Transient dynamic simulation

When the angular velocity of the hub is given,the dynamic equations of the system can be expressed by

whereT=15 s.The angular velocity of the hub reachesΩ0at timeT.The hub radiusa=0.Parameters of the rotating beam are[31]:L=8 m,S=7.2968× 10?5m2,I=8.2189 × 10?9m4,ρ=2766.7 kg/m3,andE=68.952 GPa.

The tip responses of the rotating beam with differentΩ0are shown in Figs.7–10.The simulation results of AMM,FEM,and NS-PIM are compared.AMM uses three to seven modes,FEM uses ten elements and NS-PIM uses 21 nodes.Figure 7a shows the longitudinal deformation at the tip of the beam whenΩ0=4 rad/s.It is seen that these three discrete methods almost overlap with each other.The longitudinal deformationu xcontains two parts:the axial stretch deformationw1and the coupling term of the deformationw c,which is caused by the transverse deformation.It is negative meaning that the longitudinal deformation is mainly caused by the coupling termw c.Figures 7b–c and 8 show the tip transverse bending deformation and its deformation rate which is the first derivative of deformation versus time,respectively.As shown in these figures,the results of these three discrete methods almost overlap with each other.Figure 9a shows the longitudinal deformation at the tip of the beam whenΩ0=20 rad/s.In order to illustrate the less accuracy of AMM when the deformation is large,the elasticity modulus reduces to 50 GPa.It is seen that the results of NS-PIM and FEM are almost overlap,but AMM is quite different.Figures 9b–c and 10 show the tip transverse bending deformation and its deformation velocity,respectively.It is seen that the results of AMM is quite different from the other two methods.That means AMM is not applicable when the deformation is large,it is applicable in the case of small deformation.Compared with FEM,NS-PIM is a little different.Figures 9c and 10c show that the amplitude of AMM is much smaller than FEM,the amplitude of NS-PIM is nearly the same as FEM.However,the vibration phase of NS-PIM is different from FEM.This is because of the soft characteristic of NS-PIM.In order to illustrate that AMM is not applicable in the case of large deformation,Fig.11 shows the tip transverse bending deformation whenΩ0=20 rad/s andE=20 GPa.It is seen that the largest deformation is over 5 m,belong to large deformation problem.AMM is divergent,that means AMM is not applicable in the large deformation problem.NS-PIM and FEM are convergent.But there are a little difference between NS-PIM and FEM.Still,this is because of the soft characteristic of NS-PIM.

Fig.7 Tip responses of the rotating beam for slow rotation(Ω0=4 rad/s):a The longitudinal deformation at the tip of the beam;b the transverse deformation at the tip of the beam;c the transverse deformation at the tip of the beam after 15 s

Fig.8 Tip transverse deformation rate of the beam for slow rotation(Ω0=4 rad/s):a tip transverse deformation rate from 0 to 20 s;b tip transverse deformation rate after 15 s

Fig.9 Tip responses of the rotating beam for fast rotation(Ω0=20 rad/s):a the longitudinal deformation at the tip of the beam;b the transverse deformation at the tip of the beam;c the transverse deformation at the tip of the beam after 15 s

Fig.10 Tip transverse deformation rate of the beam for fast rotation(Ω0=20 rad/s):a tip transverse deformation rate from 0 to 20 s;b tip transverse deformation rate after 15 s

Tables 4 and 5 show the relative time,the relative error and the amplitude of these three discrete methods whenΩ0=4 rad/s andΩ0=20 rad/s,respectively.AMM uses the first three to seven modes,and the relative time is com-pared with three modes AMM.The relative error is compared with ten elements FEM.NS-PIM uses 21 nodes.It is seen that three modes AMM computes the fastest,and NS-PIM is faster than FEM.WhenΩ0=4 rad/s,the deformation is small so AMM is applicable.For AMM,the error of the amplitude is small and the accuracy is not improved even if increasing the modes.Instead,the computational efficiency reduces significantly.WhenΩ0=20 rad/s,the deformation is large so AMM is not applicable and the error is big to 45.45%.Regardless of the deformation is small or large,the FEM and NS-PIM are applicable.But the amplitude of NS-PIM is larger than FEM,this is also because of the soft characteristic of NS-PIM,and its accuracy is very high even if using the linear shape function.

When there is an external torque acting on the hub,the large overall motion of the system is unknown.The dynamic responses of the beam can be calculated from Eq.(35).The parameters of the beam are the same as the situation when the angular velocity of the hub is given.Assuming the external torque of the hub is given by

Fig.11 Large deformation situation

Table 4 The relative time,the relative error and the amplitude of these three discrete methods(Ω0=4 rad/s)

Table 5 The relative time,the relative error and the amplitude of these three discrete methods(Ω0=20 rad/s)

whereT=10 s,andτ0is the amplitude of the external torque.

The rotation angle and the tip responses of the rotating beam are shown in Fig.12 whenτ0equals 10 N ·m.It is seen that these three discrete methods almost overlap with each other.When the external torque is removed,the flexible beam swings back and forth in which the rotation angle is equal to 4.5 rad.From Fig.12c,one can observe that the amplitude of NS-PIMislagerand the vibrational frequency is smaller than the other two methods.This phenomenon demonstrates again that the NS-PIM makes the system soften.Figure13 shows the tip transverse deformation of FEM and NS-PIM whenτ0equals 50 N ·m under the same calculation conditions.It is seen that the result of FEM is divergent quickly.This is because the dynamic equations of FEM are ill-conditioned.However,the result of NS-PIM is still convergent.One can conclude that the dynamic equations of NS-PIM have anti-ill solving ability compared with FEM under the same calculation conditions.

4 Conclusion

In this paper,a mesh-free method which is called NS-PIM is formulated for rigid-flexible coupled dynamic analysis of beams.We conclude the paper with the following points.

(1)Through the smoothing operation based on nodes,the requirements on the consistence of the displacement functions are further weakened.The solutions of fourth order boundary value and initial value problems can now be approximated using linear shape functions in NS-PIM.(2)For static problems,NS-PIM has high accuracy even if using linear interpolation.It provides a unique upper bound solution for the deflection.

(3)For the problem of rotating flexible beam,when the deformation is small,the results of AMM,FEM,and NS-PIM are close.Even if increasing the number of modes,the accuracy is not improved for AMM,but the computational efficiency reduces significantly.When the deformation is large,AMM is not applicable,but FEM and NS-PIM is still applicable.The soft characteristic of NS-PIM is more obvious than the situation of small deformation.What’s more,the dynamic equations of NS-PIM have anti-ill solving ability compared with FEM under the same calculation conditions.In other words,the NSPIM is more stable than FEM.

(4)Because of only displacement as independent variable at per node,NS-PIM the total DOFs is half of that of FEM model using the same number of nodes,and hence NS-PIM computes faster than FEM.

Fig.12 The rotation angle and the tip responses of the rotating beam(τ0=10 N ·m):a the rotation angle;b tip transverse deformation from 0 to 15 s;c tip transverse deformation after 10 s

Fig.13 The tip responses of the rotating beam(τ0=50 N ·m)

AcknowledgementsThe authors are grateful for the support from the National Natural Science Foundation of China(Grants 11272155,11132007,and 11502113),the Fundamental Research Funds for Central Universities(Grant 30917011103),and the China Scholarship Council for one year study at the University of Cincinnati.

1.Yoo,H.H.,Ryan,R.R.,Scott,R.A.:Dynamics of flexible beams undergoing overall motions.J.Sound Vib.181,261–278(1995)

2.Yoo,H.H.,Shin,S.H.:Vibration analysis of rotating cantilever beams.J.Sound Vib.212,807–828(1998)

3.Li,L.,Zhang,D.G.,Zhu,W.D.:Free vibration analysis of a rotating hub-functionally graded material beam system with the dynamic stiffening effect.J.Sound Vib.333,1526–1541(2014)

4.Chung,J.,Yoo,H.H.:Dynamic analysis of a rotating cantilever beam by using the finite element method.J.Sound Vib.249,147–164(2002)

5.Du,H.,Lira,M.K.,Liew,K.M.:A nonlinear finite element model for dynamics of flexible manipulators.Mech.Mach.Theory 31,1109–1119(1996)

6.Sanborn,G.G.,Shabana,A.A.:A rational finite element method based on the absolute nodal coordinate formulation.Nonlinear Dyn.58,565–572(2009)

7.Liu,G.R.,Quek,S.S.:Finite Element Method:A Practical Course,2nd edn.Butter worth-Heinemann,Burlington(2013)

8.Liu,G.R.,Gu,Y.T.:An Introduction to Meshfree Methods and Their Programming.Springer,Dordrecht(2005)

9.Sanborn,G.G.,Shabana,A.A.:On the integration of computer aided design and analysis using the finite element absolute nodal coordinate formulation.Multibody Syst.Dyn.22,181–197(2009)

10.Sugiyama,H.,Gerstmayr,J.,Shabana,A.A.:Deformation modes in the finite element absolute nodal coordinate formulation.J.Sound Vib.298,1129–1149(2006)

11.Lan,P.,Shabana,A.A.:Integration of B-spline geometry and ANCF finite element analysis.Nonlinear Dyn.61,193–206(2010)

12.Liu,Y.N.,Sun,L.,Liu,Y.H.,et al.:Multi-scale B-spline method for 2-D elastic problems.Appl.Math.Model.35,3685–3697(2011)

13.Lucy,L.B.:A numerical approach to testing of the fission hypothesis.Astron.J.8,1013–1024(1977)

14.Liu,G.R.,Liu,M.B.:Smoothed Particle Hydrodynamics:A Meshfree Practical Method.World Scientific,Singapore(2003)

15.Monaghan,J.J.:An introduction to SPH.Comput.Phys.Commun.48,89–96(1998)

16.Belytschko,Y.,Lu,Y.Y.,Gu,L.:Element-free Galerkin methods.Int.J.Numer.Methods Eng.37,229–256(1994)

17.Liu,W.K.,Jun,S.,Zhang,Y.E.:Reproducing kernel particle methods.Int.J.Numer.Methods Eng.20,1081–1106(1995)

18.Atluri,S.N.,Zhu,T.:A new meshless local Petrov–Galerkin(MLPG)approach in computational mechanics.Comput.Mech.22,117–127(1998)

19.Liu,G.R.,Gu,Y.T.:A point interpolation method for two dimensional solids.Int.J.Numer.Methods Eng.50,937–951(2001)

20.Liu,G.R.,Dai,K.Y.,Lim,K.M.,et al.:A point interpolation mesh free method for static and frequency analysis of two-dimensional piezoelectric structures.Comput.Mech.29,510–519(2002)

21.Liu,G.R.,Zhang,G.Y.,Gu,Y.T.,et al.:A meshfree radial point interpolation method(RPIM)for three-dimensional solids.Comput.Mech.36,421–430(2005)

22.Wang,J.G.,Liu,G.R.:A point interpolation meshless method based on radial basis functions.Int.J.Numer.Methods Eng.54,1623–1648(2002)

23.Liu,G.R.,Zhang,G.Y.,Dai,K.Y.,et al.:A linearly conforming point interpolation method(LC-PIM)for 2-D solid mechanics problems.Int.J.Comput.Methods 2,645–665(2005)

24.Liu,G.R.,Zhang,G.Y.:Upper bound solutions to elasticity problems:a unique property of the linearly conforming point interpolation method(LC-PIM).Int.J.Numer.Methods Eng.74,1128–1161(2008)

25.Zhang,G.Y.,Liu,G.R.,Wang,Y.Y.,et al.:A linearly conforming point interpolation method(LC-PIM)for three-dimensional elasticity problems.Int.J.Numer.Methods Eng.72,1524–1543(2007)

26.Wu,S.C.,Liu,G.R.,Zhang,H.O.,et al.:A node-based smoothed point interpolation method(NS-PIM)for thermoelastic problems with solution bounds.Int.J.Heat Mass Transf.52,1464–1471(2009)

27.Liu,G.R.:Meshfree Methods:Moving Beyond the Finite Element Method.CRC Press,Boca Raton(2002)

28.Cui,X.Y.,Liu,G.R.,Li,G.Y.,et al.:A rotation free formulation for static and free vibration analysis of thin beams using gradient smoothing technique.CMES 28,217–229(2008)

29.Liu,G.R.,Han,X.:Computational Inverse Techniques in Nondestructive Evaluation.CRC Press,Boca Raton(2003)

30.Liu,G.R.:Meshfree Methods:Moving Beyond the Finite Element Method,2nd edn.CRC Press,Boca Raton(2010)

31.Wu,S.C.,Haug,E.J.:Geometric non-linear substructuring for dynamics of flexible mechanical system.Int.J.Numer.Methods Eng.26,2211–2226(1988)