A Novel Aerodynamic Force and Flutter of the High-Aspect-Ratio Cantilever Plate in Subsonic Flow

)is paper focuses on the derivation of the aerodynamic force for the cantilever plate in subsonic flow. For the first time, a new analytical expression of the quasi-steady aerodynamic force related to the velocity and the deformation for the high-aspect-ratio cantilever plate in subsonic flow is derived by utilizing the subsonic thin airfoil theory and Kutta-Joukowski theory. Results show that aerodynamic force distribution obtained theoretically is consistent with that calculated by ANSYS FLUENT. Based on the first-order shear deformation and von Karman nonlinear geometric relationship, nonlinear partial differential dynamical equations of the high-aspect-ratio plate subjected to the aerodynamic force are established by using Hamilton’s principle. Galerkin approach is applied to discretize the governing equations to ordinary differential equations. Numerical simulation is utilized to investigate the relation between the critical flutter velocity and some parameters of the system. Results show that when the inflow velocity reaches the critical value, limit cycle oscillation occurs. )e aspect ratio, the thickness, and the air damping have significant impact on the critical flutter velocity of the thin plate.


Introduction
With the development and popularization of unmanned aerial vehicles in fields of detection, disaster prevention, and disaster reduction, researchers paid more attention on aeroelastic problems of high-aspect ratio airfoils in subsonic airflow. Flutter is self-excited oscillation of a flight vehicle under the coupling effect of the aerodynamic pressure, the elastic force, and the inertia force. Cantilever plates in the axial flow may lose stability at sufficiently high flow velocity. Analysis of the linear theory indicates that there is a critical dynamic pressure. e motion of the panel becomes unstable when the dynamic pressure is higher than the critical value. Once the instability threshold is exceeded, flutter will take place. Under the condition of flutter in the system, the energy of the surrounding airflow will be continuously pumped into the plate to sustain the flutter motion. e lightweight and high performances of the modern aircraft make the aeroelastic problems of the aircraft more prominent.
One of the key issues of aeroelasticity is flutter, which usually leads to the disaster of the aircraft [1]. us, the aeroelasticity problem of the aircraft comes into our sight. Zhang et al. [2] applied the composite multilayer plate to supersonic aircraft under the aerodynamic pressure to investigate excessive nonlinear vibration suppression of the plate. Flutter is a fluid-structure coupling problem. Chernobryvko et al. [3] discussed nonlinear dynamic stability conical shells in a supersonic gas stream. Amabili and Pellicano [4] found that flutter is very sensible to small initial imperfections of the structure. Vedeneev [5] studied flutter of an elastic thin plate and obtained the exact solution by solving the structural kinetic equation coupling with the hydrodynamic equation. When solving dynamic problems of the fluid-structure coupling system, researchers prefer adopting aerodynamic models to simulate the external flow field. Based on these aerodynamic models, aerodynamic characteristics related to motion state of the system, such as the displacement, velocity, and acceleration, are obtained. en, nonlinear dynamic equations of airfoils are derived, and aerodynamic characteristics of the system are investigated.
Aerodynamic force is a key factor in the analysis of the flutter phenomena, which directly determines the bucking form of the plate structure. So, it is necessary to select a proper aerodynamic force model in the study of flutter problem. Up to now, scholars have proposed abundant theories of the aerodynamic force. According to the dependence of aerodynamic forces on time and spatial, aerodynamic theories can be roughly divided into three types, such as the steady, quasi-steady, and unsteady aerodynamic theories. e first type is the steady aerodynamic model, namely, assuming the force acting on the lifting surface of the wing do not change with time. e steady aerodynamic model is mainly used for static aeroelastic analysis, such as the thin airfoil theory [6]. e second kind is the quasi-steady aerodynamic theory, in which we assume the aerodynamic force at any time is only related to the motion state of the wing at that time. By using this theory, the reduced frequency of motion is smaller, and vibration characteristics analysis of the airfoil subjected to the aerodynamic force can be simpler. Lin et al. [7] proposed a quasisteady piston theory to calculate the flutter problem of a cantilever plate in supersonic flow. Hu et al. [8][9][10] investigated aeroelastic vibration of the plate subjected to the aerodynamic force obtained from the first order piston theory. Brouwer and McNamara [11] optimized the piston theory and verified it. Dowell and Ganji [12] extended the piston theory to higher order terms in several expansions and analyzed the flutter of single degree of freedom panel.
Owing to the ignoring of the propagation effect of small perturbations in the subsonic airflow, the piston theory is not applicable in the subsonic airflow. us, scholars gave some new methods to compute aerodynamic. Grossman quasi-steady theory [13] can also be applied to calculate the total lift force and the total torque of high-aspect ratio wings, which can be simplified as an elastic beam without pretwist and the axial extension. Obviously, this theory is not suitable for plates with the small aspect ratio. e third one is the unsteady aerodynamic theory, which considers the influence of change of circulation and wake flow on the aerodynamic force of a moving airfoil. Xu et al. [14] deduced the solution of the unsteady aerodynamic force for a slender airfoil, which is modeled as a beam structure. However, this method can only be used to solve problems of 2D plate structures in the incompressible airflow. Based on the model constructed, by applying eodorsen unsteady aerodynamic theory, researchers investigated stability of wings. However, Cordes et al. [15] pointed out eodorsen function cannot capture the experimental transfer functions in frequency dependence when investigating the unsteady lift force of an airfoil.
Since exact solutions of unsteady potential equations are few, we can only acquire approximate analytical solutions through other methods. Peters [16] obtained a semiempirical, nonlinear, and unsteady ONERA model according to data of wind tunnel airfoils. Dunnp and Dugundjij [17] analyzed linear flutter of the composite wing using the harmonic balance method by adopting ONERA dynamic stall model. Sadr et al. [18,19] developed the ONERA model. However, the eodorsen model and the ONERA model are not suitable for the aerodynamic analysis of plate structures. Zhang and Ye [20] established the integral aerodynamic reduced order model [21] based on the Volterra series theory. However, introduction of the integral form leads to the computational complexity. Zhang and Ye [22] developed the Volterra series theory and expressed the aerodynamic force as the sum of multiple convolutions. However, the aerodynamic force obtained is implicit, so it is hard to couple the implicit expression with kinetic equations of the plate. Guo et al. [23,24] put forward the nonlinear harmonic balance method to analyze the unsteady flow in the turbomachinery and the airfoil. Based on proper orthogonal decomposition, Luo et al. [25] presented a hybrid modeling method for reconstructions of flow field and aerodynamic optimization. Nonlinear regression methods, instead of the linear regression widely used, are adopted to establish POD basis modes, which behave with good description performance in system space. Wang et al. [26] made a comprehensive review on the latest studies about the aeroelastic modeling.
e CFD method is adopted in Refs. [27][28][29]; however, this method is not conducive to perform pneumatic elastic servo analysis. Munk et al. [30] studied the limit cycle of a two-dimensional cantilever plate under subsonic flow based on the vortex lattice method. Castells Marin and Poetsch [31] used the doublet lattice method to model the lifting surface, which is more accurate than the NLRI method. Xie et al. [32] obtained the aeroelasticity deformation of the geometrically nonlinear high-aspect ratio wing, which is in great agreement with the experimental result. Pashaei et al. [33] modeled the airflow by using the vortex lattice method and studied the effect of energy harvesting properties of the metal composite on the flutter margin and limit cycle oscillation amplitudes. Ramezani et al. [34,35] established the aerodynamic model by adopting CFD/CSD coupling numerical computational method. Chen et al. [36][37][38] constructed reduced-order models of high speed vehicles based on CFD simulations. Results show that the ROM approach can significantly speedup unsteady aerodynamic calculations of a system.
Most of the results obtained from vibration of plates under subsonic airflow are related to flutter of panels. Flutter of the panel is similar to that of the wing. e main difference lies in that there is only one surface subjected to the aerodynamic force for the panel. However, for the plate, both surfaces are exposed to the airflow. Study of panel flutter began in the early 1970s. Tang et al. [39,40] theoretically and experimentally researched the aeroelastic response of a wing and a cantilever plate under subsonic airflow. Dowell et al. [41,42] proposed the linear potential flow theory, which can calculate the pressure distribution of any point on the surface of the panel. However, it is not suitable for plate or shell structures with airflow acting on both sides.
Generally, achievements have been made in the study of 2D infinite plate, panels, and beam structures in subsonic airflow. And aerodynamic models have already been applied in investigation of flutter. However, research works on explicit expressions of aerodynamic forces of cantilever plates and shells under subsonic airflow are still few. In this paper, for the first time, an analytical expression of the quasisteady aerodynamic force for the high-aspect ratio cantilever plate in subsonic flow is induced based on the subsonic thinairfoil theory and Kutta-Joukowski lift theorem. Overall, aerodynamic force theoretically calculated by using the explicit expression we derived has a good agreement with that obtained by ANSYS FLUENT. In addition, the aerodynamic model constructed based on it could be applied to flutter analysis of cantilever plates and shells with the highaspect ratio. Considering lateral vibration and deformation of the mean camber line of the cross section, nonlinear dynamic equations of transverse vibration of the high-aspect ratio cantilever plate are derived by utilizing the quasi-steady aerodynamic model. Influences of parameters of the system on the critical flutter velocity are investigated.

