Differential Quadrature Finite Element and the Differential Quadrature Hierarchical Finite Element Methods for the Dynamics Analysis of on Board Shaft

Saimi Ahmed*, Hadjoui Abdelhamid, Bensaid Ismail and Fellah Ahmed

Mechanical Systems & Materials Engineering Laboratory, University Abou Bekr Belkaid Tlemcen, Algeria

E-mail: ahmedsaimi@hotmail.fr

*Corresponding Author

Received 04 December 2020; Accepted 05 February 2021; Publication 06 May 2021

Abstract

In this paper the dynamic analysis of a shaft rotor whose support is mobile is studied. For the calculation of kinetic energy and stiffness energy, the beam theory of Euler Bernoulli was used, and the matrices of elements and systems are developed using two methods derived from the differential quadrature method (DQM). The first method is the Differential Quadrature Finite Element Method (DQFEM) systematically, as a combination of the Differential Quadrature Method (DQM) and the Standard Finite Element Method (FEM), which has a reduced computational cost for problems in dynamics. The second method is the Differential Quadrature Hierarchical Finite Element Method (DQHFEM) which is used by expressing the matrices of the hierarchical finite element method in a similar form to that of the Differential Quadrature Finite Element Method and introducing an interpolation basis on the element boundary of the hierarchical finite element method. The discretization element used for both methods is a three-dimensional beam element. In the differential quadrature finite element method (DQFEM), the mass, gyroscopic and stiffness matrices are simply calculated using the weighting coefficient matrices given by the differential quadrature (DQ) and Gauss-Lobatto quadrature rules. The sampling points are determined by the Gauss-Lobatto node method. In the Differential Quadrature Hierarchical Finite Element Method (DQHFEM) the same approaches were used, and the cubic Hermite shape functions and the special Legendre polynomial Rodrigues shape polynomial were added. The assembly of the matrices for both methods (DQFEM and DQHFEM) is similar to that of the classical finite element method. The results of the calculation are validated with the h- and hp finite element methods and also with the literature.

Keywords: Rotor dynamics, differential quadrature finite elements method, differential quadrature hierarchical finite elements method, hp-version of FEM, on-board shaft.

1 Introduction

The first to give the term finite element method (FEM) is (Clough, 1960). The finite element method has undergone several optimisations and changes over time, it has gradually become a very effective tool for numerical solutions to a wide range of engineering problems. The h-FEM version was the first method developed (O. C. Zienkiewicz, 1977), which uses elements of fixed degree and low order and whose convergence is achieved by increasing the successive mesh refinement. A new version which is the p version of the finite element method was evolved during the 1970s (Zienkiewicz et al., 1970; Szabo, 1979), which uses a single mesh element but its convergence is obtained by increasing the degree of polynomial of the element, which is a set of shape functions, the p version is called hierarchical finite element method. The p version of the finite element method has great advantages over the h version of the of the finite element method, for example the mesh of the structure only has to be generated once. Convergence is achieved by increasing the p-order (Zienkiewicz et al., 1983; Beslin and Nicolas, 1997) without mesh refinement (Babuška and Suri, 1990). Input data can be reduced to a minimum, which greatly simplifies pre-post processing (Petyt, 2010). But also it has some disadvantages, such as contact analysis problems, especially for the node at the centre of higher order quadrilateral finite elements. The rate of convergence with the p-version is based on the choice of hierarchical or non-hierarchical form functions. It only depends on the highest full polynomial you can capture with your shape features. Hierarchical form functions, however, offer the advantage of an inhernt p-refinement capability. That is, when increasing the polynomial degree of shape functions from p to p $+$ 1, only the contributions corresponding to p $+$ 1 need to be calculated. Considering non-hierarchical shape functions (based on Lagrange polynomials, for example), all shape functions must be recomputed if the polynomial degree is high. On the other hand, the degrees of freedom do not retain any physical significance for functions of hierarchical form. therefore, mass clustering techniques are not feasible, which is a distinct drawback when considering high frequency dynamics as the problems of wave propagation. These two versions of the FEM were combined to give a new version noted by the h-p version of the FEM, which was first used by Babuska (Babuska et al., 1981; Oden, 1991; Babuska and B.Q, 1992) who discovered that the finite element method converges exponentially faster when the mesh is refined using an appropriate combination of $-$h refinements (dividing the elements into smaller ones) and $-$p refinements (increasing their polynomial degree). Many scientific researchers have focused on this axis and have given rise to several quoted finite element codes (Szabo and Prob, 1985; Heuveline and Rannacher, 2003) have also worked on the h-p version of the finite element method, (Suri, 1997) used the h-p version of the finite element method for hulls, also a study on a Bernoulli beam using the h-p version of the finite element method can be found in (Bardell, 1996), in the field of rotor dynamics, work can be found using the hp version of the MEF, the dynamics of the on-board rotor by (Saimi, 2016), the study of the effect of an open transverse crack on the vibratory behaviour of rotors using the hp version of the finite element method (Fellah et al., 2019), and the numerical analysis of the dynamic behaviour of a functionally graduated shaft in a thermal environment using the hp finite element method (Hassan et al., 2020).

Over the last three decades, the Differential Quadrature Method (DQM) has emerged and gradually evolved as an efficient and accurate numerical method, and has enjoyed notable success over the last two decades (Bellman and Casti, 1971; Bert, 1988; Bert and Malik, 1996; Shu, 2000). DQM is based on the approximation of the partial derivatives of a field variable at a discrete point by a weighted linear sum of the field variable along the line passing through that point. The significant development of DQM has led to an interest in combining DQM with a variational formulation. Striz et al. (1995) took the initiative and developed the hybrid quadrature element (QEM) method for two-dimensional plane stress and plate bending problems, and plate-free vibration problems (Striz et al., 1997). Chen and New (1999) used the DQ technique to discretize the derivatives of the variable functions existing in integral statements for variational methods, GALERKIN method, etc., by deriving the finite element formulation, discretization of linear 3-D statics and the elasticity problem and the buckling problem of a plate using the principle of minimum potential energy were illustrated. This method is called the Differential Quadrature Finite Element Method (DQFEM). (Xing and Liu, 2009) presented a Differential Quadrature Finite Element Method (DQFEM) which was motivated by the complexity of imposing boundary conditions in DQM, the name is the same as that of (Chen and New, 1999), but the starting points and implementations are different. Compared to (Zhong and Yu, 2009) and (Chen and New, 1999), DQFEM (Xing and Liu, 2009) presents the following novelties: The DQ rules are reformulated and, in conjunction with the Gauss-Lobatto integral rule, are used to discretise the energy functional to derive the finite element formulation of a Euler Bernoulli beam (Yufeng et al., 2010). Lagrange interpolation functions are used as test functions, nodal shape functions as in standard FEM are not required. The DQFE element matrices are symmetrical, well-conditioned and efficiently computed by simple algebraic operations of the known weighting coefficient matrices of the reformulated DQ rules and the Gauss-Lobatto integral rule. Drawing on the polynomial hierarchical version of the Finite Element Method (p-FEM) (Cuiyun et al., 2016) combined the DQFEM method by adding shape functions as in the p version of the finite elements, DQM has had difficulty solving free vibration (Malik and Bert, 2000; Bert and M, 1996) due to its difficulty in approximating high-order derivatives, but DQFEM can solve the problem (Xing and Liu, 2009). It is clear that p-FEM can take advantage of some techniques in DQFEM, which leads (Cuiyun et al., 2016) to propose the Differential Quadrature Hierarchical Finite Element Method (DQHFEM).

Knowledge of the vibratory behaviour of rotors is of great importance for the manufacturers and operators of this equipment, many researchers focus their research on this axis and on the study of dynamic behaviour, especially the critical speeds of the rotors which differ from its natural non-rotating frequency. The main reason for this difference is known to be the gyroscopic effect. Green (1948) first studied the gyroscopic effect on the normal frequency of flexible rotors. Eshienman and Eubanks (1967) and Rao (1983) also studied the gyroscopic effect on the normal frequency of rotating shafts, as well as the effect of external loading which can change the lateral normal frequency of rotating shafts. The effect of axial force and externally applied torques on the lateral vibrations of rotating shafts has been studied by several researchers. Bokaian (1948) presented the variations of the lateral normal frequency of Euler Bernoulli beams under axial loading with various boundary conditions. Choi et al. (1995) derived the equation of motion of the rotating shaft under constant compressive axial loading by introducing gyroscopic moments. Nelson and McVaugh (1976) used the finite element method for the study of rotors on deformable bearings. The study of the dynamic behaviour of a rotor under the loading effect consists of studying the overall behaviour of a rotor whose support is subjected to any movement. This model is well adapted to understand the phenomenon of movement of rotors embedded in vehicles, aircrafts etc. This phenomenon has recently attracted the attention of researchers. Duchemin et al. (2006) made a detailed study based on the simple Rayleigh-Ritz model, on a rotor whose support is in motion. Various analytical studies have been carried out on simple movements: simple translation, sinusoidal translation, constant rotation, accelerated rotation. Dakel (2013) continued the work by adding hydrodynamic bearings by the classical finite element method. In the work of (Saimi, 2016) there is a study on the behaviour of on-board rotor by the h-p version of finite elements. Sajal et al. (2018) studied the dynamic analysis of micro-beams based on modified strain gradient theory using the differential quadrature method. Weiyan et al. (2019) to use the generalized differential quadrature method combined with Hamilton’s principle for free vibration analysis of a thin-walled rotating composite shaft. Ri et al. (2020) used DQFEM for the analysis of nonlinear forced vibrations of composite beams. Belhadj et al. (2020) also used DQM for the study of free vibrations of single-walled carbon nanotubes with rotary inertia.

In this present work, following the literature mentioned above, we study the dynamic behaviour of a rotor system with a mobile support, using two different methods, the DQFEM and the DQHFEM. The global equation of motion is obtained using Euler Bernoulli’s beam theory and Lagrange’s equation. The material of the shaft is isotropic. The natural frequencies of the system are determined using a program developed in MATLAB and the results obtained are verified with those reported in the literature. The DQFEM and DQHFEM methods have been shown to be faster than the h-FEM and hp-FEM methods in calculating convergence.

2 Model of the Rotating Shaft

In this section we will give the essential steps for modelling a rotor shaft whose support is assumed to be non-deformable and in known deterministic motion. For the analysis and prediction of rotor dynamic behaviour, differential equations are established using the approaches of the modelling used (Saimi, 2016).

The assumptions and boundary conditions taken into account in the framework of this study are:

The shaft is deformable and considered as a homogeneous beam and will be modelled using the fundamental kinematic hypothesis of Euler Bernoulli’s isotropic beam theory.

The rotor support is infinitely rigid in the first case and mobile in the second case.

The rotor rotates at a constant rotational speed $\mathrm{\Omega }$.

