Supersonic Flutter of a Spherical Shell Partially Filled with Fluid

In the present study, a hybrid finite element method is applied to investigate the dynamic behavior of a spherical shell partially filled with fluid and subjected to external supersonic airflow. The structural formulation is a combination of linear spherical shell theory and the classic finite element method. In this hybrid method, the nodal displacements are derived from exact solution of spherical shell theory rather than approximated by polynomial functions. Therefore, the number of elements is a function of the complexity of the structure and it is not necessary to take a large number of elements to get rapid convergence. Linearized first-order potential (piston) theory with the curvature correction term is coupled with the structural model to account for aerodynamic loading. It is assumed that the fluid is incompressible and has no free surface effect. Fluid is considered as a velocity potential at each node of the shell element where its motion is expressed in terms of nodal elastic displacements at the fluid-structure interface. Numerical simulation is done and vibration frequencies are obtained. The results are validated using numerical and theoretical data available in literature. The investigation is carried out for spherical shells with different boundary conditions, geometries, filling ratios, flow parameters, and radius to thickness ratios. Results show that the spherical shell loses its stability through coupled-mode flutter. This proposed hybrid finite element method can be used efficiently for analyzing the flutter of spherical shells employed in aerospace structures at less computational cost than other commercial FEM software.


Introduction
Shells of revolution, particularly spherical shells are one of the primary structural elements in high speed aircraft.
Their applications include the propellant tank or gas-deployed skirt of spacecraft.Due to the aerodynamic shape combined with thin wall thicknesses, spherical shells are more disposed to dynamic instability or flutter induced by high Mach number gas flow.It is therefore important to understand the effect of different flow parameters and loadings on their aeroelastic response.
Aeroelastic analysis of shells and plates has been studied by numerous researchers experimentally and analytically [1].Dowell gives an exhaustive study of the aeroelasticity of shells and plates in his book [2].After introducing the application of piston theory in the aeroelastic modeling presented by Ashley and Zatarian [3], a number of interesting experimental and theoretical studies were carried out to investigate supersonic flutter of cylindrical shells.In general, all of this research was concerned with the development of an analytical relation to describe the effect of shell and flow parameters on the critical flutter dynamic pressure.Aeroelastic models in combination with linear or nonlinear piston theory were coupled to the theory of shells to account for fluidstructure interaction.The resulting governing equations were treated numerically using the Galerkin method.A comprehensive experimental test was done by Fung and Olson [4].They studied the effects of shell boundary conditions and initial stress state due to internal pressure and axial load.It was observed that pressurized cylindrical shell fluttered at a lower level of freestream static pressure than predicted by theory [5].Later, Evensen and Olson [6] [7] presented a nonlinear analysis to take account of this observed effect.Dowell [8] also analyzed the behavior of a cylindrical shell in supersonic flow for different flow and shell parameters.A complete description of panel flutter modeling is given in his book [2].A study by Carter and Strearman [9] showed that agreement between the theory and experiments reported in the literature exists in cases that involve a small amount of static preload acting on the shell.Amabili and Pellicano [10] included geometric nonlinearities in their study of supersonic flutter of a circular cylindrical shell.By selecting expansion modes to discretize the aeroelastic equations, they were able to facilitate their solution, and therefore succeeded in capturing the nonlinear behavior of the shell correctly.
There are also some researchers who focused their efforts on the numerical study of this problem.The equations of virtual displacements were solved using the finite elements method.Aeroelastic governing equations were formulated by applying classical shell theory coupled with the piston theory for evaluation of aerodynamic forces.For example, Bismarck-Nasr [11] developed a FEM applied to supersonic flutter of circular shell subjected to internal pressure and axial loading.Ganapathi et al. [12] modeled an orthotropic and laminated anisotropic cylindrical shell in supersonic flow using FEM and analyzed the effect of different shell geometries on the flutter boundaries.
Aeroelasticity of conical shells has also been investigated by few researchers.The leading work in this field was conducted by Shulman [13].Ueda et al. [14] investigated theoretically and experimentally the supersonic flutter of a conical shell.Dixon and Hudson [15] studied the flutter and vibration of an orthotropic conical shell theoretically.Miserentino and Dixon [16] investigated experimentally the vibration and flutter of a pressurized truncated conical shell.Bismarck-Nasr and Costa-Savio [18] studied the supersonic flutter of conical shells using finite element method.Sunder et al. [18] successfully applied the finite element analysis to calculate the flutter of a laminated conical shell.In another study they found the optimum cone angle in aeroelastic flutter [19].Mason and Blotter [20] used a finite element technique to find the flutter boundary for a conical shell (a typical rocket nozzle element) subjected to an internal supersonic gas flow.Pidaparti and Yang Henry [21] completed a theoretical study to predict the onset of flutter instability for composites conical shells.
An analytical approach to the supersonic flutter of spherical shell becomes very complicated if one wishes to include different parameters.Therefore, the efficiency of numerical methods such as the finite element method (FEM) is an advantage for cases involving changes to all factors affecting flutter boundaries.The aim of the present study is to develop a hybrid finite element method in order to predict the aeroelastic behavior of isotropic spherical shells with different parameters as boundary conditions, geometries, flow parameters, filling ratios and radius to thickness ratios.The finite element is a spherical frustum instead of the usual rectangular shell element.Linear thin shell theory is coupled with linear piston theory.In the case of a fluid filled shell the effect of dynamic pressure acting on the wall is modeled based on a velocity potential formulation and Bernoulli's equation.It is assumed that the fluid is incompressible and has no free surface effect.The linear mass, damping and stiffness matrices are obtained.The aeroelastic equation of motion is reduced to a standard eigenvalue problem.The flutter boundary is found by analyzing the real and imaginary parts of the eigenvalues as freestream pressure is varied.

