Nonlinear analysis of dynamic stability for the thin cylindrical shells of supercavitating vehicles

Abstract

The dynamic stability of supercavitating vehicles under periodic axial loading is investigated in this article. The supercavitating vehicle is simulated as a long and thin cylindrical shell subjected to periodic axial loading and simply supported boundary conditions. The nonlinear transverse vibration differential equation is obtained in terms of nonlinear geometric equations, physical equations, and balance equations of cylindrical shells. Mathieu equation with periodic coefficients and nonlinear term is derived by employing Galerkin variational method and Bolotin method. The analytical expressions of the steady-state amplitudes of vibrations in the first- and second-order instable regions are obtained by solving nonlinear Mathieu equation derived in this article. Numerical results are presented to analyze the influence of the sailing speed, ratio of loads, the frequency of axial loads, and the mode of vibration on parametric resonance curves and to show the nonlinear parametric resonance curves incline toward the side where it is greater than the excitation frequency, which significantly extends the range of the exciting region. The presented results indicate the enlargement of the exciting region will cause shrinkage of the safe frequency range of external loads and decrease in dynamic stability of supercavitating vehicle.

Introduction

An underwater supercavitating vehicle such as the supercavitating torpedo is an ultra-high-speed underwater vehicle that rides in a cavity of water vapor initialized at the forward mounted cavitator. The liquid surrounding the body gasifies and forms a supercavity that contains the body, which results in a dramatic drag reduction. Because of this cavity, the supercavitating vehicle is only marginally in contact with the water at the cavitator and rear control fins. This minimal contact results in stability issues that are not present in conventional, fully wetted vehicles. Therefore, the structure of the vehicles must be able to withstand the tremendous propulsive forces associated with propelling a body at such high-speed underwater. The supercavitating vehicle generally has very high slenderness ratio (about 10–22), of which the dynamic buckling deformation not only influences the safety of the structure itself but also influences stability of supercavity and motion. For the long and thin supercavitating vehicle, only the cavitator, the fins, and the tail of the body have contact with the liquid. Therefore, the dynamic characteristics and the stability of the body are linked to the sailing speed and the cavitation state. The flow field characteristic of supercavity is constantly changing with the motion history and ballistic trajectory of the vehicle. ¹ When the axial forces applied to vehicle change periodically over time, as long as in certain circumstances, namely, in parameter plane meeting certain relationship between the transverse oscillation frequency, the excitation coefficient, and the frequency of external loads, parametric resonance will be excited, resulting in the amplitude of vibrations rapidly increasing to dynamic buckling of supercavitating vehicle. Hence, for long and thin supercavitating vehicle, parametric resonance is a kind of unexpected, dangerous, and secluded failure mode deserving attention.

In the early systematic studies on the dynamical stability problems of cylindrical shells, Bolotin ² pointed out that dynamic instability of cylindrical shells occurs when certain relationships between their natural frequency and the frequency of the axial load exist. Since then, dynamic stability has become an important issue for the safety design of cylindrical shells. The problem of dynamic instability of elastic structural elements induced by parametric excitation has been addressed by many researchers. Early work on this subject was reported by Evan-Iwanowski. ³ In relation to shell, Yao ⁴ examined the nonlinear elastic buckling and parametric excitation of a cylinder under axial loads. Vijayaraghavan and Evan-Iwanowski ⁵ discussed the parametric instability of circular cylindrical shells. Based on Donnell’s shell equations, Yamaki and Nagai^6,7 used Hsu’s method to examine the dynamic stability of circular cylindrical shells under both static and periodic compressive forces. Radwan and Genin ⁸ studied the stability of the steady-state response of circular cylinders subjected to harmonic excitation for the simply supported boundary condition. Using Bolotin’s method, Lam and colleagues^9–11 applied four different thin shell theories to establish differential equation of cylindrical shell under axial force and transformed it into Mathieu equation to explore the dynamic stability of cylindrical shells and analyze the parametric response. Ganapathi an colleagues^12,13 investigated the parametric instability behavior of a laminated anisotropic composite circular cylindrical shell subjected to periodic axial loading. Using a combination of the Ritz method and Bolotin’s first approximation, Liew et al.^14,15 investigated the dynamic stability of rotating cylindrical shells under static and periodic axial forces, obtained a system of Mathieu–Hill equations by application of the Ritz minimization procedure to the energy expressions. Bochkarev et al. ¹⁶ investigated the dynamic behavior of elastic coaxial cylindrical shells interacting with two flows of a perfect compressible fluid by application of the finite element method. Zhu et al. ¹⁷ presented a set of Mathieu–Hill equations for the dynamic stability analysis of piezoelectric circular cylindrical shells subjected to combined periodic axial compression and radial electric field loading. Nemoto et al. ¹⁸ used Galerkin’s method to obtain the solutions for the prebuckling motion and the perturbated motion of composite laminated cylindrical shells subjected to periodic hydrostatic pressure. Bakhtieva and Tazyukov ¹⁹ presented a new approach to building a mathematical model based on the Ostrogradsky-Hamilton principle of stationary action to solve the stability problem for a circular cylindrical shell under an axial impulsive load. Thompson ²⁰ described the static-dynamic analogy and its role in understanding the localized post-buckling of shell-like structures and showed the true significance of the Maxwell energy criterion load in predicting the sudden onset of “shock sensitivity” to lateral disturbances. Kumar et al. ²¹ used the Flügge–Luré–Byrne theory to study parametric instability of a cylinder subject to a uniform radially fluctuating pressure. Bakhtieva and Tazyukov ²² presented a modification of the algorithm to solve the task of stability for a circular cylindrical shell under axial compression.

Several studies on nonlinear models of shells for dynamic stability problems are available in the literature. Evensen ²³ presented mode expansions of shell flexural displacement involving the driven and companion modes plus an axisymmetric term. Bogdanovich ²⁴ discussed nonlinear problems of dynamic buckling for reinforced laminar cylindrical shells. Zhou and Wang ²⁵ derived the general governing equations of nonlinear dynamic stability for laminated cylindrical shells in which, factors of nonlinear large deflection, transverse shear, and longitudinal inertia force were concluded. Amabili ²⁶ investigated large-amplitude vibrations of circular cylindrical shells subjected to radial harmonic excitation in the spectral neighborhood of the lowest resonances. Using Donnell’s shallow-shell equations, Gonçalves and Del Prado ²⁷ investigated the nonlinear vibrations and instabilities of cylindrical shells under pulsating axial loads using Donnell’s shallow-shell equations. Ganapathi ²⁸ studied the dynamic stability behavior of a clamped functionally graded materials spherical shell structural element subjected to external pressure load. Darabi et al. ²⁹ presented a simple solution of the dynamic stability of functionally graded shells under periodic axial loading based on large deflection theory, obtaining the steady-state vibrations for nonlinear Mathieu equations. Paak et al. ³⁰ investigated the effect of steady viscous forces on the nonlinear behavior and stability of cantilevered shells conveying fluid. T Dey and Ramachandra ³¹ discussed the parametric vibration of a simply supported composite circular cylindrical shell under periodic partial edge loadings, obtained the characteristic features of the dynamic instability regions, the effect of dynamic load amplitude on the nonlinear response is studied. Pradyumna and Gupta ³² employed Galerkin’s method to investigate the effects of different parameters on the dynamic stability characteristic of laminated composite shells integrated with piezoelectric layers in thermal environment. Panda and Ramachandra ³³ used Hamilton’s variational principle based on higher order shear deformation theory of elastic shell theory incorporating von Kármán-type nonlinear strain displacement relations to discuss boundary conditions on the dynamic instability regions of shear deformable cross-ply laminated and composite cylindrical panels subjected to periodic nonuniform in-plane loads. Pellicano and Barbieri ³⁴ used the nonlinear Sanders-Koiter theory to analyze the non-stationary response of the circular cylindrical shells subjected to inertial axial loads. Dey and Ramachandra ³⁵ used Galerkin’s method and Incremental harmonic balance method to investigate linear and nonlinear dynamic instability behavior of cylindrical sandwich panels subjected to combined static and dynamic nonuniform in-plane loadings. M Darabi and Ganesan ³⁶ used Donnell’s shallow-shell theory and with von Karman-type of nonlinearity to investigate the dynamic instability of thin laminated composite cylindrical shells subjected to harmonic axial loading based on nonlinear analysis.