Consideration of support motion can significantly affect the shape of the equations of motion of a rotor in flexion compared to those obtained for a fixed support. In order to obtain the simplest possible modelling, the approach presented by (Duchemin et al., 2006) is used. He proposes the modelling of a rotor with a mobile support by considering the movement of the rotor with respect to the support and that of the support with respect to the ground. The transverse deflections of the mean line of the rotor shaft are studied in relation to a reference frame linked to the rigid support. The novelty of this work compared to the work of (Duchemin et al., 2006) is the modelling methods, (Duchemin et al., 2006) used the classical version of finite elements, and in this work we used the differential quadrature finite element method and the differential quadrature hierarchical finite element method, which showed a very good speed of calculation and fast convergence rate.

In order not to neglect the movement of the support, three main markers are defined:

• The Galilean fixed reference frame $R_g(x_g,y_g,z_g)$

• The frame linked to the non-deformable support $R_s(x_s,y_s,z_s)$,

• Current, rotating, shaft-linked frame $R(x,y,z)$ .

The centres of these frames are respectively O, A et C.

In order to calculate the expressions of the kinetic energy and the strain energy of the shaft, it is necessary to calculate the speed and rotation vectors of the frame R in relation to the frame $R_g$.

Since there are three frames taken into consideration, two changes in the frame system can be carried out by:

• The transformation of the frame linked to the support $R_s$ into a local frame $R$.

• The transformation of the Galilean frame $R_g$ to the frame linked to the support $R_s$.

Generally speaking, in rotor dynamics, the rotation of the frame R linked to a point in the deformable shaft relative to the frame $R_s$ linked to the support is defined by Euler’s angles ${\theta }_x,{\theta }_y$ and ${\theta }_z$, using two intermediate frames.

• A rotation angle ${\theta }_z$ around $z_s$ (intermediate frame $R1(x1,y1,z1))$

• A rotation angle ${\theta }_x$ around the new axis $x1$ (intermediate frame $R2(x2,y2,z2)$)

• A rotation angle ${\theta }_y$ around the final axis $y//y2$ (final frame $R(x,y,z)$).

Where the rotation vector of $R$ relative to $R_s$ expressed in $R$ is (1):

${\overrightarrow{\mathrm{\Omega }}}_R^{Rs}$	$={\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill {\dot{\theta }}_z\hfill \end{array}\right]}_{R_s}+{\left[\begin{array}{c}\hfill {\dot{\theta }}_x\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]}_{R_2}+{\left[\begin{array}{c}\hfill 0\hfill \\ \hfill {\dot{\theta }}_y\hfill \\ \hfill 0\hfill \end{array}\right]}_R$	(1)
${\overrightarrow{\mathrm{\Omega }}}_R^{Rs}$	$={\left[\begin{array}{c}\hfill {\dot{\theta }}_{\mathrm{x}}\cos {\theta }_{\mathrm{y}}-{\dot{\theta }}_{\mathrm{z}}\cos {\theta }_{\mathrm{x}}\sin {\theta }_{\mathrm{y}}\hfill \\ \hfill {\dot{\theta }}_{\mathrm{y}}+{\dot{\theta }}_{\mathrm{z}}\sin {\theta }_{\mathrm{x}}\hfill \\ \hfill {\dot{\theta }}_{\mathrm{x}}\sin {\theta }_{\mathrm{y}}+{\dot{\theta }}_{\mathrm{z}}\cos {\theta }_{\mathrm{x}}\cos {\theta }_{\mathrm{y}}\hfill \end{array}\right]}_R$	(2)

The movement of the support is defined by the coordinates $x_A,y_A$, and $z_A$, of the vector $\overrightarrow{OA}$ expressed in the frame $R_g$, and by the angles $\alpha ,\beta$, and $\gamma$ to switch from the frame $R_g$ to the frame $R_s$ by:

• A rotation angle $\alpha$ around $z_g$ (intermediate frame $R^{\prime }(x^{\prime },y^{\prime },z^{\prime })$)

• A rotation angle $\beta$ around the new axis $x^{\prime }$ (intermediate frame $R^{\prime \prime }(x^{\prime \prime },y^{\prime \prime },z^{\prime \prime })$)

• A rotation angle $\gamma$ around the final axis $y^{\prime \prime }=y_s$ (final frame $R_s(x_s,y_s,z_s)$)

Knowing that the rotations $\alpha ,\beta$, and $\gamma$ depends on the time t.

The rotation vector of $R_s$ relative to $R_g$ expressed in $R_s$ is written:

${\overrightarrow{\mathrm{\Omega }}}_{R_s}^{R_g}$	$={\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill \dot{\alpha }\hfill \end{array}\right]}_{R^{\prime }}+{\left[\begin{array}{c}\hfill \dot{\beta }\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]}_{R^{\prime \prime }}+{\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \dot{\gamma }\hfill \\ \hfill 0\hfill \end{array}\right]}_{R_s}$	(3)
${\overrightarrow{\mathrm{\Omega }}}_{R_s}^{R_g}$	$={\left[\begin{array}{c}\hfill \dot{\beta }\cos \gamma -\dot{\alpha }\cos \beta \sin \gamma \hfill \\ \hfill \dot{\gamma }+\dot{\alpha }\sin \beta \hfill \\ \hfill \dot{\beta }\sin \gamma +\dot{\alpha }\cos \beta \cos \gamma \hfill \end{array}\right]}_{R_s}={\left[\begin{array}{c}\hfill {\dot{\alpha }}_s\hfill \\ \hfill {\dot{\beta }}_s\hfill \\ \hfill {\dot{\gamma }}_s\hfill \end{array}\right]}_{R_s}$	(4)

To set up the translational movements of the support, the $X,Y,Z$ coordinates of the vector $\overrightarrow{OA}$ expressed in the frame $R_s$ are used.

$\overrightarrow{OA}={\left[\begin{array}{c}(x_A\cos \alpha +y_A\sin \alpha )\cos \gamma \hfill \\ -(z_A\cos \beta +\left(x_A\sin \alpha -y_A\cos \alpha \right)\sin \beta )\sin \gamma \hfill \\ z_A\sin \beta -\left(x_A\sin \alpha -y_A\cos \alpha \right)\cos \beta \hfill \\ (x_A\cos \alpha +y_A\sin \alpha )\sin \gamma \hfill \\ +(z_A\cos \beta +(x_A\sin \alpha -y_A\cos \alpha )\sin \beta )\cos \gamma \hfill \end{array}\right]}_{R_s}=\left[\begin{array}{c}\hfill X\hfill \\ \hfill Y\hfill \\ \hfill Z\hfill \end{array}\right]$

(5)

Starting from the equations of motion of the support expressed in the Galilean coordinate system $R_g$, The coordinates of the rotation vectors and the position of the frame $R$ relative to $R_g$ can be easily obtained. For the following developments, the equations are expressed as a function of ${\dot{\alpha }}_s$, ${\dot{\beta }}_s$, ${\dot{\gamma }}_s$ and $X,Y,Z$ and their derivatives in relation to time.

The rotation vector of R with respect to the frame $R_g$ is:

${\overrightarrow{\mathrm{\Omega }}}_R^{R_g}$	$={\overrightarrow{\mathrm{\Omega }}}_{R_s}^{R_g}+{\overrightarrow{\mathrm{\Omega }}}_R^{Rs}={\left[\begin{array}{c}\hfill {\dot{\alpha }}_s\hfill \\ \hfill {\dot{\beta }}_s\hfill \\ \hfill {\dot{\gamma }}_s\hfill \end{array}\right]}_{R_s}+{\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill {\dot{\theta }}_z\hfill \end{array}\right]}_{R_s}+{\left[\begin{array}{c}\hfill {\dot{\theta }}_x\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]}_{R_2}+{\left[\begin{array}{c}\hfill 0\hfill \\ \hfill {\dot{\theta }}_y\hfill \\ \hfill 0\hfill \end{array}\right]}_R$	(6)
${\overrightarrow{\mathrm{\Omega }}}_R^{R_g}$	$={\left[\begin{array}{c}({\dot{\alpha }}_s\cos {\theta }_z+{\dot{\beta }}_s\sin {\theta }_z+{\dot{\theta }}_x)\cos {\theta }_y-(({\dot{\alpha }}_s\sin {\theta }_z-{\dot{\beta }}_s\cos {\theta }_z)\sin {\theta }_x\hfill \\ +({\dot{\gamma }}_s+{\dot{\theta }}_y)\cos {\theta }_x)\sin \theta {}_y\hfill \\ -({\dot{\alpha }}_s\sin {\theta }_z-{\dot{\beta }}_s\cos {\theta }_z)\cos {\theta }_x+({\dot{\gamma }}_s+{\dot{\theta }}_y)\sin {\theta }_x+{\dot{\theta }}_y\hfill \\ ({\dot{\alpha }}_s\cos {\theta }_z+{\dot{\beta }}_s\sin {\theta }_z+{\dot{\theta }}_x)\sin {\theta }_y\hfill \\ +(({\dot{\alpha }}_s\sin {\theta }_z-{\dot{\beta }}_s\cos {\theta }_z)\sin {\theta }_x+({\dot{\gamma }}_s+{\dot{\theta }}_y)\cos {\theta }_x)\cos {\theta }_y\hfill \end{array}\right]}_R$	(7)
${\overrightarrow{\mathrm{\Omega }}}_R^{R_g}$	$=\left[\begin{array}{c}\hfill {\omega }_x\hfill \\ \hfill {\omega }_y\hfill \\ \hfill {\omega }_z\hfill \end{array}\right]$	(8)

The rotational speed of the rotor is along the y-axis; ${\theta }_z$ and ${\theta }_x$ represent the angular deformations of the shaft in the x and z directions; ${\theta }_y$ represents its angular position in relation to the support. The shaft only undergoes small deformations (elastic range), ${\theta }_z$ and ${\theta }_x$ are therefore considered to be infinitely small angles.

It is assumed that the rotor rotates at constant speed. ${\dot{\theta }}_y$:

${\dot{\theta }}_y=\mathrm{\Omega }\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{\theta }_y=\mathrm{\Omega }\mathrm{t}$

(9)

Assume $(u,y,w)$ are the displacements of a point C of the shaft in the reference frame $R_s$, $u$ and $w$ are variable whereas y is considered constant since only the bending movements of the shaft are studied, therefore:

$$\overrightarrow{AC}={\left[\begin{array}{c}\hfill u(y,t)\hfill \\ \hfill y\hfill \\ \hfill w(y,t)\hfill \end{array}\right]}_{R_s}$$

(10)

And the velocity vector of point C with respect to the Galilean frame:

${\overrightarrow{V}}^g(C)$	$={\displaystyle \frac{d^g}{dt}}\overrightarrow{OC}={\displaystyle \frac{d^s}{dt}}\overrightarrow{OC}+{\overrightarrow{\mathrm{\Omega }}}_{R_s}^{R_g}\wedge \overrightarrow{OC}\mathit{\hspace{1em}}\mathrm{with}$
$\overrightarrow{OC}$	$=\overrightarrow{OA}+\overrightarrow{AC}=\left[\begin{array}{c}\hfill X+u(y,t)\hfill \\ \hfill Y+y\hfill \\ \hfill Z+w(y,t)\hfill \end{array}\right]$	(11)
${\overrightarrow{V}}^g(C)$	$={\left[\begin{array}{c}\hfill \dot{X}+\dot{u}\hfill \\ \hfill \dot{Y}\hfill \\ \hfill \dot{Z}+\dot{w}\hfill \end{array}\right]}_{R_s}+{\left[\begin{array}{c}\hfill {\dot{\alpha }}_s\hfill \\ \hfill {\dot{\beta }}_s\hfill \\ \hfill {\dot{\gamma }}_s\hfill \end{array}\right]}_{R_s}\wedge {\left[\begin{array}{c}\hfill X+u\hfill \\ \hfill Y+y\hfill \\ \hfill Z+w\hfill \end{array}\right]}_{R_s}$	(12)
${\overrightarrow{V}}^g\left(C\right)$	$=\left[\begin{array}{c}\hfill \dot{X}+\dot{u}+{\dot{\beta }}_s(Z+w)-{\dot{\gamma }}_s(Y+y)\hfill \\ \hfill \dot{Y}+{\dot{\gamma }}_s(X+u)-{\dot{\alpha }}_s(Z+w)\hfill \\ \hfill \dot{Z}+\dot{w}+{\dot{\alpha }}_s(Y+y)-{\dot{\beta }}_s(X+u)\hfill \end{array}\right]=\left[\begin{array}{c}\hfill u_c\hfill \\ \hfill v_c\hfill \\ \hfill w_c\hfill \end{array}\right]$	(13)

The expressions of the kinetic energies of the shaft can be calculated from the expressions of the rotation vector ${\overrightarrow{\mathrm{\Omega }}}_R^{R_g}$ and the speed of the point C with respect to the frame $R_g$.

2.1 Kinetic Energy of the Shaft

The shaft is modelled by Euler Bernoulli’s beam theory, and rotates at a constant speed around its longitudinal axis. The shaft has a uniform circular cross-section.

The kinetic energy of the rotating symmetrical shaft, including the effects rotary inertia see (Saimi, 2016), can be written as

$d{E}_a$	$=({\displaystyle \frac{{\rho }_a{S}_a}{2}}(u_c{}^2+v_c{}^2+w_c{}^2)$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{2}}({\rho }_a{I}_{ax}\omega _x{}^2+{\rho }_a(I_{ax}+I_{az})\omega _y{}^2+{\rho }_a{I}_{az}\omega _z{}^2))dy$	(14)

With $S_a$, ${\rho }_a$, $I_{ax}$ and $I_{az}$ denote respectively the section (supposed to be constant), density and Section inertia.

$I_{sa}^{mo}$	$={\displaystyle \frac{I_{ax}+I_{az}}{2}}$	(15)
$E_a$	$=E_{a,1}+E_{a,2}+E_{a,3}$	(16)
$E_{a,1}$	$={\displaystyle \frac{{\rho }_a{S}_a}{2}}{\displaystyle {\int }_0^{l_a}}({\dot{u}}^2+{\dot{w}}^2+2(({\dot{x}}_0+z_0{\omega }^y-(y_0+y){\omega }^z))(\dot{u}+w{\omega }^y)$
	$\mathrm{\hspace{1em}}+2(({\dot{z}}_0+(y_0+y){\omega }^x-x_0{\omega }^y))(\dot{w}-u{\omega }^y)$
	$\mathrm{\hspace{1em}}-2(({\dot{y}}_0-z_0{\omega }^x+x_0{\omega }^z))(w{\omega }^x-u{\omega }^z)$
	$\mathrm{\hspace{1em}}+2(\dot{u}w-\dot{w}u){\omega }^y+w^2\omega ^x{}^2+(u^2+w^2)\omega ^y{}^2$
	$\mathrm{\hspace{1em}}+u^2\omega ^z{}^2-2uw{\omega }^x{\omega }^z)dy$	(17)