Equilibrium Equations
In this study the structure is modeled using hybrid finite element method which is a combination of spherical shell theory and classical finite element method.In this hybrid finite element method, the displacement functions are found from exact solution of spherical shell theory rather than approximated by polynomial functions as done in classical finite element method.In the spherical coordinate system (R, θ, ϕ) shown in Figure 1, five out of the six equations of equilibrium derived in reference [22] for spherical shells are written as follows: ( ) ( ) ( ) where N φ , N θ , N φθ are membrane stress resultants; M ϕ , M θ , M ϕθ the bending stress resultants and Q ϕ , Q θ the shear forces (Figure 2).The sixth equation, which is an identity equation for spherical shells, is not presented here.

Constitutive Relations
Strains and displacements in axial, U φ , radial, W , and circumferential, U θ directions are related as follows: Displacements U , W and V in the global Cartesian coordinate system are related to displacements The stress vector { } σ is expressed as a function of strain { } ε by: where [ ] P is the elasticity matrix for an anisotropic shell given by: Upon substitution of Equations ( 2), ( 4) and (5) into Equation (1), a system of equilibrium equations can beobtained as a function of displacements: These three linear partial differential operators 1 L , 2 L and 3 L are given in Appendix A, and ij P are ele- ments of the elasticity matrix, which for an isotropic thin shell with thickness h is given by: [ ] ( ) ( ) where is the membrane stiffness and ( ) is the bending stiffness.

Kinematic Relations
The element is a circumferential spherical frustum shown in Figure 3.It has two nodal circles with four degrees of freedom; axial, radial, circumferential and rotation at each node.This element type makes it possible to use thin shell equations easily to find the exact solution of displacement functions rather than an approximation with polynomial functions as done in classical finite element method.
For motions associated with the circumferential wave number n, we may write: The transversal displacement ( ) n w φ can be expressed as [22]: where ( ) ( ) and where ( ) Q µ φ are the associated Legendre functions of the first and second kinds respec- tively of order n and degree i µ .
The expression of the axial displacement u ϕν (ϕ) is: where the coefficient E i is given by: ( ) ( ) ( )( ) The auxiliary function ψ is given by the expression: Finally the circumferential displacement u ϕν (ϕ) can be expressed as: ( ) The degree i µ is obtained from the expression where i λ is one the roots of the cubic equation: and where ) Substituting Equations ( 18), (19) and (20) in Equations ( 9), (11) and (14) we have: ) In deriving the above relation we used the recursive relations: Using matrix formulation, the displacement functions can be expressed as follows: ( The vector { } C is given by the expression: The elements of matrix [ ] R are given in Appendix B. In the finite element method, the vector C is eliminated in favor of displacements at elements nodes.At each finite element node, the three displacements (axial, transversal and circumferential) and the rotation are applied.The displacement of node i are defined by the vector: The finite element shown in Figure 3 with two nodal lines (i and j) and eight degrees of freedom will have the following nodal displacement vector: with The terms of matrix, obtained from the values of matrix [ ] Finally, one substitutes the vector { } C into Equation ( 26) and obtains the displacement functions as fol- lows: The strain vector { } ε can be determined from the displacement functions , , and the deformation -displacement equation ( 2) as: where matrix [ ]

Mass and Stiffness Matrices
This relation can be used to find the stress vector, Equation ( 4), in terms of the nodal degrees of freedom vector: Based on the finite element formulation, the local stiffness and mass matrices are: where ρ is the density and h is the thickness of shell.
The surface element of the shell wall is 2).After integrating over θ , the pre- In the global system, the element stiffness and mass matrices are where From these equations, one can assemble the mass and stiffness matrices for each element to obtain the mass and stiffness matrices for the whole shell: [ ] M and [ ] K .Each elementary matrix is 8 × 8, therefore the final dimensions of [ ] M and [ ] K will be 4*(N + 1) where N is the number of elements of the shell.

