Effect of Thrust on the Structural Vibrations of a Nonuniform Slender Rocket

Abstract

The vibration characteristics of a nonuniform, flexible and free-flying slender rocket experiencing constant thrust is investigated. The rocket is idealized as a classic nonuniform beam with a constant one-dimensional follower force and with free-free boundary conditions. The equations of motion are derived by applying the extended Hamilton’s principle for non-conservative systems. Natural frequencies and associated mode shapes of the rocket are determined using the relatively efficient and accurate Adomian modified decomposition method (AMDM) with the solutions obtained by solving a set of algebraic equations with only three unknown parameters. The method can easily be extended to obtain approximate solutions to vibration problems for any type of nonuniform beam.

1. Introduction

The development of rockets has led to larger values of thrust-to-weight and length-to-diameter ratios as required for longer range flights. The drag force is proportional to the square of the rocket diameter, and therefore the preference is to vary the length as opposed to the diameter for a given range increase. In addition, to reduce the cost of the launching operation and initial handling and transportation, efforts are always made to reduce the rocket structural weight. These different requirements lead to a highly flexible structure, for which, dynamic response and vibrational characteristics are of great importance [1]. In general, forces acting on such a structure can be divided into conservative and non-conservative forces with the follower force (thrust) a typical example of non-conservative forces. When a given structure is both under a constant follower force and whose direction changes according to the deformation of the structure, it can undergo static divergence whereby transverse natural frequencies merge into zero and flutter, where two natural frequencies coincide resulting in the amplitude of vibration growing without bound [2].

The structural vibrations of launch vehicles have been the subject of research using various methods of calculation. The slender freely flying rocket is usually simplified as a beam with free-free boundary conditions with a follower force [3,4]. One of the first investigations of divergence and flutter was by Beal [5] where a flexible missile under constant and pulsating thrust was considered. The Galerkin method was used and it was shown that the divergence found was closely related to the first two bending moments. Wu [6,7] was the first to use a finite element technique together with the unconstrained variational approach to investigate the relation between the critical load and the eigenvalues of the system. He concluded that the optimized size and location of the concentrated mass would improve directional stability. The divergence and flutter of a spinning beam with a concentrated mass and subjected to a pulsating thrust has been investigated by Yoon and Kim [8], where they modeled the beam as a Timoshenko formulation and solved by applying a finite element method. More recently, slender rocket studies have included aerodynamic forces, gravitational forces, and internal load effects [9,10,11,12]. However, research into this area is found to be limited, with very few sustained investigations being undertaken. The present research sets a framework for future work involving more complete geometry and material discontinuities, and the ability to solve nonlinear equations, with ease, using Adomian polynomials [13].

The method of solution chosen here is the Adomian modified decomposition method [13] which is a wide-ranging method of solution for problems involving algebraic [14], differential [15], integro-differential [16], and partial differential equations [17]. Specific to this work, the Adomian decomposition and Adomian modified decomposition method have been used by several groups [18,19,20,21,22,23] for uniform and nonuniform beams, starting with either the Euler–Bernoulli or Timoshenko formulations. Mao [18] applied the Adomian modified decomposition method (AMDM) to rotating uniform beams and included a centrifugal stiffening term, Adair and Jaeger [19] applied the AMDM to rotating nonuniform beams which also included a centrifugal stiffening term, whereas Hsu et al. [20] applied the AMDM to Timoshenko beams. Lai et al. [21] and Lai et al. [22] employed both the AMDM and Adomian decomposition methods as an innovative eigenvalue solver to determine the free vibration of an Euler–Bernoulli beam under various supporting conditions. Yaman [23] used the Adomian decomposition method to investigate the influence of the orientation effect on the natural frequency of a cantilever beam carrying a tip mass, while Adair and Jaeger [24] investigated the vibrations of a beam with both a transversely and axially eccentric tip mass present using the AMDM.

The objective of the present work is to investigate the effect of the following thrust on the vibrations of a slender rocket. In order to meet this objective, the work idealizes the rocket as a nonuniform beam and the thrust as a follower force with the governing equations developed using the extended Hamilton’s principle. The AMDM is developed to calculate a solution for the dynamic analysis.

2. Structural Modeling

An assumption is made that the effect of the thrust acts only on the rocket body and does not have an influence on the fins. The body of the rocket is assumed to be a nonuniform non-extensible beam subjected to a constant follower thrust, as shown on Figure 1.

When the deformation is small,

$$u\left(x,t\right)=\underset{0}{\overset{x}{{\displaystyle \int }}}-\frac{1}{2}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\mathrm{d}x,$$

(1)

where $u\left(x,t\right)$ and $w\left(x,t\right)$ are the longitudinal and transverse flexible displacements along the beam respectively. Hence,