$E_{a,2}$	$={\displaystyle \frac{{\rho }_a{I}_{sa}^{mo}}{2}}{\displaystyle {\int }_0^{l_a}}({\left({\displaystyle \frac{\partial \dot{u}}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial \dot{w}}{\partial y}}\right)}^2+2{\displaystyle \frac{\partial \dot{w}}{\partial y}}{\omega }^x$
	$\mathrm{\hspace{1em}}+2\left({\displaystyle \frac{\partial \dot{u}}{\partial y}}{\displaystyle \frac{\partial w}{\partial y}}-{\displaystyle \frac{\partial \dot{w}}{\partial y}}{\displaystyle \frac{\partial u}{\partial y}}\right){\omega }^y-2{\displaystyle \frac{\partial \dot{u}}{\partial y}}{\omega }^z-\left({\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2-1\right)\omega ^x{}^2$
	$\mathrm{\hspace{1em}}+\left({\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\right)\omega ^y{}^2-\left({\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2-1\right)\omega ^z{}^2$
	$\mathrm{\hspace{1em}}-2{\displaystyle \frac{\partial u}{\partial y}}{\displaystyle \frac{\partial w}{\partial y}}{\omega }^x{\omega }^z-2({\displaystyle \frac{\partial u}{\partial y}}{\omega }^x+{\displaystyle \frac{\partial w}{\partial y}}{\omega }^z){\omega }^y)dy$	(18)
$E_{a,3}$	$={\rho }_a{I}_{sa}^{mo}{\displaystyle {\int }_0^{l_a}}({(\mathrm{\Omega }+{\omega }^y)}^2+{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2\omega ^x{}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\omega ^z{}^2$
	$\mathrm{\hspace{1em}}+(\mathrm{\Omega }+{\omega }^y)(-2{\displaystyle \frac{\partial \dot{u}}{\partial y}}{\displaystyle \frac{\partial w}{\partial y}}-({\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2){\omega }^y$
	$\mathrm{\hspace{1em}}+2({\displaystyle \frac{\partial u}{\partial y}}{\omega }^x+{\displaystyle \frac{\partial w}{\partial y}}{\omega }^z))+2{\displaystyle \frac{\partial u}{\partial y}}{\displaystyle \frac{\partial w}{\partial y}}{\omega }^x{\omega }^z)dy$	(19)

2.2 The Strain Energy

The strain energy is not affected by the movement of the support because it only depends on the constraints and therefore the deformation of the shaft with respect to the support. In this calculation, only the deformations due to bending is taken into account (the effects of shear are neglected).

In the case of an Euler-Bernoulli beam, shear effects are neglected and relationships between rotations ${\theta }_z$ and ${\theta }_x$ supposedly small and $u$ and $w$ displacements are expressed by (Saimi, 2016):

$U_a={\displaystyle \frac{1}{2}}{E}_a{I}_{S_a}^{mo}{\displaystyle {\int }_0^{l_a}}\left({\left({\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right)\mathit{d}y$

(20)

3 The Reformulated Differential Quadrature Rule

Known DQ rules approximate the derivative of a function at a point by a weighted linear sum of field variables along a line passing through the point. In addition to Lagrange functions, any other complete base can be used as a basis for formulating DQ rules. (Cuiyun et al., 2016) (Xing and Liu, 2009)

Thus, for a field variable $f(x)$ its derivative of order n in a discrete point $x_i$ can be expressed as:

${\displaystyle \frac{{\partial }^{n}f(x,t)}{\partial x^n}}{|}_{x_i}={\displaystyle \sum _{j=1}^N}{A}_{ij}^{(n)}f(x_j,t)\mathit{\hspace{1em}}(i=1,2,3,\mathrm{\dots },N)$

(21)

Where $A_{ij}^{(n)}$ is the weighting coefficient related to the derivative of order n, and the weighting coefficient is obtained as follows if $n=1$, so

$A_{ij}^{(1)}$	$={\displaystyle \frac{M(x_i)}{(x_i-x_j)M(x_j)}}\mathit{\hspace{1em}}i\ne j,i,j=1,2,\mathrm{\dots },N$
$A_{ii}^{(1)}$	$=-{\displaystyle \sum _{j=1,j\ne i}^n}{A}_{ij}^{(1)}\mathit{\hspace{1em}}i=1,2,\mathrm{\dots },N$	(22)

Where

$M(x_i)$	$={\displaystyle \prod _{k=1,k\ne i}^N}(x_i-x_k)$
$M(x_j)$	$={\displaystyle \prod _{k=1,k\ne i}^N}(x_j-x_k)$	(23)

If $n>1$, secondary and higher order derivatives, the weighting coefficients are determined using the following simple recurrence relationship:

$A_{ij}^{(n)}$	$=n\left(A_{ij}^{(1)}*A_{ii}^{(n-1)}-{\displaystyle \frac{A_{ij}^{(n-1)}}{(x_i-x_j)}}\right)\mathit{\hspace{1em}}i\ne j,i,j=1,2,\mathrm{\dots },N,n>1$
$A_{ii}^{(n)}$	$=-{\displaystyle \sum _{j=1,j\ne i}^N}{A}_{ij}^{(n)}\mathit{\hspace{1em}}i=1,2,\mathrm{\dots },N$	(24)

4 Gauss-Lobatto Quadrature Rule

The theory of Gauss-Lobatto quadrature rules can be found in the mathematical literature; The Gauss-Lobatto quadrature rule with a degree of accuracy $(2n-3)$ for the function $f(x)$ defined in $[-1,1]$ is:

$${\int }_{-1}^{1}f(x)\mathit{d}x=\sum _{j=1}^N{C}_{j}f(x_j)$$

(25)

With the weighting coefficient $C_j$ of the Gauss-Lobatto integration is given by:

$C_1=C_N={\displaystyle \frac{2}{N(N-1)}},C_j={\displaystyle \frac{2}{N(N-1){[P_{N-1}(x_j)]}^2}}\mathit{\hspace{1em}}(j\ne 1,N)$

(26)

$x_j$ is the $(j-1)$ zero of the first order derivative of $P_{N-1}(x)$. To solve the roots of the Legendre polynomials, we will use the recursivity formula as Equations (27) and (28), it is easy to obtain thousands of roots.

$P_{N+1}(x)={\displaystyle \frac{2N+1}{N+1}}x{P}_N(x)-{\displaystyle \frac{N}{N+1}}{P}_{N-1}(x)$

(27)

With $P_0(x)=1$, $P_1(x)=x$. The nth-order derivation of the Legendre polynomials can be determined by the following formula:

$P_{N+1}^{(n)}(x)=x{P}_N^{(n)}(x)+(N+n)P_N^{(n)}(x)$

(28)

In order to obtain a denser population near the boundaries, sampling points are selected according to the grid distribution of Gauss–Lobatto nodes.

$$x_j=-\cos \left(\frac{j-1}{N-1}\pi \right)$$

(29)

Gauss–Lobatto nodes are solved with the Newton-Raphson iteration method.

$x^{iT+1}=x^{iT}-{\mathit{F}}^{\prime }{(x^{iT})}^{-1}\mathit{F}(x^{iT}),iT=0,1,\mathrm{\dots }$

(30)

Including

$$x={[x_2,x_3,\mathrm{\dots },x_{N-1}]}^T$$	(31)
$$\mathit{F}(x)={[f(x_2),f(x_3),\mathrm{\dots },f(x_{N-1})]}^T$$	(32)
$${\mathit{F}}^{\prime }(x)={\left[\frac{\partial f(x_j)}{\partial x_i}\right]}_{(N-2)\times (N-2)}$$	(33)
$$f(x_j)=\sum _{k-1,k\ne j}^N\frac{1}{x_j-x_k}\mathit{\hspace{1em}}j=2,3,\mathrm{\dots },N-1$$	(34)
$$\frac{\partial f(x_j)}{\partial x_i}=\{\begin{array}{cc}-\sum _{k=1,k\ne j}^N\frac{1}{{(x_j-x_k)}^2},\hfill & (i=j)\hfill \\ \frac{1}{{(x_j-x_i)}^2},\hfill & (i\ne j)\hfill \end{array}$$	(35)

Where $k$ is the value of $x$ at $i{T}^{th}$ iteration step. This method is less sensitive to the initial value. The values given by Equation (29) are used as initial values.

5 The Differential Quadrature Finite Element Method

The differential quadrature finite element method was developed by (Xing and Liu, 2009), whose differential quadrature rules and Gauss-Lobatto quadrature are used to discretize the system energies.

Assuming that the deflection function is in the form:

$u(x)$	$={\displaystyle \sum _{i=1}^N}{L}_i(x)u_i$
$w(x)$	$={\displaystyle \sum _{i=1}^N}{L}_i(x)w_i$	(36)

With $L_i$ is the Lagrange polynomial, and $u_i=u(x_i)$, $w_i=w(x_i)$ are the displacements of the Gauss Lobatto quadrature points or the DQ nodal displacements of the beam finite element.

Using DQ rules and Gauss - Lobatto quadrature the expressions of kinetic energy and strain energy (14–20) can be written as follows:

$E_{a,1}$	$={\displaystyle \frac{{\rho }_a{S}_a}{2}}{\eta }_1(C,A^{(1)})$	(37)
$E_{a,2}$	$={\displaystyle \frac{{\rho }_a{I}_{sa}^{mo}}{2}}{\eta }_2(C,A^{(1)})$	(38)
$E_{a,3}$	$={\rho }_a{I}_{sa}^{mo}{\eta }_3(C,A^{(1)})$	(39)

${\eta }_1(C,A^{(1)})$, ${\eta }_2(C,A^{(1)})$, and ${\eta }_3(C,A^{(1)})$ are detailed in the appendix.

$U_a={\displaystyle \frac{1}{2}}{E}_a{I}_{S_a}^{mo}({\overline{u}}^{T}A^{(2)}{}^{T}C{A}^{(2)}\overline{u}+{\overline{w}}^{T}A^{(2)}{}^{T}C{A}^{(2)}\overline{w})$

(40)

With $A^{(1)}$ and $A^{(2)}$ indicates the matrices of the weighting coefficients of the DQ rules for the first and second order derivatives respectively calculated with Equations (22–24), with respect to the Gauss-Lobatto nodes, and

$$C=diag[C_1{C}_2\mathrm{\dots }{C}_N]$$

(41)

Where $C_j$ are the weighting coefficients of the Gauss-Lobatto integration.

${\overline{u}}^T$	$=[u_1{u}_2\mathrm{\dots }{u}_N]$
${\overline{w}}^T$	$=[w_1{w}_2\mathrm{\dots }{w}_N]$	(42)

In order to construct an element that satisfies the requirements of continuity between elements, the element displacement vectors must be:

$u^T$	$=[u_{1}u^{\prime }{}_1{u}_3\mathrm{\dots }{u}_{N-2}{u}_{N}u^{\prime }{}_N]$
$w^T$	$=[w_{1}w^{\prime }{}_1{w}_3\mathrm{\dots }{w}_{N-2}{w}_{N}w^{\prime }{}_N]$	(43)

The relation between $w$ and $\overline{w}$ is defined using the DQ rule:

$u=Q\overline{u},w=Q\overline{w}$

(44)

Where

$Q=\left[\begin{array}{cccccc}\hfill 1\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill A_{1,1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{1,2}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{1,3}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill A_{1,N-1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{1,N}^{(1)}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 1\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\ddots }\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 1\hfill \\ \hfill A_{N,1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{N,2}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{N,3}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill A_{N,N-1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{N,N}^{(1)}\hfill \end{array}\right]$

(45)

By replacing Equations (44–45) in Equations (37–40), and using the formula of the Lagrange equations, the elementary mass matrix $[M^e]$, the gyroscopic matrix $[G^e]$, $[G^e{}_{{\omega }^y}]$ and the stiffness matrices $[K^e]$, $[K^e{}_{\mathrm{\Omega }{\omega }^y}]$, $[K^e{}_{\omega ^x{}^2}]$, $[K^e{}_{\omega ^y{}^2}]$, $[K^e{}_{\omega ^z{}^2}]$, $[K^e{}_{{\omega }^x{\omega }^z}]$ are obtained.

The elementary mass matrix obtained with (DQFEM)

$[M^e]$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill Q^{-T}\overline{\mathit{C}}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill Q^{-T}\overline{\mathit{C}}{Q}^{-1}\hfill \end{array}\right]$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\hfill \end{array}\right]$	(46)

The gyroscopic matrices obtained with (DQFEM)

$[G^e]$	$=2{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\hfill \\ \hfill Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\mathrm{\Omega }$	(47)
$G^e{}_{{\omega }^y}$	$=2{\rho }_a{S}_a\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill Q^{-T}\overline{\mathit{C}}{Q}^{-1}\hfill \\ \hfill -Q^{-T}\overline{\mathit{C}}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]{\omega }^y$	(48)

The stiffness matrices obtained with (DQFEM)

$[K^e]$	$=E_a{I}_{S_a}^{mo}\left[\begin{array}{cc}\hfill Q^{-T}\overline{\mathit{A}}^{(2)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(2)}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill Q^{-T}\overline{\mathit{A}}^{(2)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(2)}{Q}^{-1}\hfill \end{array}\right]$	(49)
$[K^e{}_{\mathrm{\Omega }{\omega }^y}]$	$=2{\rho }_a{I}_{sa}^{mo}$
	$\mathrm{\hspace{1em}}\left[\begin{array}{cc}\hfill -Q^{-T}\overline{\mathit{A}}^{(1)}{}^T{\overline{\mathit{C}\mathit{A}}}^{(1)}{Q}^{-1}\hfill & \hfill [0]\hfill \\ \hfill [0]\hfill & \hfill Q^{-T}\overline{\mathit{A}}^{(1)}{}^T{\overline{\mathit{C}\mathit{A}}}^{(1)}{Q}^{-1}\hfill \end{array}\right]\mathrm{\Omega }{\omega }^y$	(50)
$[K^e{}_{\omega ^x{}^2}]$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -Q^{-T}\overline{\mathit{C}}{Q}^{-1}\hfill \end{array}\right]\omega ^x{}^2$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill -Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\omega ^x{}^2$	(51)
$[K^e{}_{\omega ^y{}^2}]$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill -Q^{-T}\overline{\mathit{C}}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -Q^{-T}\overline{\mathit{C}}{Q}^{-1}\hfill \end{array}\right]\omega ^y{}^2+{\rho }_a{I}_{sa}^{mo}$
	$\mathrm{\hspace{1em}}\left[\begin{array}{cc}\hfill Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\hfill \end{array}\right]\omega ^y{}^2$	(52)
$[K^e{}_{\omega ^z{}^2}]$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill -Q^{-T}\overline{\mathit{C}}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\omega ^z{}^2$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\hfill \end{array}\right]\omega ^z{}^2$	(53)
$[K^e{}_{{\omega }^x{\omega }^z}]$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill Q^{-T}\overline{\mathit{C}}{Q}^{-1}\hfill \\ \hfill Q^{-T}\overline{\mathit{C}}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]{\omega }^x{\omega }^z+{\rho }_a{I}_{sa}^{mo}$
	$\mathrm{\hspace{1em}}\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\hfill \\ \hfill -Q^{-T}\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{Q}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]{\omega }^x{\omega }^z$	(54)

It is easy to show that the transformation matrix Q in Equation (45) is generally well conditioned. The mass, stiffness and gyroscopic matrices of the FEM and DQFEM elements are almost the same, but Lagrange polynomials are used in Equation (36) while Hermite interpolation functions are used in FEM.

All forms of node distribution for differentiation and quadrature are $[-1,1]$. Therefore, in order to apply them in practice, the following modifications must be made to the differential and quadrature matrices,

$\overline{\mathit{C}}={\displaystyle \frac{l_e}{2}}C,{\overline{\mathit{A}}}^{(1)}={\displaystyle \frac{2}{l_e}}{A}^{(1)},{\overline{\mathit{A}}}^{(2)}={\displaystyle \frac{4}{l_e{}^2}}{A}^{(2)}$

(55)

Where $l_e$ is the length of the beam element.

The matrices for the entire system are obtained according to the MEF rules for assembling elementary matrices,

$[M]\left\{\begin{array}{c}\hfill \ddot{u}(t)\hfill \\ \hfill \ddot{w}(t)\hfill \end{array}\right\}+[G]\left\{\begin{array}{c}\hfill \dot{u}(t)\hfill \\ \hfill \dot{w}(t)\hfill \end{array}\right\}+[K]\left\{\begin{array}{c}\hfill u(t)\hfill \\ \hfill w(t)\hfill \end{array}\right\}=[0]$

(56)

6 The Differential Quadrature Hierarchical Finite Element Method

In this section we aim to illustrate the use of DQHFEM through the on-board shaft, although work has already been done (Cuiyun et al., 2016) for a uniform Euler-Bernoulli beam, we will follow the same steps to determine the differential equation of motion of the on-board shaft.

The displacement field used for the DQHFEM is the same used for hp-FEM:

$u[x(\xi )]$	$=H_1(\xi )u_1+{\displaystyle \frac{L_e}{2}}{H}_2(\xi ){\displaystyle \frac{d{u}_1}{dx}}+H_3(\xi )u_2$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{L_e}{2}}{H}_4(\xi ){\displaystyle \frac{d{u}_2}{dx}}+{\displaystyle \sum _{n=1}^M}{\psi }_n(\xi )U_n$	(57)
$w[x(\xi )]$	$=H_1(\xi )w_1+{\displaystyle \frac{L_e}{2}}{H}_2(\xi ){\displaystyle \frac{d{w}_1}{dx}}+H_3(\xi )w_2$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{L_e}{2}}{H}_4(\xi ){\displaystyle \frac{d{w}_2}{dx}}+{\displaystyle \sum _{n=1}^M}{\psi }_n(\xi )W_n$	(58)

The first four functions $H_1(\xi )$, $H_2(\xi )$, $H_3(\xi )$ and $H_4(\xi )$ are those of the finite element method necessary to describe the displacements and rotations at the nodes of the element, we use for this purpose the cubic Hermit shape functions (Bardell, 1996).

$H_1(\xi )$	$={\displaystyle \frac{1}{4}}{(1-\xi )}^2(2+\xi )\mathit{\hspace{1em}}{H}_2(\xi )={\displaystyle \frac{1}{4}}{(1-\xi )}^2(\xi +1)$
$H_3(\xi )$	$={\displaystyle \frac{1}{4}}{(1+\xi )}^2(2-\xi )\mathit{\hspace{1em}}{H}_4(\xi )={\displaystyle \frac{1}{4}}{(1+\xi )}^2(\xi -1)$	(59)

Where the local coordinates are linked to the non-dimensional coordinates by the relation:

$x={\displaystyle \frac{L_e}{2}}(\xi +1)avec-1\le \xi \le 1$

(60)

And ${\psi }_n(\xi )$ are the hierarchical functions contribute to the field of internal displacement,

${\psi }_n(\xi )={\displaystyle \frac{{({\xi }^2-1)}^2}{n(n+1)(n+2)(n+3)}}{\displaystyle \frac{d^2{P}_{n+1}}{d{\xi }^2}}$

(61)

The low order Legendre polynomial can be calculated from the Rodrigues form of special Legendre polynomials (Peano, 1976); the generating function is quoted below.

$P_n(\xi )$	$={\displaystyle \sum _{k=0}^{\frac{(n-1)}{2}}}{\displaystyle \frac{{(-1)}^k(2n-2k-7)!!}{2^{k}k!(n-2k-1)!}}{\xi }^{(n-2k-1)}\mathit{\hspace{1em}}\text{with},n>4$
$n!!$	$=n(n-2)(n-4)\mathrm{\dots }(2or\mathrm{\hspace{0.25em}1}),0!!=(-1)!!=1,$
	$\text{and }(n-1)/2\text{Refers to its own whole part}$	(62)

This Legendre polynomial expression has been used in several papers (Saimi, 2016), (Hassan et al., 2020) in the h, p and hp versions of the finite element method.

One can also use the recursivity formula Equations (27–28), the order n can reach several thousands.

The displacement vectors of the element are noted as follows:

${\mathit{u}}^{\mathit{T}}$	$=[u_{1}u^{\prime }{}_1{u}_{2}u^{\prime }{}_2{U}_1\mathrm{\dots }{U}_M]$
${\mathit{w}}^{\mathit{T}}$	$=[w_{1}w^{\prime }{}_1{w}_{2}w^{\prime }{}_2{W}_1\mathrm{\dots }{W}_M]$	(63)

So the equations of $u[x(\xi )]$ and $w[x(\xi )]$ become:

$u[x(\xi )]$	$={[N_u]}^T\mathit{u}$
$w[x(\xi )]$	$={[N_w]}^T\mathit{w}$	(64)

Where

$\mathrm{N}_u{}^T$	$=\left[H_1(\xi )-{\displaystyle \frac{L_e}{2}}{H}_2(\xi )H_3(\xi )-{\displaystyle \frac{L_e}{2}}{H}_4(\xi ){\psi }_1(\xi ){\psi }_2(\xi )\mathrm{\dots }{\psi }_M(\xi )\right]$
$\mathrm{N}_w{}^T$	$=\left[H_1(\xi ){\displaystyle \frac{L_e}{2}}{H}_2(\xi )H_3(\xi ){\displaystyle \frac{L_e}{2}}{H}_4(\xi ){\psi }_1(\xi ){\psi }_2(\xi )\mathrm{\dots }{\psi }_M(\xi )\right]$	(65)

The calculation of Gauss-Lobatto nodes ${\xi }_j$, $j=1,2,\mathrm{\dots },N$, with $N=M+4$. Defines the following motion vectors:

${\overline{\mathit{u}}}^{\mathit{T}}$	$=[u(x_1)u(x_2)\mathrm{\dots }u(x_N)]$
${\overline{\mathit{w}}}^{\mathit{T}}$	$=[w(x_1)w(x_2)\mathrm{\dots }w(x_N)]$	(66)

According to Equation (64) the relationship between (63) and (66) is defined by the following equation:

$\overline{\mathit{u}}$	$={\mathit{G}}_{\mathit{u}}\mathit{u}$
$\overline{\mathit{w}}$	$={\mathit{G}}_{\mathit{w}}\mathit{w}$	(67)

Where

${\mathit{G}}_{\mathit{u}}$	$={[[{\mathrm{N}}_u]({\xi }_1)\mathit{\hspace{1em}}[{\mathrm{N}}_u]({\xi }_2)\mathit{\hspace{1em}}\mathrm{\dots }\mathit{\hspace{1em}}[{\mathrm{N}}_u]({\xi }_N)]}^{\mathit{T}}$
${\mathit{G}}_{\mathit{w}}$	$={[[{\mathrm{N}}_w]({\xi }_1)\mathit{\hspace{1em}}[{\mathrm{N}}_w]({\xi }_2)\mathit{\hspace{1em}}\mathrm{\dots }\mathit{\hspace{1em}}[{\mathrm{N}}_w]({\xi }_N)]}^{\mathit{T}}$	(68)

By replacing Equations (22, 24, 41, 66 and 68) in Equations (37–40), and using the formula of Lagrange’s equations, the elementary mass, stiffness and gyroscopic matrices are obtained as:

The elementary mass matrix obtained with (DQHFEM)

$M^e$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill \mathit{G}_{\mathit{u}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill \mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]$	(69)

The gyroscopic matrices obtained with (DQHFEM)

$G^e$	$=2{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -\mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{w}}\hfill \\ \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\mathrm{\Omega }$	(70)
$G^e{}_{{\omega }^y}$	$=2{\rho }_a{S}_a\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{u}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{w}}\hfill \\ \hfill -\mathit{G}_{\mathit{w}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]{\omega }^y$	(71)

The stiffness matrices obtained with (DQHFEM)

$K^e$	$=E_a{I}_{S_a}^{mo}\left[\begin{array}{cc}\hfill \mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(2)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(2)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(2)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(2)}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]$	(72)
$K^e{}_{\mathrm{\Omega }{\omega }^y}$	$=2{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill -\mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]\mathrm{\Omega }{\omega }^y$	(73)
$K^e{}_{\omega ^x{}^2}$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -\mathit{G}_{\mathit{w}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]\omega ^x{}^2$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill -\mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\omega ^x{}^2$	(74)
$K^e{}_{\omega ^y{}^2}$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill -\mathit{G}_{\mathit{u}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -\mathit{G}_{\mathit{w}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]\omega ^y{}^2$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill \mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]\omega ^y{}^2$	(75)
$K^e{}_{\omega ^z{}^2}$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill -\mathit{G}_{\mathit{u}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\omega ^z{}^2$
	$\mathrm{\hspace{1em}}+{\rho }_a{I}_{sa}^{mo}\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -\mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{w}}\hfill \end{array}\right]\omega ^z{}^2$	(76)
