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

Authors: Xiyue Zhou; Lifeng Wang; Jingnong Jiang; Zhu Su
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

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.


---

## Body

## 1. 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–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 high-speed 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.


## 2. 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 Figure1. 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, kx0, kxa, ky0, and kyb represent the linear spring constants at x=0, x=a, y=0, and y=b. Kx0, Kxa, Ky0, and Kyb represent the rotational spring constants at x=0, x=a, y=0, and y=b, respectively.Figure 1
A rectangular plate with elastic constraints at all edges.The governing differential equation for free vibration of the thin plate is given by(1)D∇4wx,y−ρhω2wx,y=0,where ∇4=∂4/∂x4+2∂4/∂x2∂y2+∂4/∂y4, wx,y is the transverse displacement along z direction, ω and ρ are the natural frequency and the mass density of the plate, and D=Eh3/121−μ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:(2a)Mx=−D∂2w∂x2+μ∂2w∂y2,(2b)My=−D∂2w∂y2+μ∂2w∂x2,(2c)Mxy=−D1−μ∂2w∂x∂y,(2d)Qx=−D∂∂x∇2w+∂Mxy∂y=−D∂3w∂x3+2−μ∂3w∂x∂y2,(2e)Qy=−D∂∂y∇2w+∂Mxy∂x=−D∂3w∂y3+2−μ∂3w∂x2∂y.


### 2.2. Elastic Boundary Conditions

The boundary conditions of the rectangular thin plate with elastic constraints at all edges are given as follows:(3a)kx0w=Qx,atx=0,(3b)Kx0∂w∂x=−Mx,atx=0,(3c)kxaw=−Qx,atx=a,(3d)Kxa∂w∂x=Mx,atx=a,(3e)ky0w=Qy,aty=0,(3f)Ky0∂w∂y=−My,aty=0,(3g)kybw=−Qy,aty=b,(3h)Kyb∂w∂y=My,aty=b.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
(4a)kx0w=−D∂3w∂x3+2−μ∂3w∂x∂y2atx=0,(4b)Kx0∂w∂x=D∂2w∂x2+μ∂2w∂y2atx=0,(4c)kxaw=D∂3w∂x3+2−μ∂3w∂x∂y2atx=a,(4d)Kxa∂w∂x=−D∂2w∂x2+μ∂2w∂y2atx=a,(4e)ky0w=−D∂3w∂y3+2−μ∂3w∂x2∂yaty=0,(4f)Ky0∂w∂y=D∂2w∂y2+μ∂2w∂x2aty=0,(4g)kybw=D∂3w∂y3+2−μ∂3w∂x2∂yaty=b,(4h)Kyb∂w∂y=−D∂2w∂y2+μ∂2w∂x2aty=b.