Aerodynamic Modeling
Piston theory, introduced by Ashley and Zartarian [3], is a powerful tool for aeroelasticity modeling.In this study the fluid-structure effect due to external pressure loading can be taken into account using linearized firstorder potential theory.This pressure is expressed as: where P ∞ ,U ∞ , M and γ are the freestream static pressure, freestream velocity, Mach number and adiabatic exponent of air respectively.If the Mach number is sufficiently high ( ) 2 M ≥ , and curvature term, ( ) is neglected, the result is the so-called piston theory: where a ∞ is the free stream speed of sound.Finally, the aerodynamic pressure in terms of radial displacement is written: and the pressure loading in terms of nodal degrees of freedom is written as: where ρ ∞ the freestream air density and m R is the median radius for each element.Based on thermodynamic relations the freestream pressure and velocity can be linked together using the following relations: The matrix The matrix where matrix [ ] C is given by: [ ] The general force vector due to a pressure field is written as: The local damping matrix is given by: Finally the local stiffness matrix is given by: In the global system, the element damping and stiffness matrices are: From these equations, one can assemble the damping and stiffness matrices for each element to obtain the damping and stiffness matrices for the whole shell:

Fluid Modeling
The Laplace equation satisfied by velocity potential for inviscid, incompressible and irrotational fluid in the spherical system is written as: where the velocity components are: Using the Bernouilli equation, hydrodynamic pressure in terms of velocity potential ϕ and fluid density f ρ is found as: The impermeability condition, which ensures contact between the shell surface and the peripheral fluid, is written as: Method of separation of variables for the velocity potential solution can be done as follows: , , , , Placing this relation into the impermeability condition (50), we can find the function Hence the equation becomes With substitution of the above equation into Laplace Equation (47), the following second order equation in terms of ( ) Solution of the above differential equation yields the following: ( ) For internal flow Finally, the hydrodynamic pressure in terms of radial displacement is written: We put: And the pressure loading in terms of nodal degrees of freedom is written as: where matrix where [ ] F is expressed as: [ ] 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 The general force vector due the fluid pressure loading is given by: After substituting for pressure field vector and matrix [ ] In the global system the element stiffness and mass matrices are From these equations, one can assemble the mass for each element to obtain the mass matrix for the whole shell: The governing equation which accounts for fluid-shell interaction in the presence external supersonic airflow is derived as: where subscripts s and f refer to shells in vacuo and fluid respectively.

Eigenvalue Problem
The global fluid matrices mentioned in Equation ( 65) may be obtained, respectively, by superimposing the mass, damping and stiffness matrices for each individual fluid finite element.After applying the boundary conditions the global matrices are reduced to square matrices of order 4*(N + 1) − J, where N is the number of finite elements in the shell and J is the number of constraints applied.Finally, the eigenvalue problem is solved by means of the equation reduction technique.Equation (65) may be rewritten as follows: where δ is the global displacement vector. where and [ ] I is the identity matrix.
An in house computer code based on the finite element method was developed as part of this work to establish the structural and fluid matrices of each element based on equations developed using the theoretical approach.The calculations for each finite element are performed in two stages: the first dealing with solid shell and the second with the effect of the flowing fluid.Aeroelastic stability will be examined by studying the eigenvalues in the complex plane.When the imaginary part of ω becomes negative the amplitude of the shell motion grows exponentially with time, thus indicating dynamic instability.The flutter boundary is obtained numerically by tracing the eigenvalues to see when the sign of imaginary part just changes from positive to negative.For the fixed value of circumferential wave number n, the onset of instability is determined by varying the value of freestream static pressure.This procedure is repeated for different values of n until the minimum critical pressure is obtained.

Results and Discussion
In this section numerical results are presented and compared with existing experimental, analytical and numerical data.

Validation and Comparison
For the cases investigated in the present paper, the predicted dimensionless frequencies are expressed by the following relation: where: ω is the natural angular frequency, R is the radius of the reference surface, ρ is the density, and E is the modulus of elasticity.
Results for different boundary conditions, geometries, flow parameters and radius to thickness ratios compared to experimental, theoretical and numerical analyses are presented (see Figure 4).