$K^e{}_{{\omega }^x{\omega }^z}$	$={\rho }_a{S}_a\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \mathit{G}_{\mathit{u}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{w}}\hfill \\ \hfill \mathit{G}_{\mathit{w}}{}^T\overline{\mathit{C}}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]{\omega }^x{\omega }^z+{\rho }_a{I}_{sa}^{mo}$
	$\mathrm{\hspace{1em}}\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill -\mathit{G}_{\mathit{u}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{w}}\hfill \\ \hfill -\mathit{G}_{\mathit{w}}{}^T\overline{\mathit{A}}^{(1)}{}^T\overline{\mathit{C}}{\overline{\mathit{A}}}^{(1)}{\mathit{G}}_{\mathit{u}}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]{\omega }^x{\omega }^z$	(77)

7 Equation of Elementary Motion

By applying Lagrange’s equations to the system discretised by the DQFEM or DQHFEM methods we obtain the following system of differential equations:

$[M^e]\{\ddot{q}\}+[[G^e]\mathrm{\Omega }+[G^e{}_{{\omega }^y}]{\omega }^y]\{\dot{q}\}+[[K^e]+[K^e{}_{\mathrm{\Omega }{\omega }^y}]\mathrm{\Omega }{\omega }^y$
$\mathrm{\hspace{1em}\hspace{1em}}+[K^e{}_{\omega ^x{}^2}]+[K^e{}_{\omega ^y{}^2}]+[K^e{}_{\omega ^z{}^2}]+[K^e{}_{{\omega }^x{\omega }^z}]]\{q\}=\{0\}$	(78)

Where:

• $[M^e]$, $[[G^e]\mathrm{\Omega }+[G^e{}_{{\omega }^y}]{\omega }^y]$ , and $[[K^e]+[K^e{}_{\mathrm{\Omega }{\omega }^y}]\mathrm{\Omega }{\omega }^y+[K^e{}_{\omega ^x{}^2}]+[K^e{}_{\omega ^y{}^2}]+[K^e{}_{\omega ^z{}^2}]+[K^e{}_{{\omega }^x{\omega }^z}]]$ are respectively the elementary matrices of mass, gyroscopic and stiffness with periodic and time-varying parametric terms due to the rotations of its rigid support.

• $\{\ddot{q}\}$, $\{\dot{q}\}$, and $\{q\}$ are respectively the global acceleration, velocity and displacement vectors suitable for DQFEM or DQHFEM connectivity.

• Put $[M]$, $[G]$ and $[K]$ respectively the total matrices after assembly of mass, damping and stiffness, therefore the differential equation of motion becomes:

$$[M]\{\ddot{q}\}+[G]\{\dot{q}\}+[K]\{q\}=\{0\}$$

(79)

The assembly of the global matrices is similar to that of the classic version of the finite elements method to ensure displacement and rotational compatibility at the nodes of adjacent elements.

In the current application of DQFEM and DQHFEM, boundary conditions are applied in the same way as the hp version of the finite element method. For example, if the shaft is supported by a rigid bearing it means that the displacement at the bearing is zero.

In this work, some numerical methods for solving the equations of motion and analysis of the dynamic vibration behaviour of the rotating shaft are used. The programming language used is MATLAB, in order to program these solving methods using the DQFEM and DQHFEM methods.

The dynamic behaviour of the embedded rotary shaft is analysed using Campbell diagrams.

To study the free motion of the rotating shaft the homogeneous differential system of Equation (78) of the equation of motion of the on-board rotating shaft rotating at constant speed is solved $\mathrm{\Omega }$.

In order to transform the second order differential system to a first order system, the following expression is added to Equation (78):

$$-[M]\{\dot{q}\}+[M]\{\dot{q}\}=\{0\}$$

(80)

The result is:

$\left[\begin{array}{cc}\hfill -[M]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill ([G])\mathit{\hspace{1em}}\hfill & \hfill [M]\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \{\dot{q}\}\hfill \\ \hfill \{\ddot{q}\}\hfill \end{array}\right\}+\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill \left[M\right]\hfill \\ \hfill ([K])\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \{q\}\hfill \\ \hfill \{\dot{q}\}\hfill \end{array}\right\}=\{0\}$

(81)

Equation (69) is a first-order linear differential system in the form:

$$[B_1]\{{\dot{Q}}_r\}+[B_2]\{Q_r\}=\{0\}$$

(82)

Where:

$[B_1]$	$=\left[\begin{array}{cc}\hfill -[M]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill ([G])\mathit{\hspace{1em}}\hfill & \hfill [M]\hfill \end{array}\right]$	(83)
	$=\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [M]\hfill \\ \hfill ([K])\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]$	(84)
$\mathrm{\{}{Q}_r\}$	$=\left\{\begin{array}{c}\hfill \{q\}\hfill \\ \hfill \{\dot{q}\}\hfill \end{array}\right\}$	(85)

Either:

$$\{{\dot{Q}}_r\}=-{[B_1]}^{-1}[B_2]\{Q_r\}=[A_r]\{Q_r\}$$

(86)

With:

$[A_r]=-{[B_1]}^{-1}[B_2]=-\left[\begin{array}{cc}\hfill {[M]}^{-1}\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill {[M]}^{-2}([G])\mathit{\hspace{1em}}\hfill & \hfill {[K]}^{-1}\hfill \end{array}\right]\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [M]\hfill \\ \hfill ([K])\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \end{array}\right]$

(87)

So:

$[A_r]=\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [I]\hfill \\ \hfill -{[M]}^{-1}([K])\mathit{\hspace{1em}}\hfill & \hfill -{[M]}^{-1}([G])\hfill \end{array}\right]$

(88)

The absolute values of the eigenvalues of the matrix $[A_r]$ represent the pulsations ${\omega }_r$ of the system, they depend on the rotation speed $\mathrm{\Omega }$, the evolution of the pulsations as a function of the rotation speed is called Campbell’s diagram, the latter also represents the evolution of natural frequencies $f_r={\omega }_r/2\pi$.

Figure 1 Flowchart of the simulation process using the hierarchical differential quadrature finite element method or the differential quadrature finite element method.

8 Results and Discussion

Figure 2 Model of on-board symmetric shaft (simply supported - simply supported).

8.1 Validation

This section presents the results calculated with the DQFEM and DQHFEM methods of vibration of an on-board shaft, and compares them with the results obtained with the h and hp versions of the finite element method.

Table 1 Natural frequencies of the first three modes of vibration as a function of $\mathrm{\Omega }$

Rotating Speed
$\mathrm{\Omega }$ [rad/s]	0
Bending	(Boukhalfa
Mode	and	(René,
$\omega$ [Hz]	Hadjoui, 2014)	1988)	DQFEM	DQHFEM	hFEM	hpFEM
1	128.1837	122.7475	122.6308	121.0328	122.7300	122.6307
2	507.0806	490.9899	489.1309	489.1310	494.9093	489.1309
3	1120.8937	1104.7273	1095.3825	1095.3826	1215.6495	1095.3825
Rotating Speed	${10}^4$
$\mathrm{\Omega }$ [rad/s]	Backward mode (B)
Bending	(Boukhalfa
Mode	and	(René,
$\omega$ [Hz]	Hadjoui, 2014)	1988)	DQFEM	DQHFEM	hFEM	hpFEM
1	122.3467	119.7548	122.3145	122.3145	122.4137	122.3144
2	484.6529	479.0191	487.8730	487.8729	493.6491	487.8729
3	1073.5822	1077.7931	1092.5786	1092.5785	1212.8087	1092.5785
Rotating Speed	${10}^4$
$\mathrm{\Omega }$ [rad/s]	Forward mode (F)
Bending	(Boukhalfa
Mode	and	(René,
$\omega$ [Hz]	Hadjoui, 2014)	1988)	DQFEM	DQHFEM	hFEM	hpFEM
1	134.2918	125.8150	122.9480	122.9479	123.0472	122.9479
2	530.4377	503.2598	490.3923	490.3922	496.1727	490.3922
3	1169.8098	1132.3346	1098.1939	1098.1939	1218.4970	1098.1938

A first validation is performed on a single shaft supported in both ends with a diameter of 0.05 (m), and a length of 0.9 (m), Young’s modulus and density are respectively 2.10${}^{11}$ (N/m${}^2$), 7800 (kg/m${}^3$), 7800 (kg/m${}^3$), 7800 (kg/m${}^3$) and 7800 (kg/m${}^3$).

For a first validation the shaft is modelled by a single element, the parameters taken for the DQFEM are (h $=$ 3 elements and N $=$ 20 sampling nodes), and for the DQHFEM the same parameters were also taken to see the deference of the result between the two methods.

For the h version of the finite element method a mesh size of h $=$ 3 elements was taken, and for the hp version of the finite element method a mesh size of h $=$ 3 elements and polynomial degree p $=$ 10 were taken. The results obtained in Table 1 for the boundary condition of a shaft simply supported at both ends are compared with the results obtained analytically by (René, 1988) and (Boukhalfa and Hadjoui, 2014).

After comparing the results with other methods in the literature, it can be seen that the difference between the methods is very small.

8.2 Convergence Study of the Results

According to Figures 3–4, we notice that the hierarchical differential quadrature finite element method is the differential quadrature finite element method, have the same speed of convergence when only the number of samples is varied.

Figure 3 Convergence of the first five bending frequencies as a function of the sampling number of DQFEM nodes N (the shaft is discretised into a single element). $\mathrm{\Omega }=0$.

Figure 4 Convergence of the first five bending frequencies as a function of the sampling number of DQHFEM nodes N (the shaft is discretised into a single element). $\mathrm{\Omega }=0$.

Figure 5 Convergence of the first frequency with a combination of the mesh degree-h and the sampling points N with (Ni $=$ Ni+1 and hi $=$ hi+3). With DQHFEM.

Figure 6 Convergence of the first frequency based on the sampling number N with a fixed number of elements. With DQFEM.

Table 2 Comparison of convergence results of the first natural frequency between DQHFEM, DQFEM, and hp-FEM

Number of Sampling Points N	Number of Mesh Elements h	First Frequency (Hz), Calculated with DQHFEM	First Frequency (Hz), Calculated with DQFEM	First Frequency (Hz), Calculated with hp-FEM (p $=$ p+1, h $=$ h+3)	h-p
4	1	136.108412638979	136.108412638979	135.674253916772	1–4
5	4	122.630935036219	122.632020652930	122.239767207697	4–5
6	7	122.630797035810	122.631723912041	122.239626983926	7–6
7	10	122.630796028383	122.630935036182	122.239626978018	10–7
8	13	122.630795010945	122.630807356777	122.239626986316	13–8
9	16	122.630794801053	122.630800401490	122.239626967131	16–9
10	19	122.630794790357	122.630794789885	122.239626995141	19–10
11	22	122.630794780622	122.630794787307	122.239626913630	22–11
12	25	122.630794779322	122.630794785638	122.239627031743	25–12
13	28	122.630794742756	122.630794767887	122.239626958780	28–13
14	31	122.630794737688	122.630794762705	122.239626309820	31–14
15	34	122.630794717245	122.630794701809	122.239626718141	34–15
16	37	122.630794690141	122.630788737133	122.239626758437	37–16
17	40	122.630794140855	122.630756419200	122.239625309707	40–17
18	43	122.630793689304	122.630721166728	122.239626470715	43–18
19	46	122.630793062697	122.630230141610	122.239622160122	46–19
20	49	122.630792789800	122.626336236740	122.239619867368	49–20

In Figures 5–6 a combination of sampling number N with the number of elements h (Ni $=$ Ni+1 and hi $=$ hi+3) was used. At each iteration one sampling degree and three additional elements were added. It can be seen that the two methods DQFEM and DQHFEM converge rapidly, with N $=$ 6 and h $=$ 7 being the starting point for convergence, making them effective. Table 2 shows the details of Figures 5–6 in comparison with the hp version of the finite element method. There is a slight discrepancy between the methods. The choice of the combination of sampling number N with the number of elements h is free, which means that any combination can be chosen, the higher the value, the more convergent the results.

8.3 Influence of Boundary Conditions

In rotor dynamics, the influence of boundary conditions is a point not to be neglected; in Figures 7–8 we show the difference between several cases of boundary conditions (Simply Supported – Simply Supported; Clamped – Clamped; Clamped – Clamped; Clamped – Free). Knowing that

• Simply supported: displacement at the node is zero.

• Clamped: displacement at node zero plus rotation at node zero.

Figure 7 Influence of four cases of boundary conditions on the dynamic behaviour of the rotating shaft (DQFEM).

Figure 8 Influence of four cases of boundary conditions on the dynamic behaviour of the rotating shaft (DQHFEM).

8.4 Interpretations of Campbell Diagrams

A calculation programme for determining the Eigen-modes of the on-board rotor has been developed and to validate this programme we have the rotor shown in Figure 9.

The Table 3 shows the characteristics of the symmetrical rotor used in the literature.

In this part we have a rotor with a symmetrical shaft and a disc, in the h and hp versions of the finite element method, the disc is modelled as a point mass see (Saimi, 2016).

Figure 9 On-board symmetrical rotor.

Table 3 Characteristics of the symmetrical rotor

Density of disc material	${\rho }_d=7800$ kg/m${}^3$
Outer radius of the disc	$r_d=0.15$ m
Disc thickness	$e_d=0.03$ m
Disc position	$y_d=0.4/3$ m
Density of shaft material	${\rho }_a=7800$ kg/m${}^3$
Shaft radius	$r_a=0.01$ m
Shaft length	$l_a=0.4$ m
Shaft YOUNG module	$E_a=2*{10}^{11}$ N/m${}^2$

In the DQHFEM and DQFEM methods the disc is modelled as if it is part of the shaft, i.e. a shaft section.

The symmetrical rotor is modelled by three identical Euler Bernoulli beam finite elements. The bearings are assumed to be rigid, generating a rotor simply supported at both ends. Thus the corresponding displacements are cancelled out.

Taking into account the rotor parameters and geometry, the following Table 5 shows the rotor elements in Figure 9:

Table 4 Elements of discretization

	First Element (Shaft)	Second Element (Disc)	Third Element (Shaft)
Outside Radius (m)	0.01	0.15	0.01
Inside radius (m)	0	0	0
Length (m)	${\displaystyle \frac{0.4}{3}}-{\displaystyle \frac{0.03}{2}}$	0.03	$2{\displaystyle \frac{0.4}{3}}-{\displaystyle \frac{0.03}{2}}$

Figure 10 Campbell’s diagram calculated with DQFEM.

Figure 11 Campbell’s diagram calculated with DQHFEM.

Table 5 Natural frequencies as a function of rotor and stationary support speed (B and F represent backward and forward precision modes respectively)

		Support fixe
$\mathrm{\Omega }$ rpm	$f_r$ Hz	DQFEM	DQHFEM	hp-FEM (Saimi, 2016)	h-FEM (Saimi, 2016)
0	$f_1{}^B$	49.96	49.99	45.57	45.61
	$f_1{}^F$	49.96	49.99	45.57	45.61
	$f_2{}^B$	145.78	145.78	125.95	125.96
	$f_2{}^F$	145.78	145.78	125.95	125.96
1500	$f_1{}^B$	46.92	46.93	42.90	42.93
	$f_1{}^F$	52.58	52.58	47.68	47.72
	$f_2{}^B$	126.80	126.80	106.79	106.79
	$f_2{}^F$	169.76	169.76	150.82	150.83