### 2.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:(5)wx,y=∑m=0∞∑n=0∞Amncosλamxcosλbny+∑l=14ξbly∑m=0∞cmlcosλamx+ξalx∑n=0∞dnlcosλbny,where λam=mπ/a, λbn=nπ/b, ξalx, and ξbly must be sufficiently smooth and can satisfy the elastic boundary conditions at four edges. Here, the functions ξalx and ξbly are third-order derivable and continuous at any point of the plate. With this in mind, here, the supplementary functions are chosen as the following:
(6a)ξa1x=9a4πsinπx2a−a12πsin3πx2a,ξa2x=−9a4πcosπx2a−a12πcos3πx2a,ξa3x=a3π3sinπx2a−a33π3sin3πx2a,ξa4x=−a3π3cosπx2a−a33π3cos3πx2a,(6b)ξb1y=9b4πsinπy2b−b12πsin3πy2b,ξb2y=−9b4πcosπy2b−b12πcos3πy2b,ξb3y=b3π3sinπy2b−b33π3sin3πy2b,ξb4y=−b3π3cosπy2b−b33π3cos3πy2b.In equations (6a) and (6b), it is easy to verify that only dξa10/dx=dξa2a/dx=1, d3ξa30/dx3=d3ξa4a/dx3=1, dξb10/dy=dξb2b/dy=1, and d3ξb30/dy3=d3ξb4b/dy3=1, 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].Substituting equation (5) into equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), and (4h) results in
(7a)−kx0D∑m=0∞∑n=0∞Amncosλbny+∑l=14ξbly∑m=0∞cml+ξal0∑n=0∞dnlcosλbny=∑l=14∂3ξal0∂x3∑n=0∞dnlcosλbny+2−μ∑l=14∂ξal0∂x∑n=0∞−dnlλbn2cosλbny,atx=0,(7b)Kx0D∑l=14∂ξal0∂x∑n=0∞dnlcosλbny=∑m=0∞∑n=0∞−Amnλam2cosλbny+∑l=14ξbly∑m=0∞−cmlλam2+∂2ξal0∂x2∑n=0∞dnlcosλbny+μ∑m=0∞∑n=0∞−Amnλbn2cosλbny+∑l=14∂2ξbly∂y2∑m=0∞cml+ξal0∑n=0∞−dnlλbn2cosλbny,atx=0,(7c)kxaD∑m=0∞∑n=0∞Amn−1mcosλbny+∑l=14ξbly∑m=0∞cml−1m+ξala∑n=0∞dnlcosλbny=∑l=14∂3ξala∂x3∑n=0∞dnlcosλbny+2−μ∑l=14∂ξala∂x∑n=0∞−dnlλbn2cosλbny,atx=a,(7d)−KxaD∑l=14∂ξala∂x∑n=0∞dnlcosλbny=∑m=0∞∑n=0∞−Amnλam2−1mcosλbny+∑l=14ξbly∑m=0∞−cmlλam2−1m+∂2ξala∂x2∑n=0∞dnlcosλbny+μ∑m=0∞∑n=0∞−Amnλbn2−1mcosλbny+∑l=14∂2ξbly∂y2∑m=0∞cml−1m+ξala∑n=0∞−dnlλbn2cosλbny,atx=a,(7e)−ky0D∑m=0∞∑n=0∞Amncosλamx+∑l=14ξbl0∑m=0∞cmlcosλamx+ξalx∑n=0∞dnl=∑l=14∂3ξbl0∂y3∑m=0∞cmlcosλamx+2−μ∑l=14∂ξbl0∂y∑m=0∞−cmlλam2cosλamx,aty=0,(7f)Ky0D∑l=14∂ξbl0∂y∑m=0∞cmlcosλamx=∑m=0∞∑n=0∞−Amnλbn2cosλamx+∑l=14∂2ξbl0∂y2∑m=0∞cmlcosλamx+ξalx∑n=0∞−dnlλbn2+μ∑m=0∞∑n=0∞−Amnλam2cosλamx+∑l=14ξbl0∑m=0∞−cmlλam2cosλamx+∂2ξalx∂x2∑n=0∞dnl,aty=0,(7g)kybD∑m=0∞∑n=0∞Amn−1ncosλamx+∑l=14ξblb∑m=0∞cmlcosλamx+ξalx∑n=0∞dnl−1n=∑l=14∂3ξblb∂y3∑m=0∞cmlcosλamx+2−μ∑l=14∂ξblb∂y∑m=0∞−cmlλam2cosλamx,aty=b,(7h)−KybD∑l=14∂ξbly∂y∑m=0∞cmlcosλamx=∑m=0∞∑n=0∞−Amnλbn2−1ncosλamx+∑l=14∂2ξbly∂y2∑m=0∞cmlcosλamx+ξalx∑n=0∞−dnlλbn2−1n+μ∑m=0∞∑n=0∞−Amnλam2−1ncosλamx+∑l=14ξbly∑m=0∞−cmlλam2cosλamx+∂2ξalx∂x2∑n=0∞dnl−1n,aty=b.In equations (7a), (7b), (7c), (7d), (7e), (7f), (7g), and (7h), a total of 4M+1+4N+1 equations can be obtained by taking m=0,1,…,M, n=0,1,…,N. After simplifying, these equations can be written in matrix form as
(8)Hp=Qa.In equation (8), H and Q are constant matrices which are gotten from equations (7a), (7b), (7c), (7d), (7e), (7f), (7g), and (7h), and vectors p and a are variables which are given by
(9a)p=c01,c11,…,cM1,c02,c12,…,cM2,c03,c13,…,cM3,c04,c14,…,cM4,d01,d11,…,dN1,d02,d12,…,dN2,d03,d13,…,dN3,d04,d14,…,dN4T,(9b)a=A00,A10,…,AM0,A01,A11,…,AM1,……,A0N,A1N,…,AMNT.The relationship between the two sets of variablesp and a in equation (8) is given by means of the eight elastic boundary conditions represented by equations (7a), (7b), (7c), (7d), (7e), (7f), (7g), and (7h). Substituting the displacement expression equation (5) into the governing equation of the thin plate equation (1) leads to
(10)∑m=0∞∑n=0∞λam4+λbn4+2λam2λbn2Amncosλamxcosλbny+∑l=14∑m=0∞ξblyλam4−2∂2ξbly∂y2λam2+∂4ξbly∂y4cmlcosλamx+∑n=0∞ξalxλbn4−2∂2ξalx∂x2λbn2+∂4ξalx∂x4dnlcosλbny−ρhω2D∑m=0∞∑n=0∞Amncosλamxcosλbny+∑l=14ξbly∑m=0∞cmlcosλamx+ξalx∑n=0∞dnlcosλbny=0.Similarly, takingm=0,1,…,M, n=0,1,…,N, after simplifying, equation (10) can be written as the following matrix equation:
(11)K~a+Bp−ρhω2DM~a+Fp=0.The matricesK~, M~, B, and F are defined in the appendix.According to the relationship betweenp and a which has been defined by equation (8), equation (11) can be written as
(12)K−ρhω2DMa=0,where K=K~+BH−1Q and Μ=M~+FH−1Q.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.