$$\frac{\partial u\left(x,t\right)}{\partial x}=-\frac{1}{2}\left(\frac{\partial w\left(x,t\right)}{\partial x}\right).$$

(2)

The kinetic energy and potential energy of the system are, respectively,

$$T=\frac{1}{2}\underset{0}{\overset{l}{{\displaystyle \int }}}\rho A\left(x\right){\left(\frac{\partial w\left(x,t\right)}{\partial t}\right)}^2\mathrm{d}x,$$

(3)

$$U=\frac{1}{2}\underset{0}{\overset{l}{{\displaystyle \int }}}EI\left(x\right){\left(\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)}^2\mathrm{d}x,$$

(4)

where $l$ is length of the rocket.

The follower thrust creates virtual work and work of inertia in the lateral direction as

$$\delta W_P=P\left({\frac{\partial w\left(x,t\right)}{\partial x}\delta w\left(x,t\right)|}_0^l-\underset{0}{\overset{l}{{\displaystyle \int }}}\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\delta w\left(x,t\right)\mathrm{d}x-\frac{\partial w\left(l,t\right)}{\partial x}\delta w\left(l,t\right)\right) =-P\left({\frac{\partial w\left(x,t\right)}{\partial x}\delta w\left(x,t\right)|}_{x=0}+\underset{0}{\overset{l}{{\displaystyle \int }}}\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\delta w\left(x,t\right)\mathrm{d}x\right).$$

(5)

Therefore,

$$W=\underset{0}{\overset{l}{{\displaystyle \int }}}\rho A\left(x\right)f\left(\underset{0}{\overset{x}{{\displaystyle \int }}}-\frac{1}{2}{\left(\frac{\partial w\left(x,t\right)}{\partial x}\right)}^2\mathrm{d}x\right)\mathrm{d}x,$$

(6)

where $f$ is longitudinal acceleration.

According to Hamilton’s principle for a non-conservative system,

$$\delta \underset{t_a}{\overset{t_b}{{\displaystyle \int }}}\left(T+W-U\right)\mathrm{d}t+\underset{t_a}{\overset{t_b}{{\displaystyle \int }}}\delta W_P\mathrm{d}t=0.$$

(7)

On putting Equations (3)–(6) into Equation (7) and assuming that for any time interval $\left\{t_a,t_b\right\}$ that the virtual displacement is equal to zero, the following is obtained:

$$\frac{{\partial }^2}{\partial x^2}\left(EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)+\rho A\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}+\frac{\partial }{\partial x}\left(\left(P-f\underset{x}{\overset{l}{{\displaystyle \int }}}\rho A\left(x\right)dx\right)\frac{\partial w\left(x,t\right)}{\partial x}\right)=0.$$

(8)

According to Sugiyama and Oshima [25], when a rocket is in its positive phase, i.e., when it is driven by the thrust of the rocket, the drag as compared with the thrust is negligible and can be neglected. Therefore, in this work the acceleration can be considered as a function of thrust and mass of the rocket. So,

$$P-f\underset{x}{\overset{l}{{\displaystyle \int }}}\rho A\left(x\right)\mathrm{d}x=P-\frac{P}{{\displaystyle {{\displaystyle \int }}_0^l}\rho A\left(x\right)\mathrm{d}x}\underset{x}{\overset{l}{{\displaystyle \int }}}\rho A\left(x\right)\mathrm{d}x=\frac{P}{{\displaystyle {{\displaystyle \int }}_0^l}\rho A\left(x\right)\mathrm{d}x}\underset{0}{\overset{x}{{\displaystyle \int }}}\rho A\left(x\right)\mathrm{d}x=N\left(x\right),$$

(9)

giving

$$\frac{{\partial }^2}{\partial x^2}\left(EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)+\rho A\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}+\frac{\partial }{\partial x}\left(N\left(x\right)\frac{\partial w\left(x,t\right)}{\partial x}\right)=0.$$

(10)

It is assumed that the cross-sectional area and moment of inertia of the rocket body vary according to

$$A\left(x\right)=A_0{\left(1+\left(\alpha -1\right)\frac{x}{l}\right)}^2,$$

(11)

$$I\left(x\right)=I_0{\left(1+\left(\alpha -1\right)\frac{x}{l}\right)}^4,$$

(12)

where $\alpha =b_1/b_0$ is the ratio of the rocket cross-sectional diameters just before the nose cone and at the base. Here, $A_0=\pi b_0{}^2/4$ and $I_0=\pi b_0{}^4/64$ are the cross-sectional area and moment of inertia at $x=0$, respectively.

After setting $\beta =1-\alpha$, and setting $M={{\displaystyle \int }}_0^l\rho A\left(x\right)\mathrm{d}x$, $N\left(x\right)$ in Equation (10) can be written as