At present, few studies on parametric resonance of supercavitating vehicle have been conducted by using nonlinear models of shells for dynamic stability problems. Ruzzene, Ahn and Choi et al. from Georgia Institute of Technology have conducted researches on parametric resonance of supercavitating vehicle. On the basis of Hamilton theory, Ruzzene^37,38 established motion equation of supercavitating vehicle, to solve dynamic buckling of the structure by means of finite element discretization and further obtain instability regions of supercavitating vehicle by adopting Bolotin method. By adopting finite element method, Ahn and Ruzzene ³⁹ investigated the static and dynamic stability of supercavitating vehicles. The vehicles are modeled as thick shells, the stability performance of plain shells is compared to those of tapered and stiffened shells. Choi and Ruzzene ⁴⁰ introduced the concept of adaptive cavitator, established finite element model of supercavitating vehicle by using the beam element with shear deformation, and conducted dynamic buckling analysis of the structure to predict stability limits.

From the point of existing literature, few studies on nonlinear dynamic instability of cylindrical shell of the high-speed supercavitating vehicle based on parametric resonance have been conducted. The larger transverse deformation of motion body will directly destroy the original stability of cavitation, and the stability of cavitation determines stability of motion of supercavitating vehicle. Linear parametric resonance theory is only applicable for small deformation. With the increase of transverse amplitudes of vibrations, the influence of nonlinear factors begins to appear and becomes an indispensable factor. Hence, the instability study of the supercavitating vehicle based on parametric vibration in this article will be developed from the point of nonlinearity.

In this study, the nonlinear large deflection and curvatures are defined by Donnell’s theory, nonlinear differential equations of transverse vibration of cylindrical shell are established by adding a nonlinear term to linear strain and shearing strain. By adopting Galerkin variation method and Bolotin method, the nonlinear differential equation is transformed into Mathieu equation containing periodic coefficients and the nonlinear term. Since the steady-state amplitudes of vibrations of Mathieu equation in the first, second, and third order unstable region decrease with the increase in order rapidly, only the analytical expressions of the steady-state amplitudes of vibrations in the first- and second-order unstable regions are given by solving nonlinear Mathieu equations, so as to analyze factors affecting nonlinear parametric resonance curves of cylindrical shell. Numerical results are presented to analyze the influence of the sailing speed, ratio of loads, the frequency of axial loads, and the mode of vibration on nonlinear parametric resonance curves so as to evaluate the widths of exciting region and the amplitude of steady-state vibrations in the exciting region.

A simplified model for a supercavitating vehicle

A simplified model of the supercavitating vehicle’s configuration and of the applied forces is given in Figure 1. The body is acted upon by a system of forces corresponding to the interaction of the vehicle’s control surfaces with the cavity boundaries. The control surfaces include the immersion part at the back of the vehicle and the cavitator, whose primary function is the generation of the supercavity. The control surfaces support the vehicle in the vertical direction by providing lift and allow maneuvering according to assigned flight paths. Finally, the vehicle’s motion is sustained by a propulsion force.

Figure 1.

Configuration of the supercavitating vehicle.

The vehicle’s interactions with the fluid can be simplified as a combination of a distributed axial force, corresponding to the drag at the nose and to the propulsion at the tail, and a concentrated load representing the drag applied by the afterbody as it surfs the cavity. The drag force applied at the nose of the vehicle during its forward motion is given by ⁴¹

$$F_D=\frac{1}{2}{\rho }_w{A}_c{C}_x{V}^2$$

(1)

where ${\rho }_w$ is the density of the fluid (water), $A_c$ is the cross-sectional area of the cavitator, $C_x$ is the drag coefficient of the cavitator, and V is the velocity of forward motion of the vehicle. The drag coefficient for a cavitator with zero angle of attack can be expressed as

$$C_x=C_{x0}(1+\sigma )$$

(2)

$$C_{x0}=0.5+1.81\left(\frac{\phi }{360-0.25}\right)-2{\left(\frac{\phi }{360-0.25}\right)}^2$$

(3)

$$\sigma =\frac{2(p_{\infty }-p_c)}{{\rho }_w{V}^2}$$

(4)

where $\phi$ is the cone angle of the cavitator, and the cone angle of a disk (cavitator) is $180\circ$ . $\sigma$ is the cavitation number. $p_{\infty }$ and $p_c$ are the ambient fluid pressure and the cavity pressure, respectively

$$p_{\infty }={\rho }_w\mathrm{gH}+p$$

(5)

where H is the sailing depth, and p is standard atmospheric pressure.

The nature of supercavitating flows is extremely unstable. Cavity shape and pressure oscillations are constantly changing with the motion history and ballistic trajectory of the supercavitating vehicle, so that the cavitator drags and planing forces at the back are constantly changing over time and occur oscillation easily. Since the behavior of axial forces on supercavitating vehicles is still being evaluated and fully characterized, and due to the complexity of the considered system, the applied axial loads are simplified as the sum of a constant component and of a time varying term, which accounts for the oscillations of the force with respect to its constant, steady flow value in this article, as shown in Figure 2.

Figure 2.

The simplified configuration of supercavitating vehicle.

The simplified axial force is expressed as ³⁷

$$p\left(t\right)=P_0+P_t=p_0(1+\delta \cos \theta t)$$

(6)

$$P_0=F_D=\frac{\pi {\rho }_w{C}_x{d}_n^2{V}^2}{8}$$

(7)

where δ is the ratio of loads, $\delta =P_t/P_0$ , and $d_n$ is the diameter of cavitator. $\theta =2\pi /T$ defines the frequency of oscillations of the axial force.

Nonlinear dynamic differential equation of cylindrical shell

Simulate supercavitating vehicle as a long and thin cylindrical shell under periodic axial loading. As shown in Figure 3, consider a cylindrical shell with radius R and uniform thickness h. Curve coordinates $(\alpha ,\beta ,z)$ in the mid-surface of shell before deformation is used to describe the geometric position of finite deformation, where $\alpha ,\beta$ are the main lines of curvature in middle-plane of shell.

Figure 3.

The geometry of a thin cylindrical shell.

By introducing the influence of tangential displacement $(u,v)$ into the geometrical relationship of torsional strain $({\kappa }_{12})$ , and by introducing the nonlinear items into the linear strains $({\epsilon }_1,{\epsilon }_2)$ and the shearing strain $({\gamma }_{12})$ , the nonlinear large deflection and curvatures are defined by Donnell’s ⁴² theory as follows