3000	$f_1{}^B$	43.46	43.47	39.70	39.72
	$f_1{}^F$	54.70	54.70	49.31	49.35
	$f_2{}^B$	112.53	112.52	92.98	92.99
	$f_2{}^F$	198.50	198.50	180.97	181
4500	$f_1{}^B$	39.80	39.81	36.16	36.18
	$f_1{}^F$	56.41	56.42	50.57	50.62
	$f_2{}^B$	102.20	102.19	83.56	83.57
	$f_2{}^F$	231.37	231.37	215.53	215.57
6000	$f_1{}^B$	36.17	36.18	32.60	32.61
	$f_1{}^F$	57.81	57.81	51.55	51.61
	$f_2{}^B$	94.87	94.87	77.31	77.32
	$f_2{}^F$	267.58	267.58	253.46	253.53

The results presented in Figures 10–11 and detailed in Table 5 show that the gaps between the methods are small and do not exceed 10%.

8.5 Influence of the Angular Velocities of the Moving Support on Critical Speeds

In this section, the angular velocities of the mobile support are varied along the axis $O_x$, $O_y$, and $O_z$, by deducting their influence on the system’s own critical speeds. Critical velocities are determined from the intersection of the Campbell diagram modes with the line $fr\text{é}quances(Hz)=\mathrm{\Omega }(\mathrm{Rpm})/60$.

According to Figures 12–14, it can be seen that the constant rotation of the rigid support on the axes $O_x$, $O_y$, and $O_z$, influences the dynamic behaviour of the rotor. The angular velocities of the support along the axes $O_x$ (${\omega }^x)$, and $O_z$ (${\omega }^z)$, do not have an influence on the gyroscopic effect, which is the opposite in case of the angular velocity of the support along the axis $O_y$ (${\omega }^y)$ which explains the decrease in critical forward mode speeds and the increase in critical reverse mode speeds (Figure 11).

The angular velocities of the mobile support cause excitations that can make the rotor movement sensitive, resulting in the decrease of critical speeds.

Table 6 Critical speed of the intersection of the first backward mode with the line $f=\mathrm{\Omega }$/60

	Backward Critical Speed
Constant Rotation of the Support	Present	hp-FEM (Saimi, 2016)
${\omega }^x=25Hz$	2653 rpm	2457 rpm
${\omega }^y=25Hz$	2802 rpm	2810 rpm

The results in Figures 12–14 are obtained with the DQHFEM method since they are close to those of the DQFEM method.

Figure 12 Evolution of the critical speed of the intersection of the first backward and Forward modes with the line $f=\mathrm{\Omega }$/60 depending on the constant rotation of the support along the axis $O_x$ (${\omega }^x)$ DQHFEM.

Figure 13 Evolution of the critical speed of the intersection of the first backward and Forward modes with the line $f=\mathrm{\Omega }$/60 depending on the constant rotation of the support along the axis $O_y$ (${\omega }^y)$ DQHFEM.

Figure 14 Evolution of the critical speed of the intersection of the first backward and Forward modes with the line $f=\mathrm{\Omega }$/60 depending on the constant rotation of the support along the axis $O_z$ (${\omega }^z)$ DQHFEM.

9 Conclusion

The analysis of the vibrational behaviour of embedded rotor shafts is treated in this paper, using methods which are the differential quadrature finite element method (DQFEM) and the differential quadrature hierarchical finite element method (DQHFEM). The use of these methods is new in the context of rotor dynamics.

The calculation of the kinetic energy and stiffness energy of the on-board rotor shaft, the application of the Lagrange method, and the introduction of the characteristics of the differential quadrature finite element method (DQFEM) and the differential quadrature hierarchical finite element method (DQHFEM), gave us the equations of motion of the system. A programme for calculating natural frequencies and critical speeds was developed using the MATLAB programming language.

The results of an example tree in the literature were well validated with the DQFEM and DQHFEM methods, and the difference between the methods is small, which explains the efficiency of the methods. Convergence is obtained for a low sampling and item number compared to h-FEM, hp-FEM

Several examples were treated and this allowed us to determine the influence of different geometrical parameters of the on-board rotor as well as the influence of the movement of the support on the rotor’s behaviour.

This work allowed us to draw the following conclusions:

• The matrices of the DQFEM and DQHFEM methods are somewhat similar to those of the hp version of the finite element method.

• The calculation times using the DQFEM and DQHFEM methods are significantly faster compared to the hp version of the finite element method.

• The convergence of the results can be controlled by increasing the number of samples and the number of elements.

• The difference in results between the two methods (DQFEM and DQHFEM) is very small.

• The Legendre polynomial is used in this work, because of its very close approximation to deformed trees.

• The DQFEM and DQHFEM methods have the advantages of a simple mathematical principle, a fast convergence speed, high computational accuracy, low computation quantity and lower memory requirements, etc. According to the results obtained in this work.

• The movement of the rotor of the embedded support amplifies the gyroscopic effect caused by the coupling of the movement perpendicular to the axis of rotation, and creates an asymmetry in the movement of the rotor.

• In the case of disturbances in the vicinity of the first critical speed, the influence of the rotation of the support around the axes $O_x$ (${\omega }^x)$, $O_z$ (${\omega }^z)$, and $O_y$ (${\omega }^y)$ is very important and should not be neglected. As the angular speeds of the support increase, the critical speeds decrease.

• Following the results obtained in this work, it was concluded that the rotation of the support along an axis parallel to the axis of rotation of the rotor has a great influence on the rigidity of the system and in particular on the gyroscopic effect.

• The geometrical parameters also have an influence on the vibration behaviour of the on-board rotor. It can be seen that it causes frequency variations and therefore critical speed variations depending on the drive positions.

Appendix

a. ${\eta }_1(C,A^{(1)})$, ${\eta }_2(C,A^{(1)})$, and ${\eta }_3(C,A^{(1)})$

${\eta }_1(C,A^{(1)})$	$=({\dot{\overline{u}}}^{T}C\dot{\overline{u}}+{\dot{\overline{w}}}^{T}C\dot{\overline{w}}+2(({\dot{x}}_0+z_0{\omega }^y-(y_0+y){\omega }^z))$
	$\mathrm{\hspace{1em}\hspace{1em}}({\dot{\overline{u}}}^{T}C+{\overline{w}}^{T}C{\omega }^y)+2(({\dot{z}}_0+(y_0+y){\omega }^x-x_0{\omega }^y))$
	$\mathrm{\hspace{1em}\hspace{1em}}({\dot{\overline{w}}}^{T}C-{\overline{u}}^{T}C{\omega }^y)-2(({\dot{y}}_0-z_0{\omega }^x+x_0{\omega }^z))$
	$\mathrm{\hspace{1em}\hspace{1em}}({\overline{w}}^{T}C{\omega }^x-{\overline{u}}^{T}C{\omega }^z)+2({\dot{\overline{u}}}^{T}C\overline{w}-{\dot{\overline{w}}}^{T}C\overline{u}){\omega }^y$
	$\mathrm{\hspace{1em}\hspace{1em}}+{\overline{w}}^{T}C\overline{w}\omega ^x{}^2+({\overline{u}}^{T}C\overline{u}+{\overline{w}}^{T}C\overline{w})\omega ^y{}^2$
	$\mathrm{\hspace{1em}}+{\overline{u}}^{T}C\overline{u}\omega ^z{}^2-2{\overline{u}}^{T}C\overline{w}{\omega }^x{\omega }^z)$	(89)
