A Variable-Length Beam Element Incorporating the Effect of Spinning

Abstract

This paper proposes a novel variable-length beam element that takes into account the effect of beam spinning. This is the first such beam element of variable length based on the absolute nodal coordinate formulation. In addition to the position and slope vectors, the angles of rotation around the element axis of cross sections that contain two nodes are introduced into the element coordinates to describe the spinning of an assumed Euler-Bernoulli beam with circular cross section. The material coordinates of the two nodes are also included in the arbitrary Lagrangian-Eulerian description used for previous element to describe the varying element length caused by mass transportation at the boundaries. The proposed element facilitates convenient and effective numerical modeling of the dynamics of a circular-cross-section beam with transportation boundaries and that is spinning. Numerical examples demonstrate that the proposed element can describe the dynamic behavior of a circular-cross-section beam effectively. In the field of engineering, this novel element could be used in the dynamic analysis of drill stems, the slender workpiece of a cylindrical shaft during the turning process, and the lead screw in a ball screw mechanism.

Keywords:
Beam; circular cross section; absolute nodal coordinate formulation; arbitrary Lagrangian-Eulerian description; variable-length beam-shaft element; spinning

1 INTRODUCTION

In the late 1990s, Shabana (1996) proposed a new method for modeling the dynamics of flexible multi-body systems. Known as the absolute nodal coordinate formulation (ANCF), this method is based on the finite element method (Bonet and Wood, 1997; Khan et al., 2016) and continuum mechanics (Shabana, 2012; Asemi et al., 2014). Unlike conventional formulations (e.g., large rotation vector and floating frame of reference formulations) that use both global and local reference frames, the ANCF uses only a global reference frame. A position vector and its gradients defined in the global frame are selected as the generalized coordinates of each node. This ensures that no restriction is imposed on the amount of rotation within the finite element, thereby describing the motion of a flexible multi-body system more exactly (Shabana, 1997). As such, the ANCF has received extensive attention since it first appeared and has been developed considerably over the past 20 years. Various types of ANCF-based finite elements have been developed to apply this method to the dynamic analysis of flexible multi-body systems including structures of beams (Shabana, 1997; Sugiyama and Suda, 2007), plates (Dmitrochenko and Pogorelov, 2003; Dufva and Shabana, 2005), shells (Mikkola and Shabana, 2003), and three-dimensional (3D) solids (Kubler et al., 2003).

In the research area of ANCF-based beam elements, Shabana (1997) first developed a beam element based on the assumption of planar strain. Omar et al. (2001) then developed a two-dimensional (2D) shear-deformable beam element. Following on from 2D beam elements, Shabana and Yakoub (2001, 2001) developed a 3D beam element that uses the position vector and its three gradients as node coordinates. This element relaxes the rigid-section assumption of classical beam theory. It describes cross-sectional deformation using two position gradients, making it a more general model than the Timoshenko beam model. Unfortunately, this element is disadvantageous for the dynamic analysis of very slender beam structures because of the locking problem. Dombrowski (2002) developed a beam element based on the Euler-Bernoulli beam assumption. This element is better at describing slender beams because the number of coordinates of each node is reduced from 12 to 7 compared with Shabana’s 3D beam element. This element can be used to resolve problems with torsional beam dynamics because Tait-Brian angles are used to describe the cross-sectional spinning; however, the mass matrix of the element is no longer constant. Besides, the singularity introduced by Tait-Brian angles prevents this element from being used to describe arbitrary spatial motions. Gerstmayr and Shabana (2006) developed a 3D cable element that abandons the rotational angles of Dombrowski’s beam element and exhibits advantages in the dynamics of 3D slender beams for which the effect of spinning can be neglected.

To study the dynamics of a time-varying cable network, Lu (2007) developed a variable-length ANCF-based cable element by changing the original length of the boundary elements. This element does not strictly satisfy energy conservation, however, it still can be an engineering approximation. Based on Gerstmayr’s cable element, Tang et al. (2011) developed a new variable-length tether element that combines the Lagrangian and Eulerian descriptions. When used to resolve the dynamics of multi-body systems containing time-varying cable structures, way of length varying or mass transportation conditions of boundary elements must be defined in advance. Hence, this element is suitable only for problems involving simple transportation boundaries. Hong and Ren (2011) and Hong et al. (2011) generalized Tang’s cable element and developed a variable-length beam element (VLBE) for a one-dimensional continuum by combining the arbitrary Lagrangian-Eulerian (ALE) description with the ANCF. The material coordinates of the two nodes used to describe the mass transportation are introduced into the system’s dynamic equations as generalized coordinates. As such, this element can be used as an accurate and effective numerical modeling method for the dynamics involving complex transportation boundaries. Du et al. (2015) developed a simpler variable-length cable element that uses the position vectors of the two nodes as generalized coordinates only and uses shape functions of linear interpolating. This element maintains a straight-line configuration and therefore fails to describe the bending deformation. Hence, it is applicable to the dynamic modeling of cables that sustain only tension.

All of the aforementioned variable-length tether/beam elements neglect beam spinning; therefore, they cannot be applied to certain engineering problems such as the dynamic modeling and simulation of drill stems in the drilling process. In this work, we develop a novel variable-length ANCF-based beam element that includes the effect of spinning. Our intention is to resolve the dynamics of a circular-cross-section (CCS) beam of variable length for which the effect of spinning cannot be neglected. We refer to this novel element as a variable-length beam-shaft element (VLBSE) to reflect the introduction of both bending and spinning effects.

In addition to the position and slope vectors, we introduce the angles of rotation around the element axis of cross sections that contain two nodes into the generalized coordinates of this element to describe the spinning motion. This is done on the assumptions of an Euler-Bernoulli beam and a CCS. The material coordinates of the nodes are also included in the ALE description used for previous element to describe the variation of element length caused by mass transportation at the boundaries. The rotational angles and material coordinates appear as unknowns in the dynamic equations along with the position and slope vectors. Hence, the VLBSE provides a convenient and effective numerical modeling method for the dynamics of CCS beams with either simple or complex transportation boundaries and for which spinning cannot be neglected.

The VLBSE can be used in the dynamic analyses of (i) the drill stem of a planetary soil-drilling sampler in a deep-space exploration mission (Tian et al., 2012; Deng et al., 2014; Yang et al., 2014), (ii) the slender workpiece of a cylinder shaft in a turning process (Lv et al., 2015; Han et al., 2012), and (iii) slender rods in linear motion elements (Katz, 2001; Huang and Wang, 2014).

In the present work on the VLBSE, we begin by introducing the variables that are used to describe the motion of the element, the generalized coordinates of the element, and the corresponding shape-function matrix. We then address the axial strain, bending curvature, and change rate of rotational angle of cross section, all of which describe the deformation of the element. The mass matrix and generalized force matrix of the element are given by calculating the virtual work of the inertial and active forces. Based on these matrices, a mathematical description of the dynamical model of a CCS beam that has been discretized into elements is obtained according to the general equation of dynamics. Finally, the novel element is verified by three series of typical dynamical problems concerning an axially moving beam/string, a torsion shaft with one end fixed and the other end free, and a rotating shaft subjected to an axially moving and rotating load.