$$\{\begin{array}{c}\multicolumn{1}{c}{{\epsilon }_1=\frac{1}{R}\frac{\partial u}{\partial \alpha }+\frac{1}{2R^2}{\left(\frac{\partial w}{\partial \alpha }\right)}^2}\\ \multicolumn{1}{c}{{\epsilon }_2=\frac{1}{R}(\frac{\partial v}{\partial \beta }+w)+\frac{1}{2R^2}{\left(\frac{\partial w}{\partial \beta }\right)}^2}\\ \multicolumn{1}{c}{{\gamma }_{12}=\frac{1}{R}(\frac{\partial v}{\partial \alpha }+\frac{\partial u}{\partial \beta })+\frac{1}{R^2}\frac{\partial w}{\partial \alpha }·\frac{\partial w}{\partial \beta }}\\ \multicolumn{1}{c}{{\kappa }_1=-\frac{1}{R^2}\frac{{\partial }^{2}w}{\partial {\alpha }^2}}\\ \multicolumn{1}{c}{{\kappa }_2=-\frac{1}{R^2}(\frac{{\partial }^{2}w}{\partial {\beta }^2}+w)}\\ \multicolumn{1}{c}{{\kappa }_{12}=-\frac{1}{2R^2}(\frac{\partial u}{\partial \beta }-\frac{\partial v}{\partial \alpha }+2\frac{{\partial }^{2}w}{\partial \alpha \partial \beta })}\end{array}$$

(8)

where $u,v,w$ denote the axial, circular and radial displacement, respectively. R is the radius of cylindrical shell, ${\epsilon }_1,{\epsilon }_2$ are linear strains in the middle-plane, ${\gamma }_{12}$ is the shear strain in the middle-plane, ${\kappa }_1,{\kappa }_2$ are the bending strain, and ${\kappa }_{12}$ is the torsional strain.

In the physical equation, by introducing the influence of the bending strains $({\kappa }_1,{\kappa }_2,{\kappa }_{12})$ in the middle-plane membrane forces $(N_1,N_2,N_{12},N_{21})$ , and by introducing the influence of the middle-plane strains $({\epsilon }_1,{\epsilon }_2)$ in the middle-plane bending moment $M_1$ , we can obtain the following equation

$$\{\begin{array}{c}\multicolumn{1}{c}{N_1=\frac{\mathrm{Eh}}{1-{\nu }^2}({\epsilon }_1+\nu {\epsilon }_2+\frac{h^2}{12R}{\kappa }_1)}\\ \multicolumn{1}{c}{N_2=\frac{\mathrm{Eh}}{1-{\nu }^2}({\epsilon }_2+\nu {\epsilon }_1-\frac{h^2}{12R}{\kappa }_2)}\\ \multicolumn{1}{c}{N_{12}=\frac{\mathrm{Eh}}{2(1+\nu )}({\gamma }_{12}+\frac{h^2}{12R}{\kappa }_{12})}\\ \multicolumn{1}{c}{N_{21}=\frac{\mathrm{Eh}}{2(1+\nu )}({\gamma }_{12}-\frac{h^2}{12R}{\kappa }_{12})}\\ \multicolumn{1}{c}{M_1=-\frac{E{h}^3}{12(1-{\nu }^2)}[{\kappa }_1+\nu {\kappa }_2+\frac{1}{R}({\epsilon }_1+\nu {\epsilon }_2)]}\\ \multicolumn{1}{c}{M_2=-\frac{E{h}^3}{12(1-{\nu }^2)}({\kappa }_2+\nu {\kappa }_1)}\\ \multicolumn{1}{c}{M_{12}=\frac{E{h}^3}{24(1+\nu )}(2{\kappa }_{12}+\frac{1}{R}{\gamma }_{12})}\\ \multicolumn{1}{c}{M_{21}=\frac{E{h}^3}{12(1+\nu )}{\kappa }_{12}}\end{array}$$

(9)

where E is the elasticity modulus; h is the thickness of the cylindrical shell; $\nu$ is Poisson ratio; $N_1,N_2,N_{12},N_{21}$ are the normal forces and the shear forces in the middle-plain, respectively; and $M_1,M_2,M_{12},M_{21}$ are the middle-plain bending moments and the torques, respectively.

In the balance equation, considering the influence of the transverse shear $Q_2$ along the tangent of the circle (circumferential direction), we can obtain the following equation

$$\{\begin{array}{c}\multicolumn{1}{c}{\frac{\partial N_1}{\partial \alpha }+\frac{\partial N_{21}}{\partial \beta }+\mathrm{RX}=0}\\ \multicolumn{1}{c}{\frac{\partial N_{12}}{\partial \alpha }+\frac{\partial N_2}{\partial \beta }+Q_2+\mathrm{RY}=0}\\ \multicolumn{1}{c}{\frac{\partial Q_1}{\partial \alpha }+\frac{\partial Q_2}{\partial \beta }-N_2+\mathrm{RZ}=0}\\ \multicolumn{1}{c}{Q_1=\frac{1}{R}(\frac{\partial M_{21}}{\partial \beta }-\frac{\partial M_1}{\partial \alpha })}\\ \multicolumn{1}{c}{Q_2=\frac{1}{R}(\frac{\partial M_{12}}{\partial \alpha }-\frac{\partial M_2}{\partial \beta })}\end{array}$$

(10)

where X is the surface force along the axial direction of cylindrical shell, Y is the surface force along the circumferential direction, Z is the surface force along the normal direction, and $Q_1$ and $Q_2$ are transverse shear forces.

As a conclusion of the three equations above, the nonlinear differential equations of the long and thin cylindrical shell are as follows

$$\{\begin{array}{c}{L}_{11}\left(u\right)+L_{12}\left(v\right)+L_{13}\left(w\right)+{\Delta }_1+\frac{R^2(1-{\nu }^2)}{\mathrm{Eh}}X=0\hfill \\ \multicolumn{1}{c}{L_{21}\left(u\right)+L_{22}\left(v\right)+L_{23}\left(w\right)+{\Delta }_2+\frac{R^2(1-{\nu }^2)}{\mathrm{Eh}}Y=0}\\ \multicolumn{1}{c}{L_{31}\left(u\right)+L_{32}\left(v\right)+L_{33}\left(w\right)+{\Delta }_3-\frac{R^2(1-{\nu }^2)}{\mathrm{Eh}}Z=0}\end{array}$$

(11)

Equation (11) is a three differential equations about the displacement $u,v,w$ , with neglecting of relevant minor terms, and ${\Delta }_1,{\Delta }_2$ and ${\Delta }_3$ are nonlinear terms, only related to the displacement w, where

$$\{\begin{array}{c}{L}_{11}=\frac{{\partial }^2}{\partial {\alpha }^2}+\frac{1-\nu }{2}\frac{{\partial }^2}{\partial {\beta }^2},L_{12}=\frac{1+\nu }{2}\frac{{\partial }^2}{\partial \alpha \partial \beta }\hfill \\ \multicolumn{1}{c}{L_{13}=\nu \frac{\partial }{\partial \alpha }-c^2\frac{{\partial }^3}{\partial {\alpha }^3}+\frac{1-\nu }{2}{c}^2\frac{{\partial }^3}{\partial \alpha \partial {\beta }^2}}\end{array}$$

(12)

$$\{\begin{array}{c}\multicolumn{1}{c}{L_{21}=\frac{1+\nu }{2}\frac{{\partial }^2}{\partial \alpha \partial \beta }}\\ \multicolumn{1}{c}{L_{22}=\frac{1-\nu }{2}\frac{{\partial }^2}{\partial {\alpha }^2}+\frac{{\partial }^{2}v}{\partial {\beta }^2}}\\ \multicolumn{1}{c}{L_{23}=\frac{\partial }{\partial \beta }-\frac{3-\nu }{2}{c}^2\frac{{\partial }^3}{\partial {\alpha }^2\partial \beta }}\end{array}$$

(13)

$$\{\begin{array}{c}\multicolumn{1}{c}{L_{31}=\nu \frac{\partial }{\partial \alpha }-c^2(\frac{{\partial }^3}{\partial {\alpha }^3}-\frac{1-\nu }{2}\frac{{\partial }^3}{\partial \alpha \partial {\beta }^2})}\\ \multicolumn{1}{c}{L_{32}=\frac{\partial }{\partial \beta }-c^2\frac{3-\nu }{2}\frac{{\partial }^3}{\partial {\alpha }^2\partial \beta }}\\ \multicolumn{1}{c}{L_{33}=c^2({\nabla }^2{\nabla }^2+2\frac{{\partial }^2}{\partial {\beta }^2}+1)+1}\end{array}$$

(14)

where

$${\nabla }^2{\nabla }^2=\frac{{\partial }^4}{\partial {\alpha }^4}+2\frac{{\partial }^4}{\partial {\alpha }^2\partial {\beta }^2}+\frac{{\partial }^4}{\partial {\beta }^4}$$

(15)

$$c^2=\frac{h^2}{12R^2}$$

(16)

In order to get the approximate solutions of equation set (11), assume the form of the transverse displacement function of the cylindrical shell subjected to periodic axial loads is as follows

$$w(\alpha ,\beta ,t)=f(t)\sin n\alpha \cos k\beta$$

(17)

where $n=i\pi R/L$ , i and k are the axial and circumferential half wave number. L is the length of the cylindrical shell. Equation (17) has given the exact solutions for the linear problem of cylindrical shell subjected to periodic axial uniform loads. Substituting equation (17) into ${\Delta }_1,{\Delta }_2,{\Delta }_3$ , we can obtain the following equation

$${\Delta }_1=(\nu k^2-n^2)\frac{n{f}^2}{2R}\sin 2n\alpha {\cos }^{2}k\beta -\frac{(1+\nu )k^2}{2}\frac{n{f}^2}{2R}\sin 2n\alpha \cos 2k\beta$$

(18)

$${\Delta }_2=(k^2-\nu n^2)\frac{k{f}^2}{2R}\sin 2k\beta {\sin }^{2}n\alpha -\frac{(1+\nu )n^2}{2}\frac{k{f}^2}{2R}\sin 2k\beta \cos 2n\alpha$$

(19)

$${\Delta }_3=\frac{c^2{n}^2{f}^2}{R}\cos 2n\alpha [n^2{\cos }^{2}k\beta -k^2{\sin }^{2}k\beta +\frac{(1-\nu )k^2}{2}]+\frac{f^2{k}^2}{2R}{\sin }^{2}n\alpha {\sin }^{2}k\beta +\frac{\nu n^2{f}^2}{2R}{\cos }^{2}n\alpha {\cos }^{2}k\beta$$

(20)

Assume the form of solution of the first two formulas in the equation set (11) as follows

$$\{\begin{array}{c}u=A_1\sin 2n\alpha {\cos }^{2}k\beta +B_1\sin 2n\alpha \cos 2k\beta +u_0\hfill \\ \multicolumn{1}{c}{v=A_2\sin 2k\beta {\sin }^{2}n\alpha +B_2\sin 2k\beta \cos 2n\alpha +v_0}\\ \multicolumn{1}{c}{w=w_0=f(t)\sin n\alpha \cos k\beta }\end{array}$$

(21)

where $u_0,v_0$ , and $w_0$ satisfy the system of linear equations

$$\{\begin{array}{c}{L}_{11}\left(u\right)+L_{12}\left(v\right)+L_{13}\left(w\right)=0\hfill \\ \multicolumn{1}{c}{L_{21}\left(u\right)+L_{22}\left(v\right)+L_{23}\left(w\right)=0}\end{array}$$

(22)

And they have to satisfy the conditions of boundary. Substituting equation (21) into the first two formulas in the equation set (11) (neglect the tangential surface forces, that is, $X=Y=0$ ), the expressions of the coefficient $A_1,B_1,A_2,B_2$ can be obtained as follows

$$\{\begin{array}{c}\multicolumn{1}{c}{A_1=\frac{f^2(\nu k^2-n^2)}{8\mathrm{nR}},B_1=-\frac{\nu k^2{f}^2}{16\mathrm{nR}}}\\ \multicolumn{1}{c}{A_2=\frac{f^2(k^2-\nu n^2)}{8\mathrm{kR}},B_2=-\frac{\nu n^2{f}^2}{16\mathrm{kR}}}\end{array}$$

(23)

Combining equations (8), (9), (21), and (23), the axial force $N_1$ and the circumferential force $N_2$ can be obtained as follows

$$\{\begin{array}{c}{N}_1=\frac{Eh\nu k^2{f}^2}{8\left(1-{\nu }^2\right)R^2}\left[1+\frac{n^2}{\nu k^2}+\frac{\left(1-{\nu }^2\right)n^2}{\nu k^2}\cos 2k\beta \right]+N_1^0\\ N_2=\frac{Eh\nu n^2{f}^2}{8\left(1-{\nu }^2\right)R^2}\left[1+\frac{k^2}{\nu n^2}-\frac{\left(1-{\nu }^2\right)k^2}{\nu n^2}\cos 2n\alpha \right]+N_2^0\end{array}$$

(24)

where

$$\{\begin{array}{c}\multicolumn{1}{c}{N_1^0=\frac{\mathrm{Eh}}{(1-{\nu }^2)R}[\frac{\partial u_0}{\partial \alpha }+\nu \frac{\partial v_0}{\partial \beta }+(\nu +c^2{n}^2)w_0]}\\ \multicolumn{1}{c}{N_2^0=\frac{\mathrm{Eh}}{(1-{\nu }^2)R}[\nu \frac{\partial u_0}{\partial \alpha }+\frac{\partial v_0}{\partial \beta }+(1-k^2{c}^2)w_0]}\end{array}$$

(25)

Assume the form of the expressions of $u_0,v_0$ , and $w_0$ as follows

$$\{\begin{array}{c}{u}_0(\alpha ,\beta ,t)=U\cos n\alpha \cos k\beta \hfill \\ v_0(\alpha ,\beta ,t)=V\sin n\alpha \sin k\beta \hfill \\ w_0(\alpha ,\beta ,t)=f\sin n\alpha \cos k\beta \hfill \end{array}$$

(26)

where they can satisfy both boundary conditions (the radial and the circumferential displacement are 0, and the axial displacement is not 0) and the equation set (22). Substituting equation (26) into equation (22), we can obtain the following equation

$$\{\begin{array}{c}U=\frac{n[\nu n^2-k^2+c^2(n^4-k^4)]f}{{(n^2+k^2)}^2}\hfill \\ \multicolumn{1}{c}{V=-\frac{k[(\nu +2)n^2+k^2+2c^2{n}^2(n^2+k^2)]f}{{(n^2+k^2)}^2}}\end{array}$$

(27)

For the cylindrical shell subjected to the periodic axial loads $P_0+P_t\cos \theta t$ , the balance condition along the axial direction is as follows

$$\underset{0}{\overset{2\pi }{\int }}{N}_1(\alpha ,\beta )d\beta =-\frac{(P_0+P_t\cos \theta t)}{R}$$

(28)

where $p_0$ is the static component of pulsating axial load, $p_t$ is the harmonic component of pulsating axial load, and $\theta$ is the angular frequency of oscillations of the axial load.

Combining equations (24)–(28), $N_1$ and $N_2$ can be obtained

$$\{\begin{array}{c}\multicolumn{1}{c}{N_1=-\frac{P_0+P_t\cos \theta t}{2\pi R}+\frac{\mathrm{Eh}{f}^2}{8R^2}{n}^2\cos 2k\beta +\frac{\mathrm{Eh}{n}^2{k}^2}{(1-{\nu }^2)R}\left[\frac{(1-\nu )(1+\nu +2c^2{k}^2+2c^2{n}^2)}{{(n^2+k^2)}^2}\right]·f\sin n\alpha \cos k\beta }\\ \multicolumn{1}{c}{N_2=\frac{\mathrm{Eh}(\nu n^2+k^2)f^2}{8(1-{\nu }^2)R^2}-\frac{\mathrm{Eh}{f}^2}{8R^2}{k}^2\cos 2n\alpha +\frac{\mathrm{Eh}}{(1-{\nu }^2)R}[\frac{(1-{\nu }^2)n^4+(\nu -4)c^2{n}^2{k}^4}{{(n^2+k^2)}^2}-\frac{c^2\nu n^6+3c^2{n}^4{k}^2+c^2{k}^6}{{(n^2+k^2)}^2}]·}\\ \multicolumn{1}{c}{f\sin n\alpha \cos k\beta }\end{array}$$

(29)

In the third formula in the differential equations (11), the surface force along the normal direction of cylindrical shell is ²

$$Z=\frac{1}{R^2}[N_1\frac{{\partial }^{2}w}{\partial {\alpha }^2}+N_2(\frac{{\partial }^{2}w}{\partial {\beta }^2}+w)]-m\frac{{\partial }^{2}w}{\partial t^2}$$

(30)

where m is mass of per unit area in the middle-plane of cylindrical shell, $m=\rho h$ .

By combining the third formula of the differential equations (11), equation (20), (21), (23), (26), (27), (29), and (30), the nonlinear differential equations of transverse vibration of the cylindrical shell can be obtained as follows

$$\begin{array}{c}[f″+{\left({\omega \prime }_{n,k}\right)}^{2}f-n^2\frac{(P_0+P_t\cos \theta t)}{2\pi m{R}^3}f]\sin n\alpha \cos k\beta \hfill \\ \multicolumn{1}{c}{+\frac{\mathrm{EhB}}{m{R}^2(1-{\nu }^2)}{f}^2+[\frac{(1-{\nu }^2)n^4}{8R^2}\cos 2k\beta +}\\ \frac{(1-{\nu }^2)k^2(1-k^2)}{8R^2}\cos 2n\alpha -\frac{(1-k^2)(\nu n^2+k^2)}{8R^2}]·\hfill \\ \multicolumn{1}{c}{f^3\sin n\alpha \cos k\beta =0}\end{array}$$

(31)

where

$$\begin{array}{l}B=\frac{{\nu }^2{k}^2}{8R}\cos 2n\alpha +\frac{k^2}{4R}{\sin }^{2}n\alpha +\frac{\nu n^2}{8R}+\frac{n^2}{R}·\\ \frac{\left(1-\nu \right)n^2{k}^2\left(1+\nu +2c^2{k}^2+2c^2{n}^2\right)}{{\left(n^2+k^2\right)}^2}{\sin }^{2}n\alpha {\cos }^{2}k\beta \\ -\frac{\left(1-{\nu }^2\right)n^4-c^2\nu n^6+\left(\nu -4\right)c^2{n}^2{k}^4-3c^2{n}^4{k}^2-c^2{k}^6}{{\left(n^2+k^2\right)}^2}·\\ \frac{\left(1-k^2\right)}{R}·{\sin }^{2}n\alpha {\cos }^{2}k\beta \end{array}$$

(32)

where $\omega {\prime }_{n,k}$ is the inherent frequency of free vibration of cylindrical shell with the mode of vibration $(n,k)$

$${\left({\omega \prime }_{n,k}\right)}^2=\frac{\mathrm{Eh}{c}^{2}g\prime (n,k)}{m{R}^2(1-{\nu }^2)}$$

(33)

where

$$g\prime (n,k)=\frac{{(n^2+k^2-1)}^2{(n^2+k^2)}^2}{{(n^2+k^2)}^2}+\frac{2(1-\nu )(n^6-n^2{k}^4)+\frac{1-{\nu }^2}{c^2}{n}^4}{{(n^2+k^2)}^2}$$

(34)

Conducting Galerkin’s variation on the differential equation (31), the differential equation of nonlinear transverse vibration can be derived as follows

$$f″+{\left({\omega \prime }_{n,k}\right)}^{2}f-n^2\frac{(P_0+P_t\cos \theta t)}{2\pi m{R}^3}f+\frac{[n^4+k^2(k^2-1)-2(1-k^2)(\nu n^2+k^2)]\mathrm{Eh}}{16m{R}^4}{f}^3=0$$

(35)

Introducing the following notations of the critical load $P_{*}$ and the nonlinear coefficient ${\gamma }_{n,k}$

$$P_{*n,k}=\frac{2\pi m{R}^3{\left({\omega \prime }_{n,k}\right)}^2}{n^2}=\frac{2\pi \mathrm{Dg}\prime (n,k)}{n^{2}R}$$

(36)

$${\gamma }_{n,k}=\frac{[n^4+k^2(k^2-1)-2(1-k^2)(\nu n^2+k^2)]\mathrm{Eh}}{16m{R}^4}$$

(37)

Thus, the differential equation of transverse vibration can be written as follows

$$f″+{\left({\omega \prime }_{n,k}\right)}^2(1-\frac{P_0+P_t\cos \theta t}{P_{*n,k}})f+{\gamma }_{n,k}{f}^3=0$$

(38)

The differential equation of equation (38) is the second-order homogeneous differential equation with periodic coefficients and nonlinear terms, which is used for describing parametric vibration of cylindrical shell subjected to periodic axial loads. According to Bolotin method, equation (38) is written as the form of Mathieu equation as follows

$$f″+{\Omega }_{n,k}^2(1-2{\mu }_{n,k}\cos \theta t)f+{\gamma }_{n,k}{f}^3=0$$

(39)

where ${\Omega }_{n,k}$ is the inherent frequency of transverse vibration acted upon by the static component $P_0$ of pulsating axial loads and ${\mu }_{n,k}$ is the excitation coefficient as

$${\Omega }_{n,k}=\omega {\prime }_{n,k}\sqrt{1-\frac{P_0}{P_{*n,k}}},{\mu }_{n,k}=\frac{P_t}{2(P_{*n,k}-P_0)}$$

(40)

Amplitude of vibrations at the parametric resonance of cylindrical shell

As mentioned above, equation (39) is nonlinear Mathieu equation where the nonlinear term ${\gamma }_{n,k}{f}^3$ represents the effect of large deflection. According to Liapunov principle, dynamic unstable region is determined by the linear parts of equation (39). ² So, here the discussion is made only for the parametric resonance of the system. The basic solutions of Mathieu equation include two periodical solutions, that is, periodical solutions of T and 2T with $T=2\pi /\theta$ .

Steady-state amplitude of vibrations in the first order unstable region

In order to obtain the steady-state amplitude of vibrations in the first unstable region, Bolotin’s method ² is employed for parametric vibration. Assume the form of solution of period 2T as follows

$$f\left(t\right)=\sum _{j=1,3,5}^{\infty }(a_j\sin \frac{j\theta t}{2}+b_j\cos \frac{j\theta t}{2})$$

(41)

where $a_j$ and $b_j$ are arbitrary vectors.

Substituting equation (41) into equation (39), we obtain

$$\begin{array}{c}\sum _{j=1,3,5}^{\infty }({\Omega }_{n,k}^2-\frac{j^2{\theta }^2}{4})(a_j\sin \frac{j\theta t}{2}+b_j\cos \frac{j\theta t}{2})\hfill \\ \multicolumn{1}{c}{-{\Omega }_{n,k}^2{\mu }_{n,k}\sum _{j=1,3,5}^{\infty }{a}_j[\sin \frac{(j+2)\theta t}{2}+\sin \frac{(j-2)\theta t}{2}]}\\ -{\Omega }_{n,k}^2{\mu }_{n,k}\sum _{j=1,3,5}^{\infty }{b}_j[\frac{\cos (j+2)\theta t}{2}+\cos \frac{(j-2)\theta t}{2}]\hfill \\ \multicolumn{1}{c}{+{\psi }^{*}\left(f\right)=0}\end{array}$$

(42)

where ${\psi }^{*}(f)={\gamma }_{n,k}{f}^3$ is the nonlinear term, which is expanded by Fourier Series as follows

$${\psi }^{*}\left(f\right)=\sum _{j=1,3,5}^{\infty }({\Phi }_j\sin \frac{j\theta t}{2}+{\Psi }_j\cos \frac{j\theta t}{2})$$

(43)

where

$$\{\begin{array}{c}{\Phi }_j=\frac{\theta }{2\pi }\underset{0}{\overset{\frac{4\pi }{\theta }}{\int }}{\gamma }_{n,k}{f}^3\sin \frac{j\theta t}{2}\mathrm{dt}\hfill \\ \multicolumn{1}{c}{{\Psi }_j=\frac{\theta }{2\pi }\underset{0}{\overset{\frac{4\pi }{\theta }}{\int }}{\gamma }_{n,k}{f}^3\cos \frac{j\theta t}{2}\mathrm{dt}}\end{array}(j=1,2,3,\cdots )$$

(44)

Substituting equation (43) into equation (42), equating the coefficient of $\sin (j\theta t/2)$ and $\cos (j\theta t/2)$ terms as follows

$$\{\begin{array}{c}\multicolumn{1}{c}{[{\Omega }_{n,k}^2(1+{\mu }_{n,k})-\frac{{\theta }^2}{4}]a_1-{\Omega }_{n,k}^2{\mu }_{\mathrm{nk}}{a}_3+{\Phi }_1=0}\\ \multicolumn{1}{c}{[{\Omega }_{n,k}^2(1-{\mu }_{n,k})-\frac{{\theta }^2}{4}]b_1-{\Omega }_{n,k}^2{\mu }_{n,k}{b}_3+{\Psi }_1=0}\\ \multicolumn{1}{c}{({\Omega }_{n,k}^2-\frac{j^2{\theta }^2}{4})a_j-{\Omega }_{n,k}^2{\mu }_{n,k}(a_{j-2}+a_{j+2})+{\Phi }_j=0}\\ \multicolumn{1}{c}{({\Omega }_{n,k}^2-\frac{j^2{\theta }^2}{4})b_j-{\Omega }_{n,k}^2{\mu }_{n,k}(b_{j-2}+b_{j+2})+{\Psi }_j=0}\end{array}$$

(45)

If we investigate the vibration for the principal resonance $\theta =2\Omega$ , we can neglect the influence of higher harmonics in the expansion of equation (41) and can assume

$$f\left(t\right)=a_1\sin \frac{\theta t}{2}+b_1\cos \frac{\theta t}{2}$$

(46)

Neglecting terms containing higher harmonics, the system of equation (45) is simplified as follows

$$\{\begin{array}{c}[{\Omega }_{n,k}^2(1+{\mu }_{n,k})-\frac{{\theta }^2}{4}]a_1+{\Phi }_1=0\hfill \\ \multicolumn{1}{c}{[{\Omega }_{n,k}^2(1+{\mu }_{n,k})-\frac{{\theta }^2}{4}]b_1+{\Psi }_1=0}\end{array}$$

(47)

Substituting the solution of equation (46) into equation (44) and neglecting terms containing higher harmonics, we obtain the expression of two nonlinear terms ${\Phi }_1$ and ${\Psi }_1$ as follows

$${\Phi }_1=\frac{3{\gamma }_{n,k}{A}_1^2}{4}{a}_1,{\Psi }_1=\frac{3{\gamma }_{n,k}{A}_1^2}{4}{b}_1$$

(48)

So, equation (47) can be written as follows

$$\{\begin{array}{c}\multicolumn{1}{c}{[(1+{\mu }_{n,k})-\frac{{\theta }^2}{4{\Omega }_{n,k}^2}]a_1+\frac{3{\gamma }_{n,k}{A}_1^2}{4{\Omega }_{n,k}^2}{a}_1=0}\\ \multicolumn{1}{c}{[(1-{\mu }_{n,k})-\frac{{\theta }^2}{4{\Omega }_{n,k}^2}]b_1+\frac{3{\gamma }_{n,k}{A}_1^2}{4{\Omega }_{n,k}^2}{b}_1=0}\end{array}$$

(49)

where $A_1=\sqrt{a_1^2+b_1^2}$ is the amplitude of steady-state vibrations in the first resonance region. Obviously, when $a_1=b_1=A_1=0$ , the equation set (49) is true, corresponding to the situation where the cylindrical shell subjected to axial loading has no transverse oscillation.

In order to obtain the non-zero solution for the equation set (49), regard it as homogeneous linear equations with respect to $a_1$ and $b_1$ . This system has solutions that differ from zero only in the case where the determinant composed of the coefficient disappears

$$\left|\begin{array}{cc}1+\mu -{\left(\frac{\theta }{2{\Omega }_{n,k}}\right)}^2+\frac{3\gamma A_1^2}{4{\Omega }_{n,k}^2}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1-\mu -{\left(\frac{\theta }{2{\Omega }_{n,k}}\right)}^2+\frac{3\gamma A_1^2}{4{\Omega }_{n,k}^2}\hfill \end{array}\right|=0$$

(50)

Expanding the determinant and solving the resulting equation with respect to the amplitude, $A_1$ , of the steady-state vibrations in the first order unstable region, the following equation is obtained

$$\{\begin{array}{c}\multicolumn{1}{c}{A_{11}=\frac{2{\Omega }_{n,k}}{\sqrt{3\gamma }}\sqrt{{\left(\frac{\theta }{2{\Omega }_{n,k}}\right)}^2-1+\mu }}\\ \multicolumn{1}{c}{A_{12}=\frac{2{\Omega }_{n,k}}{\sqrt{3\gamma }}\sqrt{{\left(\frac{\theta }{2{\Omega }_{n,k}}\right)}^2-1-\mu }}\end{array}$$

(51)

Steady-state amplitude of vibrations in the second-order unstable region

In order to obtain the steady-state amplitude of vibrations in the second-order unstable region, assume the form of solution of period T is as follows

$$f\left(t\right)=b_0+a_2\sin \theta t+b_2\cos \theta t$$

(52)

Substituting equation (52) into equation (39), we obtain the equation set containing the coefficient $a_2$ and $b_0,b_2$ as follows

$$\{\begin{array}{c}{\Omega }_{n,k}^2(b_0-{\mu }_{n,k}{b}_2)+{\Psi }_0=0\hfill \\ \multicolumn{1}{c}{({\Omega }_{n,k}^2-{\theta }^2)a_2+{\Phi }_2=0}\\ \multicolumn{1}{c}{({\Omega }_{n,k}^2-{\theta }^2)b_2-2{\mu }_{n,k}{\Omega }_{n,k}^2{b}_0+{\Psi }_2=0}\end{array}$$

(53)

where ${\Phi }_2$ and ${\Psi }_2$ are calculated based on equation (44). The coefficient ${\Psi }_0$ is

$${\Psi }_0=\frac{\theta }{4\pi }\underset{0}{\overset{\frac{4\pi }{\theta }}{\int }}{\gamma }_{n,k}{f}^3\mathrm{dt}$$

(54)

Substituting equation (52) into equation (44) and (54) and neglecting terms containing higher harmonics, we obtain the expression of three nonlinear terms ${\Psi }_0$ , ${\Phi }_2$ and ${\Psi }_2$ as follows

$$\{\begin{array}{c}\multicolumn{1}{c}{{\Psi }_0=3{\gamma }_{n,k}{b}_0(\frac{2}{3}{b}_{0.}^2+A_2^2)}\\ \multicolumn{1}{c}{{\Phi }_2=\frac{3}{4}{\gamma }_{n,k}{a}_2(\frac{1}{4}{b}_{0.}^2+A_2^2)}\\ \multicolumn{1}{c}{{\Psi }_2=\frac{3}{4}{\gamma }_{n,k}{b}_2(\frac{1}{4}{b}_{0.}^2+A_2^2)}\end{array}$$

(55)

So, equation (53) can be written as follows

$$\{\begin{array}{c}(b_0-{\mu }_{n,k}{b}_2)+\frac{3{\gamma }_{n,k}{b}_0}{{\Omega }_{n,k}^2}(\frac{2}{3}{b}_0^2+A_2^2)=0\hfill \\ \multicolumn{1}{c}{(1-\frac{{\theta }^2}{{\Omega }_{n,k}^2})a_2+\frac{3{\gamma }_{n,k}{a}_2}{4{\Omega }_{n,k}^2}(\frac{1}{4}{b}_0^2+A_2^2)=0}\\ \multicolumn{1}{c}{(1-\frac{{\theta }^2}{4{\Omega }_{n,k}^2})b_2-2{\mu }_{n,k}{b}_0+\frac{3{\gamma }_{n,k}{b}_2}{4{\Omega }_{n,k}^2}(\frac{1}{4}{b}_0^2+A_2^2)=0}\end{array}$$

(56)

where $A_2=\sqrt{a_2^2+b_2^2}$ is the amplitude of steady-state vibrations in the second resonance region.

In order to obtain the non-zero solution for the equation set (56), regard it as homogeneous linear equations with respect to $a_2$ , $b_0$ and $b_2$ . This system has solutions that differ from zero only in the case where the determinant composed of the coefficient disappears

$$(1-\frac{{\theta }^2}{{\Omega }_{n,k}^2}+\frac{3{\gamma }_{n,k}}{4{\Omega }_{n,k}^2}{A}_2^2)·\left|\begin{array}{cc}1+\frac{3{\gamma }_{n,k}}{{\Omega }_{n,k}^2}{A}_2^2\hfill & -{\mu }_{n,k}\hfill \\ \multicolumn{1}{c}{-2{\mu }_{n,k}}& 1-\frac{{\theta }^2}{{\Omega }_{n,k}^2}+\frac{3{\gamma }_{n,k}}{4{\Omega }_{n,k}^2}{A}_2^2\hfill \end{array}\right|=0$$

(57)

Expanding the determinant and solving the resulting equation with respect to the amplitude, $A_2$ , of the steady-state vibrations in the second-order unstable region, the following equation is obtained

$$\{\begin{array}{c}\multicolumn{1}{c}{A_{21}=\frac{2{\Omega }_{n,k}}{\sqrt{3{\gamma }_{n.k}}}\sqrt{{\left(\frac{\theta }{{\Omega }_{n,k}}\right)}^2-1}}\\ \multicolumn{1}{c}{\begin{array}{c}\multicolumn{1}{c}{A_{22}=\frac{{\Omega }_{n,k}}{\sqrt{6{\gamma }_{n,k}}}\sqrt{\frac{4{\theta }^2}{{\Omega }_{n,k}^2}-5+\sqrt{{(\frac{4{\theta }^2}{{\Omega }_{n,k}^2}-3)}^2+32{\mu }_{n,k}^2}}}\\ \multicolumn{1}{c}{A_{23}=\frac{{\Omega }_{n,k}}{\sqrt{6{\gamma }_{n,k}}}\sqrt{\frac{4{\theta }^2}{{\Omega }_{n,k}^2}-5-\sqrt{{(\frac{4{\theta }^2}{{\Omega }_{n,k}^2}-3)}^2+32{\mu }_{n,k}^2}}}\end{array}}\end{array}$$

(58)

It can be proved that, within the equation set (58), $A_{23}$ is always imaginary number, so the steady-state amplitudes of vibrations in the second-order unstable region are $A_{21}$ and $A_{22}$ . Hence, the second-order unstable region is surrounded by two periodic solutions with the amplitudes $A_{21}$ , $A_{22}$ and period T.

Above all, for nonlinear analysis of dynamic stability for cylindrical thin shells, it can be known from equation (51) and (58) that, considering geometric nonlinearity factors, the amplitude $A_{\mathrm{ij}}$ of steady-state vibrations of the solution of period at the boundary of the dynamic unstable region is determined with the change of frequency $\theta$ of external loading. On the contrary, in linear analysis of dynamic stability for cylindrical thin shells, it is believed that infinitely increasing solution exists at the boundary of the dynamic unstable region, which is contradicted to the fact obviously. Therefore, nonlinear analysis of dynamic stability for cylindrical thin shells can quantitatively provide the amplitude of parametric vibration of the supercavitating vehicle under the different frequency of external loading, so as to lay foundation for determination of wetted area of the tail of the vehicle and analysis of hydrodynamic force.

Numerical results and nonlinear analysis of dynamic stability for supercavitating vehicle

The depth of sailing of supercavitating vehicle is $H=10\mathrm{m}$ , and the density of the flow field is ${\rho }_w=1000\mathrm{kg}/{\mathrm{m}}^3$ . For natural supercavity, the saturated vapor pressure is $p_c=2350\mathrm{Pa}$ $(20\circ \mathrm{C})$ , standard atmospheric pressure is $p=101,325\mathrm{Pa}$ . As for geometrical parameters of the cabin, the radius is $R=0.2\mathrm{m}$ , the length is $L=4\mathrm{m}$ , the thickness is $h=3\mathrm{mm}$ , and the diameter of cavitator is $d_n=0.12\mathrm{m}$ . About the physical parameter of material, elasticity modulus is $E=209\mathrm{GPa}$ , the material density is $\rho =7850\mathrm{kg}/{\mathrm{m}}^3$ and Poisson ratio is $\mu =0.3$ .

It can be concluded from section “Amplitude of vibrations at the parametric resonance of cylindrical shell” that, the amplitude $A_1(A_2)$ of the steady-state vibrations in the dynamic unstable region is related to both the mode $(n,k)$ of vibration of transverse vibrations and the frequency $\theta$ (or $\theta /2\Omega$ ) of axial loads, as well as the ratio $\delta$ of loads and the sailing speed V. Therefore, the influence of these factors on the curve of nonlinear parametric resonance will be analyzed in following part.

Analysis for the curves of nonlinear parametric resonance in the dynamic unstable region

Nonlinear parametric resonance curves under different modes of vibration with given speed and the ratio of loads

The first nonlinear parametric resonance curves under different modes of vibration with the given speed $(V=400\mathrm{m}/\mathrm{s})$ and the ratio $(\delta =0.4)$ of loads are shown in Figure 4 in which the stable and unstable solutions are specified with solid and dashed lines respectively.

Figure 4.

The first parametric resonance curves. (a) The mode $(i,1)$ of vibration cases, (b) The mode $(i,2)$ of vibration cases, (c) The mode $(i,3)$ of vibration cases.

It can be known from equation (48) that, if $A_1=0$ , then ${\theta }_{11}/2{\Omega }_{n,k}=\sqrt{1+\mu }$ and ${\theta }_{12}/2{\Omega }_{n,k}=\sqrt{1-\mu }$ . Take a mode $(i=1,k=1)$ as an example, ${\theta }_{11}/2{\Omega }_{n,k}$ corresponds to point G in Figure 1, and ${\theta }_{12}/2{\Omega }_{n,k}$ corresponds to point C in Figure 1. So, the field of abscissa axis surrounded by two curves (a solid line and a dashed line) stands for the exciting region of this mode.

Conclusions can be drawn from the three figures above as follows:

1. When the frequency of external loads applied to cylindrical shell is less than the excitation frequency (i.e. the line BC in Figure 1 (B is a point on the left side of C)), the cylindrical shell has no the transverse vibration.

2. The parametric resonance curves of any order all incline towards the side where it is greater than the frequency of exciting region. So, the frequency of external loads gradually increases from the side where it’s less than the lower bound of exciting region (B → C → G → F in Figure 1 (F is a point on the right side of G)), and the amplitude of transverse vibration of cylindrical shell increases along BCHD. For the part FG outside region of resonance, if external disturbance is not enough to make the transverse vibration of shell go beyond unstable solutions (curve EG), the amplitude of transverse vibration of shell is still 0. Until the frequency reaches the point G, the shell may have the unexpected transverse vibration with the amplitude at the point H. If disturbance can make the transverse vibration of shell go beyond unstable solutions to reach stable solution (curve DH), the amplitude of transverse vibration of cylindrical shell suddenly increases to the curve DH from zero, which will result in the part GF becoming the potential exciting region, and the exciting region causing danger will be enlarged.

3. It can be known from comparison between the vertical coordinates in Figure 4(a)–(c) that, when $k=2$ or 3, the steady-state amplitudes of vibrations in the unstable region is the same order of magnitude of the thickness of shell $(h=3\mathrm{mm})$ , but when $k=1$ , the steady-state amplitudes of vibrations are two orders of magnitudes more than that when $k=2$ or 3.

4. Comparing the widths of unstable regions in different modes of vibration cases in Figure 4, it is found that the unstable regions have a maximum width for the mode $(i=k)$ . Especially, the steady-state amplitude of vibration for the mode(1, 1) is much larger than those associated with other modes in the unstable regions, which indicates that instability for the mode(1, 1) is most likely to occur.

Nonlinear parametric resonance curves under the mode(1,1) of vibration with the ratio $(\delta =0.4)$ of loads

The nonlinear parametric resonance curves for the mode(1,1) of vibration with the ratio $(\delta =0.4)$ of loads are shown in Figure 5. Figure 5 shows the relations between the amplitude $A_1$ of the steady-state vibrations, the frequency $\theta /2\Omega$ of axial loads and the sailing speed V in the exciting region, which can make an a priori prediction of the sailing speed and the amplitude of the steady-state vibrations for the problems of dynamic stability of the supercavitating vehicle. The amplitude of the steady-state vibrations increases with the sailing speed under the same frequency (see $A(\theta /2\Omega =1.002,V=160\mathrm{m}/\mathrm{s},A_1=0.07145)\to B(\theta /2\Omega =1.002,V=285\mathrm{m}/\mathrm{s},A_1=0.08563)$ in Figure 5) and with the frequency under the same sailing speed in the dynamic unstable region (see $B(V=285\mathrm{m}/\mathrm{s},\theta /2\Omega =1.002,A_1=0.08563)\to C(V=285\mathrm{m}/\mathrm{s},\theta /2\Omega =1.002,A_1=0.1655)$ in Figure 5).

Figure 5.

The nonlinear parametric resonance curves for the mode(1,1) of vibration with the ratio $(\delta =0.4)$ of loads.

The critical sailing speed can be determined by the amplitude of the steady-state vibrations under the given ratio of loads and frequency, as shown in Figure 6.

Figure 6.

Amplitude of the steady-state vibrations and critical sailing speed for the mode(1,1) of vibration with $\delta =0.4$ and $\theta /2\Omega =1$ .

Comparison between linear and nonlinear theory of parametric resonance curves

By taking the mode $(i=1,k=1)$ as an example, comparing linear and nonlinear theory of parametric vibration, conclusions can be drawn as follows:

1. In linear theory of parametric vibration, the inclination of parametric resonance curves does not exist, and the only exciting region is the line CG in Figure 4(a). But in nonlinear theory of parametric vibration considering geometric nonlinearity, due to the inclination of parametric vibration curve towards the side where it is more than the excitation frequency, the exciting region includes the line CG and GF. Hence, by considering geometric nonlinearity, the enlargement of the exciting region will cause shrinkage of the safe frequency range of external loads and decrease in dynamic stability of the supercavitating vehicle.

2. The linear theory of parametric vibration assumes the amplitude of vibrations in the exciting region is infinitely increasing, which does not conform to engineering practice. Actually, because of geometric nonlinearity, the amplitude of vibrations in the exciting region is a finite value. For the supercavitating vehicle, the deformation of structure will directly affect the wetted area which determines stability of cavitation.

Influence of the speed and the ratio of loads on nonlinear parametric resonance curves with the given mode of vibration $(i=2,k=2)$

The first nonlinear parametric resonance curves under different sailing speeds V with the given the ratio $(\delta =0.4)$ of loads are shown in Figure 5. The first nonlinear parametric resonance curves under different ratios $\delta$ of loads with the given the sailing speed $(V=400\mathrm{m}/\mathrm{s})$ of loads are shown in Figure 6.

We can draw a conclusion from Figures 7 and 8 that the widths of exciting region and the amplitude of steady-state vibrations in the exciting region increase with the increase in the sailing speed V and ratios δ of loads.

Figure 7.

The first nonlinear parametric resonance curves under different sailing speeds.

Figure 8.

The first nonlinear parametric resonance curves under different ratios of loads.

The same conclusion as that in the first order can be drawn as well by analyzing the second nonlinear parametric resonance curves as shown in Figure 8 which provides the comparison between the nonlinear resonance curves in the first and second order ( $\delta =0.4$ and $V=400\mathrm{m}/\mathrm{s}$ ).

It can be known from Figure 9 that the shape of nonlinear resonance curve in the second order is similar to that in the first order, both inclining towards the side of greater frequency; the widths of exciting region and the amplitude of steady-state vibrations in the second order are much smaller than that in the first order.

Figure 9.

Comparison of nonlinear parametric resonance between the first and the second order.

Conclusion

The nonlinear dynamic stability for the thin cylindrical shell of the supercavitating vehicle subjected to periodic axial loading has been investigated in this article. A numerical example is given to analyze factors affecting nonlinear parametric resonance curves of the cylindrical shell. The following conclusions have been drawn:

1. The nonlinear parametric resonance curves incline towards the side where it is greater than the excitation frequency, which enlarges the width of the exciting region, so as to result in the shrinkage of the safe frequency range of external loads and decrease in the dynamic stability of the vehicle.

2. The amplitude of steady-state vibrations in the exciting region can be obtained by analysis of nonlinear parametric resonance, which is the basis of fluid–structure interaction analysis of the supercavitating vehicle.

3. The width of the exciting region increases with the increase in the sailing speed V and the ratios $\delta$ of loads of the supercavitating vehicle; the amplitude of steady-state vibrations in the exciting region increases with the increase in the sailing speed V and the ratios $\delta$ of loads of the supercavitating vehicle.

4. when $k=2$ or 3, the steady-state amplitudes of vibrations in the unstable region is the same order of magnitude of the thickness of shell $(h=3\mathrm{mm})$ , but when $k=1$ , the steady-state amplitudes of vibrations is two orders of magnitudes more than that when $k=2$ or 3.

The supercavitating vehicle is here modeled as the thin cylindrical shell. The dynamic stability analysis of cylindrical shells is performed in order to evaluate the widths of exciting region and the amplitude of steady-state vibrations in the exciting region. The analysis demonstrates the complex characteristics of nonlinear parametric resonance with respect to thin cylindrical shell structures. The presented results suggest that optimal structural designs of the supercavitating vehicles may be identified to achieve stability at given operating speeds and under specified periodic axial loading. The analysis should then be extended to evaluate optimal design configurations for supercavitating vehicles with superior buckling characteristics. The optimization procedure will consider the amplitude of steady-state vibrations of transverse oscillation as an objective function, so that an improvement in the thickness of the vehicles can be achieved.