## 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 Figure1. 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, kx0, kxa, ky0, and kyb represent the linear spring constants at x=0, x=a, y=0, and y=b. Kx0, Kxa, Ky0, and Kyb represent the rotational spring constants at x=0, x=a, y=0, and y=b, respectively.Figure 1
A rectangular plate with elastic constraints at all edges.The governing differential equation for free vibration of the thin plate is given by(1)D∇4wx,y−ρhω2wx,y=0,where ∇4=∂4/∂x4+2∂4/∂x2∂y2+∂4/∂y4, wx,y is the transverse displacement along z direction, ω and ρ are the natural frequency and the mass density of the plate, and D=Eh3/121−μ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:(2a)Mx=−D∂2w∂x2+μ∂2w∂y2,(2b)My=−D∂2w∂y2+μ∂2w∂x2,(2c)Mxy=−D1−μ∂2w∂x∂y,(2d)Qx=−D∂∂x∇2w+∂Mxy∂y=−D∂3w∂x3+2−μ∂3w∂x∂y2,(2e)Qy=−D∂∂y∇2w+∂Mxy∂x=−D∂3w∂y3+2−μ∂3w∂x2∂y.


## 2.2. Elastic Boundary Conditions

The boundary conditions of the rectangular thin plate with elastic constraints at all edges are given as follows:(3a)kx0w=Qx,atx=0,(3b)Kx0∂w∂x=−Mx,atx=0,(3c)kxaw=−Qx,atx=a,(3d)Kxa∂w∂x=Mx,atx=a,(3e)ky0w=Qy,aty=0,(3f)Ky0∂w∂y=−My,aty=0,(3g)kybw=−Qy,aty=b,(3h)Kyb∂w∂y=My,aty=b.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
(4a)kx0w=−D∂3w∂x3+2−μ∂3w∂x∂y2atx=0,(4b)Kx0∂w∂x=D∂2w∂x2+μ∂2w∂y2atx=0,(4c)kxaw=D∂3w∂x3+2−μ∂3w∂x∂y2atx=a,(4d)Kxa∂w∂x=−D∂2w∂x2+μ∂2w∂y2atx=a,(4e)ky0w=−D∂3w∂y3+2−μ∂3w∂x2∂yaty=0,(4f)Ky0∂w∂y=D∂2w∂y2+μ∂2w∂x2aty=0,(4g)kybw=D∂3w∂y3+2−μ∂3w∂x2∂yaty=b,(4h)Kyb∂w∂y=−D∂2w∂y2+μ∂2w∂x2aty=b.