Analysis of the Aerodynamic Force for the High-Aspect-Ratio Cantilever Plate in Subsonic Flow
Field. e schematic diagram of the cantilever plate considered is shown in Figure 1. e wing is simplified as a high-aspect ratio cantilever plate. e span length, chord length, and thickness of the plate are a, b, and h, respectively. e velocity of the subsonic airflow along the chordwise direction is denated as U ∞ . Cross section of the cantilever plate is marked as A. (X, Y, Z) is the inertial coordinate system, and the origin of it is in point O. e ⇀ 0 x is the spanwise direction, e ⇀ 0 y is the chordwise direction, and e ⇀ 0 z is the thickness direction, respectively. Based on the strip assumption, Kutta-Joukowski lift theorem, and linearized small perturbation theory, we induce an aerodynamic force model of a high-aspect ratio cantilever plate in subsonic airflow. e strip assumption of the plate with high-aspect ratio can be briefly introduced as follows. Actually, the airflow on wings is a three-dimensional fluid. However, if then the geometric dimension of the cantilever plate does not change along the chordwise direction, the aspect ratio is high and the inflow velocity is not change along the spanwise. e velocity of the fluid can be considered as a component in 2D plane (YOZ plane) and the component in the spanwise axis (X-axis) is zero in the most part of the spanwise region. us, every chord section can be analyzed as a 2D airfoil with an infinite spanwise length. Diagrams of each section along the chordwise direction, such as section A, can be shown in Figure 1.
At the beginning of the 20th century, Joukowski proposed Kutta-Joukowsi lift theorem, which established the relation of the lift and circular rector of moving objects in the air, as shown in equation (1) . In the incompressible, lowvelocity, inviscid, and straight uniform flow field, the force distribution on unit length of the spanwise of the closed 2D wing is perpendicular to the direction of the airflow. Its value can be expressed by product of the density of the fluid ρ, flow velocity U ∞ , and the vortex strength c(Y) of unit arc length of the wing: (1) e vortex strength c(Y) on unit arc length is positive in a clockwise direction. When using Kutta-Joukowski lift theorem (inviscid potential flow theory) to solve the lift of the airfoil, the chief problem among the issues is how to obtain the local vortex strength c(Y). e airfoil, whose ratio of the maximum thickness and the chord length is less than 12%, is defined as a thin airfoil. For a thin airfoil, we can use linearized small perturbation theory of the low speed flow around a thin airfoil in aerodynamics to calculate the local vortex strength c(Y). e problem about flow around a thin airfoil means the small attack angle and the small bending. In addition, a thin airfoil problem means the thin thickness. us, boundary conditions and pressure coefficient of the airfoil can be linearized. erefore, based on principle of superposition, the attack angle, the bending, and the thickness can be considered separately and then superposed. e potential flow around a thin airfoil can be decomposed into three simple linear potential flows, which include the flow around a curved plate without an attack angle, a symmetrical airfoil without an attack angle, and a flat plate with a small attack angle. e flow around a symmetrical airfoil without an attack angle cannot generate the lift force, so we only need to consider the lift forces generated by the small attack angel and the bending of the airfoil. When the straight uniform airflow flows across the mean camber line of 2D airfoil with a small attack angle, we can use the surface vortex on the mean camber line to simulate the distribution of the local vortex strength c(Y), which has a trigonometric series solution. Lateral displacement of the mean camber line of the chordwise cross section is set up as W Y . e chordwise position of an arbitrary point in the middle arc line can be written as When the plate vibrates slightly, a very small displacement d W Y appears. According to the mathematical definition of limit, the tangent value d W Y /d Y is equal to the corresponding angle value. us, we express the attack angle as d W Y /d Y , that is, the tangent value of the attack angle. Figure 1: Model of the high-aspect-ratio cantilever plate.