Spherical Shells in Vacuo
Case 1: clamped spherical shell with 0 φ = 10˚ Narassihan and Alwar [23] investigated the case of an axisymmetric clamped spherical shell.The analysis is based on the application of the Chebyshev-Galerkin spectral method for the evaluation of free vibration frequencies and mode shapes.Sai Ram and Sreedhar Babu [24] analyzed the same case with the classical finite element method using 80 elements.Each element is an eight nodded degenerated is oparametric shell element with nine degrees of freedom at each node.With our model and using 6 finite elements, the natural frequencies were computed; the results are shown in Table 1.
Case 2: clamped spherical shell with 0 φ = 30˚ This case was investigated analytically by Kalnins [25] using classical theory and transverse vibration theory.With our theory, we used 8 finite elements to study the spherical shell with the results shown in Table 2.The frequencies we obtained with our model are very comparable to Kalinin's values.
Case 4: spherical shell with 0 φ = 90˚ Kraus [22] investigated the case of simply supported spherical shell using a general theory, which included the effects of transverse shear stress and rotational inertia.For cases both with and without these effects, he determined the natural frequencies for the shell motion that was independent of θ for circumferential mode number 0 n = .Tessler and Spiridigliozzi [31], Gautham and Ganesan [34] analyzed the case of clamped he- mispherical shell.Ventsel et al. [35] studied the case of simply supported spherical shell using the boundary elements method for various circumferential mode numbers ( ) . With our model and using 12 finite elements, the natural frequencies were computed for clamped and simply supported shells.The results are shown respectively in Table 5 and Table 6.

Flutter of Spherical Shells
The problem treated for validation is the flutter boundary of a simply-supported spherical shell subjected to external supersonic airflow.As there is no information available for flutter of spherical shells, this case has been compared with simply-supported cone studied by various authors.The conical shell has the following data: Young's Modulus, E = 6.5 106 lb-in-2, Poisson's ratio, ν = 0.29, material mass density, ρ = 8.33 10-4 lb-s2-in-4, shell thickness, h = 0.051 in, cone semi-vertex angle α = 5˚.The supersonic airflow has freestream Mach number, M∞ = 3, stagnation temperature, T∞ = 288.15K.The results are shown in Table 7 where Λ is the dynamic pressure parameter defined as: where ( ) is the bending stiffness.
When results are summarized and compared with other finite element and analytical solutions, this method shows good convergence using only 15 elements with small disagreements.It should be noted that the previous analytical methods [13] [15] use Donnel-Mushtari simplified shell theory while [14] uses Novozhilov's thin shell with the different method of application of finite element solution.On the other hand, a complete form of the linear piston theory is used by [21] as in the present study and the results are very close; but the expression used by Dixon and Hudson [15], Ueda et al. [14] for the piston theory does not have a curvature term which has caused greater differences in the results.