2 MOTION AND DEFORMATION OF THE ELEMENT

The following assumptions are made about the beam before establishing the VLBSE.

• (a) The beam cross sections are rigid and remain perpendicular to the beam axis during deformation because of the slenderness of the beam.

• (b) Any bending-induced rotational inertia effects relative to the neutral axis of the beam cross sections are neglected.

Together, assumptions (a) and (b) constitute the Euler-Bernoulli beam assumption.

• (c) The beam cross sections are circular and uniform.

Assumption (c) is consistent with the Euler-Bernoulli beam assumption because a non-CCS would distort if the beam were subjected to torsion.

The motion of an element can be described by two features: translation of an arbitrary point on the element axis, and rotation of the corresponding cross section. Translation of a point on the element axis can be described by the position vector r of the point. Rotation of a cross section can be described using only the third Tait-Brian angle instead of the direction cosine matrix of the element local reference frame with respect to the global reference frame in the general case with the Euler-Bernoulli beam assumption according to Dombrowski’s work (Dombrowski, 2002). This angle is written as ϑ, the physical meaning of which is the rotational angle of the cross section around the element axis.

Hence, the motion of the VLBSE is described by four variables that can be expressed as a generalized position vector of an arbitrary point on the element axis:

$$\eta ={\left[\begin{array}{cc}{r}_{}^T& \theta \end{array}\right]}_{}^T={\left[\begin{array}{cccc}{r}_1^{}& r_2^{}& r_3^{}& \theta \end{array}\right]}_{}^T$$ (1)

A schematic of the VLBSE is shown in Figure 1. In addition to r and ϑ, r ' is the slope vector of the arbitrary point on the element axis. The derivative is with respect to p, i.e., r ' = ∂r/∂p, where p is the material coordinate of the point, which is the arc length along the unstrained beam from the reference point. Terms r ₁, r ' ₁, ϑ ₁, p ₁, r ₂, r ' ₂, ϑ ₂, and p ₂ are the corresponding parameters of the two nodes. Function f(p, t) is the distributed external load per unit length on the element, and F _C is the concentrated force on the element. The concentrated torque T is not marked in the figure.

Figure 1:
Schematic of the VLBSE.

A control volume for the element can be expressed as [p ₁(t), p ₂ (t)] using an interval representation. In the ALE framework, both p ₁ and p ₂ are time-dependent functions, which means that there is mass transportation across both element boundaries and that an element can slide and stretch arbitrarily on the beam. However, if p ₁ and p ₂ remain constant, there is no mass transportation across the element boundaries; the VLBSE becomes a fixed-length beam-shaft element in the Lagrangian framework. The material coordinate p can be normalized into a local coordinate s that satisfies

$$s=\frac{2p-p_1^{}(t)-p_2^{}(t)}{p_2^{}(t)-p_1^{}(t)}$$ (2)

From equation (2), the local coordinates of the first and second nodes are −1 and 1, respectively. The length of the unstrained element is

$$L_e^{}(t)=p_2^{}(t)-p_1^{}(t)$$ (3)

Of the four components of the generalized position vector, the first three corresponding to the position vector are interpolated using third order polynomials to ensure continuity of the position and slope vectors at the nodes. The position vector can then be written as

$$r=N_1^{}{r}_1^{}+N_2^{}r{\text{'}}_1^{}+N_3^{}{r}_2^{}+N_4^{}r{\text{'}}_2^{}$$ (4)

where

$$\begin{array}{l}{N}_1^{}\text{=}\frac{1}{4}{\left(s-1\right)}_{}^2\left(2+s\right),N_2^{}\text{=}\frac{L_e^{}}{8}{\left(s-1\right)}_{}^2\left(s+1\right),\\ N_3^{}\text{=}\frac{1}{4}{\left(s+1\right)}_{}^2\left(2-s\right),N_4^{}\text{=}\frac{L_e^{}}{8}{\left(s+1\right)}_{}^2\left(s-1\right).\end{array}$$ (5)

The fourth component of the generalized position vector corresponding to the rotational angle is interpolated using a linear polynomial:

$$\theta =H_1^{}{\theta }_1^{}+H_2^{}{\theta }_2^{},$$ (6)

where

$$H_1^{}=\frac{1}{2}\left(1-s\right),H_2^{}=\frac{1}{2}\left(1+s\right).$$ (7)

Equations (4) and (6) can be unified as

$$\left[\begin{array}{c}r\\ \theta \end{array}\right]=\left[\begin{array}{cccccc}{N}_1^{}{I}_{3\times 3}^{}& N_2^{}{I}_{3\times 3}^{}& 0_{3\times 1}^{}& N_3^{}{I}_{3\times 3}^{}& N_4^{}{I}_{3\times 3}^{}& 0_{3\times 1}^{}\\ 0_{1\times 3}^{}& 0_{1\times 3}^{}& H_1^{}& 0_{1\times 3}^{}& 0_{1\times 3}^{}& H_2^{}\end{array}\right]\left[\begin{array}{c}{r}_1^{}\\ r{\text{'}}_1^{}\\ {\theta }_1^{}\\ r_2^{}\\ r{\text{'}}_2^{}\\ {\theta }_2^{}\end{array}\right],$$ (8)

which can be simplified to

$$\eta \left(p,t\right)=S_e^{}\left(p,p_1^{}(t),p_2^{}(t)\right)q_e^{}(t)=\left[S_{e1}^{}{S}_{e2}^{}\right]\left[\begin{array}{l}{q}_{e1}^{}\\ q_{e2}^{}\end{array}\right]=S_{e1}^{}{q}_{e1}^{}+S_{e2}^{}{q}_{e2}^{}.$$ (9)

To partition the matrix S _e, two submatrices

$$S_{et}^{}=\left[S_{et1}^{}{S}_{et2}^{}\right]=\left[\begin{array}{cccccc}{N}_1^{}{I}_{3\times 3}^{}& N_2^{}{I}_{3\times 3}^{}& 0_{3\times 1}^{}& N_3^{}{I}_{3\times 3}^{}& N_4^{}{I}_{3\times 3}^{}& 0_{3\times 1}^{}\end{array}\right],$$ (10)

$$S_{er}^{}=\left[S_{e\text{r}1}^{}{S}_{er2}^{}\right]=\left[\begin{array}{cccccc}{0}_{1\times 3}^{}& 0_{1\times 3}^{}& H_1^{}& 0_{1\times 3}^{}& 0_{1\times 3}^{}& H_2^{}\end{array}\right],$$ (11)

are defined to represent the first three lines and the fourth line, respectively, leading to the following equations:

$$S_e^{}=\left[\begin{array}{c}{S}_{et}^{}\\ S_{er}^{}\end{array}\right]r=S_{et}^{}{q}_e^{}\theta =S_{er}^{}{q}_e^{}.$$ (12)

The first and second derivatives of the generalized position vector with respect to time (i.e., the generalized velocity and acceleration, respectively) can be written as

$$\dot{\eta }=S\dot{q},$$ (13)

$$\ddot{\eta }=S\ddot{q}+{\ddot{\eta }}_p^{},$$ (14)

where