Shock and Vibration
Substitute equation (2) into d W Y /d Y , then the expression of the attack angle is changed into function K θ , which is related to W Y and θ. erefore, expression of W Y is derived as after mentioned. If the attack angle K θ is integrable, as given in equation (4), the local vortex strength c(Y) can be expressed as c(θ): where According to equation (3), when the analytic expressions of the attack angle and the mean camber line are given, there is a unique trigonometric series solution of the local vortex strength c(θ). Coefficients of the solution can be confirmed by equations (4a) and (4b).
Linearized small perturbation theory of the flow around the thin airfoil is a steady theory, so it can only solve the steady local vortex strength of the thin airfoil without deformation and vibration. However, when a plate vibrates, additional attack angle caused by vibration velocity and trailing vortex caused by changing circular rector will impact local vortex strength, as shown in Figure 2. Lateral vibration velocity V w and inflow velocity U ∞ will generate a new additional attack angle θ 2 , which is approximately V w /U ∞ , as shown in Figure 3. e additional attack angle of the thin plate θ 2 is defined as V w /U ∞ , so flow theory of the thin airfoil can be used to calculate the local vortex strength caused by θ 2 . Substitute equation (2) into V w /U ∞ , expression of the additional attack angle is marked as Q θ , which is related to V w and θ. us, expression of V w is focused on. If the additional attack angle Q θ is integrable, the local vortex strength caused by vibration velocity can be expressed as a function c ′ (θ) as follows: where e total circular rector in the nonviscous flow field must be conserved. us, the integrable total circular rector of every point in that field must be zero. When a plate vibrates, equivalent wake vortex is generated by change of the circular rector. e wake vortex generated has a strong influence on the local vortex strength of airfoils. erefore, the wake vortex will also affect deformation of the plate. Deformation will also generate new wake vortex. Both the wake vortex and deformation of the plate have influence on the local vortex. Comparatively speaking, wake vortex has a smaller influence on the local vortex strength, so we introduce the quasi-steady hypothesis, in which the influences of the wake vortex on the local vortex strength are neglected. To sum up, under the quasi-steady hypothesis, we only need to consider the influence of the mean camber line's deformation and the additional angle of attack on local vortex strength.

