Hypersonic Aeroelastic Response of Elastic Boundary Panel Based on a Modified Fourier Series Method

In this paper, an analytical method is proposed to directly obtain the aeroelastic time domain response of the elastic boundary panel. Based on a modified Fourier series method (MFSM), the vibration analysis of elastic boundary panels is carried out, after the dynamics equation of the panel is obtained. Then, the vibrational functions are combined with the supersonic piston theory to establish the aeroelastic equation. The aeroelastic time domain response of the panel is obtained to analyze the flutter speed of the panel more intuitively. Finally, the flutter speeds of panels with different length-width ratios, thicknesses, and elastic boundary conditions are discussed in detail.


Introduction
Hypersonic aircraft has gradually become a new-generation aircraft developed by various countries.Hypersonic vehicle is a frontier research field in aerospace engineering.Up to now, there are still many problems that have to be overcome.One of the most common problems of the hypersonic vehicle is flutter.As early as 1970, Dowell summarized the stability analysis of the flutter problems for plates and shells [1].So far, a large number of scholars have adopted various methods for flutter analysis.
The piston theory is commonly used in supersonic aerodynamics.Flutter analysis based on local flow piston theory was developed by Yang and Song [2], which was used to calculate a supersonic wing with attack angle.Then, Zhang et al. [3] summarized and utilized the method of supersonic flutter analysis based on local piston theory.The local piston theory was improved by Yang et al. [4] to analyze aeroelastic behaviors of curved panels.With the local flow field parameters obtained by CFD technique, the modified local piston theory is more reliable and enlarges the application range of piston theory for curved panels.In recent years, the first-order piston theory was used by Yao et al. [5] for the high-speed rotating cantilever rectangular plate.Besides, the third-order piston theory was used by Shao et al. [6] for a smart laminated panel.The piston theory has been applied well in supersonic aerodynamic analysis.
In terms of structure, the finite element method (FEM) is widely used for panel flutter analysis.Mei [7] developed a FEM to analyze the flutter behavior of nonlinear plates.Later, the method was applied to the large-amplitude panel flutter of thin laminates [8] and composite panels [9].Then, the method was used to analyze curved panel flutter [10], and the use of FEM is gradually well-developed [11][12][13].
In recent years, orthogonal decomposition methods have been used for flutter analysis.Xie and Xu [14] first applied the method for von Karman plate under supersonic flow.The method was improved and used for a cantilever plate in supersonic flow [15].The reduced-order model was proved to be suitable for nonlinear aeroelastic oscillations [16].Flutter analysis of a nonlinear panel was carried out based on proper orthogonal decomposition method [17].
In addition, the differential quadrature method (DQM) [18] and other methods have also been used for flutter analysis of plates.However, the structural models used in the flutter analysis of plates are usually plates with simple classical boundaries, and the study on elastic boundary plates focuses on free vibration instead of flutter.
By now, there have been many researches on the vibration analysis of plates but most of them limited to classical boundary conditions, i.e., fixed boundary, simply supported boundary, and free boundary.A lot of the studies were devoted to the vibration of rectangular plates with general elastic constraints along the edge.Li [19] analyzed the vibration of rectangular plates with general elastic boundary supports.Then, Du et al. [20] used a modified Fourier series method to obtain an analytical solution for the in-plane vibrations of a rectangular plate with elastically restrained edges.Li et al. [21] summarized the modified Fourier series method (MFSM) and proposed a complete set of analytical solutions for the transverse vibration of rectangular plates with general elastic boundary supports.Later on, the MFSM has been applied to solve many problems of plates with elastic boundary restraints, such as free vibration of two elastically coupled rectangular plates [22], modal analysis of general plate structures [23], and modeling analysis of elastically restrained panel [24].This method is also used well in triangular plates [25], blades [26], circular plates [27], confocal annular elliptic plates [28], and so on.
In engineering, the boundary conditions of the highspeed aircraft are much more complicated, which is difficult to be described accurately with classical boundary conditions.Elastic boundary conditions are more universal and flexible and can be reduced to classical boundary conditions.At present, there are just few works directly combining the vibration analysis and aerodynamic force of the elastic boundary panel.Zhou et al. [29] presented a method for supersonic flutter analysis for a Mindlin orthotropic plate with general boundary conditions and analyzed the factors that affect flutter characteristics.
In this paper, an analytical method is proposed to obtain the aerodynamic elastic response of the elastic boundary panel in time domain directly, so as to study the flutter problem more intuitively.