$$S=\left[S_{e1}^{}\frac{\partial S_e^{}}{\partial p_1^{}}{q}_e^{}{S}_{e2}^{}\frac{\partial S_e^{}}{\partial p_2^{}}{q}_e^{}\right]$$ (15)

is the shape-function matrix of the element,

$$q={\left[\begin{array}{cccc}{q}_{e1}^T& p_1^{}& q_{e2}^T& p_2^{}\end{array}\right]}_{}^T={\left[\begin{array}{cccccccc}{r}_1^T& r{\text{'}}_1^T& {\theta }_1^{}& p_1^{}& r_2^T& r{\text{'}}_2^T& {\theta }_2^{}& p_2^{}\end{array}\right]}_{}^T$$ (16)

represents the generalized coordinates of the element, and

$${\ddot{\eta }}_p^{}=2\left(\frac{\partial S_e^{}}{\partial p_1^{}}{\dot{p}}_1^{}+\frac{\partial S_e^{}}{\partial p_2^{}}{\dot{p}}_2^{}\right){\dot{q}}_e^{}+\left(\frac{{\partial }^2{S}_e^{}}{\partial p_1^2}{\dot{p}}_1^2+2\frac{{\partial }^2{S}_e^{}}{\partial p_1^{}\partial p_2^{}}{\dot{p}}_1^{}{\dot{p}}_2^{}+\frac{{\partial }^2{S}_e^{}}{\partial p_2^2}{\dot{p}}_2^2\right)q_e^{}$$ (17)

is an additional generalized acceleration due to mass transportation across the element boundaries. This term vanishes if there is no mass transportation, in which case both p ₁ and p ₂ are constants.

Matrix S is partitioned in the same way as S _e (i.e., two submatrices S _t and S _r represent the first three lines and the fourth line, respectively), leading to the following equations:

$$S=\left[\begin{array}{c}{S}_t^{}\\ S_r^{}\end{array}\right]\dot{r}=S_t^{}\dot{q}\dot{\theta }=S_r^{}\dot{q},$$ (18)

where

$$S_t^{}=\left[S_{et1}^{}\frac{\partial S_{et}^{}}{\partial p_1^{}}{q}_e^{}{S}_{et2}^{}\frac{\partial S_{et}^{}}{\partial p_2^{}}{q}_e^{}\right],$$ (19)

$$S_r^{}=\left[S_{er1}^{}\frac{\partial S_{er}^{}}{\partial p_1^{}}{q}_e^{}{S}_{er2}^{}\frac{\partial S_{er}^{}}{\partial p_2^{}}{q}_e^{}\right].$$ (20)

Because the VLBSE is part of the beam, the deformation description of the beam also applies to the element. Due to the Euler-Bernoulli beam assumption, only one the nine components of the Green-Lagrange strain tensor is nonzero, which is the axial strain, besides, the bending curvature of the beam axis and the change rate of the rotational angle of the cross section are also taken into account to describe the beam deformation. These three components are given by

$$\epsilon =\frac{1}{2}\left(r{\text{'}}_{}^{T}r\text{'}-1\right),$$ (21)

$$\kappa =\frac{\left|r\text{'}\times r\text{'}\text{'}\right|}{{\left|r\text{'}\right|}_{}^3},$$ (22)

$$\phi =\frac{\partial \theta }{\partial p}=\theta \text{'},$$ (23)

where r " is the second derivative of r with respect to p, i.e., r " = ∂ ² r / ∂p ². The axial strain measures the extension of the beam axis, the bending curvature measures the change of orientation of the tangent vector of beam axis, and the change rate of the rotational angle measures the twist of beam. The detailed expressions of r ', r " and ϑ' can be written as:

$$r\text{'}=\frac{\partial r}{\partial p}=\frac{\partial S_{et}^{}}{\partial p}{q}_e^{}=\frac{\partial s}{\partial p}\frac{\partial S_{et}^{}}{\partial s}{q}_e^{},$$ (24)

$$r\text{'}\text{'}=\frac{{\partial }^{2}r}{\partial p_{}^2}=\frac{{\partial }^2{S}_{et}^{}}{\partial p_{}^2}{q}_e^{}={\left(\frac{\partial s}{\partial p}\right)}_{}^2\frac{{\partial }^2{S}_{et}^{}}{\partial s_{}^2}{q}_e^{},$$ (25)

$$\theta \text{'}=\frac{\partial \theta }{\partial p}=\frac{\partial S_{er}^{}}{\partial p}{q}_e^{}=\frac{\partial s}{\partial p}\frac{\partial S_{er}^{}}{\partial s}{q}_e^{}.$$ (26)

3 MASS MATRIX AND GENERALIZED FORCE MATRIX OF THE ELEMENT

The mass matrix and generalized force matrix of the VLBSE need to be given first to obtain a mathematical description of the dynamical model of the CCS beam. These can be obtained by calculating the virtual work of the inertial and active forces on the element. The active forces can be classified into viscoelastic and external forces: the latter include distributed external loads, concentrated forces, and concentrated torques.

The virtual work due to inertial force arises from both translational acceleration of material particles on the element axis and rotational angular acceleration of the cross sections:

$$\delta W_i^{}={\displaystyle {\int }_{p_1^{}}^{p_2^{}}\left[\delta r_{}^T\left(-\rho A\ddot{r}\right)dp+\delta \theta \left(-\rho I\ddot{\theta }\right)dp\right]}=-\delta q_{}^T\left(M_{ele}^{}\ddot{q}+Q_p^{}\right),$$ (27)

where ϱ, A, and I are material density, area, and polar moment of inertia, respectively, of the element cross section. Terms M _ele and Q _p are the mass matrix and the matrix of the additional generalized inertial force caused by mass transportation, respectively. Under the assumption of uniform material density, these can be written as

$$M_{ele}^{}=\rho \frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{S}_{}^T\left[\begin{array}{cc}A{I}_{3\times 3}^{}& 0_{3\times 1}^{}\\ 0_{1\times 3}^{}& I\end{array}\right]S}ds,$$ (28)

$$Q_p^{}=\rho \frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{S}_{}^T\left[\begin{array}{cc}A{I}_{3\times 3}^{}& 0_{3\times 1}^{}\\ 0_{1\times 3}^{}& I\end{array}\right]{\ddot{\eta }}_p^{}}ds.$$ (29)

A block diagonal matrix appears in the integrands of these two matrices, which is unlike existing elements in the literature. The virtual work of viscoelastic forces arises from tension, bending, and torsion of the element:

$$\delta W_s^{}=\delta W_{s\epsilon }^{}+\delta W_{s\kappa }^{}+\delta W_{s\phi }^{},$$ (30)

where _^δWsε , _^δWsκ , and _^δWsφ are the virtual-work contributions of tension, bending, and torsion, respectively. Specifically,

$$\begin{array}{c}\delta W_{s\epsilon }^{}=-{\displaystyle {\int }_{p_1^{}}^{p_2^{}}EA\left(\epsilon +c\dot{\epsilon }\right)\delta \epsilon dp}=-{\displaystyle {\int }_{-1}^{1}EA\left(\epsilon +c\dot{\epsilon }\right)\frac{\partial \epsilon }{\partial q}\delta q\frac{\partial p}{\partial s}ds}\\ =-\delta q_{}^T\frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{\left(\frac{\partial \epsilon }{\partial q}\right)}_{}^{T}EA\left(\epsilon +c\dot{\epsilon }\right)ds}=-\delta q_{}^T{Q}_{s\epsilon }^{}\end{array}$$ , (31)

