1 Introduction

Military aircrafts usually carry external missiles, fuel tanks, and other objects in the fuselage or wings. The flutter problem of wings and plates is an important part in the process of the aeroelastic design of a aircraft. At present, the location of a store for aircrafts is mainly distributed at two positions, wings and belly parts. The trends of aircraft design are leading to the placement of heavy lumped mass on the panels of the wings or the fuselage. The installation of the store tends to change the vibration characteristics and flutter boundary of the structures. In aircraft design, the plate of the belly is usually connected by riveting and regarded as the simply supported boundary, and the wing or control surface can be regarded as a cantilever plate.

For the research about the wings with an external store, e.g., masses, missile, and tank, Tang et al.^[1] studied a combination system of a delta plate with a store, where the flutter velocity was sensitive to the store span location and store attachment stiffness. Fazelzadeh et al.^[2] investigated the static and dynamic aeroelastic stability of a cantilever wing/store system. Mazidi et al.^[3] considered the aeroelastic response of a wing with an engine subjected to time-dependent thrust excitations. Firouz-Abadi et al.^[4] studied the aeroelastic instability of a wing with a fuel tank store. Lei and Ye^[5] simulated the process of store separation with the wing deformation numerically. Gern^[6] and Golparvar et al.^[7] investigated the vibration characteristics of a wing with a store and the effects of different parameters such as the weights, the inertia of stores, and the thickness and store center of gravity on the flutter. Runyan and Sewall^[8] showed the effect of the concentrated weights in a wing model at various positions on the flutter characteristics. The aforementioned studies show the effects of the store on the flutter behavior, but they often focus on the subsonic aerodynamic stability of the wing plate. With the increase in the speed of an aircraft, it is necessary to pursue an approach to study the flutter problem of both the wings and the plates in a supersonic flow. Since more and more composite laminated structures are used in real aircrafts, it is important to carry out research on the flutter characteristics of composite laminated plates with external stores in the supersonic regions.

For the flutter problem of plates and shells, theoretical results have provided a good understanding. However, the studies about the flutter characteristics of plates with complex boundaries or elastic supports are rare. Dowell and Weiliang^[9] investigated the flutter problem of a cantilever isotropic plate with a Rayleigh-Ritz approach, and showed that the length-to-width ratio of the plate had a significant effect on the flutter boundary. Kouchakzadeh et al.^[10] studied the aeroelasticity problem of the laminated composite panel in a supersonic airflow with a simply supported boundary with the Galerkin method. Pourtakdoust and Fazelzadeh^[11] investigated the chaotic behavior of the nonlinear viscoelastic plate in a supersonic region with a Galerkin approach. Singha and Ganapathi^[12], Kuo^[13], and Marques et al.^[14] investigated the flutter behavior of composite plates in the supersonic flow with a linear piston theory through the finite element method (FEM). However, there are difficulties in the mode shape construction for the structure with a complex boundary due to the lack of a precise analytical mode function.

For complex boundaries, Bhat^[15] used the Gram-Schmidt process to generate the orthogonal polynomials so as to solve the natural frequency of the rectangular plates with various boundaries, and obtained superior results for lower modes. Amabili and Garziera^[16] used the artificial springs to investigate the vibration problems of the cylindrical shells with complex boundary constraints, where the mode of the less-restrained condition was used to expand the solutions of displacement. The distributions of spring stiffness were systematically represented by the cosine series. Jin et al.^[17] studied the free vibration of composite laminated shells under general elastic boundary with the Rayleigh-Ritz procedure, expressed the displacements by a Fourier cosine series, and illustrated the effects of boundary restraining stiffness on frequency. Mahi et al.^[18] investigated the free vibration of different kinds of composite panels with different boundaries by use of a series of characteristic orthogonal polynomials. Zhang and Li^[19] and Wang et al.^[20] solved the vibration problem of rectangular plates with elastic edge restraints, where the translational and rotational restraints were arbitrarily applied to an edge and their stiffness distributions were described by a set of cosine functions. Li and Narita^[21] studied the critical aerodynamic pressure for different types of laminated plates with general boundary conditions with the consideration of the thermal effects. The modes of plates or shells with elastic support are rarely used to solve the complicated problems such as the analysis of response or flutter. In this paper, artificial springs are used to solve the flutter boundary of the simply supported and cantilever plates, respectively.