Structure Model of the Panel with Elastic Boundary
2.1.Governing Equation of the Plate.The panel is an essential structural element widely used in aircrafts, missiles and launch vehicles, and so on.The thickness of the panel is much smaller than that of the other two directions in practice, and the panel is usually considered as a thin plate.Consider a rectangular plate with elastic constraints at any edge(s), as shown in Figure 1.The middle surface of the plate is the xy plane, and the z-axis is perpendicular to the plate.The length of the plate in the x and y directions are denoted by a and b, and the thickness is h.As shown in Figure 1, k x0 , k xa , k y0 , and k yb represent the linear spring  International Journal of Aerospace Engineering constants at x = 0, x = a, y = 0, and y = b.K x0 , K xa , K y0 , and K yb represent the rotational spring constants at x = 0, x = a, y = 0, and y = b, respectively.The governing differential equation for free vibration of the thin plate is given by D∇ 4 w x, y − ρhω 2 w x, y = 0, 1 where ∇ 4 = ∂ 4 /∂x 4 + 2∂ 4 /∂x 2 ∂y 2 + ∂ 4 /∂y 4 , w x, y is the transverse displacement along z direction, ω and ρ are the natural frequency and the mass density of the plate, and D = Eh 3 /12 1 − μ 2 is the flexural rigidity, where E and μ are the Young modulus and the Poisson ratio, respectively.In terms of transverse displacements, the bending moments, twisting moment, and transverse shearing forces can be expressed as follows: The general elastic boundary conditions expressed by equations (3a), (3b), (3c), (3d), (3e), (3f), (3g), and (3h) can be reduced to classical homogeneous boundary conditions by setting the corresponding spring constants.For example, if all the linear spring constants are set to be extremely large, while the constants of all rotational springs are set as 0, then it can be reduced to the case of simply supported ones.
Substituting equations (2a), (2b), (2c), (2d), and (2e) into equations (3a), (3b), (3c), (3d), (3e), (3f), (3g), and (3h) leads to 3. Analytical Solution for the Plate with Elastic Boundary.By using the Fourier series method, the transverse displacement of the plate can be expanded into the following form: A mn cos λ am x cos λ bn y where λ am = mπ/a, λ bn = nπ/b, ξ l a x , and ξ l b y must be sufficiently smooth and can satisfy the elastic boundary conditions at four edges.Here, the functions ξ l a x and ξ l b y are third-order derivable and continuous at any point 3 International Journal of Aerospace Engineering of the plate.With this in mind, here, the supplementary functions are chosen as the following: In equations (6a) and (6b), it is easy to verify that only , and all the other first and third derivatives at the edges are equal to zero.By choosing the supplementary functions in such a way, the two-dimensional Fourier series expansions in equation ( 5) can represent a residual displacement field which is third-order derivable and continuous over the entire x − y domain.Such supplementary functions make the solution expressed in equation ( 5) to satisfy the governing equation at every field point and the boundary conditions at every boundary point.More importantly, it is now guaranteed to converge uniformly at a substantially improved speed for any boundary condition [21].
According to the relationship between p and a which has been defined by equation (8), equation ( 11) can be written as By solving the corresponding characteristic equation of equation (12), the eigenvalues λ 2 = ρhω 2 /D and eigenvectors a can be obtained as well as a series of natural frequencies ω and corresponding vectors a and p.
where p ∞ , a ∞ , and κ are constants, indicating the atmospheric pressure of the incoming flow, the speed of sound of the incoming flow, and the specific heat capacity ratio, respectively.The model in equation (13) describes the relationship between the lateral velocity v z of the panel and the local pressure p.
Given that the partial derivative of transverse displacements w x, y, t to x is equal to the velocity of flow and the partial derivative to y is 0, the lateral velocity v z of the panel can be expressed as follows: where U ∞ is the flow velocity.
On the basis of equations ( 13) and ( 14), the pressure difference between the upper and lower surfaces of the panel can be obtained as By separating the variables of the transverse displacement and cut at Nth modes, the form of the solution is assumed as follows: Substituting equations ( 15) and ( 17) into equation ( 16) leads to 7 International Journal of Aerospace Engineering where q i t represents the second derivative of q i t versus time, and W i x, y represents the ith vibration mode function of the panel.
According to the modal orthogonality, after multiplying both ends of equation ( 18) by W i and integrating the new equation, the ith order equation can be obtained All of N equations can be rewritten as a form of state-space equation as follows: By the MFSM described before, the solution for the mode function can be obtained Without loss of generality, W x, y is used to represent W i x, y in equation (21).When given a and p, W x, y can be reduced to the corresponding W i x, y .Substituting equation (21) into equations (20a), (20b), and (20c) results in q 1 , q 2 , … , q N , then the aeroelastic response of the plate in the time domain can be obtained.