$$\frac{P\rho A_0}{M}\underset{0}{\overset{x}{{\displaystyle \int }}}\left(1-2\beta \frac{x}{l}+{\beta }^2{\left(\frac{x}{l}\right)}^2\right)\mathrm{d}x,$$

(13)

which leads to

$$N\left(x\right)=\frac{P\rho A_0}{3M}\left[3x-\frac{3\beta }{l}{x}^2+\frac{{\beta }^2}{l^2}{x}^3\right].$$

(14)

In addition, according to modal analysis for harmonic free vibration, $w\left(x,t\right)$ can be separated in space and time as

$$w\left(x,t\right)=\phi \left(x\right)h\left(t\right),$$

(15)

where $\phi \left(x\right)$ is the modal deflection and $h\left(t\right)$ is the harmonic function of time. If $\omega$ denotes the circular frequency of $h\left(t\right)$, then ${\partial }^{2}w\left(x,t\right)/\partial t^2=-{\omega }^2\phi \left(x\right)h\left(t\right)$ and the eigenvalue of the problem of Equation (10) reduces to the differential equation

$$\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\left[{\left(1-\beta \frac{x}{l}\right)}^4\frac{{\mathrm{d}}^2\phi \left(x\right)}{\mathrm{d}{x}^2}\right]+\frac{P\rho A_{0}l}{3E{I}_{0}M}\frac{\mathrm{d}}{\mathrm{d}x}\left[\left[3\frac{x}{l}-3\beta {\left(\frac{x}{l}\right)}^2+{\beta }^2{\left(\frac{x}{l}\right)}^3\right]\frac{\mathrm{d}\phi \left(x\right)}{\mathrm{d}x}\right]-\frac{\rho A_0{\omega }^2}{E{I}_0}{\left(1-\beta \frac{x}{l}\right)}^2\phi \left(x\right)=0.$$

(16)

Without loss of generality, the following dimensionless quantities are introduced:

$$\eta =\frac{x}{l}, \phi \left(\eta \right)=\frac{\phi \left(x\right)}{l}, {\lambda }_1=\frac{\rho A_0{\omega }^2{l}^4}{E{I}_0}, {\lambda }_2=\frac{P\rho A_0}{3E{I}_{0}M}.$$

Equation (16) can now be written as

$$\frac{{\mathrm{d}}^2}{\mathrm{d}{\eta }^2}\left[{\left(1-\beta \eta \right)}^4\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}\right]+{\lambda }_2\frac{\mathrm{d}}{\mathrm{d}\eta }\left[\left[3\eta -3\beta {\eta }^2+{\beta }^2{\eta }^3\right]\frac{\mathrm{d}\phi \left(\eta \right)}{\mathrm{d}\eta }\right]-{\lambda }_1{\left(1-\beta \eta \right)}^2\phi \left(\eta \right)=0.$$

(17)

Equation (17) is expanded as

$$\frac{{\mathrm{d}}^4\phi \left(\eta \right)}{\mathrm{d}{\eta }^4}-\frac{8\beta }{\left(1-\beta \eta \right)}\frac{{\mathrm{d}}^3\phi \left(\eta \right)}{\mathrm{d}{\eta }^3}+\left[\frac{12{\beta }^2}{{\left(1-\beta \eta \right)}^2}+{\lambda }_2\frac{\left(3\eta -3\beta {\eta }^2+{\beta }^2{\eta }^3\right)}{{\left(1-\beta \eta \right)}^4}\right]\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}-\frac{{\lambda }_1}{{\left(1-\beta \eta \right)}^2}\phi \left(\eta \right)=0.$$

(18)

3. Boundary Conditions

Generally, the boundary conditions can be considered as

$${EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}|}_{x=0}={EI\left(x\right)\frac{{\partial }^{3}w\left(x,t\right)}{\partial x^3}|}_{x=l}=0,$$

(19)

$${\frac{\partial }{\partial x}\left(EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)|}_{x=0}={\frac{\partial }{\partial x}\left(EI\left(x\right)\frac{{\partial }^{3}w\left(x,t\right)}{\partial x^3}\right)|}_{x=l}=0.$$

(20)

With reference to Figure 2, the boundary conditions, which do not vary with time and considered in dimensionless form, can be written as follows:

$$\begin{array}{l}{\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}|}_{\eta =0}-K_{RL}{\frac{\mathrm{d}\phi \left(\eta \right)}{\mathrm{d}\eta }|}_{\eta =0}=0,\\ {\frac{{\mathrm{d}}^3\phi \left(\eta \right)}{\mathrm{d}{\eta }^3}|}_{\eta =0}-4\beta {\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}|}_{\eta =0}+K_{TL}{\phi \left(\eta \right)|}_{\eta =0}=0,\end{array}$$