Interpolation Functions of the Mean Camber Line's Deformation and Vibration
Velocity. In order to use equations (3)-(6a) and (6b) to calculate the local vortex strength caused by the additional attack angle and deformation of the mean camber line, it is necessary to know analytic expressions of the mean camber line's deformation of the chord section W Y and the vibration velocity distribution function V w . Moreover, K θ and Q θ should be integrable on θ. However, when cantilever plate and shell structures vibrate, deformation and velocity are different at different time, which may not satisfy the condition mentioned above. us, a fitting method of the interpolation function is applied to express the deformation function W Y and the vibration velocity distribution function V w . erefore, the sum of several interpolation functions can satisfy the condition. If these interpolation functions are integrable on θ after substituting Y into the deformation function W Y and the distribution functionV w , W Y and V w of the vibration velocity can be approximately expressed as integrable functions about θ. Nonequidistant Lagrange interpolation method satisfies the condition abovementioned, so this method is applied to fit W Y and V w . Interpolation functions are several polynomial functions of curves that pass through points given on the 2D plane.
Since the first mode of the cantilever plate vibration is the bending deformation, the first interpolation function WX 1 is set up as a constant, which is the average value of deformation of the mean camber line at Y � 0 and Y � b, as shown in Figure 4. WX 1 is expressed as After the first fitting of the mean camber line's deformation, difference between the deformation curve CC 0 of the mean camber line and the first interpolation function WX 1 is CC 1 , as shown in Figure 5.
Since the second order mode of the cantilever plate vibration is the torsion deformation, the second interpolation function WX 2 is set up as a linear function, which is shown in Figure 6.
Let values of the interpolation function at Y � 0 and Y � b be equal to values of residual CC 1 at Y � 0 and Y � b, respectively. us, slope of the first interpolation function can be calculated as follows: e second interpolation function WX 2 can be written as follows: Difference between residual CC 1 after the first fitting and the second interpolation function WX 2 is marked as CC 2 , which can be expressed as e third interpolation function WX 3 is set up as a quadratic function, as shown in Figure 7. Values of CC 2 at Y � 0 and Y � b are 0. In order not to increase the residual at Y � 0 and Y � b in the third fitting, the third interpolation function is taken as in order to get a minimal difference between the third interpolation function and the residual CC 2 . According to the equation, unknown coefficient K 2 and the third interpolation function WX 3 can be calculated. WX 3 is written as follows: Difference between the residual CC 2 in the second fitting and the third interpolation function WX 3 is marked as CC 3 . e difference can be expressed as e fourth interpolation function WX 4 is chosen as a cubic function, which is shown in Figure 8. Values of CC 3 at  Figure 5: e first interpolation function on the plane (Y, Z). Figure 6: e second interpolation function on the plane (Y, Z).