## 2.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:(5)wx,y=∑m=0∞∑n=0∞Amncosλamxcosλbny+∑l=14ξbly∑m=0∞cmlcosλamx+ξalx∑n=0∞dnlcosλbny,where λam=mπ/a, λbn=nπ/b, ξalx, and ξbly must be sufficiently smooth and can satisfy the elastic boundary conditions at four edges. Here, the functions ξalx and ξbly are third-order derivable and continuous at any point of the plate. With this in mind, here, the supplementary functions are chosen as the following:
(6a)ξa1x=9a4πsinπx2a−a12πsin3πx2a,ξa2x=−9a4πcosπx2a−a12πcos3πx2a,ξa3x=a3π3sinπx2a−a33π3sin3πx2a,ξa4x=−a3π3cosπx2a−a33π3cos3πx2a,(6b)ξb1y=9b4πsinπy2b−b12πsin3πy2b,ξb2y=−9b4πcosπy2b−b12πcos3πy2b,ξb3y=b3π3sinπy2b−b33π3sin3πy2b,ξb4y=−b3π3cosπy2b−b33π3cos3πy2b.In equations (6a) and (6b), it is easy to verify that only dξa10/dx=dξa2a/dx=1, d3ξa30/dx3=d3ξa4a/dx3=1, dξb10/dy=dξb2b/dy=1, and d3ξb30/dy3=d3ξb4b/dy3=1, 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].Substituting equation (5) into equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), and (4h) results in
(7a)−kx0D∑m=0∞∑n=0∞Amncosλbny+∑l=14ξbly∑m=0∞cml+ξal0∑n=0∞dnlcosλbny=∑l=14∂3ξal0∂x3∑n=0∞dnlcosλbny+2−μ∑l=14∂ξal0∂x∑n=0∞−dnlλbn2cosλbny,atx=0,(7b)Kx0D∑l=14∂ξal0∂x∑n=0∞dnlcosλbny=∑m=0∞∑n=0∞−Amnλam2cosλbny+∑l=14ξbly∑m=0∞−cmlλam2+∂2ξal0∂x2∑n=0∞dnlcosλbny+μ∑m=0∞∑n=0∞−Amnλbn2cosλbny+∑l=14∂2ξbly∂y2∑m=0∞cml+ξal0∑n=0∞−dnlλbn2cosλbny,atx=0,(7c)kxaD∑m=0∞∑n=0∞Amn−1mcosλbny+∑l=14ξbly∑m=0∞cml−1m+ξala∑n=0∞dnlcosλbny=∑l=14∂3ξala∂x3∑n=0∞dnlcosλbny+2−μ∑l=14∂ξala∂x∑n=0∞−dnlλbn2cosλbny,atx=a,(7d)−KxaD∑l=14∂ξala∂x∑n=0∞dnlcosλbny=∑m=0∞∑n=0∞−Amnλam2−1mcosλbny+∑l=14ξbly∑m=0∞−cmlλam2−1m+∂2ξala∂x2∑n=0∞dnlcosλbny+μ∑m=0∞∑n=0∞−Amnλbn2−1mcosλbny+∑l=14∂2ξbly∂y2∑m=0∞cml−1m+ξala∑n=0∞−dnlλbn2cosλbny,atx=a,(7e)−ky0D∑m=0∞∑n=0∞Amncosλamx+∑l=14ξbl0∑m=0∞cmlcosλamx+ξalx∑n=0∞dnl=∑l=14∂3ξbl0∂y3∑m=0∞cmlcosλamx+2−μ∑l=14∂ξbl0∂y∑m=0∞−cmlλam2cosλamx,aty=0,(7f)Ky0D∑l=14∂ξbl0∂y∑m=0∞cmlcosλamx=∑m=0∞∑n=0∞−Amnλbn2cosλamx+∑l=14∂2ξbl0∂y2∑m=0∞cmlcosλamx+ξalx∑n=0∞−dnlλbn2+μ∑m=0∞∑n=0∞−Amnλam2cosλamx+∑l=14ξbl0∑m=0∞−cmlλam2cosλamx+∂2ξalx∂x2∑n=0∞dnl,aty=0,(7g)kybD∑m=0∞∑n=0∞Amn−1ncosλamx+∑l=14ξblb∑m=0∞cmlcosλamx+ξalx∑n=0∞dnl−1n=∑l=14∂3ξblb∂y3∑m=0∞cmlcosλamx+2−μ∑l=14∂ξblb∂y∑m=0∞−cmlλam2cosλamx,aty=b,(7h)−KybD∑l=14∂ξbly∂y∑m=0∞cmlcosλamx=∑m=0∞∑n=0∞−Amnλbn2−1ncosλamx+∑l=14∂2ξbly∂y2∑m=0∞cmlcosλamx+ξalx∑n=0∞−dnlλbn2−1n+μ∑m=0∞∑n=0∞−Amnλam2−1ncosλamx+∑l=14ξbly∑m=0∞−cmlλam2cosλamx+∂2ξalx∂x2∑n=0∞dnl−1n,aty=b.In equations (7a), (7b), (7c), (7d), (7e), (7f), (7g), and (7h), a total of 4M+1+4N+1 equations can be obtained by taking m=0,1,…,M, n=0,1,…,N. After simplifying, these equations can be written in matrix form as
(8)Hp=Qa.In equation (8), H and Q are constant matrices which are gotten from equations (7a), (7b), (7c), (7d), (7e), (7f), (7g), and (7h), and vectors p and a are variables which are given by
(9a)p=c01,c11,…,cM1,c02,c12,…,cM2,c03,c13,…,cM3,c04,c14,…,cM4,d01,d11,…,dN1,d02,d12,…,dN2,d03,d13,…,dN3,d04,d14,…,dN4T,(9b)a=A00,A10,…,AM0,A01,A11,…,AM1,……,A0N,A1N,…,AMNT.The relationship between the two sets of variablesp and a in equation (8) is given by means of the eight elastic boundary conditions represented by equations (7a), (7b), (7c), (7d), (7e), (7f), (7g), and (7h). Substituting the displacement expression equation (5) into the governing equation of the thin plate equation (1) leads to
(10)∑m=0∞∑n=0∞λam4+λbn4+2λam2λbn2Amncosλamxcosλbny+∑l=14∑m=0∞ξblyλam4−2∂2ξbly∂y2λam2+∂4ξbly∂y4cmlcosλamx+∑n=0∞ξalxλbn4−2∂2ξalx∂x2λbn2+∂4ξalx∂x4dnlcosλbny−ρhω2D∑m=0∞∑n=0∞Amncosλamxcosλbny+∑l=14ξbly∑m=0∞cmlcosλamx+ξalx∑n=0∞dnlcosλbny=0.Similarly, takingm=0,1,…,M, n=0,1,…,N, after simplifying, equation (10) can be written as the following matrix equation:
(11)K~a+Bp−ρhω2DM~a+Fp=0.The matricesK~, M~, B, and F are defined in the appendix.According to the relationship betweenp and a which has been defined by equation (8), equation (11) can be written as
(12)K−ρhω2DMa=0,where K=K~+BH−1Q and Μ=M~+FH−1Q.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.