(21)

$$\begin{array}{l}{\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}|}_{\eta =1}+K_{RR}{\frac{\mathrm{d}\phi \left(\eta \right)}{\mathrm{d}\eta }|}_{\eta =1}=0,\\ {\frac{{\mathrm{d}}^3\phi \left(\eta \right)}{\mathrm{d}{\eta }^3}|}_{\eta =1}-\frac{4\beta }{1-\beta }{\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}|}_{\eta =1}-K_{TR}{\phi \left(\eta \right)|}_{\eta =1}=0,\end{array}$$

(22)

where the coefficients are in nondimensionalized form with

$$K_{TL}=\frac{k_{TL}{l}^3}{E{I}_0}, K_{TR}=\frac{k_{TR}{l}^3}{E{I}_1}, K_{RL}=\frac{k_{RL}l}{E{I}_0}, K_{RR}=\frac{k_{RR}l}{E{I}_1},$$

where $I_1$ is the moment of inertia at $x=l.$

4. Application of AMDM

According to the AMDM, $\phi \left(\eta \right)$ in Equation (18) can be expressed as an infinite series, i.e.,

$$\phi \left(\eta \right)={\displaystyle \sum }_{m=0}^{\infty }{C}_m{\eta }^m,$$

(23)

where the unknown coefficients $C_m$ are determined recurrently. If a linear operator $G\equiv {\mathrm{d}}^4/\mathrm{d}{\eta }^4$ is used, then the inverse operator of $G$ is a four-fold integral operator defined as

$$G^{-1}=\underset{0}{\overset{\eta }{{\displaystyle \int }}}\underset{0}{\overset{\eta }{{\displaystyle \int }}}\underset{0}{\overset{\eta }{{\displaystyle \int }}}\underset{0}{\overset{\eta }{{\displaystyle \int }}}(\cdots )\mathrm{d}\eta \mathrm{d}\eta \mathrm{d}\eta \mathrm{d}\eta .$$

(24)

Equation (18) can now be written as

$$\phi \left(\eta \right)=\mathsf{\Phi }\left(\eta \right)+G^{-1}\{\frac{8\beta }{\left(1-\beta \eta \right)}\frac{{\mathrm{d}}^3\phi \left(\eta \right)}{\mathrm{d}{\eta }^3}-\left[\frac{12{\beta }^2}{{\left(1-\beta \eta \right)}^2}+{\lambda }_2\frac{\left(3\eta -3\beta {\eta }^2+{\beta }^2{\eta }^3\right)}{{\left(1-\beta \eta \right)}^4}\right]\frac{{\mathrm{d}}^2\phi \left(\eta \right)}{\mathrm{d}{\eta }^2}+\frac{{\lambda }_1}{{\left(1-\beta \eta \right)}^2}\phi \left(\eta \right)\}.$$

(25)

By substituting the decomposition $\phi \left(\eta \right)={{\displaystyle \sum }}_{m=0}^{\infty }{C}_m{\eta }^m$ into Equation (25), the following is obtained:

$$\phi \left(\eta \right)=\mathsf{\Phi }\left(\eta \right)+G^{-1}\{\frac{8\beta }{\left(1-\beta \eta \right)}{\displaystyle \sum }_{m=0}^{\infty }\left(m+1\right)\left(m+2\right)\left(m+3\right)C_{m+3}{\eta }^m-\left[\frac{12{\beta }^2}{{\left(1-\beta \eta \right)}^2}+{\lambda }_2\frac{\left(3\eta -3\beta {\eta }^2+{\beta }^2{\eta }^3\right)}{{\left(1-\beta \eta \right)}^4}\right]{\displaystyle \sum }_{m=0}^{\infty }\left(m+1\right)\left(m+2\right)C_{m+2}{\eta }^m+\frac{{\lambda }_1}{{\left(1-\beta \eta \right)}^2}{\displaystyle \sum }_{m=0}^{\infty }{C}_m{\eta }^m\},$$

(26)

where $\mathsf{\Phi }\left(\eta \right)={{\displaystyle \sum }}_{m=0}^3{C}_m{\eta }^m=\phi \left(0\right)+{\phi }^{\prime }\left(0\right)\eta +{\phi }^{″}\left(0\right){\eta }^2/2+{\phi }^{‴}\left(0\right)/6$ is the initial term.

The following are a list of power series:

$$\begin{array}{l}{\left(1-\beta \eta \right)}^{-1}={\displaystyle {\displaystyle \sum }_{j=0}^{\infty }}{\left(\beta \eta \right)}^j, \beta \ne 0,\\ {\left(1-\beta \eta \right)}^{-2}={\displaystyle {\displaystyle \sum }_{j=0}^{\infty }}\left(j+1\right){\left(\beta \eta \right)}^j, \beta \ne 0,\\ {\left(1-\beta \eta \right)}^{-4}={\displaystyle {\displaystyle \sum }_{j=0}^{\infty }}\frac{\left(j+1\right)\left(j+2\right)\left(j+3\right)}{6}{\left(\beta \eta \right)}^j, \beta \ne 0.\end{array}$$

(27)

Using these relationships, the following are obtained:

$${\left(1-\beta \eta \right)}^{-1}{\displaystyle \sum }_{m=0}^{\infty }\left(m+1\right)\left(m+2\right)\left(m+3\right)C_{m+3}{\eta }^m={\displaystyle \sum }_{m=0}^{\infty }{\eta }^m{\displaystyle \sum }_{j=0}^m{\beta }^{m-j}\left(j+1\right)\left(j+2\right)\left(j+3\right)C_{j+3}.$$

(28)

$${\left(1-\beta \eta \right)}^{-2}{\displaystyle \sum }_{m=0}^{\infty }\left(m+1\right)\left(m+2\right)C_{m+2}{\eta }^m={\displaystyle \sum }_{m=0}^{\infty }{\eta }^m{\displaystyle \sum }_{j=0}^m{\beta }^{m-j}\left(m-j+1\right)\left(j+1\right)\left(j+2\right)C_{j+2}.$$

(29)

$${\left(1-\beta \eta \right)}^{-4}{\displaystyle \sum }_{m=0}^{\infty }\left(m+1\right)\left(m+2\right)C_{m+2}{\eta }^m={\displaystyle \sum }_{m=0}^{\infty }{\displaystyle \sum }_{j=0}^m{\beta }^{m-j}\frac{\left(m-j+1\right)\left(m-j+2\right)\left(m-j+3\right)}{6}\left(j+1\right)\left(j+2\right)C_{j+2}.$$

(30)

On integrating Equation (26) and incorporating Equations (28)–(30)

$$\begin{array}{l}{\displaystyle {\displaystyle \sum }_{m=0}^{\infty }}{C}_m{\eta }^m=\phi \left(0\right)+{\phi }^{\prime }\left(0\right)\eta +\frac{\phi ″\left(0\right)}{2}{\eta }^2+\frac{\phi ″\prime \left(0\right)}{6}{\eta }^3\\ +{\displaystyle {\displaystyle \sum }_{m=0}^{\infty }}\{\frac{{\eta }^{m+4}}{\left(m+1\right)\left(m+2\right)\left(m+3\right)\left(m+4\right)}{\displaystyle {\displaystyle \sum }_{j=0}^m}8{\beta }^{m-j+1}\left(j+1\right)\left(j+2\right)(j\\ +3)C_{j+3}-\left(j+1\right)\left(j+2\right)[12{\beta }^{m-j+2}(m-j\\ +1)+{\lambda }_2(3\eta -3\beta {\eta }^2\\ +{\beta }^2{\eta }^3\left){\beta }^{m-j}\frac{\left(m-j+1\right)\left(m-j+2\right)\left(m-j+3\right)}{6}\right]C_{j+2}+{\lambda }_1{\beta }^{m-j}(m-j\\ +1\left)C_j\right\}.\end{array}$$

(31)

The recurrence relation for the coefficients $C_k$ can now be stated as

$$C_0=\phi \left(0\right), C_1={\phi }^{\prime }\left(0\right), C_2=\frac{{\phi }^{″}\left(0\right)}{2}, C_3=\frac{{\phi }^{‴}\left(0\right)}{6},$$

(32)

and for $m\ge 4$ as

$$\begin{array}{l}{C}_m\\ =\frac{1}{m\left(m-1\right)\left(m-2\right)\left(m-3\right)}{\displaystyle {\displaystyle \sum }_{j=0}^{m-4}}\left(j+1\right)\left(j+2\right)\left(j+3\right)8{\beta }^{m-j-3}{C}_{j+3}-(j\\ +1\left)\right(j\\ +2\left)\right[{\beta }^{m-j-2}(m-j\\ -3)+{\lambda }_2(3\eta -3\beta {\eta }^2\\ +{\beta }^2{\eta }^3\left){\beta }^{m-j-4}\frac{\left(m-j-3\right)\left(m-j-2\right)\left(m-j-1\right)}{6}\right]C_{j+2}+{\lambda }_1{\beta }^{m-j-4}(m\\ -j-3\left)C_j\right\}.\end{array}$$

(33)