Shock and Vibration
In order not to increase the residual at Y � 0, Y � b, and Y � b/2 in the fourth fitting, the fourth interpolation function is taken as ) in order to make the fourth interpolation function fit the residual CC 3 to the greatest extent. According to the equation, unknown coefficient K 3 and the expression of WX 4 can be calculated. WX 4 is shown as follows: Similarly, difference between the residual CC 3 in the third fitting and the fourth interpolation function WX 4 is marked as CC 4 , which can be expressed as e fifth interpolation function WX 5 is set up as the form of a quartic function, as shown in Figure 9. Values of In order not to increase the residual at Y � 0, Y � b, and Y � b/2 in the fourth fitting, the fifth interpolation function is taken as . Difference at extreme points of the interpolation function 6 in order to make the fifth interpolation function fit the residual CC 4 to the greatest extent. According to the equation, unknown coefficient K 4 and the fifth interpolation function WX 5 can be calculated. WX 5 is shown as follows: Difference between the residual CC 4 in the fourth fitting and the fifth interpolation function WX 5 is marked as CC 5 . e difference can be expressed as After five times of deformation fitting abovementioned, the total residual CC 5 brought by the mean camber line's deformation CC 0 of five interpolation functions of the thin plate have the tendency to gradually converge to zero, as shown in Figure 10. If the mean camber line is more complicated, it is necessary to conduct more interpolation functions. For cantilever plates and shells with high-aspect ratio, chord deformation is relatively simple, so it is enough to take top four interpolation functions to fit the deformation.
Velocity is the derivative of displacement, which is continuous while the plate is vibrating, so vibration velocity distribution of the plate is also continuous. If we represent the mean camber line's vibration velocity of the chord section as a curve, a method that is totally similar to the above can be used to fit the curve. Form of the interpolation function of velocity distribution is exactly the same as that of the interpolation function of deformation. erefore, it only needs to replace items about deformation in equations (7) and (9)-(12) with a corresponding item about velocity. en, five interpolation functions for fitting velocity distribution are obtained as follows: Figure 9: e fifth interpolation function on the plane (Y, Z).

Lift Force Calculation. Precondition of calculating lift
force of the high-aspect-ratio cantilever plate is to obtain the local vortex strength, as shown in equation (1). Local vortex strength of the cantilever thin plate with high-aspect ratio is mainly caused by two factors, deformation of the mean camber line and the vibration velocity. Based on linearized small perturbation theory, the total local vortex strength caused by these two factors can be obtained by using linear superposition, as shown in equation (14). In order to calculate the local vortex caused by the vibration velocity, local vortex strength caused by each interpolation function of the velocity is calculated, respectively. Linear superposition is conducted to obtain total local vortex strength caused by the vibration velocity. Because vibration velocity of cantilever plate is relatively simple, it is accurate enough to take four interpolation functions to fit velocity distribution of the plate. us, top four interpolation functions VX 1 , VX 2 , VX 3 , and VX 4 are selected to calculate the local vortex. Firstly, angle of attack functions VX 1 /U ∞ , VX 2 /U ∞ , VX 3 /U ∞ , and VX 4 /U ∞ of four interpolation functions are calculated, respectively. And substitute equation (2) into VX 1 /U ∞ , VX 2 /U ∞ , VX 3 /U ∞ , and VX 4 /U ∞ . en, angle of attack functions are translated into functions related to θ. ese four functions are set up as Q 1 , Q 2 , Q 3 , and Q 4 , which are substituted into equations (6a) and (6b), respectively. A 0 ′ and A n ′ calculated by equations (6a) and (6b) are substituted into equation (5) to solve the corresponding local vortex: cVX 1 , cVX 2 , cVX 3 , and cVX 4 of top four interpolation functions. To sum up, the total local vortex strength caused by the mean camber line's deformation and the lateral vibration velocity is c z , which can be expressed as linear superposition of the local vortex strength calculated by abovementioned eight interpolation functions. c z is written as follows: where Shock and Vibration 7 According to equations (1) and (14), the aerodynamic force Δp can be expressed as Since the aerodynamic force expression is analytic, it is convenient to use analytic and semianalytic method to study the flutter problem of the cantilever plate.

