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.

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

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 [711]. 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 Figure 1. The initial shape of the beam can be described by the undeformed configuration , and the configuration of the deformed beam is described by . is the reference configuration, which is the common configuration of the strait beam generally.

Let and denote the centerline of and , respectively, and define as the running length coordinate of the material point along the center-line of the beam. (for ) at material point are orthogonal unit basis vectors. Among these vectors, is tangent to the center-line of , while and are fixed in the cross section of the beam. at material point are another orthogonal unit basis vectors. Unlike , because of the deformation of the beam, 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 , respectively. The curvature vector contains three components describing twist and bending relative to ; thus, it can be expressed by the following equation:

The strain vector is shown as follows:where is the axial strain measure of the beam and and are the transverse shear strain measures of the beam.

Considering small deformations, the spatial rate of change of the centerline position can be expressed by [8]where .

And the spatial rate of change of the basis vectors can be expressed by

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:where ; and are the linear velocity and angular velocity of the cross section, respectively; and are the linear momentum and angular momentum per unit length, which is associated with and ; and and are the internal force vector and moment vector, respectively. The unit vector is measured relative to .

The above two equations can describe the dynamic response of the beam’s cross section located at any arbitrary material point while undergoing external distributed force and moment vector per unit length.

The constitutive law associating internal force vector and moment vector with curvature vector and strain vector can be expressed by the following equation:where is the constitutive matrix and and are the initial curvature and strain vectors, respectively.

The intrinsic formula takes , , , and as the variables. There are another three equations implemented as follows to get a close form:where and (for ); is the mass per unit length; and are offsets from the reference line of the cross-sectional mass centroid; and , , and 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 and . Then, the expression of equation (9) can be simplified as follows:where , .

Substituting equation (10) into equation (3) will lead to

Then substituting equation (10) into equation (4) will give

Since , equation (12) can be expressed by a simple form

Then, we have the following four equations:

For isotropic materials, the constitutive law can be decoupled as

For simple expression, express by the following form:then the constitutive law can be expressed in the form

By substituting equation (17) into the four equations, the following equations are obtained:

Now we re-express equation (18) in another form as follows:where

3. Discretization for the Intrinsic Beam Formula

We use a series of basis functions to interpolating , , , and , , where , , , and . For simplicity of the expressions, we use to replace . Then, the expressions of the interpolation will be

To each interior point located at , which is also the th influence domain, we have

As to , if each influence domain contains 3 grid points and , for example, then, the number of the interior point will be . The following expressions will be satisfied as

Writing equation (27) in a matrix form, the following equation is obtained:

Since the expression of is known to us, the unknown can be determined by

Substituting equation (29) into equation (22) will lead to

S(x) is denoted by the following equation:

Then, within the th influence domain can be expressed by

Following the same process, we will have

Align all of the variables in the following form:where denotes generalized grid point displacement for the th influence domain, , , and are the generalized grid point displacements of the th, th, and th grid point, respectively, andwhere the superscript denotes the th grid point. Or, equation (35) can be expressed bywhere , , , and .

Then within to can be expressed bywhere the shape function is a matrix. , , and are all diagonal matrices, and the corresponding expressions arewhere the matrix is a identity matrix.

Substituting equation (37) into equation (20) will result in

As to the interpolation meshless method, equation (39) must be satisfied at each grid point ; thus, the following equations are obtained:where

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 function as follows:then we will havewhere , , , and . Denote , , , and .

Substituting , into equation (41), the following is obtained:

To each tip point, we consider that each influence domain contains 2 grid points and . 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:

DenotingAll of the variables at can be expressed by

And all of the variables at can be expressed by

Then, the generalized response vector at each tip can be expressed bywhere , , , and are all diagonal matrices; thus, and . The expressions of and are listed in the following:

Substituting equation (49) into equation (20) and replacing as and , respectively, will givewhere has the same form as shown in equation (41).

By the same steps, we depart and as follows:where , , , and , . Denote , , , and .

Substituting , , into equation (41) will give the same form as shown in equation (42).


To equations (40) and (51), we expand , 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:where

The expression of has been illustrated as equation (44), and 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

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 of , 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.

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 Figure 3. In this figure, x1 = 0 to x1 = 75 represent the coordinates of the selected point in the natural coordinate system of the beam.

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 Figure 3.

Figure 3 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 Figure 3, 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.

Figure 4 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 and of the displacement components are symmetrical due to the symmetry of the load. This is consistent with physical phenomenon.

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.

Data Availability

The geometrically nonlinear analysis for elastic beam.zip data used to support the findings of this study have been deposited in the “figshare” repository (https://doi.org/10.6084/m9.figshare.7558421.v1).

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This study was funded by the Natural Science Foundation of Jiangsu Province, China (Grant no. BK20160782), and Fundamental Research Funds for Central University (Grant no. NS2016098).