${\eta }_2(C,A^{(1)})$	$=({\dot{\overline{u}}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\dot{\overline{u}}+{\dot{\overline{w}}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\dot{\overline{w}}$
	$\mathrm{\hspace{1em}}+2C{A}^{(1)}\dot{\overline{w}}{\omega }^x+2({\dot{\overline{u}}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}$
	$\mathrm{\hspace{1em}}-{\dot{\overline{w}}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{u})\omega {}^y-2CA{}^{(1)}\dot{\overline{u}}\omega {}^z$
	$\mathrm{\hspace{1em}}-({\overline{u}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{u}-1)\omega ^x{}^2$
	$\mathrm{\hspace{1em}}+({\overline{u}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{u}+{\overline{w}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w})\omega ^y{}^2$
	$\mathrm{\hspace{1em}}-({\overline{w}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}-1)\omega ^z{}^2-2{\overline{u}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}{\omega }^x{\omega }^z$
	$\mathrm{\hspace{1em}}-2(C{A}^{(1)}\overline{u}{\omega }^x+C{A}^{(1)}\overline{w}{\omega }^z){\omega }^y)$	(90)
${\eta }_3(C,A^{(1)})$	$=({(\mathrm{\Omega }+{\omega }^y)}^2+{\overline{u}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{u}\omega ^x{}^2$
	$\mathrm{\hspace{1em}}+{\overline{w}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}\omega ^z{}^2+(\mathrm{\Omega }+{\omega }^y)$
	$\mathrm{\hspace{1em}}(-2{\dot{\overline{u}}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}$
	$\mathrm{\hspace{1em}}-({\overline{u}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{u}+{\overline{w}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}){\omega }^y$
	$\mathrm{\hspace{1em}}+2(C{A}^{(1)}\overline{u}{\omega }^x+C{A}^{(1)}\overline{w}{\omega }^z))$
	$\mathrm{\hspace{1em}}+2{\overline{u}}^{T}A^{(1)}{}^{T}C{A}^{(1)}\overline{w}{\omega }^x{\omega }^z)$	(91)

b. Transformation matrix Q (DQFEM)

Figure 15 Sampling points of a beam element.

Figure 15 Represents sampling points of order N of a beam element. In order to construct an element which satisfies the requirements for continuity between elements, the sampling points must be in the following form:

Figure 16 Sampling points of a beam element.

The transformation between Figures 15 and 16 is made using the DQ rule,

$u=\left(\begin{array}{c}\hfill u_1\hfill \\ \hfill u^{\prime }{}_1=A_{1,1}^{(1)}{u}_1+A_{1,2}^{(1)}{u}_2+A_{1,3}^{(1)}{u}_3\hfill \\ \hfill +\mathrm{\cdots }+A_{1,N-1}^{(1)}{u}_{N-1}+A_{1,N}^{(1)}{u}_N\hfill \\ \hfill u_3\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill u_{N-2}\hfill \\ \hfill u_N\hfill \\ \hfill u^{\prime }{}_N=A_{N,1}^{(1)}{u}_1+A_{N,2}^{(1)}{u}_2+A_{N,3}^{(1)}{u}_3\hfill \\ \hfill +\mathrm{\cdots }+A_{N,N-1}^{(1)}{u}_{N-1}+A_{N,N}^{(1)}{u}_N\hfill \end{array}\right)$

(92)

We can write this vector as a product of a matrix by a vector

$u=\left[\begin{array}{cccccc}\hfill 1\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill A_{1,1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{1,2}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{1,3}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill A_{1,N-1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{1,N}^{(1)}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 1\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\ddots }\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 1\hfill \\ \hfill A_{N,1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{N,2}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{N,3}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill A_{N,N-1}^{(1)}\mathit{\hspace{1em}}\hfill & \hfill A_{N,N}^{(1)}\hfill \end{array}\right]\left(\begin{array}{c}\hfill u_1\hfill \\ \hfill u^{\prime }{}_1\hfill \\ \hfill u_3\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill u_{N-2}\hfill \\ \hfill u_N\hfill \\ \hfill u^{\prime }{}_N\hfill \end{array}\right)$

(93)

Acknowledgements

We acknowledge with grateful thanks the support by the laboratory of mechanical and material systems engineering, as well as the General Directorate of Scientific Research and Technological Development of the Ministry of Higher Education of Algeria.

We also want to thank the effort made by European Journal of Computational Mechanics team, and the reviewers.

References

Babuska and B.Q, G., 1992. The h, p and h-p version of the finite element method: basis theory and applications. Advances in Engineering Software, 15(3–4), pp. 159–174.

Babuska, I., B A, S. and I N, K., 1981. The p-Version of the Finite Element Method. SIAM Journal on Numerical Analysis, 18(3), pp. 515–545.

Babuška, I. and Suri, M., 1990. The p- and h-p versions of the finite element method, an overview. Computer Methods in Applied Mechanics and Engineering, 80(1–3), pp. 5–26.

Bardell, N. S., 1996. An engineering application of the hp version of the finite element method to the static analysis of a Euler-Bernoulli beam. Computers & structures, 59(2), pp. 195–211.

Belhadj, A., Boukhalfa, A. and Belalia, S. A., 2020. Free vibration investigation of single walled carbon nanotubes with rotary inertia. Nanomaterials Science & Engineering, 2(3), pp. 103–112.

Bellman, R. and Casti, J., 1971. Differential quadrature and long term integration. Journal of Mathematical Analysis and Applications, Volume 34, pp. 235–238.

Bert, C. and M, M., 1996. The differential quadrature method for irregular domains and application to plate vibration. International Journal of Mechanical Sciences, 38(6), pp. 589–606.

Bert, C. W. J. S. K. a. S. A. G., 1988. Two new approximate methods for analyzing free vibration of structural components. AIAA Journal, Volume 26, pp. 612–618.

Bert, C. W. and Malik, M., 1996. Differential quadrature method in computational mechanics: A review. Applied Mechanics Reviews, Volume 49, pp. 1–28.

Beslin, O. and Nicolas, J., 1997. A hierarchical functions set for predicting very high order plate bending modes with any boundary conditions. Journal of Sound and Vibration, 202(5), pp. 633–655.

Bokaian, A., 1988. Natural frequencies of beams under compressive axial load. Journal of Sound and Vibration, Volume 126, pp. 49–65.

Boukhalfa, A. and Hadjoui, A., 2014. Dynamic analysis of a spinning functionally graded material shaft by the p–version of the finite element method. Latin American Journal of Solids and Structures.

Chen and New, C., 1999. A differential quadrature finite element method. Applied mechanics in the Americas. Proceedings of the 6th Pan-American Congress of Applied Mechanics and 8th International Conference on Dynamic Problems in Mechanics PACAM VI, Rio de Janeiro, Brazil; US, pp. 305–308.

Choi, I. S., Pierre, C. and Ulsoy, A. G., 1992. Consistent modeling of rotating Timoshenko shaft subject to axial loads. Journal of Vibration and Acoustics, Volume 114, pp. 249–259.

Clough, R. W., 1960. The finite element method in plane stress analysis Conference on Electronic Computation, Pittsburgh, PA. Proceedings of the Second American Society of Civil Engineers, pp. 345–378.

Cuiyun, L. et al., 2016. A differential quadrature hierarchical finite element method and its applications to vibration and bending of Mindlin plates with curvilinear domains. International Journal for Numerical Methods in Engineering, 109(2), pp. 174–197.

Dakel, M., 2013. Steady-state dynamic behaviour of an on-board rotor under combined base motions. Journal of Vibration and Control, Volume 20, pp. 2254–2287.

Duchemin, M., Berlioz, A. and Ferraris, G., 2006. Dynamic behavior and stability of a rotor under base excitation. Journal of Vibration and Acoustics, Volume 128, pp. 576–585.

Eshlenman, R. L. and Eubanks, R. A., 1967. On the critical speeds of a continuous shaft-disk system. Journal of Engineering for Industry, Volume 80, pp. 645–652.

Fellah, A., Hadjoui, A. and Bekhaled, B. S. A., 2019. Study of the Effect of an Open Transverse Crack on the Vibratory Behavior of Rotors Using the hp Version of the Finite Element Method. Journal of Solid Mechanics, 11(1), pp. 181–200.

Green, R. B., 1948. Gyroscopic effects on the critical speeds of flexible rotors,. Transaction of the American society of Mechanic Engineers., Volume 70, pp. 309–376.

Hassan, A., Abdelhamid, H. and Saimi, A., 2020. Numerical analysis on the dynamics behavior of FGM rotor in thermal environment using hp finite element method. Mechanics Based Design of Structures and Machines, pp. 1–24.

Heuveline, V. and Rannacher, R., 2003. Duality-Based Adaptivity in the Hp-Finite Element Method. Journal of Numerical Mathematics, 11(2), pp. 95–113.

Malik, M. and Bert, C., 2000. Vibration analysis of plates with curvilinear quadrilateral planforms by DQM using blending functions. Journal of Sound and Vibration, 230(4), pp. 949–954.

Nelson, H. D. and McVaugh, J. M., 1976. The dynamics of rotor-bearing systems using finite elements. ASME Journal of Engineering for Industry, Volume 98, pp. 593–600.

O. C. Zienkiewicz, M. L., 1977. The finite element method, 3rd edn. Wiley Online Library, pp. 1054–1054.

Oden, J. T. a. D. L., 1991. h-p adaptive finite element methods in computational fluid dynamics. Computer Methods in Applied Mechanics and Engineering, Volume 89, pp. 11–40.

Peano, A., 1976. Hierarchies of conforming finite elements for plane elasticity and plate bending.. Computers & Mathematics with Applications, Volume 2, pp. 211–224.

Petyt, M., 2010. Introduction to Finite Element Vibration Analysis (2nd edn). New York: Cambridge University Press.

Rao, S. S., 1983. Rotor dynamics. New York: NY: Wiley.

René, J. G., 1988. Vibrations des structures : interactions avec les fluides, sources d’excitation aléatoires. Paris: Eyrolles.

Ri, K. et al., 2020. Nonlinear forced vibration analysis of composite beam combined with DQFEM and IHB. AIP Advances, 10(8).

Saimi, A. H. A., 2016. An engineering application of the h-p version of the finite elements method to the dynamics analysis of a symmetrical on-board rotor. European Journal of Computational Mechanics, pp. 388–416.

Sajal, S. S. et al., 2018. Dynamic analysis of microbeams based on modified strain gradient theory using differential quadrature method. European Journal of Computational Mechanics , 27(3), pp. 187–203.

Shu, C., 2000. Differential Quadrature and its Application in Engineering. Springer-Verlag d. London: s.n.

Striz, A. G., Chen, W. L. and Bert, C. W., 1995. High accuracy plane stress and plate elements in the quadrature element method. Proceedings of the 36th AIAA/ASME/ASCE/AHS/ASC, pp. 957–965.

Striz, A. G., Chen, W. L. and Bert, C. W., 1997. Free vibration of plates by the high accuracy quadrature element method. Journal of Sound and Vibration, Volume 202, pp. 689–702.

Suri, M., 1997. A reduced constraint h$\mathrm{?}$$\mathrm{?}$ finite element method for shell problems. Mathematics of computation, 66(217), pp. 15–29.

Szabo, B., 1979. Some recent developments in finite element analysis. Computers & Mathematics with Applications , 5(2), pp. 99–115.

Szabo, B. and Prob, 1985. Theoretical manual release 1. St Louis, Missouri: Noetic Technologies Corporation.

Weiyan, Z., Feng, G. and Yongsheng, R., 2019. Generalized Differential Quadrature Method for Free Vibration Analysis of a Rotating Composite Thin-Walled Shaft. Mathematical Problems in Engineering, Volume 2019.

Xing, Y. F. and Liu, B., 2009. High-accuracy differential quadrature finite element method and its application to free vibrations of thin plate with curvilinear domain. International Journal for Numirical methods in engineering. Vol 80. issue 13. pp. 1718–1742.

Yufeng, X., Bo, L. and Guang, L., 2010. A Differential Quadrature Finite Element Method. International Journal of Applied Mechanics, 2(1), pp. 207–227.

Zhong, H. and Yu, T., 2009. A weak form quadrature element method for plane elasticity problems. Applied Mathematical Modelling, 33(10), pp. 3801–3814.

Zienkiewicz, O., De, S., Gago, J. and Kelly, D., 1983. The hierarchical concept in finite element analysis. Computers & Structures, 16(1–4), pp. 53–65.

Zienkiewicz, O., Irons, B., Scott, F. and Campbell, J., 1970. Three-dimensional stress analysis. Proceedings of IUTAM Symposium on High Speed Computing of Elastic Structures, Liege, pp. 413–431.

Biographies

Saimi Ahmed obtained his Ph.D in Mechanics of Materials and Structures from the University of Tlemcen, Algeria, in 2017. He is currently a Senior Lecturer at the National High School of Hydraulics Blida, Algeria. A researcher member of Mechanical Systems and Structural Engineering Laboratory, IS2M/UABT. His research interests are: Finite element methods, Structural vibration, Structural dynamics, Dynamics of rotors, Dynamics of rotating machines, computational mechanics, FG materials, Composite materials.

Hadjoui Abdelhamid obtained his Ph.D in mechanical engineering from the University of Tlemcen, Algeria. He is currently a professor at the University of Tlemcen, Algeria (UABT). Research Director in Mechanical Systems and Structural Engineering Laboratory, IS2M/UABT. His research interests are: Materials Engineering, Structural Engineering, Mechanical Engineering, Structural Analysis, Finite elements Modeling, Structural Dynamics, Simulation, Dynamic Analysis, Modal Analysis, Structural Vibration, Vibration Analysis.

Bensaid Ismail received his B.Sc, M.Sc and Ph.D degrees in Mechanical Engineering from Abou Beckr Belkaid University Tlemcen, Algeria. He is currently working in the level of the Mechanical engineering department at the same University. Dr. Bensaid does research in Mechanical and structural Engineering, Materials, Composite, Maintenance, Nanostructures and Dynamical Systems. He, as an author/co-author, has published more than 18 articles in various journals.

Fellah Ahmed obtained a Ph. D in rehabilitation and reliability of structures and mechanical equipment from the University of Tlemcen, Algeria, in 2019. He is currently a research associate in the laboratory of engineering of mechanical systems and structures, IS2M/UABT. His research interests are: Finite element methods, Structural vibration, Structural dynamics, Dynamics of rotors, Dynamics of rotating machines, computational mechanics, FG materials, Composite materials.

European Journal of Computational Mechanics, Vol. 29_4-6, 303–344.
doi: 10.13052/ejcm1779-7179.29461
© 2021 River Publishers