Aerodynamic Correction and Error Analysis. Value of item
�������� (b/Y) − 1in equations (15a) and (15h) at Y � 0 is infinite, which leads to an infinite leading edge lift force. As a matter of fact, leading edge lift force of the wing cannot be infinite. Appearance of such a singularity at the leading edge of the wing which is attributed to the basic solution of the thin-wall theory gives no consideration to flow around the leading edge, namely, when air flows past the leading edge of the thin plate, part of air will pass through the upper panel from the lower panel. Neglecting thickness of the plate, the thin-airfoil theory leads to an infinite streaming velocity and an infinite lift force at the leading edge. As a result, it is necessary to correct this problem.
According to Ref. [43], although there is a singular point at the leading edge of the plate, the pressure distribution on 95% chord length range near the trailing edge has a good consistency with that of actual measurement. us, it is necessary to add a correlation coefficient in e infinity value of this function at Y � 0 is corrected to be equal to the value of the original curve at Y � 0.95b. After trial, the item 95b. e value of these two functions at the trailing edge portion of the plate changes little, as shown in Figure 11. e corrected aerodynamic expression is denoted as Δp′.
If the air on the plate flows at a speed greater than 0.3 times the speed of sound, influences of the compressibility of air on aerodynamic force cannot be neglected. us, it is essential to modify the impact of compression. Von Karman-Chandra Formula is used to estimate the influence of air compressibility on aerodynamic force, and equation (17) is the relationship between the two: aerodynamic pressure Δp p on the plate surface in nonsticky, steady, and subsonic velocity and 2D compressible flow field and the corresponding pressure Δp ′ in the incompressible flow.
Ma ∞ is the ratio of the flow velocity U ∞ to the local velocity of sound.
To sum up, after correction and considering the compressibility of air, the aerodynamic force expression Δp p is as follows: Aerodynamic force Δp p is the linear superposition of aerodynamic forces calculated by several interpolation functions. Moreover, inflow air must satisfy hypotheses of irrotational and nonviscous. us, the aerodynamic force Δp p calculated by equation (17) is an approximate result.
ere is an error between it and the real aerodynamic force. Effect of this approach is evaluated by estimating magnitude of the error between the two.
Mean camber line's deformation and the lateral vibration velocity are mainly considered in theoretical calculation of the aerodynamic force. e essential reason why the lift force can be generated is to change the attack angle of the wing, which indirectly affects the aerodynamic force. us, what is need is to make a comparison between the lift force generated by deformation of the mean camber line of the plate and that of the corresponding finite element model to illustrate effectiveness of the aerodynamic force theoretically calculated.
ANSYS FLUNT finite element software is applied to calculate the aerodynamic force distribution. A spline curve of definite shape whose chord length is 1 meter is drew in Computer Aided Design (CAD). Coordinates of control points of the spline are shown in Figure 12. A thin shell model of 0.01 meter thickness is constructed by stretching the spline curve we drew to 10 meters along the spanwise 8 Shock and Vibration direction. Assuming the elasticity modulus of the shell structure in subsonic, compressible, and nonsticky turbulence flow field is infinite. e velocity of the airflow is U ∞ � 200 m/s. After numerical simulation, pressure cloud pictures of the upper and lower surface can be obtained, as shown in Figures 13 and 14.
e lift force can be calculated by the pressure difference between the upper and lower surface of the airfoil, as shown in Figure 15.
If we take the spline curve in Figure 12 as a mean camber line of the cross section at a certain moment of the deformed high-aspect ratio cantilever plate structure, the aerodynamic force distribution on the section can be calculated according to equation (17). As shown in Figure 12, in calculation of the aerodynamic force, these items are introduced: , t � 0.581, In calculation of the lift force generated only by deformation of the mean camber line, velocities in the aerodynamic force expression Δp p are taken as zero, namely, Substituting abovementioned data into equation (17), aerodynamic force distribution calculated is shown in Figure 15. Trend of value theoretical calculated has a good agreement with that obtained by ANSYS FLUENT. In contrast with value obtained by ANSYS FLUENT, the result theoretical calculated is relatively small. In particular, the error is nearly about 20% at the trailing edge, as shown in Figure 15.
e comparison above indicates that the aerodynamic force fitted by utilizing interpolation functions can predict the influence of deformation on aerodynamic distribution approximately. Value of the aerodynamic force obtained by ANSYS is relatively large. Source of the error in the theoretical value is mainly due to the hypothesis of nonsticky and the neglecting of influence brought by the trailing vortex.

Establishment and Dispersion of Aerodynamic Equations for the Cantilever Plate
e schematic diagram of the wing is shown in Figure 1, which is simplified as a cantilever plate in subsonic flow. e rectangular plate is characterized by a × b × h, where a is the span length, b is the chord length, and h is the thickness. A is the cross section of the wing. (X, Y, Z) is the inertial coordinate system and the origin of it is located at the leading edge of the fixed end of the cantilever plate, which is marked by point O.
Considering the first-order shear deformation, Kirchhoff hypothesis [44], and scale effect, the displacement filed of the plate is established. Displacement of any point along x, y, and z direction can be expressed by that of the neutral plane of the plate as follows: v(x, y, z, t) � v 0 (x, y, t) + zφ y (x, y, t), w(x, y, z, t) � w 0 (x, y, t), 37.6 cm 58.10cm 8.3cm -45 cm -30cm Figure 12: e coordinate graph of the spline curve.