## 3. Analysis of Aeroelastic Response

### 3.1. Aerodynamic Model

The aerodynamic model is based on the second-order piston theory which can be expressed as follows:(13)pp∞=1+κvza∞+κκ+14vza∞2,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 vz of the panel and the local pressure p.Given that the partial derivative of transverse displacementswx,y,t to x is equal to the velocity of flow and the partial derivative to y is 0, the lateral velocity vz of the panel can be expressed as follows:
(14)vz=∂wx,y,t∂t=∂w∂x∂x∂t+∂w∂y∂y∂t+∂w∂t=∂w∂xU∞+∂w∂t,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
(15)Δp=−2κp∞a∞∂w∂xU∞+∂w∂t.


### 3.2. Aeroelastic Response in the Time Domain

After adding the aerodynamic force, the vibration equation of the rectangular plate should be rewritten as(16)ρh∂2wx,y,t∂t2+D∇4wx,y,t=Δp.By separating the variables of the transverse displacement and cut atNth modes, the form of the solution is assumed as follows:
(17)wx,y,t=∑i=1NWix,yqit.Substituting equations (15) and (17) into equation (16) leads to
(18)ρh∑i=1NWix,yq¨it+D∑i=1N∇4Wix,yqit=−2κp∞a∞U∞∑i=1N∂Wix,y∂xqit+∑i=1NWix,yq̇it,where q¨it represents the second derivative of qit versus time, and Wix,y represents the ith vibration mode function of the panel.According to the modal orthogonality, after multiplying both ends of equation (18) by Wi and integrating the new equation, the ith order equation can be obtained
(19)Miq¨i+Kiqi=−4qdMa∞∑j=1Nqj⋅∬∂Wj∂xWidxdy+q̇jU∞∬WjWidxdy.All ofN equations can be rewritten as a form of state-space equation as follows:
(20a)q̇1⋮q̇Nq¨1⋮q¨N=0IA21A22q1⋮qNq̇1⋮q̇N,(20b)A21=−4qd/Ma∞∬∂W1/∂xW1dxdy−K1M1⋯−4qd/Ma∞∬∂WN/∂xW1dxdyM1⋮⋱⋮−4qd/Ma∞∬∂W1/∂xWNdxdyMN⋯−4qd/Ma∞∬∂WN/∂xWNdxdy−KNMN,(20c)A22=−4qdMa∞U∞∬W1W1dxdy⋯−4qdMa∞U∞∬WNW1dxdy⋮⋱⋮−4qdMa∞U∞∬WNW1dxdy⋯−4qdMa∞U∞∬WNWNdxdy.By the MFSM described before, the solution for the mode function can be obtained(21)Wx,y=∑m=0∞∑n=0∞Amncosλamxcosλbny+∑l=14ξbly∑m=0∞cmlcosλamx+ξalx∑n=0∞dnlcosλbny.Without loss of generality,Wx,y is used to represent Wix,y in equation (21). When given a and p, Wx,y can be reduced to the corresponding Wix,y. Substituting equation (21) into equations (20a), (20b), and (20c) results in q1,q2,…,qN, then the aeroelastic response of the plate in the time domain can be obtained.


