A Combined Energy Method for Flutter Instability Analysis of Weakly Damped Panels in Supersonic Airﬂow

: A combined energy method is proposed to investigate the ﬂutter instability characteristics of weakly damped panels in the supersonic airﬂow. Based on the small damping assumption, the motion governing partial differential equation (PDE) of the panel aeroelastic system, is built by adopting the ﬁrst-order piston theory and von Karman large deﬂection plate theory. Then by applying the Galerkin procedure, the PDE is discretized into a set of coupled ordinary differential equations, and the system reduced order model (ROM) with two degrees of freedom is obtained. Considering that the panel aeroelastic system is non-conservative in the physical nature, and assuming that the panel exhibits a single period oscillation on the ﬂutter occurrence, the non-conservative energy balance principle is applied to the linearized ROM within one single oscillation period. The obtained result shows that the ROM modal coordinate amplitudes ratio is regulated by the modal damping coefﬁcients ratio, though each modal damping coefﬁcient is small. Furthermore, as the total damping dissipation energy can be eliminated due to its smallness, the He’s energy balance method is applied to the undamped ROM, therefore the critical non-dimensional dynamic pressure on the ﬂutter instability occurrence, and the oscillation circular frequency amplitude relationship (linear and nonlinear form) are derived. In addition, the damping destabilization paradoxical inﬂuence on the system ﬂutter instability is investigated. The accuracy and efﬁciency of the proposed method are validated by comparing the results with that obtained by using Routh Hurwitz criteria.


Introduction
Skin panels of high speed vehicles in supersonic airflow may experience the flutter oscillations due to the aeroelastic interaction mechanism. The oscillations can cause severe damage to the skin panel structures, and threaten the vehicle's integrity. Therefore, the panel aeroelastic behaviors, including the instability and the nonlinear flutter oscillations characteristics, have drawn great attentions to relevant researchers [1][2][3][4][5][6][7][8][9][10][11]. For the panel in supersonic airflow, there exist two types of panel flutter behavior, i.e., the single mode type and the coupled modes type [3][4][5]. These two types panel flutter instability characteristics are both related to the system damping. The single mode type flutter behavior can occur in low supersonic airflow, while the coupled modes type flutter behavior mainly occurs in high supersonic airflow. In this study, we mainly investigate the coupled modes type panel flutter instability characteristics using a combined energy method.
In the panel flutter instability research field, several sophisticated methods have been developed, such as, the classical Routh Hurwitz criteria [3,10], the eigenvalue analysis method [4,5], and the asymptotic method [6]. Though the above methods are feasible and effective, their evaluation procedures are complicated on some level, and lack of definite physical meanings from the energy transfer perspective. Recently, a specific critical parameters evaluation method based on non-conservative energy balance principle was proposed to predict the system flutter instability characteristics, but the procedure is also complicated [11].
For the panel aeroelastic system in supersonic airflow, the piston theory is usually adopted to simulate the aeroelastic interaction mechanism. This piston aerodynamics can play a role as the circulatory force exerting on the aeroelastic system, through the aerodynamic stiffness term ∂w/∂x, where w and x are the panel deflection and the chordwise coordinate respectively. If the continuous panel aeroelastic system is discretized and transformed into the modal coordinates, each order modal oscillator is coupled through the aerodynamic stiffness, and the circulatory system can experience the coupled modes type flutter instability. Based on this consideration, Hamilton energy conservation law can be applied to investigate the flutter instability, and the energy transfer characteristics of the undamped panel aeroelastic system.
In addition to the system stiffness, the system damping also exhibit significant influence on the system flutter instability characteristics, especially the damping paradox [11][12][13][14][15][16][17][18][19]. In this field, the damping paradox was firstly noted by Dowell [1]. Then, the complex role played by the system damping was investigated by Bolotin using Routh Hurwitz criteria [3]. Recently, a non-conservative energy balance principle was adopted to evaluate the damping paradoxical influence [11]. Actually, the damping paradox has been analyzed based on some simple non-conservative models, such as Ziegler pendulum and Pflüger column. Based on these models, Kirillov concluded that the damping paradox is related to the system symmetry breaking from the PT-symmetry perspective [13][14][15]. Bigoni supplied both the theoretical and experimental results to validate the damping destabilization paradox [16][17][18]. Luongo gave an interpretation that the damping destabilization paradox is related to the sign of the projection of the damping force on the eigenvector of the dual basis with the concept of "modal damping" [19]. Based on these studies, it can be concluded that the damping destabilization paradox is more related to the damping distribution characteristics in the modal coordinates, rather than the system damping dissipated energy.
Therefore, considering the reality that the damping within the panel aeroelastic system is small, thus we try to propose a combined energy method to investigate the flutter instability characteristics of weakly damped panels in this study. The remainder of this theoretical study is organized as follows. In Section 2, the system motion governing equation is built. In Section 3, the incremental non-conservative energy balance equation within a closed trajectory is derived. In Section 4, based on the small damping assumption, He's energy balance method [20][21][22] is adopted to investigate the system Hamilton energy balance relationship, and to derive the oscillation circular frequency amplitude relationship (the linear form). Additionally, the panel flutter instability characteristics and the damping paradox can be investigated. In Section 5, the nonlinear frequency-amplitude relationship is obtained. Parallel with the above studies, the proposed combined method is validated by comparing the results in this study with that derived by Routh Hurwitz criteria. Finally in Section 6, some conclusions are drawn.