The coefficients $C_m$ can be found from the recurrent Equations (32) and (33) and the solution for $\phi \left(\eta \right)$ subsequently calculated from Equation (26). The series solution is $\phi \left(\eta \right)={{\displaystyle \sum }}_{m=0}^{\infty }{C}_m{\eta }^m$, although all of the coefficients $C_m$ cannot be determined, and thus the solutions must be approximated by the truncated series ${{\displaystyle \sum }}_{m=0}^{n-1}{C}_m{\eta }^m.$ The successive approximations are ${\phi }^{\left[n\right]}\left(\eta \right)={{\displaystyle \sum }}_{m=0}^{n-1}{C}_m{\eta }^m$, as $n$ increases and the boundary conditions are met.

Thus, ${\phi }^{\left[1\right]}\left(\eta \right)=C_0, {\phi }^{\left[2\right]}\left(\eta \right)={\phi }^{\left[1\right]}\left(\eta \right)+C_1\eta , \mathrm{and} {\phi }^{\left[3\right]}\left(\eta \right)={\phi }^{\left[2\right]}+C_2{\eta }^2$ serve as approximate solutions with increasing accuracy as $n\to \infty$. The four coefficients $C_m\left(m=0, 1, 2, 3\right)$ in Equation (32) depend on the boundary conditions of Equations (21) and (22). As in this case, the boundary conditions are free-free, so $K_{RR}\to 0,K_{TR}\to 0,K_{RL}\to 0,K_{TL}\to 0$, and so it is most convenient to choose the two coefficients $C_0$ and $C_1$ as zero, and the other two coefficients $C_2$ and $C_3$ as arbitrary constants with the coefficients $C_m$ obtained using Equation (33). By substituting ${\phi }^{\left[n\right]}\left(\eta \right)={{\displaystyle \sum }}_{m=0}^{n-1}{C}_m{\eta }^m$ into the last two boundary conditions of Equation (22) two algebraic equations involving $C_2$ and $C_3$ are obtained

$${\displaystyle \sum }_{m=0}^{n-3}\left(m+1\right)\left(m+2\right)C_{m+2}=f_{12}^{\left[n\right]}\left({\lambda }_1\right)C_2+f_{13}^{\left[n\right]}\left({\lambda }_1\right)C_3=0,$$

(34)

$${\displaystyle \sum }_{m=0}^{n-4}\left(m+1\right)\left(m+2\right)\left(m+3\right)C_{m+3}=f_{22}^{\left[n\right]}\left({\lambda }_1\right)C_2+f_{23}^{\left[n\right]}\left({\lambda }_1\right)C_3=0,$$

(35)

and for non-trivial solutions of $C_2$ and $C_3$ the frequency equation can be written as

$$f_{12}^{\left[n\right]}\left({\lambda }_1\right)f_{23}^{\left[n\right]}\left({\lambda }_1\right)-f_{22}^{\left[n\right]}\left({\lambda }_1\right)f_{13}^{\left[n\right]}\left({\lambda }_1\right)=0.$$

(36)

The ith estimated eigenvalue ${\lambda }_{1\left(i\right)}^{\left[n\right]}$ corresponding to $n$ is obtained from Equation (36), i.e., the ith estimated dimensionless natural frequency ${\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}=\sqrt{{\lambda }_{1\left(i\right)}^{\left[n\right]}}$ is also obtained and $n$ is determined by the following equation:

$$\left|{\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}-{\mathsf{\Omega }}_{n\left(i\right)}^{\left[n-1\right]}\right|\le \epsilon ,$$

(37)

where ${\mathsf{\Omega }}_{n\left(i\right)}^{\left[n-1\right]}$ is the ith estimated dimensionless natural frequency corresponding to $n-1$ and $\epsilon$ is a preset sufficiently small value. If Equation (37) is satisfied, then ${\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}$ is the ith dimensionless natural frequency ${\mathsf{\Omega }}_{n\left(i\right)}$. By substituting ${\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}$ into Equation (36), we obtain

$$C_3=-\frac{f_{r2}^{\left[n\right]}\left({\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}\right)}{f_{r3}^{\left[n\right]}\left({\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}\right)}{C}_2, r=1,2,$$

(38)

and all of the other coefficients $C_m$ can be obtained from Equations (32) and (33). In addition, the ith mode shape ${\phi }_i^{\left[n\right]}$ corresponding to the ith eigenvalue ${\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}$ is obtained by

$${\phi }_i^{\left[n\right]}\left(\eta \right)={\displaystyle \sum }_{m=0}^{n-1}{C}_m^{\left[i\right]}{\eta }^m,$$

(39)

where $C_m^{\left[i\right]}$ is $C_m\left(\eta \right)$ in which ${\lambda }_1$ is substitutes by ${\lambda }_{1\left(i\right)}$ and ${\phi }_i^{\left[n\right]}$ is the ith eigenfunction corresponding to the ith eigenvalue ${\lambda }_{1\left(i\right)}$. By normalizing Equation (39), the ith normalized eigenfunction is defined as