## 3.1. Aerodynamic Model

The aerodynamic model is based on the second-order piston theory which can be expressed as follows:(13)pp∞=1+κvza∞+κκ+14vza∞2,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 vz of the panel and the local pressure p.Given that the partial derivative of transverse displacementswx,y,t to x is equal to the velocity of flow and the partial derivative to y is 0, the lateral velocity vz of the panel can be expressed as follows:
(14)vz=∂wx,y,t∂t=∂w∂x∂x∂t+∂w∂y∂y∂t+∂w∂t=∂w∂xU∞+∂w∂t,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
(15)Δp=−2κp∞a∞∂w∂xU∞+∂w∂t.


## 3.2. Aeroelastic Response in the Time Domain

After adding the aerodynamic force, the vibration equation of the rectangular plate should be rewritten as(16)ρh∂2wx,y,t∂t2+D∇4wx,y,t=Δp.By separating the variables of the transverse displacement and cut atNth modes, the form of the solution is assumed as follows:
(17)wx,y,t=∑i=1NWix,yqit.Substituting equations (15) and (17) into equation (16) leads to
(18)ρh∑i=1NWix,yq¨it+D∑i=1N∇4Wix,yqit=−2κp∞a∞U∞∑i=1N∂Wix,y∂xqit+∑i=1NWix,yq̇it,where q¨it represents the second derivative of qit versus time, and Wix,y represents the ith vibration mode function of the panel.According to the modal orthogonality, after multiplying both ends of equation (18) by Wi and integrating the new equation, the ith order equation can be obtained
(19)Miq¨i+Kiqi=−4qdMa∞∑j=1Nqj⋅∬∂Wj∂xWidxdy+q̇jU∞∬WjWidxdy.All ofN equations can be rewritten as a form of state-space equation as follows:
(20a)q̇1⋮q̇Nq¨1⋮q¨N=0IA21A22q1⋮qNq̇1⋮q̇N,(20b)A21=−4qd/Ma∞∬∂W1/∂xW1dxdy−K1M1⋯−4qd/Ma∞∬∂WN/∂xW1dxdyM1⋮⋱⋮−4qd/Ma∞∬∂W1/∂xWNdxdyMN⋯−4qd/Ma∞∬∂WN/∂xWNdxdy−KNMN,(20c)A22=−4qdMa∞U∞∬W1W1dxdy⋯−4qdMa∞U∞∬WNW1dxdy⋮⋱⋮−4qdMa∞U∞∬WNW1dxdy⋯−4qdMa∞U∞∬WNWNdxdy.By the MFSM described before, the solution for the mode function can be obtained(21)Wx,y=∑m=0∞∑n=0∞Amncosλamxcosλbny+∑l=14ξbly∑m=0∞cmlcosλamx+ξalx∑n=0∞dnlcosλbny.Without loss of generality,Wx,y is used to represent Wix,y in equation (21). When given a and p, Wx,y can be reduced to the corresponding Wix,y. Substituting equation (21) into equations (20a), (20b), and (20c) results in q1,q2,…,qN, then the aeroelastic response of the plate in the time domain can be obtained.