We use the Dirac function to describe the interaction between the plate and store, and establish a dynamic model of composite laminated plate with a store with the Hamilton principle. The solution procedure is divided into two steps. Firstly, the weight coefficients of the mode functions and frequencies of composite laminated plate are obtained by use of the Rayleigh-Ritz method, and the boundary constraint is modeled as artificial tension and torsion springs. The simulation of different boundaries can be realized by setting different stiffness values of the artificial springs. According to the deformation coordination relation between the plate and the store, the mode shape of the plate sub-system is brought into the motion equation of the store subsystem with the consideration of pitching and plunging, and the expressions of the interaction force between the plate and the external store are derived. Secondly, the mode shape and generalized coordinates are taken into the Hamilton equation to obtain the governing equations. The convergence of the number of the orthogonal polynomials on the frequency is studied. The validity of the proposed approach is given by comparing the frequency and mode shape with the FEM results. The effects of the ply angle of composite plates on the frequency and flutter boundary are highlighted. Three parameters of the external store, i.e., the mounting position, the center of gravity position, and the two mounting points spacing, are chosen to discuss the flutter boundary.

2 Theoretical formulation

A schematic of the rectangle-plate/store geometry is shown in Fig. 1, where P₁ and P₂ are attachment points between the plate and the store. There are a global coordinate system o-xyz, and a local coordinate system O-ξηζ. Here, e₁ denotes the distance between the mounting points, and e₂ denotes the distance between the mounting position and the edge of the plate for the simply supported structure under the belly parts of aircraft (see Fig. 1(a)) or for a cantilever supported plate (wing) (see Fig. 1(b)). ξ₁ and ξ₂ are used to determine the center of gravity (ξ₁/ξ₂). The airflow is along the x-direction. Two cases are considered, the first one is simply supported boundary, and the second one is the cantilever boundary. The installation layout between the plate and the store is shown in Fig. 1(c). The form of composite plates is shown in Fig. 1(d).

2.1 Energy equation of the composite laminated plate

According to the Kirchhoff theory, the displacement of the laminated plate is given as follows:

(1)

where u, v, and w are the displacements of the plate in the x-, y-, and z-directions, respectively, and w₀(x, y, t) is the transverse displacement of the middle plane of the plate.

The strain-displacement relations can be expressed as follows:

(2)

where ε is the strain vector of the arbitrary point, and κ is the variation of the curvature vector in the middle plane expressed by

(3)

The constitutive equation of an orthotropic lamina is

(4)

where σ₁, σ₂, and τ₁₂ are the normal and shear stresses, and ε₁, ε₂, and γ₁₂ are the strains. Q₁₁=E₁/(1-υ₁₂υ₂₁), Q₂₂=E₂/(1-υ₁₂υ₂₁), Q₁₂=υ₂₁E₂/(1-υ₁₂υ₂₁), and Q₆₆=G₁₂ are the stiffness coefficients, in which E₁, E₂, and G₁₂ are the elastic and shear moduli, and υ₁₂ and υ₂₁ are Poisson's ratios.

The constitutive equation of the kth lamina of the laminated plate is

(5)

where σ_x, σ_y, and τ_xy are the stresses in the x-and y-directions and the xy-plane, respectively, and Q_k is the transformed reduced stiffness matrix of the kth lamina expressed by

(6)

In the above equation, T_k is the transformation matrix expressed by

(7)

where θ_k is the ply angle.

The membrane stress resultant vector N and the bending stress resultant vector M can be written as follows:

(8)

(9)

where n is the number of the layers, z_k and z_k-1 are the coordinates of the upper and lower surfaces of the kth lamina in the z-direction, respectively, and B_ij and D_ij are defined by

(10)

The kinetic energy of the composite laminated plate is

(11)

where "." represents the first derivative of time, ρ is the density of the laminated plate, and V is the volume of the laminated plate. The laminated plate is symmetrically paved, and thus B_ij=0. The strain energy due to the vibration stresses is

(12)

where A is the area of the plate.

2.2 Simulation of the artificial springs

We use the artificial springs to describe the connections between the plate and the supporting members. The artificial springs are evenly distributed on the edges of the plate, the constraints are realized by the tension springs and torsion springs, and the corresponding stiffness are k^w and k^ϑ (see Fig. 2). We set different stiffness values. Then, the artificial springs can be used to describe arbitrary boundary conditions. Because the actual plate structure may be bolted, riveted, or elastically supported, the aim of using the artificial springs is to simulate not only the classical boundary but also the elastic boundary.