$${\overline{\phi }}_i^{\left[n\right]}\left(\eta \right)=\frac{{\phi }_i^{\left[n\right]}\left(\eta \right)}{\sqrt{{{\displaystyle \int }}_0^1{\left[{\phi }_i^{\left[n\right]}\right]}^2\mathrm{d}\eta }},$$

(40)

where ${\overline{\phi }}_i^{\left[n\right]}\left(\eta \right)$ is the ith mode shape function of the beam corresponding to the ith natural frequency ${\omega }_i^{\left[n\right]}=\sqrt{{\lambda }_{1\left(i\right)}^{\left[n\right]}}\sqrt{\frac{E{I}_0}{\rho A_0{l}^4}}={\mathsf{\Omega }}_{n\left(i\right)}^{\left[n\right]}\sqrt{\frac{E{I}_0}{\rho A_0{l}^4}}.$

5. Numerical Results

The convergence plots for the first three dimensionless natural frequencies, with the truncation factor, $\beta$ set at 0.25 and thrust set at zero, are shown in Figure 3. By solving Equation (36) for $n$ and taking the real root for ${\lambda }_{1\left(1\right)}={\mathsf{\Omega }}_{n\left(1\right)}^2$, it can be shown that for $n=42$

$$\left|{\mathsf{\Omega }}_{n\left(1\right)}^{\left[42\right]}-{\mathsf{\Omega }}_{n\left(1\right)}^{\left[41\right]}\le \epsilon =0.00001\right|,$$

(41)

with

$${\mathsf{\Omega }}_{n\left(1\right)}\approx {\mathsf{\Omega }}_{n\left(1\right)}^{\left[42\right]}=15.1278.$$

(42)

For this particular case, it was found that convergence occurred for ${\mathsf{\Omega }}_{n\left(1\right)}$ after 42 iterations, for ${\mathsf{\Omega }}_{n\left(2\right)}$ after 52 iterations, and for ${\mathsf{\Omega }}_{n\left(3\right)}$ after 60 iterations. It is considered that the convergence is quite rapid for each of the natural frequencies.

The development method of solution was in part validated with published results. In

Table 1. Frequency parameters of a uniform beam with circular cross-section with free-free boundary conditions.