Shock and Vibration
where u 0 , v 0 , and w 0 are displacements along X, Y, and Z e constitutive relation of the isotropic material is expressed as follows: where E is the elasticity modulus of the material and υ is Poisson's ratio.   In order to simplify the calculation, relation between the internal forces and stresses is expressed as follows: Substituting nonlinear strains of equations (21a)-(21d) obtained by von Karman nonlinear geometric relationship into equations (22a)-(22e) and combining the result with equations (23a)-(23c), the potential energy of the plate are obtained as follows: e kinetic energy and the virtual work done by external forces are as follows: where ρ c is material density and C is the air damping. Nonlinear dynamic equations of motion for the rotating blade are established by using Hamilton's principle: e lateral displacement of the thin plate is more obvious compared with others. us, displacements u 0 and v 0 on the middle surface and rotation angle φ x and φ y are neglected. Substituting equations (24)-(26) into equation (27), nonlinear dynamic equation of the plate in the lateral direction is obtained as follows: 1 12 π 2 E (z/zy)ϕy(x, y, t) + z 2 /zy 2 w(x, y, t) h 2 + 2υ + Eυh(z/zx)w(x, y, t)((z/zy)w(x, y, t)) z 2 /zy zx w(x, y, t) − υ 2 + 1 + Eh(z/zx)w(x, y, t)((z/zx)w(x, y, t)) z 2 /zx 2 w(x, y, t) − υ 2 + 1 + (1/2υ)((z/zy)w(x, y, t)) 2 hE z 2 /zx 2 w(x, y, t) − υ 2 + 1 + (1/2)((z/zx)w(x, y, t)) 2 hE z 2 /zx 2 w(x, y, t) − υ 2 + 1 + π 2 E((z/zx)w(x, y, t)) 2 z 2 /zy 2 w(x, y, t) h 12 ×(2 + 2υ) + υ((z/zx)w(x, y, t)) z 2 /zy zx w(x, y, t)Eh(z/zy)w(x, y, t) − υ 2 + 1 + ((z/zy)w(x, y, t)) z 2 /zy 2 w(x, y, t)Eh(z/zy)w(x, y, t) I i is calculated by e boundary conditions of the fix end of the plate are e boundary conditions of the free end of the plate are Avramov and Mikhlin [45] have performed a lot of review of theoretical developments of nonlinear normal modes for continuum mechanical systems. Mode function that we choose must satisfy the first two order modes of lateral nonlinear vibration and the boundary conditions given in equations (30a), (30b), (31a), and (31b). According to mode approximation functions for the cantilever plate given in Ref. [46], the following mode functions are chosen where k 1 , k 2 , β 1 , and β 2 are calculated by e lateral displacement of the plate can be expressed as follows: Substituting equations (33) and (34) and parameters in Table 1 into equation (28), ordinary differential equations of transverse vibration for the first two modes of the cantilever plate are obtained by using the Galerkin method [47] as follows: In equations (35a) and (35b), x 1 and x 3 are generalized coordinates of the first-order and the second-order modals. e coefficient of U ∞ · _ x 1 is moved to the left end of the equation and is neglected, which shows that the aerodynamic force acts as a negative damping. us, amplitude of the vibrating structure increases continuously by absorbing energy from the airflow, which may lead to the occurrence of flutter.