## 4. 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 Table1.Table 1
Physical parameters of the panel.
Parameter
Meaning
Case 1
Case 2
a
Panel length inx-direction
1 m
0.4 m
b
Panel width iny-direction
1, 1.5, 2, 2.5, 3 m
0.5 m
h
Panel thickness
0.002 m
0.0015 m
ρ
Mass density
2700 kg/m^3
2700 kg/m^3
E
Young’s modulus
69 GPa
69 GPa
μ
Poisson’s ratio
0.33
0.33The 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 as1.0×108 N/m, while the rotational springs set to be 0.Table2 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.Table 2
Frequencies for S-S-S-S square panel with different truncate numbers.
M=N
The first 6 natural frequencies (Hz)
1
2
3
4
5
6
5
9.66
24.23
24.23
38.73
48.78
48.99
10
9.70
24.27
24.27
38.81
48.56
48.59
15
9.71
24.28
24.28
38.83
48.55
48.57
16
9.71
24.28
24.28
38.83
48.55
48.57
17
9.71
24.28
24.28
38.83
48.55
48.57
18
9.71
24.28
24.28
38.83
48.55
48.57
19
9.71
24.28
24.28
38.84
48.55
48.57
20
9.71
24.28
24.28
38.84
48.55
48.56
25
9.71
24.28
24.28
38.84
48.55
48.56
30
9.71
24.28
24.28
38.84
48.55
48.56In Table3, 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 MFSM are compared with results from FEM and Navier’s analytical solutions. A unified grid cell of size 0.01m×0.01m 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.Table 3
Frequencies for S-S-S-S panel with different aspect ratios.
r=b/a
The first 6 natural frequencies (Hz)
1
2
3
4
5
6
1.0
9.710MFSM
24.277
24.277
38.836
48.554
48.566
9.672FEM
24.147
24.147
38.347
48.290
48.290
9.713NAS
24.283
24.283
38.853
48.567
48.567
1.5
7.013MFSM
13.487
21.582
24.280
28.053
38.840
6.996FEM
13.424
21.523
24.148
27.825
38.348
7.015NAS
13.491
21.585
24.283
28.061
38.853
2.0
6.070MFSM
9.711
15.782
20.639
24.279
24.282
6.060FEM
9.674
15.703
20.606
24.148
24.148
6.071NAS
9.713
15.784
20.641
24.283
24.283
2.5
5.633MFSM
7.964
11.849
17.289
20.203
22.533
5.627FEM
7.939
11.796
17.199
20.181
22.447
5.634NAS
7.965
11.850
17.290
20.204
22.535
3.0
5.396MFSM
7.015
9.713
13.491
18.351
19.966
5.391FEM
6.996
9.674
13.425
18.250
19.950
5.396NAS
7.015
9.713
13.491
18.347
19.966
MFSM for results from modified Fourier series method. FEM for results from finite element method. NAS for results from Navier’s analytical solutionωmn=π/2m2/a2+n2/b2D/ρh.Figure 2
The structural finite element mesh with unified grid cells.Figure3 shows the 1st, 2nd, 3rd, and 5th natural frequencies of the panel with four clamped edges (C-C-C-C). The 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×108 N/m and 1.0×108 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.Figure 3
The 1st, 2nd, 3rd, and 5th frequencies of the C-C-C-C panel with different widths.Figure 4
The 1st, 2nd, 3rd, and 5th frequencies of the panel with different widths.
(a)
S-F-S-F panel
(b)
C-F-C-F panel


### 4.2. 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 Table1.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/m3, 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.05m×0.05m is adopted in this paper, shown in Figure 7.Figure 5
Vibration response of the S-S-S-S panel.
(a)
U∞=850m/s
(b)
U∞=868m/s
(c)
U∞=870m/s
(d)
U∞=890m/sFigure 6
V-g figure of the panel (867.5 m/s).Figure 7
The aerodynamic mesh of the lifting surface with unified aerodynamic grid cells.Figure8 presents the flutter speeds of the panels with one side keeps 1.0 m and the thickness are h=0.003 m and h=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 8
Flutter speeds of panels in different widths with different thicknesses.
(a)
h=0.003 m
(b)
h=0.0035 mFigure9 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 9
Flutter speeds of panels in different thicknesses with widthb=1.3 m.Figure10 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×1010 N/m, while the rotational springs Kx0, Kxa, Ky0, and Kyb change from 0 to 1.0×1010 N/rad, shown in Table 4. When the stiffness of torsional springs reaches 1.0×106 N/rad, the boundary condition can be considered as clamped, and the flutter speed results equal to 1668 m/s.Figure 10
Flutter speed of different boundary conditions.Table 4
Boundary conditions.
Case
1
2
3
4
5
6
7
8
Linear springs (N/m)kx0,kxa,ky0,kyb
1010
1010
1010
1010
1010
1010
1010
1010
Rotational springs (N/rad)Kx0,Kxa,Ky0,Kyb
10
102
103
104
105
106
108
1010