Uniform Cylinder, $\mathit{\beta }=0.0$
	${\mathbf{\Omega }}_{\mathit{n}\left(1\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(2\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(3\right)}$
Current	22.3745	61.6912	121.0452
Ref. [26]	22.3733	61.6728	120.9030

, a comparison is made with the results of Krynicki and Mazurkiewicz [26] for a solid cylindrical beam with free-free boundary conditions and with the thrust set at zero. It can be found that there is reasonable agreement between the two sets of results, although, generally, the current method gives slightly higher solutions for all three natural frequencies.

Calculations were also made and compared with published data for a “complete” cone, which means the cone was not truncated but sharp when $\eta =1$. The boundary conditions were again free-free and the thrust was set to zero. As can be seen from

Table 2. Frequency parameters of a “complete” rocket as compared with the frequency parameters of a free-free complete cone.

‘Complete’ Rocket, $\mathit{\beta }=1.0$
	${\mathbf{\Omega }}_{\mathit{n}\left(1\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(2\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(3\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(4\right)}$
Current	17.4195	35.0436	57.4867	84.8549
Ref. [27]	17.4164	35.0270	57.3796	84.6028

, there was fairly good agreement with those reported by Naguleswaran [27], although, again, the current method gave higher results for each of the natural frequencies.

By substituting ${\mathsf{\Omega }}_{n\left(1\right)}^{\left[42\right]}$ as given by Equation (42) into Equations (33) and (39) (with ${\lambda }_2=0)$ and normalizing using Equation (40), a polynomial can be obtained to describe the first mode shape function. The same procedure was employed for other natural frequencies and mode shapes. The mode shapes related to the first four natural frequencies are shown on Figure 4.

The variation of the natural frequencies with the truncation factor, again without thrust, is given in

Table 3. First three dimensionless natural frequencies of a free-free truncated cone for different truncation factors, $\beta$ with free-free boundary conditions.

Truncated Cone
$\mathit{\beta }$		${\mathbf{\Omega }}_{\mathit{n}\left(1\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(2\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(3\right)}$
0.05	Current	15.9551	33.1259	56.5495
	Ref. [27]	15.9543	33.1251	56.5484
0.15	Current	14.9399	34.8542	63.8621
	Ref. [27]	14.9381	34.8535	63.8434
0.25	Current	15.1278	37.8976	71.5367
	Ref. [27]	15.1266	37.8972	71.5171
0.35	Current	15.7681	41.1205	78.8673
	Ref. [27]	15.7674	41.1193	78.8259
0.45	Current	16.6135	44.3543	86.0330
	Ref. [27]	16.6130	44.3537	85.8162
0.55	Current	17.5662	47.5667	90.2381
	Ref. [27]	17.5653	47.5650	-
0.65	Current	18.5803	50.6751	93.7655
	Ref. [27]	18.5794	-	-

. It can be seen that the values are in reasonable agreement with those reported by Naguleswaran [27]. Additionally, it is notable that the first natural frequency does not vary monotonically with $\beta$, whereas the second and third natural frequencies do.

As may be intuitively obvious, the increase in thrust force will reduce the natural frequencies. The calculations made for the first three modes for increasing nondimensional thrust values and for different truncation factors, $\beta$, are shown in

Table 4. Effect of thrust magnitude on the natural frequency of the first three modes.

Truncated Cone
$\mathit{\beta }$		${\mathbf{\Omega }}_{\mathit{n}\left(1\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(2\right)}$	${\mathbf{\Omega }}_{\mathit{n}\left(3\right)}$
0.05	${\lambda }_2=0$	15.9551	33.1259	56.5495
	${\lambda }_2=1$	11.9663	22.7243	35.3434
	${\lambda }_2=2$	7.6451	14.2773	22.9732
	${\lambda }_2=3$	5.9561	10.2790	16.7838
0.25	${\lambda }_2=0$	15.1278	37.8976	71.5367
	${\lambda }_2=1$	11.0306	25.7700	44.7104
	${\lambda }_2=2$	6.9336	16.3336	29.0618
	${\lambda }_2=3$	5.9881	11.7861	21.2321
0.65	${\lambda }_2=0$	18.5803	50.6751	93.7655
	${\lambda }_2=1$	13.1610	34.9608	60.9476
	${\lambda }_2=2$	8.1289	22.2970	39.3815
	${\lambda }_2=3$	6.5804	16.2160	28.1300

. As can be seen, the increases in thrust have a considerable effect on the reduction of all three natural frequencies for each of the $\beta$ values. For example, for the first natural frequency when $\beta =0.05$, there is a large reduction in natural frequency of approximately 63%.

To find the effect of the increase in thrust on mode shape, calculations were made with increasing nondimensional thrust values, ${\lambda }_2$ for the fundamental mode shape. The increase in thrust is seen in Figure 5 to reduce the magnitude of the mode shape in the vicinity of the aft of the rocket, i.e., close to the application of the thrust, whilst close to the rocket nose the magnitude of the mode shape increases. This phenomenon, according to Pourtakdoust and Assadian [1] motivates inertial measuring units (IMU) induced vibrations on considering the forward location of the IMU. In Pourtakdoust and Assadian [28] it was found that the IMU location in the aft of the vehicle could lead to dynamic instability. Increasing thrust caused a backward shift in the extremum point of the fundamental mode shape curve which allowed for a larger margin of safety for the IMU induced vibration.

Finally, to show the ease of changing boundary conditions, and for further validation, a comparative study with calculations obtained by a finite-element method (FEM) for various boundary conditions was undertaken. The analytical results were compared with the results obtained using the commercial FEM software ANSYS [29] and previously reported by Avcar [30].

In this study, the beam had a constant square cross-section and homogeneous material properties. The analytical calculations were converted to dimensional quantities for the purpose of comparison, with $E=70\times {10}^9$ N/m², $\rho =2700$ kg/m³, $l=3$ m, and $A=0.03$ m². The boundary conditions used were, clamped-clamped and clamped-free, respectively. For the clamped-clamped case, the nondimensional coefficients for translation and rotation, as found in Equations (21) and (22), become $K_{RL}\to \infty , K_{RR}\to \infty , K_{TL}\to \infty , K_{TR}\to \infty$. For the clamped-free case, $K_{RL}\to \infty , K_{RR}\to 0, K_{TL}\to \infty , K_{TR}\to 0.$ The comparison between the current method and the finite-element method is shown in Figure 6, where $f_n$ and $n$ are the natural frequency and mode number, respectively.

As can be seen, the two methods of calculation give close results, although the AMDM calculations tend to be consistently slightly higher. The maximum discrepancy found between the two calculation methods was 0.132%.

6. Conclusions

The proposed method has been effective in obtaining the closed-form of the free vibration of a slender rocket with free-free boundary conditions. The ith natural frequency and mode shape function can be derived using this method, with the results of this relatively simple approach obtaining results similar to more complex methods. Two specific advantages of using the AMDM approach are the fast convergence rates noted and the high accuracy of the calculations. A practical advantage of using AMDM is that boundary conditions can be easily changed for different scenarios.