Flutter Boundary
Flutter which is observed in all the papers using piston theory is a coupled-mode flutter.Indeed, let us consider motion of the shell eigenvalues in the complex ω plane.If the freestream pressure is not very high, and the shell is stable, all complex frequencies are located in the top ω half-plane.Let us now increase freestream pressure.The first and the second complex frequencies move toward each other, almost merge, and then go away from each other in vertical directions (Figure 5).Thus, interaction of two modes occurs.Physically this interaction of the shell modes happens through the effect of the air flow.
A simply supported spherical shell with 0 φ = 30˚ is treated here.The complex frequencies only for the first and second modes versus freestream dynamic pressure are plotted in Figure 6.Aerodynamic pressure is evaluated using Equation (36).In Figure 6(a) the real part of the complex frequency increases for the first mode while for the second mode it decreases as the freestream dynamic pressure parameter Λ increases.For higher  Pidaparti and Yang Henri [21] Shulman [13] Bismark-Nasr [17] 520(5) a 590( 5) 609( 5) 576( 5) 669( 6) 702( 6) values of dynamic pressure these real parts, representing the oscillation frequency, eventually coalesce into a single mode.Further increasing the dynamic pressure of the flow causes the shell to lose its stability at Λ cr = 410.This instability is due to coupled-mode flutter where the imaginary part of complex frequency (representing the damping term of the aeroelastic system) becomes zero for certain critical pressure (Figure 6(b)).
The same behaviour is observed by real and imaginary parts of complex frequencies as the static pressure increases (Figure 7) but the onset of flutter is at Λ cr = 410 if the freestream static pressure is evaluated using Equation (37).Prediction of the critical freestream static pressure using Equation ( 36) provides approximately the same results when evaluating the pressure field using Equation (37).As expected, using the piston theory with the correction term to account for shell curvature produces a better approximation for the pressure loading acting on a curved shell exposed to supersonic flow.
In Figure 8 the onset of flutter for different angles is plotted.By increasing the angle 0 φ , flutter instability occurs at lower pressure.This decrease in Λ cr with 0 φ is attributed to the fact that the natural frequencies al- ways decrease as the angle 0 φ is increased.The effect of radius to thickness ratio R/h is presented in Figure 9.This figure shows an increase of Λ cr with an increase of radius to thickness ratio.This increase in Λ cr is attributed to the fact that the mass of shell is greater when the shell is thick, and the effect of pressure is less important for a thick shell than for a thin shell.On the other hand, when the shell is thin it becomes unstable at higher dynamical pressure levels due to an increase in stiffness because of a decrease in thickness.The same conclusion is reported in [12] for the flutter of  In order to study the effect of filling ratio, Figure 10 shows the critical value of freestream static pressure for different filling ratios, H/R.Shell geometry and flow parameters are the same as the previous case study with liquid filled density ρf =9.355 × 10-5lb s2 in-4.It is seen that the value of critical dynamic pressure parameter decreases as the filling ratio increases from a low value.This rapid change in critical dynamic pressure at low filling ratios and its almost steady behaviour at large filling ratios indicates that the fluid near the bottom of the shell is largely influenced by elastic deformation when a shell is subjected to supersonic flow.
The effect of boundary conditions on the flutter onset is presented in Table 8.It is seen that for freely simply supported ends, v = w = 0, flutter onset occurs at Λ = 510.5 which indicates more flutter resistance compared to simply supported or clamped ends.It is indicated that there is no difference for flutter onset when the shell is either clamped or simply supported.We obtained the same results in the conical shells subjected to supersonic flow.

Conclusion
An efficient hybrid finite element method is presented to investigate the aeroelastic stability of an empty or partially liquid filled spherical shell subjected to external supersonic flow.Linear shell theory is coupled with first order piston theory to account for aerodynamic pressure.The effect of curvature correction in piston theory was   analyzed.Fluid structure interaction due to hydrodynamic pressure of internal fluid is also taken into account.
The study has been done for shells with various geometries, radius to thickness ratios, filling ratios and boundary conditions.all study cases one type of instability is found; coupled-mode flutter in the first and second mode.Increasing the radius to thickness ratio leads the onset of flutter to occur at higher dynamic pressure.Decreasing the angle 0 φ of the spherical shell causes the flutter boundary to occur at lower dynamic pressure.A lower filling ratio has more flutter resistance than a higher filling ratio.The proposed hybrid finite element formulation can give reliable results at less computational cost compared to commercial software since the latter imposes some restrictions when such analysis is done.

References
sin sin sin        ( )  ( ) ( )( ) In deriving the above relation we used the recursive relations:

Figure 1 .
Figure 1.Geometry of the spherical shell.
matrix of the constant C i as a function of the degree of freedom:

N
in the above equation, the local matrix from the following: added fluid mass.The eigenvalue problem is given by:

Figure 5 .
Figure 5. Trajectories of the complex frequencies loci in the complex ω plane during the changing of the dynamic pressure.

Figure 6 .
Figure 6.(a) Real part and (b) imaginary part of the complex frequencies versus the freestream static pressure parameter; static pressure evaluated by Equation (36).

Figure 7 .
Figure 7. (a) Real part and (b) imaginary part of the complex frequencies versus the freestream static pressure parameter; static pressure evaluated by Equation (37).

Figure 8 .Table 8 .Figure 9 .
Figure 8. Variation of the critical freestream static pressure parameter with angle 0 φ for simply supported shell.

Figure 10 .
Figure 10.Variation of the critical freestream static pressure parameter with R/H for simply supported shell.

Table 1 .
Normalized natural frequencies for 10˚ clamped spherical shell with R/h = 200.

Table 2 .
Normalized natural frequencies for 30˚ clamped spherical shell with R/h = 20.

Table 3 .
Normalized natural frequencies for 60˚ clamped spherical shell with R/h = 20.

Table 4 .
Normalized natural frequencies for 60˚ simply supported spherical shell with R/h = 20.

Table 5 .
Normalized natural frequencies for 90˚ clamped spherical shell with R/h = 10.

Table 6 .
Normalized natural frequencies for 90˚ simply supported spherical shell.

Table 7 .
Comparison of critical dynamical pressure parameter (simply supported case). ) )