where E and c are the elasticity modulus and damping coefficient, respectively. Term Q _sε is the generalized viscoelastic force matrix resulting from the tension of the element:

$$Q_{s\epsilon }^{}=\frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{\left(\frac{\partial \epsilon }{\partial q}\right)}_{}^{T}EA\left(\epsilon +c\dot{\epsilon }\right)ds},$$ (32)

where,

$$\dot{\epsilon }=r{\text{'}}_{}^T\dot{r}\text{'},$$ (33)

$$\dot{r}\text{'}=\frac{d}{dt}\left(\frac{\partial r}{\partial p}\right)=\frac{\partial }{\partial p}\left(\frac{dr}{dt}\right)=\frac{\partial S_t^{}}{\partial p}\dot{q},$$ (34)

and

$$\frac{\partial \epsilon }{\partial q}=r{\text{'}}_{}^T\frac{\partial r\text{'}}{\partial q},$$ (35)

$$\frac{\partial r\text{'}}{\partial q}=\frac{{\partial }^{2}r}{\partial p\partial q}=\frac{{\partial }^2\left(S_{et}^{}{q}_e^{}\right)}{\partial p\partial q}.$$ (36)

Similarly,

$$\delta W_{s\kappa }^{}=-{\displaystyle {\int }_{p_1^{}}^{p_2^{}}EJ\left(\kappa +c\dot{\kappa }\right)\delta \kappa dp}=-\delta q_{}^T{Q}_{s\kappa }^{},$$ (37)

$$\delta W_{s\phi }^{}=-{\displaystyle {\int }_{p_1^{}}^{p_2^{}}GI\left(\phi +c\dot{\phi }\right)\delta \phi dp}=-\delta q_{}^T{Q}_{s\phi }^{},$$ (38)

where J is the moment of inertia of the element cross section and G is the shear modulus. Terms Q _sκ and Q _sφ are the generalized viscoelastic force matrices of element resulting from bending and torsion, respectively:

$$Q_{s\kappa }^{}=\frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{\left(\frac{\partial \kappa }{\partial q}\right)}_{}^{T}EJ\left(\kappa +c\dot{\kappa }\right)ds},$$ (39)

$$Q_{s\phi }^{}=\frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{\left(\frac{\partial \phi }{\partial q}\right)}_{}^{T}GI\left(\phi +c\dot{\phi }\right)ds}.$$ (40)

The virtual work of the viscoelastic forces can then be written as

$$\delta W_s^{}=-\delta q_{}^T{Q}_{s\epsilon }^{}-\delta q_{}^T{Q}_{s\kappa }^{}-\delta q_{}^T{Q}_{s\phi }^{}=-\delta q_{}^T{Q}_s^{},$$ (41)

where Q _s is the generalized viscoelastic force matrix and can be written as

$$Q_s^{}=Q_{s\epsilon }^{}+Q_{s\kappa }^{}+Q_{s\phi }^{}.$$ (42)

In particular, when beam torsion or bending is negligible, Q _sφ or Q _sκ can be omitted accordingly. The virtual work of a distributed external load can be written as

$$\delta W_f^{}={\displaystyle {\int }_{p_1^{}}^{p_2^{}}\delta r_{}^{T}fdp}=\delta q_{}^T{Q}_f^{},$$ (43)

where Q _f is the generalized external force matrix of the distributed external load:

$$Q_f^{}=\frac{\partial p}{\partial s}{\displaystyle {\int }_{-1}^1{S}_t^{T}fds}.$$ (44)

The generalized external force matrix of a concentrated force and a concentrated torque can be obtained similarly:

$$Q_F^{}=S_t^T{F}_C^{},$$ (45)

$$Q_T^{}=S_r^{T}T.$$ (46)

The total virtual work of an element can then be obtained as

$$\delta W_{ele}^{}=-\delta q_{}^T\left(M_{ele}^{}\ddot{q}+Q_{ele}^{}\right),$$ (47)

where Q _ele is the generalized force matrix of the element:

$$Q_{ele}^{}=Q_p^{}+Q_s^{}-Q_f^{}-Q_F^{}-Q_T^{}.$$ (48)

To resolve the dynamics of a CCS beam with mass transportation, we divide the beam into several VLBSEs that are connected by nodes. Suppose that the beam consists of N elements; therefore, there are N+1 nodes. The generalized coordinates of the beam are

$$q={\left[r_1^{T}r{\text{'}}_1^T{\theta }_1^{}{p}_1^{}{r}_2^{T}r{\text{'}}_2^T{\theta }_2^{}{p}_2^{}\dots r_{N+1}^{T}r{\text{'}}_{N+1}^T{\theta }_{N+1}^{}{p}_{N+1}^{}\right]}_{}^T.$$ (49)

In general, these coordinates are not completely independent, for example: in the case that the beam is connected at node i with a fixed point C in the inertial reference frame by means of a spherical joint, the constraint equation is given as

$$r_{\text{i}}^{}=r_C^{},$$ (50)

in another case that all of the elements of the beam are of equal length, the constraint equations are given as

$$p_i^{}=\frac{N-i+1}{N}{p}_1^{}+\frac{i-1}{N}{p}_{N+1}^{}\left(i=2,\mathrm{3,...,}N\right).$$ (51)

More constraint equations of different joints are given by García-Vallejo et al. (2003). These constraint equations can be expressed uniformly in the form of an algebraic equation as

$$C\left(q,t\right)=0.$$ (52)

By introducing a Lagrange multiplier λ, the dynamical equation of the beam can be obtained as

$$M\ddot{q}+Q+C_q^T\lambda =0,$$ (53)

where M and Q are the mass matrix and generalized force matrix, respectively, of the beam. These derive from M _ele and Q _ele of each element, with C _q being the Jacobian matrix of the constraint.

The dynamical equation (53) and the constraint equation (52) constitute a set of differential algebraic equations (DAEs) and should be solved synchronously. In this paper, such DAEs are integrated using the classical differential algebraic system solver DASSL (Brenan et al., 1996), whose core algorithm is the backward differential formulation, i.e., BDF.

4 NUMERICAL EXAMPLES AND DISCUSSION

The VLBSE is verified in this section. The first example is used to show that this element retains the ability of the VLBE to resolve the dynamics of a CCS beam with mass transportation in the case without spinning. The second example is used to demonstrate the ability of this novel element to resolve purely torsional dynamics. Synthetic verification of the VLBSE is carried out in the third example. Finally, potential engineering applications of the VLBSE are discussed.

4.1 Free Flexural Vibration Analysis of an Axially Moving Beam

When a beam moves continuously along its own axis, as shown in Figure 2, its transverse vibrational characteristics are affected significantly by the global axial motion. In Figure 2, OXZ represents the inertial reference frame, and x is the coordinate in the X direction. The target object in this example is the beam section between the two supports. The beam is in continuous motion; hence, the deflection of cross section should be seen as a function of the field coordinate x and time t in the framework of Eulerian description (Liu et al., 2011).

Figure 2:
Axially moving beam.

Suppose that the beam cross section is uniform and circular. The length of the beam section between the supports is l, its cross-sectional area is A, its material density is ϱ, its flexural stiffness is EJ, and its deflection is expressed as w(x, t). The beam travels at a constant speed v ₀ in the X direction under the action of tension F. The dynamic equation of free flexural vibration can be written as (Liu, 2012; Wickert and Mote, 1990)

$$\rho A\frac{{\partial }^{2}w}{\partial t_{}^2}+2\rho A{v}_0^{}\frac{{\partial }^{2}w}{\partial x\partial t}+EJ\frac{{\partial }^{4}w}{\partial x_{}^4}-\left(F-\rho A{v}_0^2\right)\frac{{\partial }^{2}w}{\partial x_{}^2}=0.$$ (54)

To study the natural frequency characteristics more conveniently, the flexural stiffness of the beam can be neglected. In that case, the axially moving beam degenerates to an axially moving string, and the dynamic equation can be simplified to

$$\rho A\frac{{\partial }^{2}w}{\partial t_{}^2}+2\rho A{v}_0^{}\frac{{\partial }^{2}w}{\partial x\partial t}+\left(\rho A{v}_0^2-F\right)\frac{{\partial }^{2}w}{\partial x_{}^2}=0.$$ (55)

The first-order natural frequency of the free transverse vibration of the string can be written as

$$f_1^{}=\frac{1}{2lc}\left(c_{}^2-v_0^2\right),$$ (56)

where $$c=\sqrt{F/(\rho A)}$$ is the axial propagation speed of the string’s transverse vibration.

Conventional ANCF elements in the Lagrangian framework are not suitable for modeling the beam section that flows between the supports in this example, because their target object is the continuum in initial configuration throughout. If the target object is the continuum in a certain fixed region in space, the ALE description degenerates into the Eulerian description. Therefore, the transverse vibrational characteristics of an axially moving beam/string are studied in the following using Hong’s VLBE (Hong and Ren, 2011; Hong et al., 2011) and the VLBSE proposed in this paper, which are both established in the ALE framework.

The geometric, material, and kinematic parameters of the axially moving beam/string are given in Table 1, where r is the cross-sectional radius.

Table1:
Simulation parameters of an axially moving beam/string.

Parameter	l	r	ϱ	E	F	v 0
Value	1 m	0.005 m	2700 kg/m³	70 GPa	1000 N	40 m/s

• (1) ANCF solutions of free transverse vibration of an axially moving string

The string section between the two supports is divided into four elements of equal length. The string begins vibrating freely because of the initial transverse velocities of three intermediary nodes. Transverse displacements of the material point halfway between the supports are compared in Figure 3 for two types of element.

Figure 3:
Comparison of transverse displacements of the material point halfway between supports.

The transverse-displacement-time curves obtained using the two different elements are nearly coincident. According to equation (56), the first-order natural frequency of the transverse vibration of the string should be 22.6857 Hz; the numerical solutions give 22.6723 Hz (VLBSE) and 22.6786 Hz (VLBE) with relative errors that are both less than 0.1%. These results demonstrate that, in the dynamic analysis of a CCS beam subjected to only tension, the VLBSE can describe the continuum dynamics effectively.

• (2) ANCF solutions of free flexural vibration of an axially moving beam

Mesh generation for a beam section between two supports is the same as that for a string section, but now bending is taken into account in the dynamic equations. The beam begins vibrating freely because of the initial transverse velocities of three intermediary nodes. The deflections of the material point halfway between the supports are compared in Figure 4 for two types of element.

Figure 4:
Comparison of deflections of the material point halfway between supports.

The deflection-time curves obtained using the two different elements are nearly coincident. The numerical solutions for the first-order natural frequency of flexural vibration of the beam are 32.2097 Hz (VLBSE) and 32.2425 Hz (VLBE), with a relative error of approximately 0.1%. Furthermore, the first-order natural frequency of transverse vibration of an axially moving beam is clearly greater than that of an axially moving string as the result of flexural stiffness. In the case described by the simulation parameters given in Table 1, the natural frequency increases from 22.7 Hz to 32.2 Hz. The above results demonstrate that, in the dynamic analysis of a CCS beam subjected to both tension and bending, the VLBSE can describe the continuum dynamics effectively.

From the above analysis of the free transverse vibration of an axially moving beam/string, we see that the VLBSE retains the ability of the VLBE to resolve the dynamics of a CCS beam with mass transportation but without spinning. In fact, the VLBSE degenerates to the VLBE if the rotational angles ϑ₁ and ϑ₂ both remain zero.

4.2 Free Torsional Vibration Analysis of a Cylindrical Shaft

The model for free torsional vibration of a cylindrical shaft is shown in Figure 5. The length of the cylindrical shaft is l, the material density is ϱ, and the shear modulus is G. The X axis is set to coincide with the central axis of the shaft, with x the coordinate in the X direction. The rotational angle of the cross section is expressed as ϑ(x,t), and the torque acting on the cross section is T(x, t).

Figure 5:
Model for free torsional vibration of a cylindrical shaft.

This torsional vibration model can be expressed mathematically as a one-dimensional wave equation (Gatti and Ferrari, 2003):

$$\frac{{\partial }_{}^2\theta }{\partial t_{}^2}=c_{}^2\frac{{\partial }_{}^2\theta }{\partial x_{}^2},$$ (57)

where $$c=\sqrt{G/\rho }$$ is the axial propagation speed of the torsional vibration. The torsional natural frequencies _^fi of a shaft that is fixed at x = 0 and free at x = l are

$$f_i^{}=\left(\frac{2i-1}{4}\right)\frac{c}{l}=\frac{2i-1}{4l}\sqrt{\frac{G}{\rho }}\left(i=1,\mathrm{2,...}\right)$$ (58)

The VLBSE in the ALE framework degenerates into a conventional element in the Lagrangian framework if the material coordinates of both nodes are invariant. Therefore, it applies to resolve problems of torsional vibration of the shaft in this example. Free torsional vibration of the shaft are analyzed below using the ANCF. The shaft is divided into different numbers (from 1 to 10) of VLBSEs of equal length. The geometric and material parameters of the shaft are given in Table 2, where r is the cross-sectional radius, and E and μ are the elastic modulus and Poisson’s ratio, respectively, of the material.

Table2:
Simulation parameters of a cylindrical shaft.

Parameter	l	r	ϱ	E	μ
Value	1 m	0.005 m	2700 kg/m³	70 GPa	0.3

Table 3 gives the numerical solutions and their relative errors for the first five orders of torsional natural frequencies of the cylindrical shaft for different numbers of elements; “--” means that the corresponding order of natural frequencies does not exist.

Table3:
Numerical solutions and their relative errors for torsional natural frequencies (Hz) of a shaft.

Order	1st	2nd	3rd	4th	5th
Element number
1	869.9	--	--	--	--
	10.20%
2	809.9	2829.7	--	--	--
	2.60%	19.50%
3	799.9	2609.7	4739.5	--	--
	1.33%	10.20%	20.07%
4	789.9	2509.7	4549.5	6589.3	--
	0.06%	5.97%	15.26%	19.24%
5	789.9	2459.8	4349.6	6529.3	8399.2
	0.06%	3.86%	10.20%	18.15%	18.22%
6	789.9	2429.8	4229.6	6279.4	8489.2
	0.06%	2.60%	7.15%	13.63%	19.48%
7	789.9	2409.8	4159.6	6089.4	8249.2
	0.06%	1.75%	5.38%	10.19%	16.10%
8	789.9	2399.8	4109.6	5959.4	8019.2
	0.06%	1.33%	4.11%	7.84%	12.87%
9	789.9	2389.8	4069.6	5869.4	7839.2
	0.06%	0.91%	3.10%	6.21%	10.33%
10	789.9	2389.8	4049.6	5809.4	7699.2
	0.06%	0.91%	2.60%	5.12%	8.36%
Analytical solution	789.4	2368.3	3947.2	5526.1	7105.0

Table 3 reveals a decrease in shaft torsional stiffness with element number, in that the numerically determined torsional natural frequencies decrease accordingly and converge toward the analytical solutions. The relative error of the first-order torsional natural frequency drops below 1% at element number 4. The relative error of the second-order torsional natural frequency drops below 1% at element number 9. Eventually, the relative error of even the fifth-order torsional natural frequency drops to 8.36% when the shaft is divided into 10 elements.

This example demonstrates that in the dynamic analysis of a CCS beam subjected to only torsion, the VLBSE can describe the continuum dynamics effectively if the mesh is fine enough.

4.3 Dynamic Response of a Spinning Shaft Subjected to an Axially Moving and Rotating Load

A simply supported cylindrical shaft spinning at constant angular velocity is subjected to a constant-magnitude load that remains perpendicular to the shaft axis. The load travels at constant speed along the shaft axis and also rotates at a constant angular speed around the shaft axis, as shown in Figure 6 (a) (Katz, 2001). The length of the shaft is l, the cross-sectional radius is r, the elastic modulus is E, the material density is ϱ, and the angular speed of the shaft is Ω. We establish an inertial frame of reference whose origin is located at one end of the shaft. The X ₃ axis is along the shaft central axis, and the X ₁ and X ₂ axes form a right-handed system together with the X ₃ axis. The load P coincides with the X ₁ axis initially; its axial speed is v and its angular speed is ω. The bending of the shaft is described by the deflections of its central axis in the X ₁ and X ₂ directions, which are expressed as U ₁ and U ₂, respectively, as shown in Figure 6 (b). The deflections U ₁ and U ₂ and their derivatives with respect to time are all zero initially.

Figure 6:
Dynamic model of a spinning shaft subjected to an axially moving and rotating load.

Using complex notation, the dynamic equation of the spinning shaft (Katz, 2001; Katz et al., 1988) can be written as

$$EJ\frac{{\partial }^{4}U}{\partial z_{}^4}+\rho A\frac{{\partial }^{2}U}{\partial t_{}^2}-\rho J\left[\frac{{\partial }_{}^{4}U}{\partial z_{}^2\partial t_{}^2}-i2\Omega \frac{{\partial }_{}^{3}U}{\partial z_{}^2\partial t}\right]=P(z,t)$$ (59)

where J and A are the moment of inertia and the cross-sectional area, respectively, z is the coordinate along the X ₃ axis, and i is the imaginary unit. The load and deflection terms are

$$P(z,t)=P\delta (z-vt)e_{}^{i\omega t},$$ (60)

$$U=U_1^{}+i{U}_2^{}.$$ (61)

The analytical solutions of equation (59) in the general case (ω ≠ Ω) can be given using modal analysis and integral transformation methods (Katz, 2001) as

$\begin{array}{l}{U}_1\left(z,t\right)=\frac{2P}{{\varrho }^{Al}}\sum _{n=1}^{\infty }\frac{{\vartheta }_n}{1+{\beta }_n^2}\\ \times \frac{1}{\left({\omega }_{n1}-{\omega }_{n2}\right)\left(r_1+r_2\right)\left({\omega }_{n1}-r_1\right)\left({\omega }_{n1}+r_2\right)\left({\omega }_{n2}-r_1\right)\left({\omega }_{n2}+r_2\right)}\\ \times \left\{\begin{array}{c}-({\omega }_{n2}-r_1)({\omega }_{n2}+r_2)(r_1+r_2)\sin {\omega }_{n1}t+({\omega }_{n1}-r_1)({\omega }_{n1}+r_2)(r_1+r_2)\sin {\omega }_{n2}t\\ -({\omega }_{n1}-{\omega }_{n2})({\omega }_{n1}+r_2)({\omega }_{n2}+r_2)\sin r_{1}t-({\omega }_{n1}-{\omega }_{n2})({\omega }_{n1}-r_1)({\omega }_{n2}-r_1)\sin r_{2}t\end{array}\right\}\\ \times \sin \frac{n\pi z}{l}\end{array}$ (62)

$\begin{array}{l}{U}_2\left(z,t\right)=\frac{2P}{{\varrho }^{Al}}\sum _{n=1}^{\infty }\frac{{\vartheta }_n}{1+{\beta }_n^2}\\ \times \frac{1}{\left({\omega }_{n1}-{\omega }_{n2}\right)\left(r_1+r_2\right)\left({\omega }_{n1}-r_1\right)\left({\omega }_{n1}+r_2\right)\left({\omega }_{n2}-r_1\right)\left({\omega }_{n2}+r_2\right)}\\ \times \left\{\begin{array}{c}({\omega }_{n2}-r_1)({\omega }_{n2}+r_2)(r_1+r_2)\mathrm{co}\mathrm{s}{\omega }_{n1}t-({\omega }_{n1}-r_1)({\omega }_{n1}+r_2)(r_1+r_2)\mathrm{co}\mathrm{s}{\omega }_{n2}t\\ +({\omega }_{n1}-{\omega }_{n2})({\omega }_{n1}+r_2)({\omega }_{n2}+r_2)\mathrm{co}\mathrm{s}{r}_{1}t-({\omega }_{n1}-{\omega }_{n2})({\omega }_{n1}-r_1)({\omega }_{n2}-r_1)\mathrm{co}\mathrm{s}{r}_{2}t\end{array}\right\}\\ \times \sin \frac{n\pi z}{l}\end{array}$ (63)

where ϑ_n is the frequency that represents the axial motion of the load, ^_βn is the Rayleigh coefficient, _{^ωn 1} and _{^ωn 2} are the positive and negative values, respectively, of the natural frequencies, and r ₁ and r ₂ are the frequencies which describe the moving and rotating load.

To verify the VLBSE further, the deflections of the shaft are calculated here using the ANCF and are compared with the analytical solutions. The shaft is divided into four VLBSEs, with the third node placed at the point at which the load acts. The shaft section on either side of the load is divided into two elements of equal length, thereby causing the length of every element to vary as the load travels along the shaft.

To study the nature of the transient response intuitively, the deflections in the two directions are normalized with respect to U _s into two non-dimensional parameters U ₁ / U _s and U ₂ / U _s, where U _s is the static midpoint deflection of a simply supported beam of the same dimensions subjected to a static load P at its midspan. The geometric, material, and kinematic parameters of the shaft are given in Table 4.

Table4:
Simulation parameters of a sppining shaft.

Parameter	l	β 1	E	ϱ	Ω
Value	1 m	0.02	206 GPa	7700 kg/m³	1 rad/s

Two more non-dimensional parameters a and _^kω , which represent the load speed ratio and the number of rotations that the load completes while traversing the shaft length, respectively, are introduced to describe the axial speed v and the angular speed ω of the load. Their expressions are also given by Katz (2001). In the following, the normalized deflections of the shaft are analyzed for two different sets of simulation parameters (a, _^kω ).

• (1) a=0.5, _^kω =1

Figure 7 shows the configurations of the shaft and its projections on the X ₁ X ₃ and X ₂ X ₃ planes obtained using two methods for loading point at 1/4, 2/4, 3/4, and 4/4 of the shaft length away from the origin. Figure 8 compares the ANCF and analytical solutions for the normalized deflections at the loading point and the midspan of the shaft.

Figure 7:
Configurations of the shaft and its projections for a load at (a) 1/4, (b) 2/4, (c) 3/4, and (d) 4/4 of the shaft length away from the origin (a=0.5, _^kω =1).

Figure 8:
Comparison of normalized deflections at (a) the loading point and (b) the midspan of the shaft (a=0.5, _^kω =1).

When the load is at 85% of the shaft length from the origin, the ANCF solution for the normalized deflection at the loading point in the X ₁ direction reaches a maximum absolute error of −0.0085, whereas the absolute error of the ANCF normalized deflection at the loading point in the X ₂ direction reaches a maximum of −0.0033 at the 78% position. When the load is at 80% from the origin, the ANCF solution for the normalized deflection at the midspan in the X ₁ direction reaches a maximum absolute error of −0.0049, whereas the absolute error of the ANCF normalized deflection at the midspan in the X ₂ direction reaches a maximum of −0.0053 at the 72% position.

• (2) a=1.5, _^kω =0.5

For this set of simulation parameters, the axial speed v has been increased twofold, and the angular speed ω has been increased by half in comparison with the first set of simulation parameters. Figure 9 shows the configurations of the shaft and its projections on the X ₁ X ₃ and X ₂ X ₃ planes obtained using two methods for loading point at 1/4, 2/4, 3/4, and 4/4 of the shaft length away from the origin. Figure 10 compares the ANCF and analytical solutions for the normalized deflections at the loading point and the midspan of the shaft.

Figure 9:
Configurations of the shaft and its projections for a load at (a) 1/4, (b) 2/4, (c) 3/4, and (d) 4/4 of the shaft length away from the origin (a=1.5, _^kω =0.5).

Figure 10:
Comparison of normalized deflections at (a) the loading point and (b) the midspan of the shaft (a=1.5, _^kω =0.5).

When the load is at 61% of the shaft length away from the origin, the absolute error of the ANCF solution for the normalized deflection at the loading point in the X ₁ direction reaches a maximum of 0.0090, whereas the absolute error of the ANCF normalized deflection at the loading point in the X ₂ direction reaches a maximum of 0.0055 at the 84% position. When the load is at 99% of the shaft length, the ANCF solution for the normalized deflection at the midspan in the X ₁ direction reaches a maximum absolute error of 0.0176, whereas the absolute error of the ANCF normalized deflection at the midspan in the X ₂ direction reaches a maximum of 0.0079 at the 100% position.

We see from the above results that in the problem of the dynamic response of a spinning shaft subjected to an axially moving and rotating load, the numerical solutions for the normalized deflections of the shaft obtained using the ANCF are in good agreement with the analytical solutions.

This example demonstrates that in the dynamic analysis of a CCS beam with mass transportation and that is spinning, the VLBSE can describe the continuum dynamics effectively.

4.4 Potential Engineering Applications of the VLBSE

Strictly speaking, a CCS beam with mass transportation at its boundaries and that is spinning is rarely used in engineering. However, in the following, we discuss several special cases in which a slender rod can be treated approximately using the proposed VLBSE.

In the drilling process of a soil-drilling sampler for planetary surface exploration (Tian et al., 2012; Deng et al., 2014; Yang et al., 2014), part of the slender drill stem is in the soil while the residual is outside the soil, and the outer part is likely constrained at certain positions to avoid buckling. To study the dynamic behavior of the drill stem, it could be divided into several VLBSEs using the constraint points and the intersection of the soil interface and drill stem as nodes to establish the dynamic model.

During the turning process of a slender shaft (Lv et al., 2015; Han et al., 2012), complex interactions between the rotating workpiece and the translating cutter may cause excessive vibrations that not only reduce machining accuracy and lead to poor surface quality of the workpiece but also impose restrictions on production efficiency. To study the dynamic characteristics of the shaft, it could be divided into several VLBSEs using the intersection of the shaft and the cutter as one of the nodes. This dynamic model could also reflect the difference in shaft diameter on either side of the cutter due to material removal.

The VLBSE could also be used in the dynamic modeling of slender rods in linear motion elements such as spline shafts in rotary spline screws and rotating ball splines, as well as the lead screw in a ball screw mechanism (Katz, 2001; Huang and Wang, 2014).

5 CONCLUSIONS

In this paper, we proposed a VLBSE based on the ANCF in the ALE framework. This novel element takes spinning into account, which is an effect that is neglected in most variable-length tether/beam elements that exist in the literature.

In addition to the position and slope vectors, the angles of rotation around the element axis of cross sections that contain two nodes are cast as generalized coordinates of the element to describe the spinning of an assumed Euler-Bernoulli beam. The material coordinates of the two nodes are also included in the ALE description to deal with mass transportation at the boundaries. Thus, there are eight coordinates in each node: the position vector, slope vector, rotational angle, and material coordinate. This novel element could therefore be used to realize the dynamic modeling of a CCS beam whose length varies because of either simple or complex transportation boundaries and for which the effect of spinning is not negligible.

Several typical numerical examples were considered based on the ANCF. These were (i) free transverse vibration analysis of an axially moving string and free flexural vibration analysis of an axially moving beam, (ii) free torsional vibration analysis of a cylindrical shaft with one end fixed and the other end free, and (iii) dynamic response analysis of a spinning shaft subjected to an axially moving and rotating load. It is shown that the VLBSE could describe the dynamic behavior of a spinning CCS beam with mass transportation effectively.

Our VLBSE enriches the element library of the ANCF and allows this modeling method to be used in more types of dynamical problems in the field of flexible multi-body systems.

Acknowledgements

The authors are particularly grateful for the technical assistance from Suzhou Tongyuan Soft. & Ctrl. Tech. Co., Ltd., Associate Professor Wei Cheng and Doctor Guo Jiawen of Harbin Institute of Technology, Professor Ren Gexue of Tsinghua University, Associate Professor Tian Qiang of Beijing Institute of Technology.

Funding

This paper was supported by the National Defense Special Field Science and Technology Project (No. TY3Q20110005), the National Key Basic Research Program of China (Nos. 2013CB733006 and 2013CB733004) and the Self-Planned Task of State Key Laboratory of Robotics and System (Harbin Institute of Technology) (No. SKLRS201616B).

References

• Asemi, S.R., Farajpour, A., Borghei, M. and Hassani, A.H., (2014). Thermal effects on the stability of circular graphene sheets via nonlocal continuum mechanics. Latin American Journal of Solids and Structures 11: 704-724.
• Bonet, J. and Wood, R.D., (1997). Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press (Cambridge, UK).
• Brenan, K.E., Campbell, S.L. and Petzold, L.R., (1996). Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, SIAM(Philadelphia).
• Deng, Z.Q., Yang, S., Sun, J., Jiang, S.Y., Quan, Q.Q., Tang, J.Y., (2014). Dynamic modeling and analysis of rotating mechanism in planetary soil drilling sampler. In: 2014 IEEE International Conference on Mechatronics and Automation 663-668.
• Dmitrochenko, O.N. and Pogorelov, D.Y., (2003). Generalization of plate finite elements for absolute nodal coordinate formulation. Multibody System Dynamics 10: 17-43.
• Dombrowski, S.V., (2002). Analysis of large flexible body deformation in multibody systems using absolute coordinates. Multibody System Dynamics 8: 409-432.
• Du, J.L., Cui, C.Z., Bao, H. and Qiu, Y.Y., (2015). Dynamic analysis of cable-driven parallel manipulators using a variable length finite element. Journal of Computational and Nonlinear Dynamics 10: 011013:1-011013:7.
• Dufva, K. and Shabana, A.A., (2005). Analysis of thin plate structures using the absolute nodal coordinate formulation. Proceedings of the Institution of Mechanical Engineers Part K-Journal of Multi-Body Dynamics 219: 345-355.
• García-Vallejo, D., Escalona, J.L., Mayo, J. and Domínguez, J., (2003). Describing Rigid-Flexible Multibody Systems Using Absolute Coordinates. Nonlinear Dynamics 34: 75-94.
• Gatti, P.L. and Ferrari, V., (2003). Applied Structural and Mechanical Vibrations-Theory, methods and measuring instrumentation, Taylor & Francis Group (Abingdon, UK).
• Gerstmayr, J. and Shabana, A.A. (2006). Analysis of thin beams and cables using the absolute nodal co-ordinate formulation. Nonlinear Dynamics 45: 109-130.
• Han, X.G., Wang, M.J. and Ouyang, H.J., (2012). Establishment of vibration model and experimental analyses of workpiece with variable mass and diameter in turning process. Journal of Dalian University of Technology 52: 514-521.
• Hong, D.F. and Ren, G.X., (2011). A modeling of sliding joint on one-dimensional flexible medium. Multibody System Dynamics 26: 91-106.
• Hong, D.F., Tang, J.L. and Ren, G.X., (2011). Dynamic modeling of mass-flowing linear medium with large amplitude displacement and rotation. Journal of Fluids and Structures 27: 1137-1148.
• Huang, L.L. and Wang, K.S., (2014). Research on reliability of dynamic performance and precision of screw drive system. Journal of Beijing Information Science and Technology University 29: 75-79.
• Katz, R., (2001). The dynamic response of a rotating shaft subject to an axially moving and rotating load. Journal of Sound and Vibration 246: 757-775.
• Katz, R., Lee, C.W., Ulsoy, A.G. and Scott, R.A., (1988). The dynamic response of a rotating shaft subject to a moving load. Journal of Sound and Vibration 122: 131-148.
• Khan, A.A., Alam, M.N., Rahman, N.U. and Wajid, M., (2016). Finite Element Modelling for Static and Free Vibration Response of Functionally Graded Beam. Latin American Journal of Solids and Structures 13: 690-714.
• Kubler, L., Eberhard, P. and Geisler, J. (2003) Flexible multibody systems with large deformations and nonlinear structural damping using absolute nodal coordinates. Nonlinear Dynamics 34: 31-52.
• Liu, Y.Z., (2012). Exact dynamical model of axially moving beam with large deformation. Chinese Journal of Theoretical and Applied Mechanics 44: 832-838.
• Liu, Y.Z., Chen, L.Q. and Chen, W.L., (2011). Mechanics of Vibration (2nd ed.), Higher Education Press (Beijing).
• Lu, Y.J. (2007). Positioning and Orientating Control Study on the Feed Support System of a Large Radio Telescope, Ph.D. Thesis (in Chinese), Tsinghua University, China.
• Lv, B.L., Li, W.Y. and Ouyang, H.J., (2015). Moving force-induced vibration of a rotating beam with elastic boundary conditions. International Journal of Structural Stability and Dynamics 15: 1450035:1-1450035:24.
• Mikkola, A.M. and Shabana, A.A., (2003). A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications. Multibody System Dynamics 9: 283-309.
• Omar M.A. and Shabana, A.A., (2001). A two-dimensional shear deformable beam for large rotation and deformation problems. Journal of Sound and Vibration 243: 565-576.
• Shabana, A.A. and Yakoub, R.Y., (2001). Three dimensional absolute nodal coordinate formulation for beam elements: theory. Journal of Mechanical Design 123: 606-613.
• Shabana, A.A., (1996). An absolute nodal coordinates formulation for the large rotation and deformation analysis of flexible bodies, Technical Report (no. MBS96-1-UIC), University of Illinois at Chicago, USA.
• Shabana, A.A., (1997). Definition of the Slopes and the Finite Element Absolute Nodal Coordinate Formulation. Multibody System Dynamics 1: 339-348.
• Shabana, A.A., (2012). Computational Continuum Mechanics (2nd ed.), Cambridge University Press (New York).
• Sugiyama, H. and Suda, Y., (2007). A curved beam element in the analysis of flexible multi-body systems using the absolute nodal coordinates. Proceedings of the Institution of Mechanical Engineers Part K-Journal of Multi-Body Dynamics 221: 219-231.
• Tang, J.L., Ren, G.X., Zhu, W.D. and Ren, H., (2011). Dynamics of variable-length tethers with application to tethered satellite deployment. Communications in Nonlinear Science and Numerical Simulation 16: 3411-3424.
• Tian, Y., Deng, Z.Q., Tang, D.W., Jiang, S.Y. and Hou, X.Y., (2012). Structure parameters optimization and simulation experiment of auger in lunar soil drill-sampling device. Journal of Mechanical Engineering 48: 10-15.
• Wickert, J.A. and Mote, C.D., (1990). Classical vibration analysis of axially moving continua. Journal of Applied Mechanics-Transactions of the ASME 57: 738-744.
• Yakoub, R.Y. and Shabana, A.A., (2001). Three dimensional absolute nodal coordinate formulation for beam elements: implementation and applications. Journal of Mechanical Design 123: 614-621.
• Yang, S., Deng, Z.Q., Sun, J., Jiang, S.Y., Quan, Q.Q., Li, P., Hou, X.Y., (2014). Dynamic analysis of feeding drive mechanism in planetary soil automatic drilling sampler. In: 2014IEEE International Conference on Information and Automation 763-768.

Publication Dates

• Publication in this collection
  16 Aug 2017

History

• Received
  04 Apr 2017
• Reviewed
  12 May 2017
• Accepted
  05 June 2017