Limit Cycle
Based on equations (35a) and (35b), Runge-Kutta algorithm is utilized to construct numerical simulation of flutter phenomenon for the cantilever plate subjected to the aerodynamic force in subsonic air flow. Initial values of the ordinary differential equations (35a) and (35b) are given as x 1 � 0.001, x 3 � 0.001, _ x 1 � 0, and _ x 3 � 0. When the inflow velocity U ∞ is chosen as 72 m/s, the system is in the stable state, as shown in Figure 16, in which vibration amplitude of the plate at X � 5 and Y � 0.8 decreases as time increases. Transverse displacement of the system is convergent when the inflow velocity fails to reach the critical flutter velocity. Figure 16(a) gives the waveform on the plane (t, w), where w is the transverse vibration displacement. Figure 16(b) represents the phase portrait on where v is the transverse vibration velocity. e adjacent trajectory gradually becomes a closed loop, as shown in Figure 16(b).
When the inflow speed U ∞ is set as 80.43 m/s, the system keeps in a critical state between stability and instability, as shown in Figure 17(a). e amplitude of the plate at X � 5 and Y � 0.8converges to a constant. e phase portrait of the transverse displacement w and the transverse velocity v of the vibrating plate at X � 5 and Y � 0.8 is shown in Figure 17(b). e phase portrait is an isolated closed loop, which do not have closed rails nearby.
is is called a limit cycle. At the moment, the critical flutter velocity is reached.
When the inflow velocity U ∞ is increased from 80.43 m/s to 89 m/s gradually, the system appears to diverge and becomes instable, as shown in Figure 18(a). e phase portrait of transverse displacement w and transverse velocity v of the vibrating plate at X � 5 and Y � 0.8 is shown in Figure 18(b). e phase portrait is far from a closed loop.
In this paper, the high-aspect-ratio subsonic cantilever plate system is nonconservative, which means the system will perform steady periodic vibration under specific conditions. As shown in equations (35a) and (35b), the aerodynamic force play the role of negative damping in vibration of the plate. us, when the inflow velocity reaches a certain value, there will be periodic vibration of equal amplitudes, as shown in Figure 17. is indicates that energy supplied by the aerodynamic force and that dissipated by damping are in balance. In fact, since there is no periodic external force in the system, the system can only absorb energy from surroundings by the negative damping caused by the aerodynamic force.
Periodic motions in nonconservative system belong to self-excited vibrations, which corresponds to stable limit cycles of the autonomous system, that is, when the initial conditions are disturbed, the original vibration state can still be restored. If the system subjects to some small disturbance, the original periodic motion will not be changed. For example, when the inflow velocity is chosen as 80.43 m/s, the system vibrates periodically with equal amplitudes. When the inflow velocity increases to a value smaller than 89 m/s, the system still keeps periodic motion with greater amplitude. Limit cycles of the linear system is often unstable and linear theory is not suitable for the system with large inflow velocity; thus, in this paper, nonlinear theory is more appropriate for limit cycles analysis of the cantilever plate.
Stability of the linear system can be judged by calculating eigenvalues, and the critical flutter velocity can be calculated. It is not easy to solve eigenvalues for nonlinear systems, so the critical flutter velocity cannot be solved by this method. We can assume the inflow velocity to be a certain value and substitute it into the ordinary differential equations of motion for the high-aspect-ratio cantilever plate. Based on the equations, phase portrait of the system are numerically obtained. When limit cycle appears in the phase space, the inflow velocity is exactly the critical flutter velocity.
In order to study the influence of system parameters on critical flutter velocity, different inflow velocity U ∞ is substituted into equations (35a) and (35b). When the phase portrait turns to be a limit cycle, the inflow velocity U ∞ reaches the critical flutter velocity. Parameters in Table 1 except thickness are brought into the equation. Corresponding critical flutter velocities of plates with different thickness are calculated, as shown in Figure 19(a). e critical flutter velocity increases as the thickness of the plate increases.
Change the length of the plate only to detect the influence of the aspect ratio on the critical flutter velocity. Other parameters are shown in Table 1. e critical flutter velocity is decreased with the increase of the aspect ratio, as shown in Figure 19(b). Keep other parameters the same as that in Table 1, and only change the damping coefficient to calculate the corresponding critical flutter velocity. As shown in Figure 19(c), when the damping coefficient increases, the critical flutter velocity increases.

Conclusions
In order to use analytical or semianalytical method to analyze flutter of the high-aspect ratio cantilever plate, it is necessary to obtain an aerodynamic analytical expression of the plate under subsonic airflow. Giving the assumption of irrotationality, nonsticky, incompressibility, and strip theory of high-aspect ratio and quasi-steady state, considering the skeleton line's deformation of airfoils' chord section and influences of lateral vibration velocity on quasi-steady aerodynamic force based on the thin-airfoil theory under subsonic flow and Kutta-Joukowski lift theorem, the authors induced a new quasi-steady aerodynamic analytical expression high-aspect ratio cantilever plate under subsonic airflow. It is confirmed that the analytical expression is consistent with the lift distribution tendency of finite element calculation. According to Reddy's first-order shearing plate theory and geometric equations of von Karman's large deformation, Hamilton principle was applied to establish nonlinear partial differential kinetic equation of cantilever plate, which suffers from aerodynamic force. e Galerkin method was used to obtain second-order ordinary differential governing equations of motion. Numerical simulation is carried out on the discrete nonlinear dynamic equations of ordinary differential equations to study the relation between critical flutter velocity and system parameters. e results showed that when inflow velocity reached the critical value, the limit cycle occurred. e increment of the aspect ratio or the thickness can both result in the decrement of the critical velocity of flutter. On the contrary, increment of air damping make the critical flutter velocity increase obviously.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.

Authors' Contributions
Li Ma and Minghui Yao contributed equally to this work.