## 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 Table1.Table 1
Physical parameters of the panel.
Parameter
Meaning
Case 1
Case 2
a
Panel length inx-direction
1 m
0.4 m
b
Panel width iny-direction
1, 1.5, 2, 2.5, 3 m
0.5 m
h
Panel thickness
0.002 m
0.0015 m
ρ
Mass density
2700 kg/m^3
2700 kg/m^3
E
Young’s modulus
69 GPa
69 GPa
μ
Poisson’s ratio
0.33
0.33The 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 as1.0×108 N/m, while the rotational springs set to be 0.Table2 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.Table 2
Frequencies for S-S-S-S square panel with different truncate numbers.
M=N
The first 6 natural frequencies (Hz)
1
2
3
4
5
6
5
9.66
24.23
24.23
38.73
48.78
48.99
10
9.70
24.27
24.27
38.81
48.56
48.59
15
9.71
24.28
24.28
38.83
48.55
48.57
16
9.71
24.28
24.28
38.83
48.55
48.57
17
9.71
24.28
24.28
38.83
48.55
48.57
18
9.71
24.28
24.28
38.83
48.55
48.57
19
9.71
24.28
24.28
38.84
48.55
48.57
20
9.71
24.28
24.28
38.84
48.55
48.56
25
9.71
24.28
24.28
38.84
48.55
48.56
30
9.71
24.28
24.28
38.84
48.55
48.56In Table3, 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 MFSM are compared with results from FEM and Navier’s analytical solutions. A unified grid cell of size 0.01m×0.01m 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.Table 3
Frequencies for S-S-S-S panel with different aspect ratios.
r=b/a
The first 6 natural frequencies (Hz)
1
2
3
4
5
6
1.0
9.710MFSM
24.277
24.277
38.836
48.554
48.566
9.672FEM
24.147
24.147
38.347
48.290
48.290
9.713NAS
24.283
24.283
38.853
48.567
48.567
1.5
7.013MFSM
13.487
21.582
24.280
28.053
38.840
6.996FEM
13.424
21.523
24.148
27.825
38.348
7.015NAS
13.491
21.585
24.283
28.061
38.853
2.0
6.070MFSM
9.711
15.782
20.639
24.279
24.282
6.060FEM
9.674
15.703
20.606
24.148
24.148
6.071NAS
9.713
15.784
20.641
24.283
24.283
2.5
5.633MFSM
7.964
11.849
17.289
20.203
22.533
5.627FEM
7.939
11.796
17.199
20.181
22.447
5.634NAS
7.965
11.850
17.290
20.204
22.535
3.0
5.396MFSM
7.015
9.713
13.491
18.351
19.966
5.391FEM
6.996
9.674
13.425
18.250
19.950
5.396NAS
7.015
9.713
13.491
18.347
19.966
MFSM for results from modified Fourier series method. FEM for results from finite element method. NAS for results from Navier’s analytical solutionωmn=π/2m2/a2+n2/b2D/ρh.Figure 2
The structural finite element mesh with unified grid cells.Figure3 shows the 1st, 2nd, 3rd, and 5th natural frequencies of the panel with four clamped edges (C-C-C-C). The 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×108 N/m and 1.0×108 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.Figure 3
The 1st, 2nd, 3rd, and 5th frequencies of the C-C-C-C panel with different widths.Figure 4
The 1st, 2nd, 3rd, and 5th frequencies of the panel with different widths.
(a)
S-F-S-F panel
(b)
C-F-C-F panel


## 4.2. 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 Table1.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/m3, 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.05m×0.05m is adopted in this paper, shown in Figure 7.Figure 5
Vibration response of the S-S-S-S panel.
(a)
U∞=850m/s
(b)
U∞=868m/s
(c)
U∞=870m/s
(d)
U∞=890m/sFigure 6
V-g figure of the panel (867.5 m/s).Figure 7
The aerodynamic mesh of the lifting surface with unified aerodynamic grid cells.Figure8 presents the flutter speeds of the panels with one side keeps 1.0 m and the thickness are h=0.003 m and h=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 8
Flutter speeds of panels in different widths with different thicknesses.
(a)
h=0.003 m
(b)
h=0.0035 mFigure9 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 9
Flutter speeds of panels in different thicknesses with widthb=1.3 m.Figure10 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×1010 N/m, while the rotational springs Kx0, Kxa, Ky0, and Kyb change from 0 to 1.0×1010 N/rad, shown in Table 4. When the stiffness of torsional springs reaches 1.0×106 N/rad, the boundary condition can be considered as clamped, and the flutter speed results equal to 1668 m/s.Figure 10
Flutter speed of different boundary conditions.Table 4
Boundary conditions.
Case
1
2
3
4
5
6
7
8
Linear springs (N/m)kx0,kxa,ky0,kyb
1010
1010
1010
1010
1010
1010
1010
1010
Rotational springs (N/rad)Kx0,Kxa,Ky0,Kyb
10
102
103
104
105
106
108
1010


## 5. 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.



---

*Source: 5164026-2019-03-24.xml*