Since the artificial springs are assumed to be massless, only the potential energy is considered, which can be expressed as follows:

(13)

where k_0(a)^w, k_a(a)^w, k_0(b)^w, k_b(b)^w and k_0(a)^ϑ, k_a(a)^ϑ, k_0(b)^ϑ, k_b(b)^ϑ are the stiffness values of the tension spring and the torsion spring at x=0, x=a, y=0, and y=b, respectively, the letter in parentheses (a) represents the x-axis, and the letter in parentheses (b) represents the y-axis.

2.3 Rayleigh-Ritz method

The Rayleigh-Ritz method is used to obtain the weight coefficients of the mode function and frequency. The transverse displacement of the plate can be expressed as follows:

(14)

where ω represents the circular frequency, and W(x, y) is the mode shape that can be expressed as the Gram-Schmidt orthogonal polynomials, i.e.,

(15)

in which a_mn is the unknown coefficient, ϕ_m (x) and φ_n (y) are, respectively, the orthogonal polynomials in the x-and y-directions, m_t and n_t are the truncation numbers of the polynomials.

Substituting Eqs. (14) and (15) into Eqs. (11)-(13) yields the Rayleigh quotient to the corresponding factor as follows:

(16)

The minimizing condition with the energy function yields the generalized eigenvalue problem as follows:

(17)

where K_ε, K_s, and M are m_t×n_t matrices, and a=(a₁₁, a₁₂, ..., a_mtnt)^T. The matrices K_ε, K_s, and M are

where

The construction of the Gram-Schmidt orthogonal polynomials is as follows. The first term of the orthogonal polynomials cluster φ₀ (ξ) is given so that the certain boundary condition can be satisfied. According to the steps of the orthogonal polynomial cluster configuration, the remaining items can be obtained as follows:

(18)

(19)

where

Each polynomial is normalized, i.e.,

(20)

If the first term satisfies the geometry boundary conditions, it is easy to verify that other terms will also satisfy the geometry boundary conditions. The relationship between any two orthogonal polynomials is

(21)

The first item of the orthogonal polynomials with different boundary conditions is shown in Table 1. The letter F represents freedom, C represents clamped, and S represents simply supported.

2.4 Motion equation of the store

The displacement matching relationship between the plate and the store is shown in Fig. 3. The store is regarded as a rigid body, including two degrees of freedom, i.e., plunging and pitching. The motion equation of the store is

(22)

where ξ₁ and ξ₂ are the coordinates of the attachment point for the center of the gravity position in the O-ξηζ plane, θ is the pitch angle, m is the mass of the store, J is the moment of inertia, F is the equivalent aerodynamic force of the store acting on the center of gravity, M_ΔP is the equivalent aerodynamic force moment, and F₁ and F₂ are the interaction forces between the plate and the store.

According to the displacement matching relationship between the plate and the store, the transverse displacements w₁ and w₂ are obtained as follows:

(23)

Combining Eq. (22) and Eq. (23), we have

(24)

2.5 Aerodynamic load

The piston theory proposed by Lighthill and Ashley is the most widely used aerodynamic theory in supersonic region. The aerodynamic pressure caused by the supersonic flow can be described by the first-order quasi-stable piston theory as follows:

(25)

where q=ρ_av_∞²/2 is the dynamic pressure, ρ_a is the incoming density, v_∞ is the airflow velocity, and Ma_∞ is the Mach number.

2.6 Governing equations

With the mode shapes W(x, y) in Subsection 2.2 to obtain the governing equation by use of the assumed mode method, we can express the displacement of the plate as follows:

(26)

where g_ij (t) is the generalized coordinate, W_ij is the corresponding modal shapes obtained from the Rayleigh-Ritz approach which satisfies the geometry boundary conditions, m and n are the modal orders, and

The transverse displacements of the mount points on the composite plate are

(27)

where (x₁, y₁) and (x₂, y₂) are the position coordinates of the mount points in the Global coordinate system o-xyz, and y₁=y₂.

It is acceptable to calculate the flutter boundary by neglecting the aerodynamics pressure on the store^[22]. Therefore, the effect of the store aerodynamic load is ignored. Only the mass and moment of inertia are accounted for in the modeling of the store. The concentrated force in Eq. (24) can be written as follows:

(28)

The governing equations of the composite plate with an external store is established by the Hamilton principle. The expression is

(29)

where the virtual work of the kinetic energy δT and the virtual work of the potential energy δU are

(30)

(31)

The virtual work of the external load δV is

(32)

Substituting Eqs. (30)-(32) into Eq. (29) yields the governing equation as follows:

(33)

where

3 Solving method

Substituting the vibration displacement of the plate r(t)=r₀e^λt into Eq. (33) yields

(34)

The condition of Eq. (34) having nonzero solutions is

(35)

Then, the complex eigenvalues can be obtained as follows:

(36)

where the frequency and the damping ratio of the system are

(37)

The flutter is a kind of dynamic instability phenomena. The system equilibrium position loses its stability when its flight velocity reaches the critical flutter velocity. The critical velocity is a state parameter between convergence and instability. With the increase in the flow velocity, when the real part of the eigenvalue values is from negative to positive or the damping ratio is from positive to negative, zero equilibrium point of the system becomes unstable. One can find out a linear critical velocity by inspecting the maximal real part of the eigenvalues. Further calculation finds that, without the consideration of the damping, the critical flutter speed coincides with the velocity when two consecutive low-level frequencies of the system are equal. Neglecting the aerodynamic damping, the frequency coincidence theory is taken as the basis for the occurrence of the flutter in this paper.

4 Numerical results

A composite laminated plate is taken as a numerical example to demonstrate the solution procedure. Referring to Tang et al.^[1], we choose the parameters of the plate and store as listed in Table 2, where a, b, and h are the geometric parameters of the simply supported plate for the belly structure, and a', b' and h' are the geometric parameters of the cantilever plate for the wing. In this paper, under the wing condition, the spring stiffness value of the constraint edge is set to be infinity so as to simulate the cantilever plate. For simply supported plate, the stiffness value of the tension spring is set to be infinity, and the stiffness value of the torsion spring is set to be zero.

4.1 Vibration characteristics analysis

In order to ensure the validity of the calculation, the convergence analysis about the orthogonal polynomial orders on the accuracy of the frequency and mode is discussed. The FEM results of ANSYS are regarded as a reference value. The error is determined as follows:

(38)

Taking the symmetric ply [0°/90°/90°/0°] plates with four-side fixed boundary as an example, the digitals 3×3, 4×4, 5×5, and 6×6 represent the truncation orders of the orthogonal polynomials in the x-and y-directions. With the increase in the number of the polynomials, the frequency converges to the FEM result in Table 3. When six orders of polynomials are selected in both the x-and the y-directions, the error between the FEM and Rayleigh-Ritz method is no more than 1%, and sufficient precision can be satisfied.

The natural frequencies of the [0°/90°/90°/0°] composite laminated plate with the classical boundary are compared with the FEM results in Table 4, where the letter F represents freedom, C is for clamped, and S is for simply supported. The stiffness value of the artificial springs for the fixed support is set to be infinity, whether it is a torsion spring or a tension spring. The number of the orthogonal polynomials is eight in both the x-and the y-directions, and the maximum error is less than 2%, which fully satisfies the engineering requirements. With different stiffness values, the artificial springs can be used to describe all the classical boundaries.

The mode shapes obtained from the FEM and Rayleigh-Ritz method for the cantilever plate are shown in Fig. 4, where the upper corresponds to the finite element results and the lower is for the present method. The results show that the mode shapes are in good agreement with the FEM. This ensures the accuracy of the proposed approach with the artificial springs to simulate the boundary constraints.

In order to continuously demonstrate the effect of spring stiffness on the natural frequency, (Condition 1) the tension spring value of k^w is fixed, i.e., k^w=10⁹ N/m, and the torsion stiffness value continuously changes, which results in the change from the simple supported to the fixed boundary, and the corresponding first two frequencies are shown in Fig. 5. When the stiffness value of each spring is 10⁹, the natural frequency of the plate converges to a fixed value, which is the frequency of the plate with fixed boundary, and the torsion spring in the x-direction has a more significant effect on the first frequency (comparing the convergent frequency when the torsion spring stiffness is 10⁹ N·m/rad in the x-and y-directions) because it can provide a larger constraint stiffness. In Condition 2, the value of k^ϑ is fixed to be k^ϑ=10⁹ N·m/rad. When k^w is small, the first natural frequency of the plate is close to zero. The tension spring in the x-direction has a more significant effect on the frequency because of a longer distance in the length direction (see Fig. 6). When the stiffness of the artificial spring is greater than 10⁷, the frequency does not change with the increase in the spring stiffness value.