Results and Discussion
4.1.Vibration Analysis.The vibration of plates with various boundary conditions and different aspect ratios is discussed here.First, consider a panel fully simply supported along four edges, and the parameters of the panel are shown in Case 1 of Table 1.
The simply supported boundary can be viewed as a special case when the linear spring constants become extremely large, which is set to be a very large number as 1 0 × 10 8 N/m, while the rotational springs set to be 0.
Table 2 presents the first six natural frequencies (for a = b = 1 m) calculated with different truncate numbers to verify the convergence of the solution.Since the results show no difference for M = N larger than 20, the series expansion will be truncated to M = N = 20 in all the subsequent cases.
In Table 3, the first six natural frequencies of the panel are given for the simply supported panels (S-S-S-S) with different aspect ratios.The frequencies obtained from International Journal of Aerospace Engineering MFSM are compared with results from FEM and Navier's analytical solutions.A unified grid cell of size 0 01 m × 0 01 m is used to obtain the converged FEM results, as Figure 2. Obviously, Table 3 shows that the results of MFSM are very close to the exact solutions.The analytical method based on MFSM in this article is precise enough and even much closer to Navier's analytical solution than the results based on FEM. Figure 3 shows the 1st, 2nd, 3rd, and 5th natural frequencies of the panel with four clamped edges (C-C-C-C).The  9 International Journal of Aerospace Engineering clamped boundary can be viewed as another special case when both the linear and rotational spring constants become extremely large, which are set to be very large numbers as 1 0 × 10 8 N/m and 1 0 × 10 8 N/rad.In Figure 3, the 2nd and 3rd modes of the panel have the same natural frequency when b = a = 1 0 m, and the frequencies calculated by MFSM are compared with results obtained by FEM.Furthermore, Figure 4 shows the 1st, 2nd, 3rd, and 5th natural frequencies of the panel with S-F-S-F (which are simply supported at x = 0 and x = a and free at y = 0 and y = b) and C-F-C-F (which are clamped at x = 0 and x = a and free at y = 0 and y = b) boundary conditions which are obtained by MFSM and FEM.

Aeroelastic Analysis.
A specific example involving various boundary conditions is discussed here.The boundary condition is gradually shifted from four simply supported edges (S-S-S-S) to four clamped edges (C-C-C-C).The parameters of the panel are shown in Case 2 of Table 1.
First, considering the simplest case, S-S-S-S, the method is suitable for supersonic aeroelastic calculation, when Mach number is much larger than 1 (Ma ∞ ≫ 1).In this example, the air density at high altitude is 0.3 kg/m 3 , and the specific heat capacity ratio κ is 1.4 typicality.From Figure 5, it can be clearly seen the panel vibration response of convergence (Figures 5(a) and 5(b)) and divergence (Figures 5(c) and 5(d)).Naturally, the critical velocity between convergence and divergence is the flutter speed.The flutter speed in this example is 869 m/s, which is very close to the result calculated by the V-g method, 867.5 m/s, seen in Figure 6.The V-g method introduces artificial damping g into the system.When the artificial damping g is negative, no flutter will occur, and the point where g equals 0 is considered as the critical flutter speed.Figure 7 shows the aerodynamic mesh of the lifting surface employed in the V-g method.The airflow is forward along the x-axis, which is parallel to the chord.When meshing, the chordwise edge of the grid must be parallel to the airflow.Considering that the shape of the panel is rectangular, a unified aerodynamic grid cell of size 0 05 m × 0 05 m is adopted in this paper, shown in Figure 7. 10 International Journal of Aerospace Engineering 0 0035 m.In Figure 8, the overall trend is that the flutter speed decreases slowly as the width increases.Moreover, as the width-length ratio increases, the flutter speeds predicted by these two methods match better.Figure 9 shows flutter speeds of panels, which length is 1.0 m and the width is 1.3 m.The flutter speed increases with the increasing of thickness.In addition, the results of MFSM together with the piston theory are slightly lower than those of reference values computed by the V-g method in flutter computation, which is based on FEM for the simulation of the structural behavior and on lifting surface theory for the aerodynamics.
Figure 10 shows that the flutter speed of the panel increases gradually when the simply supported boundary shifts to the clamped boundary.In this case, the simply supported boundary shifts to the clamped boundary, which means the linear spring constants keep extremely large which are set as very large numbers such as 1 0 × 10 10 N/m, while the rotational springs K x0 , K xa , K y0 , and K yb change from 0 to 1 0 × 10 10 N/rad, shown in Table 4.When the stiffness of torsional springs reaches 1 0 × 10 6 N/rad, the boundary condition can be considered as clamped, and the flutter speed results equal to 1668 m/s.

Conclusion
An analytical method based on the MFSM is proposed to obtain the aeroelastic time domain response of the elastic boundary panel.Compared with the results from FEM and Navier's analytical solution, the results based on the MFSM prove to be accurate and reliable.Based on the supersonic piston theory, the aeroelastic time domain response of the panel is obtained to analyze the flutter speed.The flutter speeds of panels with different length-width ratios, thicknesses, and elastic boundary conditions are analyzed in detail.Obviously, this proposed method based on MFSM and supersonic piston theory is reliable in flutter analysis of the panel with elastic boundary conditions.A 3f

Figure 1 :
Figure 1: A rectangular plate with elastic constraints at all edges.

Figure 2 :Figure 3 :
Figure 2: The structural finite element mesh with unified grid cells.

Figure 5 :
Figure 5: Vibration response of the S-S-S-S panel.

Figure 9 :
Figure 9: Flutter speeds of panels in different thicknesses with width b = 1 3 m.

Table 2 :
Frequencies for S-S-S-S square panel with different truncate numbers.

Table 3 :
Frequencies for S-S-S-S panel with different aspect ratios.
MFSM for results from modified Fourier series method.FEM for results from finite element method.NAS for results from Navier's analytical solution ω mn = π/2 m 2 /a 2 + n 2 /b 2 D/ρh .

Table 1 :
Physical parameters of the panel.