Formulation of Motion Equation
Considering a two-dimensional, isotropic flat panel with simply supported boundary conditions in supersonic airflow (see Figure 1) [1], the panel motion governing partial differential equation (PDE) can be obtained using the piston theory and von Karman large deflection plate theory, as, where, the panel transverse deflection is w(x, t), the large deflections induced panel in- dx, the panel bending stiffness is D = Eh 3 12(1−ν 2 ) , the panel density is ρ m , the panel length is l (along the chord-wise x direction, see Figure 1), the panel thickness is h, the panel material elastic modulus is E, Poisson's ratio is v, the air density is ρ a , the airflow velocity is V ∞ , and Mach number is Ma. Ad- indicates the aerodynamic loading by using the first-order piston theory, q = ρ a V 2 ∞ /2 is the dynamic pressure of the supersonic airflow, β = √ Ma 2 − 1 is the Prandtl Glauert factor.
plane force is , the panel density is m ρ , the panel length is l (along the chord-wise x direction, see Figure 1), the panel thickness is h , the panel material elastic modulus is E , Poisson's ratio is v , the air density is a ρ , the airflow velocity is ∞ V , and Mach number is Ma . Additionally, indicates the aerodynamic loading by using the first-order piston theory, 2 2 ∞ = V ρ q a is the dynamic pressure of the supersonic airflow, is the Prandtl Glauert factor. Following with Dowell's procedure [1], and introducing the following non-dimensional parameters: , Equation (1) can be rewritten as,  is the corresponding mode shape function. Using this expression, the panel non-dimensional in-plane force can be rewritten as Substituting the above equations into the governing equation (Equation (2)), then multiplying each term by ) sin( ξ mπ and integrating each term over the non-dimensional length ξ , a set of coupled nonlinear ordinary differential equations (ODEs) can be obtained as, Following with Dowell's procedure [1], and introducing the following non-dimensional (1) can be rewritten as, the panel simply supported boundary conditions are, Equations (2) and (3) constitute a boundary value problem. To solve this problem, the Galerkin method is applied. Considering the boundary conditions described as Equation (3), the panel non-dimensional deflection W(ξ, τ) can be expressed as, where, q m (τ) is the mth modal coordinate and sin(mπξ) is the corresponding mode shape function. Using this expression, the panel non-dimensional in-plane force R x = N x l 2 /D can be rewritten as R x = 3 N ∑ m=1 q 2 m (mπ) 2 . Substituting the above equations into the governing equation (Equation (2)), then multiplying each term by sin(mπξ) and integrating each term over the non-dimensional length ξ, a set of coupled nonlinear ordinary differential equations (ODEs) can be obtained as, where N is the selected modes number. It is noted that in Equation (5), only the aerodynamic damping is considered, as each order modal damping coefficient is λR M /π 4 . Since that there may exist other damping models (such as the structural viscoelastic damping) with arbitrary modal distribution characteristics, so that each order modal damping Mathematics 2021, 9, 1090 4 of 11 coefficient λR M /π 4 is replaced by ε s=1,2 to consider the arbitrary damping models. Additionally, the following equations are obtained, where, the sth order modal damping coefficient ε s is assumed having small values (ε 1,2 << 1), i.e., the system is weakly damped. As the system coupled modes type flutter instability characteristics of the panel initial flat equilibrium can be analyzed using the panel first two modes (N = 2), a weakly damped reduced order model (ROM) with two degrees of freedom (2-DOF) can be built, ..
This ROM can be treated as two weakly damped (ε 1,2 << 1) nonlinear oscillators with coupled aerodynamic stiffness terms ( 8λ 3π 4 q 1,2 ). Additionally, taking the difference between the non-conservative energy and the conservative energy into account, a combined energy method, consisting of non-conservative energy balance principle and Hamilton energy conservation law (He's energy balance method), is proposed to investigate the panel flutter instability characteristics in the following sections.

Non-Conservative Energy Balance Principle
In this section, a non-conservative energy balance principle is applied to the ROM described as Equation (7). It should be noted that a complicated parameters evaluation procedure based on the non-conservative energy balance principle has been proposed by the author in Reference [11]. While in this study, the system critical parameters are evaluated based on Hamilton energy conservation law, while the non-conservative energy balance principle is adopted to derive the relationship between the modal coordinate amplitudes and the modal damping coefficients. Therefore, the proposed combined energy method in this study is different to the previous one [11].
The power flow equation for the ROM is derived by multiplying each term of the first equation and the second equation in Equation (7) with . q 1 and . q 2 respectively, as, For the panel flutter oscillation on the flutter occurrence, since the modal coordinate q 1 and q 2 are very small, the above power flow equation can be linearized by eliminating the nonlinear terms, as, ..
Considering that the system non-conservative energy is path dependent, it is necessary to select a specific time interval [τ 1 , τ 2 ] to guarantee a closed trajectory for removing the influence of the system conservative energy. Therefore, the relationship between the system stiffness related circulatory type energy and the damping related dissipation type energy can be derived. On the flutter occurrence, assuming that the panel is experiencing a neutrally stable single period oscillation with the circular frequency ω, thus the time interval [τ 1 , τ 2 ] can be specified as one single oscillation period 0, 2π ω , and the non-conservative energy balance equation for the linearized ROM can be derived as, where, ∆E 1 and ∆E 2 denote the incremental non-conservative energy of the first and second order mode within one single oscillation period, respectively. For the panel's neutrally stable oscillation on the flutter occurrence, there exists the following non-conservative energy balance equation, Additionally, the selected two modal coordinate, q 1 (τ) and q 2 (τ) are assumed as, where, a and b are the two modal coordinate amplitudes, θ is the positive modal phase shift induced by the system damping, and ω is the oscillation circular frequency. Substituting Equation (12) into Equation (10) and integrating each term over 0, 2π ω , obtains, or in the form, As on the flutter occurrence, there exists a condition, a b T = 0 0 T , which demands the coefficients matrix determinant being zero, thus we have, yields the phase shift, Based on the small damping assumptions, the phase shift can be very small. In addition, the energy dissipated by the system damping, should be balanced by the energy conserved by the system aerodynamic stiffness, With considering the energy balance relationship described as Equation (11), the following energy balance equations can be obtained as, 8 3 2λ where η = ε 2 /ε 1 denotes the ratio of the modal damping coefficients. From Equation (20), it can be observed that the modal coordinate amplitudes ratio can be regulated by the modal damping coefficients ratio η. Additionally, it is noted that this relationship is only dependent on the ratio of two modal damping coefficients, but independent on their values. Correspondingly, the mode configurations on the system flutter instability occurrence can be regulated by η as shown in Figure 2. For fixed b = 1, the first order modal coordinate amplitude a can be increased with increasing η. This variation tendency demonstrates a side effect of the damping paradoxical influence (discussed in Section 4). Furthermore, the ROM fluttering configurations can be affected by η. Therefore, the panel geometric dependent in-plane stress distribution is altered, so that the panel nonlinear flutter behavior is also affected by the ratio η. Based on the above discussions, the maximum deflection of the fluttering configuration on the system flutter instability occurrence can be altered as shown in Figure 3.
With considering the energy balance relationship described as Equation (11), the following energy balance equations can be obtained as, Substituting Equation (16) into Equation (19) obtains, denotes the ratio of the modal damping coefficients. From Equation (20), it can be observed that the modal coordinate amplitudes ratio can be regulated by the modal damping coefficients ratio η . Additionally, it is noted that this relationship is only dependent on the ratio of two modal damping coefficients, but independent on their values. Correspondingly, the mode configurations on the system flutter instability occurrence can be regulated by η as shown in Figure 2. For fixed , the first order modal coordinate amplitude a can be increased with increasing η . This variation tendency demonstrates a side effect of the damping paradoxical influence (discussed in Section 4).
Furthermore, the ROM fluttering configurations can be affected by η . Therefore, the panel geometric dependent in-plane stress distribution is altered, so that the panel nonlinear flutter behavior is also affected by the ratio η . Based on the above discussions, the maximum deflection of the fluttering configuration on the system flutter instability occurrence can be altered as shown in Figure 3.

He's Energy Balance Method
In this section, He's energy balance method, based on Hamilton energy conservation law, is applied to investigate the panel flutter instability characteristics. The reason for using this method to deal with the weakly damped panel's flutter instability, is that the system damping dissipated energy can be eliminated due to its smallness. Even though, the modal coordinate amplitudes ratio can be regulated by the modal damping coefficients ratio (Equation (20)). After eliminating the damping dissipation energy, there exists Hamilton energy to sustain the system flutter oscillation. Based on the panel aeroelastic system's reduced order model, this Hamilton energy is mainly related to the system coupled aerodynamic stiffness, as the aerodynamic stiffness achieve the flow energy transfer process. Therefore, the He's energy balance method and He's frequency amplitude formulation are practicable to theoretically derive the system non-dimensional critical dynamic pressure, and the circular frequency.
Based on the small damping assumption, the phase shift θ can also be neglected due to its smallness, and

He's Energy Balance Method
In this section, He's energy balance method, based on Hamilton energy conservation law, is applied to investigate the panel flutter instability characteristics. The reason for using this method to deal with the weakly damped panel's flutter instability, is that the system damping dissipated energy can be eliminated due to its smallness. Even though, the modal coordinate amplitudes ratio can be regulated by the modal damping coefficients ratio (Equation (20)). After eliminating the damping dissipation energy, there exists Hamilton energy to sustain the system flutter oscillation. Based on the panel aeroelastic system's reduced order model, this Hamilton energy is mainly related to the system coupled aerodynamic stiffness, as the aerodynamic stiffness achieve the flow energy transfer process. Therefore, the He's energy balance method and He's frequency amplitude formulation are practicable to theoretically derive the system non-dimensional critical dynamic pressure, and the circular frequency.
Based on the small damping assumption, the phase shift θ can also be neglected due to its smallness, and q 1 (τ) and q 2 (τ) in Equation (12), are approximated as, Following with the He's energy balance method [20], a time scale parameter κ = ωτ is introduced into Equation (7). Then, eliminating the nonlinear stiffness and damping terms, and denoting the differentiation with respect to κ by the primes, the following equation of the linearized undamped 2-DOF ROM can be obtained, Its Hamiltonian equation is, Using Equation (21), the initial conditions (κ = 0) for Equation (23) are, and the Hamilton equation is (κ = 0), In addition, considering the condition κ = π 4 , and using the Hamilton energy conservation law, the following equation is derived, Solving this homogeneous linear equation, the oscillation circular frequency amplitude relationship for the linearized ROM can be obtained as, It is observed from Equation (27) that the system critical oscillation circular frequency is dependent on the two modal coordinate amplitudes, while the modal coordinate amplitudes ratio is regulated by the modal damping coefficients ratio (Equation (20)). Therefore, substituting Equation (20) into Equation (27) yields, Based on Equation (28), it is concluded that the flutter oscillation circular frequency is related with the modal damping coefficients ratio η, and this equation meets well with the one obtained by the non-conservative energy balance principle based parameters Mathematics 2021, 9, 1090 8 of 11 evaluation procedure [11]. Additionally, for η = 1, we have ω = √ 17/2, this equation is identical with the one obtained in studies [4,10,11]. The variation of ω with changing ε 1 and ε 2 is shown in Figure 4. While the variation of ω with changing η is shown in Figure 5. It is shown that with increasing η, ω decreased gradually.
Based on Equation (28), it is concluded that the flutter oscillation circular frequency is related with the modal damping coefficients ratio η , and this equation meets well with the one obtained by the non-conservative energy balance principle based parameters evaluation procedure [11]. Additionally, for , this equation is identical with the one obtained in studies [4,10,11]. The variation of ω with changing 1 ε and 2 ε is shown in Figure 4. While the variation of ω with changing η is shown in Figure 5. It is shown that with increasing η , ω decreased gradually.
This equation is identical with the one in References [3,4,10,11], if the system damping is assumed with a small value. From Equation (29), it is observed that the system critical non-dimensional dynamic pressure is regulated by the ratio , which is depending on the modal damping coefficients ratio η . As for the ratio η η + 1 2 , it is easily concluded that for This equation is identical with the one in References [3,4,10,11], if the system damping is assumed with a small value. From Equation (29), it is observed that the system critical non-dimensional dynamic pressure is regulated by the ratio 2 √ η 1+η , which is depending on the modal damping coefficients ratio η. As for the ratio 2 √ η 1+η , it is easily concluded that for η > 0, there always exists 2 √ η 1+η ≤ 1, i.e., the system damping "Ziegler" paradoxical influence, and only when η = 1 or the undamped situation, 2λ π 4 can reach its max value, 2λ π 4 η=1 = 45 8 . The variation of 2λ π 4 with changing ε 1 and ε 2 is shown in Figure 6. It is shown that this symmetric surface separates the stable domain (below the surface) from the unstable one (above the surface). In mathematics, this surface is well known as the "Whitney's umbrella" surface. It should be noted that there exists a ridge line on the surface, dividing the surface into two symmetric surface, and along this ridge line, there exits η = 1. Additionally, the variation 2λ π 4 with changing η is shown in Figure 7. It is shown that with increasing η, 2λ π 4 firstly increased sharply to its maximum value 2λ = 45 8 , then decreased gradually.

Nonlinear Frequency Amplitude Relationship
In addition to the critical parameters (circular frequency ω and non-dimensional dynamic pressure 4 2 π λ ) on the system instability occurrence, the nonlinear circular frequency amplitude relationship for the system post flutter oscillation is derived theoretically by applying He's energy balance method in this section.
Substituting Equation (21) and the time scale parameter κ into Equation (7)

Nonlinear Frequency Amplitude Relationship
In addition to the critical parameters (circular frequency ω and non-dimensional dynamic pressure 4 2 π λ ) on the system instability occurrence, the nonlinear circular frequency amplitude relationship for the system post flutter oscillation is derived theoretically by applying He's energy balance method in this section.
Substituting Equation (21) and the time scale parameter κ into Equation (7) π κ = , the following equation can be obtained,

Nonlinear Frequency Amplitude Relationship
In addition to the critical parameters (circular frequency ω and non-dimensional dynamic pressure 2λ π 4 ) on the system instability occurrence, the nonlinear circular frequency amplitude relationship for the system post flutter oscillation is derived theoretically by applying He's energy balance method in this section.

Conclusions
In this paper, a combined energy method was proposed to theoretically investigate the flutter instability characteristics of weakly damped panels in the supersonic airflow. With this method, both the critical non-dimensional dynamic pressure and the oscillation circular frequency on the system flutter occurrence were derived theoretically. The system damping paradoxical effect was investigated. The accuracy and efficiency of the proposed techniques were validated by comparing the derived results with that obtained by other methods. Additionally, the nonlinear circular frequency amplitude relationship of the panel post flutter oscillation was derived theoretically for the first time, which can enhance our understanding of the panel post-flutter behavior.
Author Contributions: Conceptualization, X.W.; writing-original draft preparation, X.W. and X.X.; writing-review and editing, X.W. and G.Z.; supervision, Z.Y.; project administration, X.W.; funding acquisition, X.W. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data included in this study are available upon request by contactwith the corresponding author.