4.2 Effects of composite layers and different boundaries on the flutter boundary

The laminated composite plate in the ply condition [90°/0°/90°/0°]_s investigated by Zhao^[23] is adopted to compare the flutter boundary (see Table 5). The results show that the proposed method and formula derivation for the flutter boundary by use of the artificial springs are valid.

Two kinds of ply conditions, i.e., symmetrical paving [0°/θ/θ/0°] and asymmetrical paving [0°/θ/0°/θ], are chosen to study the effect of the composite layer on the frequency and flutter boundary, where θ changes from -90° to 90° (see Fig. 7). For the symmetric ply condition, when θ is positive, the natural frequency gradually increases with the increase in the ply angle, and the first frequency changes very small, but the second frequency obviously increases. When the ply angle is close to 90°, the growth rate slows down. For the asymmetrical ply condition, with the increase in the ply angle, the natural frequency increases and there is a peak at about 50°. The ply angle and ply arrangement have significant effects on the frequency of the plate, which will change the critical flutter boundary of the structure.

The flutter boundary without taking into account the store effect is shown in Fig. 8, where the left coordinate is frequency and the right coordinate is the critical flutter velocity. For the positive angle, the critical flutter velocity gradually becomes small with the increase in the ply angle, and there is a maximum point at 0 degree due to the maximum stiffness, which indicates that the laminated plates with smaller ply angle have better aerodynamic stability. The premise is to ensure the sufficient strength. Figure 9 shows that, with the increase in the flow velocity, the second-order mode and the third-order mode coincide at the point v=342 m·s^-1, and then separate. Until the first-mode and the third-mode coincide, there is no separation. The mode jumping phenomenon appears at v=342 m·s^-1, because of the accuracy problem of the numerical calculation, it is difficult to observe the phenomenon from the frequency chart. In order to better display the phenomenon, the mode shapes at v=330 m·s^-1 (1, 2, 3, and 4) and v=350 m·s^-1 (Ⅰ , Ⅱ, Ⅲ, and Ⅳ) considering the aerodynamic loads are shown in Fig. 10. It is found that, the second mode and the third mode exchange, and a mode jumping can be clearly found, which shows that one single coincidence point cannot be used as the basis of flutter occurrence. Under the action of aerodynamic loads, the mode shape of the plate has changed significantly.

The flutter boundary under different boundary conditions without store is studied in Fig. 11. With the increase in the flow velocity, the first two frequencies of the plate are gradually approaching, and finally a curve is formed when the flutter occurs. It is found that the flutter velocity of the fixed plate is obviously higher than that of the simply supported plate, which indicates that appropriately enhancing the boundary constraints through increasing the stiffness value of artificial spring can improve the critical flutter velocity. The mode jumping phenomenon can be observed in this condition.

For the plate with store, the distance between the fore and aft attachment point is e₁=0.2 m, the locations in the plate are (0.05, 0.2) and (0.25, 0.2), and |ξ₁/ξ₂|=1/3 for the simply supported boundary in Fig. 12. With the increase in the flow velocity, the first two frequencies of the structure coincide at v=416 m·s^-1, and flutter occurs. The flutter velocity with store is 416 m·s^-1, which is higher than 382 m·s^-1 without store. The adding of store increases the flutter velocity of the system. Therefore, the store can also be used as a method of flutter control for the simply supported boundary.

4.3 Effects of the external store on the flutter boundary

In this section, the ply condition [0°/90°/90°/0°] is chosen to study the effect of the parameters of the store on the flutter characteristics. The effects of the store, e.g., the mounting position, the center of gravity position, and the two mounting points spacing on the flutter are discussed.

The frequency of the cantilever plate with store is shown in Fig. 13, where e₂ is the mounting distance between the store and the fixed edge. The first frequency changes little, but the second frequency increases with the increase in e₂. The effect of the mounting position on the flutter boundary is shown in Fig. 14, and e₂ changes from 0.1 m to 0.5 m (the length of the cantilever plate is 0.6 m, 0 < e₂ < 0.6). The frequency coalescence between the first mode and the second mode leads to flutter instability. When e₂ increases, the flutter velocity becomes smaller. The flutter velocity of the cantilever plate without store is 901 m·s^-1, which is greater than 884 m·s^-1 for the case e₂=0.1 m·s^-1. Adding a store reduces the critical flutter velocity, which is contrary to the simply supported boundary. When the store is near the constraint edge, the system can obtain a better aerodynamic stability.

Because the simply supported plate is symmetrical, the mounting position e₂ changes from 0.05 m to 0.2 m (the width of the simply supported plate is 0.4 m, 0 < e₂ < 0.2), and only one side is selected. In Fig. 15, the results show that the structure has a high flutter velocity near the center, but it does not continuously increase with the increase in the mounting position, and the peak will appear in the middle e₂=0.1 m.

The center of the gravity position |ξ₁/ξ₂| changes from 0.15 to 5 (ξ₁ +ξ₂=0.2 m) (see Fig. 16). The change of the center of gravity will affect the flutter velocity of the structure, and the flutter occurs in the coincidence of the second and third order modes rather than the first two modes when the center of gravity is backward, which shows that the critical flutter velocity has a jump. When the flow velocity continues to increase, from |ξ₁/ξ₂| > 0.66, a complex frequency coincidence phenomenon takes place, and several coincidence curves appear. At this time, the flutter has already occurred. When |ξ₁/ξ₂| is greater than a value, with the center of gravity position backing, the flutter velocity gradually becomes small, but the effect is not significant. When the center of the gravity position is located at 1/3 of the distance between two attaching points, a larger flutter velocity can be obtained.

The two mounting point spacing of the store e₁ changes from 0.06 m to 0.27 m (see Fig. 17). With the increase in the flow velocity, the flutter velocity is not positively correlated with the two mounting point spacing. When the two mounting point spacing is small, e.g., e₁=0.15 m, which is about 1/3 of the length of the plate, it gets a high flutter velocity. The same as the case of the effect of the center of the gravity position, on the flutter, the frequency coalescence between the second mode and the third mode leads to flutter instability rather than the first mode and the second mode when the two mounting point spacing becomes big. The alterations of the coupling mode orders are observed. When e₁ increases, the flutter velocity gradually decreases.

For the simply supported plate, the location of the store will significantly affect the flutter boundary of the system, and when the store is near the middle position of the composite laminated plate, the gravity center is located at 1/3 of the distance of the two attaching points, the spacing of the mounting points is reduced appropriately, and the critical flutter velocity of the structure will improve, which lays the foundation for the actual design. Carrying the store may increase the critical flutter velocity of the plate with simply supported boundary. For the cantilever plate, when the store is close to the fixed edge, the system can get a big flutter velocity. The store reduces the aerodynamic stability of the structure, which is contrary to the simply supported boundary.

5 Conclusions

The vibration characteristics of the composite plate structure under different boundary conditions are investigated with the Rayleigh-Ritz method, and the boundary constraint is modeled as the artificial springs. The effects of a store on the flutter boundary are discussed.

(ⅰ) The frequency and mode of the plate with different boundaries are obtained by setting different stiffness values of the artificial springs and increasing the numbers of the orthogonal polynomials, the precision of the result will be further improved. The ply angle of composite material can significantly affect the frequency of the structure, which in turn affects the flutter characteristics.

(ⅱ) For the simply supported plate, the flutter velocity can be raised when the store is located at the center of the plate. When the gravity position is located at 1/3 of the distance of the two attaching points, the flutter velocity gets the peak value, and the change of the center of the gravity position affects the coincidence frequency, which shows that the critical flutter velocity has a jump. After the flutter occurs, the superposed modes can be re-separated, and gradually coincide with other modes, and a complex frequency coincidence phenomenon occurs. For the cantilever plate, when the store is near the constraint edge, the system can obtain a better aerodynamic stability. Adding a store reduces the critical flutter velocity, which is contrary to the simply supported boundary.

(ⅲ) With the increase in the flow velocity, the mode jumping phenomenon and alterations of the coupling mode orders are observed, and only one frequency coincidence point cannot be used as the basis of the flutter occurrence. The methods in this paper break through the boundary constraints, which can be used to solve the flutter boundary and response of the plate with the general elastic support, e.g., the bolting or riveting in the practical engineering.