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

Authors: Li Ma; Minghui Yao; Wei Zhang; Kai Lou; Dongxing Cao; Yan Niu
Journal: Shock and Vibration
 (2020)

---

## Abstract

This 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. The aspect ratio, the thickness, and the air damping have significant impact on the critical flutter velocity of the thin plate.


---

## Body

## 1. 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. The 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. The 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]. Thus, 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. Then, 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. The first type is the steady aerodynamic model, namely, assuming the force acting on the lifting surface of the wing do not change with time. The steady aerodynamic model is mainly used for static aeroelastic analysis, such as the thin airfoil theory [6]. The 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 quasi-steady piston theory to calculate the flutter problem of a cantilever plate in supersonic flow. Hu et al. [8–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. Thus, 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. The 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 Theodorsen unsteady aerodynamic theory, researchers investigated stability of wings. However, Cordes et al. [15] pointed out Theodorsen 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 Theodorsen 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. The CFD method is adopted in Refs. [27–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–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. The 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 quasi-steady aerodynamic force for the high-aspect ratio cantilever plate in subsonic flow is induced based on the subsonic thin-airfoil 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 high-aspect 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.


## 2. Aerodynamic Force Derivation for the High-Aspect-Ratio Cantilever Plate in Subsonic Flow

### 2.1. Analysis of the Aerodynamic Force for the High-Aspect-Ratio Cantilever Plate in Subsonic Flow Field

The schematic diagram of the cantilever plate considered is shown in Figure1. The wing is simplified as a high-aspect ratio cantilever plate. The span length, chord length, and thickness of the plate are a, b, and h, respectively. The 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⇀x0 is the spanwise direction, e⇀y0 is the chordwise direction, and e⇀z0 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.Figure 1
Model of the high-aspect-ratio cantilever plate.The 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. The 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. Thus, 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, low-velocity, 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 γY of unit arc length of the wing:(1)dL=ρU∞γYdY.The vortex strengthγ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 γY.The 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γY. The 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. Thus, boundary conditions and pressure coefficient of the airfoil can be linearized. Therefore, based on principle of superposition, the attack angle, the bending, and the thickness can be considered separately and then superposed. The 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. The 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 γY, which has a trigonometric series solution. Lateral displacement of the mean camber line of the chordwise cross section is set up as WY. The chordwise position of an arbitrary point in the middle arc line can be written as(2)Y=b21−cosθ,0≤θ≤π.When the plate vibrates slightly, a very small displacementdWY appears. According to the mathematical definition of limit, the tangent value dWY/dY is equal to the corresponding angle value. Thus, we express the attack angle as dWY/dY, that is, the tangent value of the attack angle. Substitute equation (2) into dWY/dY, then the expression of the attack angle is changed into function Kθ, which is related to WY and θ. Therefore, expression of WY is derived as after mentioned. If the attack angle Kθis integrable, as given in equation (4), the local vortex strength γY can be expressed as γθ:(3)γθ=2U∞A0cotθ2+∑n=1∞Ansinnθ,where(4a)A0=α−1π∫0πKθdθ,(4b)An=2π∫0πKθ⋅cosnθdθ,n=1,2,3,4,….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 γθ. 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 Figure2. Lateral vibration velocity Vw and inflow velocity U∞ will generate a new additional attack angle θ2, which is approximately Vw/U∞, as shown in Figure 3.Figure 2
Additional attack angle of the plate.Figure 3
Equivalent diagram of the additional attack angle of the wing.The additional attack angle of the thin plateθ2 is defined as Vw/U∞, so flow theory of the thin airfoil can be used to calculate the local vortex strength caused by θ2. Substitute equation (2) into Vw/U∞, expression of the additional attack angle is marked as Qθ, which is related to Vw and θ. Thus, expression of Vw 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 γ′θ as follows:(5)γ′θ=2U∞A0′cotθ2+∑n=1∞An′sinnθ,where(6a)A0′=α−1π∫0πQθdθ,(6b)An′=2π∫0πQθ⋅cosnθdθ,n=1,2,3,4,….The total circular rector in the nonviscous flow field must be conserved. Thus, 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. The wake vortex generated has a strong influence on the local vortex strength of airfoils. Therefore, 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.


### 2.2. 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 WY and the vibration velocity distribution function Vw. 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. Thus, a fitting method of the interpolation function is applied to express the deformation function WYand the vibration velocity distribution function Vw. Therefore, the sum of several interpolation functions can satisfy the condition. If these interpolation functions are integrable on θ after substituting Y into the deformation function WY and the distribution functionVw, WY and Vw 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 WY and Vw. 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 functionWX1 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. WX1is expressed as(7)WX1=12wx,b,t+wx,0,t.Figure 4
Deformation of the mean camber line.After the first fitting of the mean camber line’s deformation, difference between the deformation curveCC0 of the mean camber line and the first interpolation function WX1 is CC1, as shown in Figure 5.Figure 5
The first interpolation function on the planeY,Z.Since the second order mode of the cantilever plate vibration is the torsion deformation, the second interpolation functionWX2 is set up as a linear function, which is shown in Figure 6.Figure 6
The second interpolation function on the planeY,Z.Let values of the interpolation function atY=0 and Y=b be equal to values of residual CC1 at Y=0 and Y=b, respectively. Thus, slope of the first interpolation function can be calculated as follows:(8)1bwx,b,t−WX1−wx,0,t−WX1.The second interpolation functionWX2 can be written as follows:(9)WX2=1bwx,b,t−wx,0,t⋅Y+wx,0,t−WX1.Difference between residualCC1 after the first fitting and the second interpolation function WX2 is marked as CC2, which can be expressed as CC2=W−WX1−WX2.The third interpolation functionWX3 is set up as a quadratic function, as shown in Figure 7. Values of CC2 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 K2⋅Yb−Y. The extreme point of the interpolation function appears at Y=b/2; thus, it is necessary to make K2×Yb−YY=b/2=CC2Y=b/2 in order to get a minimal difference between the third interpolation function and the residual CC2. According to the equation, unknown coefficient K2 and the third interpolation function WX3 can be calculated. WX3 is written as follows:(10)WX3=wx,b2,t−Wx14b2⋅Yb−Y.Figure 7
The third interpolation function on the planeY,Z.Difference between the residualCC2 in the second fitting and the third interpolation function WX3 is marked as CC3. The difference can be expressed as CC3=W−WX1−WX2−WX3.The fourth interpolation functionWX4 is chosen as a cubic function, which is shown in Figure 8. Values of CC3 at Y=0, Y=b, and Y=b/2 are 0. 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 K3×Y×b−Y×b/2−Y. Extreme points of the interpolation function occur at Y=b/2−3b/6 and Y=b/2+3b/6, so it is essential to establish the equation K3×Y×b−Y×b/2−YY=b/2−3b/6=CC3Y=b/2−3b/6 in order to make the fourth interpolation function fit the residual CC3to the greatest extent. According to the equation, unknown coefficient K3 and the expression of WX4 can be calculated. WX4 is shown as follows:(11)WX4=wx,y,t−WX1−WX2−WX3Y=b/2−3b/6b/2−3b/6b/2+3b/63b/6×Yb−Yb2−Y.Figure 8
The fourth interpolation function on the planeY,Z.Similarly, difference between the residual CC3 in the third fitting and the fourth interpolation function WX4 is marked as CC4, which can be expressed as CC4=W−WX1−WX2−WX3−WX4.The fifth interpolation functionWX5 is set up as the form of a quartic function, as shown in Figure 9. Values of CC4 at Y=0, Y=b, Y=b/2, and Y=b/2−3b/6 are zero. 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 K4×Y×b−Y×b/2−Y×b/2−3b/6−Y. Difference at extreme points of the interpolation function Y=b/2+3b/6 is still exist, so it is necessary to let K4×Y×b−Y×b/2−Y×b/2−3b/6−YY=b/2−3b/6=CC4Y=b/2−3b/6=CC4Y=b/2+3b/6 in order to make the fifth interpolation function fit the residual CC4 to the greatest extent. According to the equation, unknown coefficient K4 and the fifth interpolation function WX5 can be calculated. WX5 is shown as follows:(12)WX5=wx,y,t−WX1−WX2−WX3−WX4Y=b/2+3b/6b/2+3b/6b/2−3b/6−3b/6b×Yb−Yb2−Yb2−36b−Y.Figure 9
The fifth interpolation function on the planeY,Z.Difference between the residualCC4 in the fourth fitting and the fifth interpolation function WX5 is marked as CC5. The difference can be expressed as CC5=W−WX1−WX2−WX3−WX4−WX5. Values of CC5 at Y=0, Y=b, Y=b/2, Y=b/2−3b/6, and Y=b/2+3b/6 are 0.After five times of deformation fitting abovementioned, the total residual CC5 brought by the mean camber line’s deformation CC0 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.Figure 10
Residual error of the deformation in fitting.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. Therefore, it only needs to replace items about deformation in equations (7) and (9)–(12) with a corresponding item about velocity. Then, five interpolation functions for fitting velocity distribution are obtained as follows:(13)VX1=∂wx,b,t/∂t+∂wx,0,t/∂t2,VX2=1b∂wx,b,t∂t−∂wx,0,t∂tY+∂wx,0,t∂t−VX1,VX3=∂wx,b/2,t∂t−VX14b2b−YY,VX4=∂wx,b,t/∂t−VX1−VX2−VX3Y=b/2−3b/6b/2−3b/6b/2+3b/63b/6×Yb−Yb2−Y,VX5=∂w/∂t−VX1−VX2−VX3−VX4Y=b/2+3b/6b/2+3b/6b/2−3b/6−3b/6b×Yb−Yb2−Yb2−36b−Y.


### 2.3. 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 the previous section, we give the interpolation functions of the two variables. In order to calculate the local vortex strength caused by deformation of the mean camber line, local vortex strength caused by every interpolation function can be calculated, respectively, and then the linear superposition can be carried out. Chordwise deformation of the plate structure is relatively simple, so top four interpolation functions are enough to fit the deformation of the cantilever plate accurately. WX1, WX2, WX3, and WX4 are taken to calculate local vortex. First of all, attack angle caused by four interpolation functions dWX1/dY, dWX2/dY, dWX3/dY, and dWX4/dY are calculated, respectively. Substitute equation (2) into dWX1/dY, dWX2/dY, dWX3/dY, and dWX4/dY; then, functions of the attack angle become related to θ. These four functions are set up as K1, K2, K3, and K4, which are substituted into equations (4a) and (4b), respectively. A0 and A calculated by equations (4a) and (4b) are substituted into equation (3) to solve the corresponding local vortex of top four interpolation functions: γWX1, γWX2, γWX3, and γWX4.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. Thus, top four interpolation functionsVX1, VX2, VX3, and VX4 are selected to calculate the local vortex. Firstly, angle of attack functions VX1/U∞, VX2/U∞, VX3/U∞, and VX4/U∞of four interpolation functions are calculated, respectively. And substitute equation (2) into VX1/U∞, VX2/U∞, VX3/U∞, and VX4/U∞. Then, angle of attack functions are translated into functions related to θ. These four functions are set up as Q1, Q2, Q3, and Q4, which are substituted into equations (6a) and (6b), respectively. A0′ and An′ calculated by equations (6a) and (6b) are substituted into equation (5) to solve the corresponding local vortex: γVX1, γVX2, γVX3, and γVX4 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 γz, which can be expressed as linear superposition of the local vortex strength calculated by abovementioned eight interpolation functions. γz is written as follows:(14)γz=γWX1+γWX2+γWX3+γWX4+γVX1+γVX2+γVX3+γVX4=pw1⋅wx,b,t+pw2⋅wx,0,t+pw3⋅wx,b2,t+pw4⋅wx,b2−3b6,t+pv1⋅∂wx,b,t∂t+pv2⋅∂wx,0,t∂t+pv3⋅∂wx,b/2,t∂t+pv4⋅∂wx,b/2−3b/6,t∂t,where(15a)pw1=123−36Yb+10−63U∞b×Yb−Y2b2+3−72U∞bbY−1,(15b)pw2=−26−63+123+36YbU∞bYb−Y2b2+7+32U∞bbY−1,(15c)pw3=483Yb16−243U∞bYb−Y2b2+23U∞bbY−1,(15d)pw4=363−723YbU∞bYb−Y2b2−33U∞bbY−1,(15e)pv1=12−43Yb−Y2b23/2+3−32−4Yb×Yb−Y2b2−12bY−1,(15f)pv2=−43−12Yb−Y2b23/2−12bY−1+11+32−4YbYb−Y2b2,(15g)pv3=−163Yb−Y2b23/2+8Yb+23−4×Yb−Y2b2−bY−1,(15h)pv4=24Yb−Y2b23/23−3Yb−Y2b23.According to equations (1) and (14), the aerodynamic force Δp can be expressed as(16)Δp=ρU∞γz.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.


### 2.4. Aerodynamic Correction and Error Analysis

Value of itemb/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. Thus, it is necessary to add a correlation coefficient inb/Y−1. The 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 itemb/Y−1 in equation (15a) is corrected asb−Y/Y+0.05b. At the moment, the value of b−Y/Y+0.05b at Y=0 is equal to the value of b/Y−1 at Y=0.95b. The value of these two functions at the trailing edge portion of the plate changes little, as shown in Figure 11. The corrected aerodynamic expression is denoted as Δp′.Figure 11
Correction of the coefficient of the aerodynamic force.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. Thus, 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 Δpp 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Δpp is as follows:(17)Δpp=Δp′1−Ma∞2+1/2Δp′⋅1−1−Ma∞2.Aerodynamic forceΔpp is the linear superposition of aerodynamic forces calculated by several interpolation functions. Moreover, inflow air must satisfy hypotheses of irrotational and nonviscous. Thus, the aerodynamic force Δpp calculated by equation (17) is an approximate result. There 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. The essential reason why the lift force can be generated is to change the attack angle of the wing, which indirectly affects the aerodynamic force. Thus, 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 Figure12. A thin shell model of 0.01 meter thickness is constructed by stretching the spline curve we drew to 10 meters along the spanwise direction. Assuming the elasticity modulus of the shell structure in subsonic, compressible, and nonsticky turbulence flow field is infinite. The velocity of the airflow is U∞=200m/s. After numerical simulation, pressure cloud pictures of the upper and lower surface can be obtained, as shown in Figures 13 and 14.Figure 12
The coordinate graph of the spline curve.Figure 13
Pressure of the lower surface of the airfoil.Figure 14
Pressure of the upper surface of the airfoil.The lift force can be calculated by the pressure difference between the upper and lower surface of the airfoil, as shown in Figure15.Figure 15
Comparison between the theoretical result and the ANSYS result.If we take the spline curve in Figure12 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:(18)wx,0,t=0.376,wx,b2,t=0.083,wx,b2−3b6,t=0.581,wx,b,t=−0.3,wx,b2+3b6,t=−0.45.In calculation of the lift force generated only by deformation of the mean camber line, velocities in the aerodynamic force expressionΔppare taken as zero, namely,(19)∂wx,0,t∂t=0,∂wx,b/2,t∂t=0,∂wx,b/2−3b/6,t∂t=0,∂wx,b,t∂t=0,∂wx,b/2+3b/6,t∂t=0.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.The 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.


## 2.1. Analysis of the Aerodynamic Force for the High-Aspect-Ratio Cantilever Plate in Subsonic Flow Field

The schematic diagram of the cantilever plate considered is shown in Figure1. The wing is simplified as a high-aspect ratio cantilever plate. The span length, chord length, and thickness of the plate are a, b, and h, respectively. The 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⇀x0 is the spanwise direction, e⇀y0 is the chordwise direction, and e⇀z0 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.Figure 1
Model of the high-aspect-ratio cantilever plate.The 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. The 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. Thus, 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, low-velocity, 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 γY of unit arc length of the wing:(1)dL=ρU∞γYdY.The vortex strengthγ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 γY.The 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γY. The 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. Thus, boundary conditions and pressure coefficient of the airfoil can be linearized. Therefore, based on principle of superposition, the attack angle, the bending, and the thickness can be considered separately and then superposed. The 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. The 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 γY, which has a trigonometric series solution. Lateral displacement of the mean camber line of the chordwise cross section is set up as WY. The chordwise position of an arbitrary point in the middle arc line can be written as(2)Y=b21−cosθ,0≤θ≤π.When the plate vibrates slightly, a very small displacementdWY appears. According to the mathematical definition of limit, the tangent value dWY/dY is equal to the corresponding angle value. Thus, we express the attack angle as dWY/dY, that is, the tangent value of the attack angle. Substitute equation (2) into dWY/dY, then the expression of the attack angle is changed into function Kθ, which is related to WY and θ. Therefore, expression of WY is derived as after mentioned. If the attack angle Kθis integrable, as given in equation (4), the local vortex strength γY can be expressed as γθ:(3)γθ=2U∞A0cotθ2+∑n=1∞Ansinnθ,where(4a)A0=α−1π∫0πKθdθ,(4b)An=2π∫0πKθ⋅cosnθdθ,n=1,2,3,4,….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 γθ. 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 Figure2. Lateral vibration velocity Vw and inflow velocity U∞ will generate a new additional attack angle θ2, which is approximately Vw/U∞, as shown in Figure 3.Figure 2
Additional attack angle of the plate.Figure 3
Equivalent diagram of the additional attack angle of the wing.The additional attack angle of the thin plateθ2 is defined as Vw/U∞, so flow theory of the thin airfoil can be used to calculate the local vortex strength caused by θ2. Substitute equation (2) into Vw/U∞, expression of the additional attack angle is marked as Qθ, which is related to Vw and θ. Thus, expression of Vw 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 γ′θ as follows:(5)γ′θ=2U∞A0′cotθ2+∑n=1∞An′sinnθ,where(6a)A0′=α−1π∫0πQθdθ,(6b)An′=2π∫0πQθ⋅cosnθdθ,n=1,2,3,4,….The total circular rector in the nonviscous flow field must be conserved. Thus, 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. The wake vortex generated has a strong influence on the local vortex strength of airfoils. Therefore, 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.


## 2.2. 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 WY and the vibration velocity distribution function Vw. 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. Thus, a fitting method of the interpolation function is applied to express the deformation function WYand the vibration velocity distribution function Vw. Therefore, the sum of several interpolation functions can satisfy the condition. If these interpolation functions are integrable on θ after substituting Y into the deformation function WY and the distribution functionVw, WY and Vw 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 WY and Vw. 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 functionWX1 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. WX1is expressed as(7)WX1=12wx,b,t+wx,0,t.Figure 4
Deformation of the mean camber line.After the first fitting of the mean camber line’s deformation, difference between the deformation curveCC0 of the mean camber line and the first interpolation function WX1 is CC1, as shown in Figure 5.Figure 5
The first interpolation function on the planeY,Z.Since the second order mode of the cantilever plate vibration is the torsion deformation, the second interpolation functionWX2 is set up as a linear function, which is shown in Figure 6.Figure 6
The second interpolation function on the planeY,Z.Let values of the interpolation function atY=0 and Y=b be equal to values of residual CC1 at Y=0 and Y=b, respectively. Thus, slope of the first interpolation function can be calculated as follows:(8)1bwx,b,t−WX1−wx,0,t−WX1.The second interpolation functionWX2 can be written as follows:(9)WX2=1bwx,b,t−wx,0,t⋅Y+wx,0,t−WX1.Difference between residualCC1 after the first fitting and the second interpolation function WX2 is marked as CC2, which can be expressed as CC2=W−WX1−WX2.The third interpolation functionWX3 is set up as a quadratic function, as shown in Figure 7. Values of CC2 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 K2⋅Yb−Y. The extreme point of the interpolation function appears at Y=b/2; thus, it is necessary to make K2×Yb−YY=b/2=CC2Y=b/2 in order to get a minimal difference between the third interpolation function and the residual CC2. According to the equation, unknown coefficient K2 and the third interpolation function WX3 can be calculated. WX3 is written as follows:(10)WX3=wx,b2,t−Wx14b2⋅Yb−Y.Figure 7
The third interpolation function on the planeY,Z.Difference between the residualCC2 in the second fitting and the third interpolation function WX3 is marked as CC3. The difference can be expressed as CC3=W−WX1−WX2−WX3.The fourth interpolation functionWX4 is chosen as a cubic function, which is shown in Figure 8. Values of CC3 at Y=0, Y=b, and Y=b/2 are 0. 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 K3×Y×b−Y×b/2−Y. Extreme points of the interpolation function occur at Y=b/2−3b/6 and Y=b/2+3b/6, so it is essential to establish the equation K3×Y×b−Y×b/2−YY=b/2−3b/6=CC3Y=b/2−3b/6 in order to make the fourth interpolation function fit the residual CC3to the greatest extent. According to the equation, unknown coefficient K3 and the expression of WX4 can be calculated. WX4 is shown as follows:(11)WX4=wx,y,t−WX1−WX2−WX3Y=b/2−3b/6b/2−3b/6b/2+3b/63b/6×Yb−Yb2−Y.Figure 8
The fourth interpolation function on the planeY,Z.Similarly, difference between the residual CC3 in the third fitting and the fourth interpolation function WX4 is marked as CC4, which can be expressed as CC4=W−WX1−WX2−WX3−WX4.The fifth interpolation functionWX5 is set up as the form of a quartic function, as shown in Figure 9. Values of CC4 at Y=0, Y=b, Y=b/2, and Y=b/2−3b/6 are zero. 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 K4×Y×b−Y×b/2−Y×b/2−3b/6−Y. Difference at extreme points of the interpolation function Y=b/2+3b/6 is still exist, so it is necessary to let K4×Y×b−Y×b/2−Y×b/2−3b/6−YY=b/2−3b/6=CC4Y=b/2−3b/6=CC4Y=b/2+3b/6 in order to make the fifth interpolation function fit the residual CC4 to the greatest extent. According to the equation, unknown coefficient K4 and the fifth interpolation function WX5 can be calculated. WX5 is shown as follows:(12)WX5=wx,y,t−WX1−WX2−WX3−WX4Y=b/2+3b/6b/2+3b/6b/2−3b/6−3b/6b×Yb−Yb2−Yb2−36b−Y.Figure 9
The fifth interpolation function on the planeY,Z.Difference between the residualCC4 in the fourth fitting and the fifth interpolation function WX5 is marked as CC5. The difference can be expressed as CC5=W−WX1−WX2−WX3−WX4−WX5. Values of CC5 at Y=0, Y=b, Y=b/2, Y=b/2−3b/6, and Y=b/2+3b/6 are 0.After five times of deformation fitting abovementioned, the total residual CC5 brought by the mean camber line’s deformation CC0 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.Figure 10
Residual error of the deformation in fitting.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. Therefore, it only needs to replace items about deformation in equations (7) and (9)–(12) with a corresponding item about velocity. Then, five interpolation functions for fitting velocity distribution are obtained as follows:(13)VX1=∂wx,b,t/∂t+∂wx,0,t/∂t2,VX2=1b∂wx,b,t∂t−∂wx,0,t∂tY+∂wx,0,t∂t−VX1,VX3=∂wx,b/2,t∂t−VX14b2b−YY,VX4=∂wx,b,t/∂t−VX1−VX2−VX3Y=b/2−3b/6b/2−3b/6b/2+3b/63b/6×Yb−Yb2−Y,VX5=∂w/∂t−VX1−VX2−VX3−VX4Y=b/2+3b/6b/2+3b/6b/2−3b/6−3b/6b×Yb−Yb2−Yb2−36b−Y.


## 2.3. 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 the previous section, we give the interpolation functions of the two variables. In order to calculate the local vortex strength caused by deformation of the mean camber line, local vortex strength caused by every interpolation function can be calculated, respectively, and then the linear superposition can be carried out. Chordwise deformation of the plate structure is relatively simple, so top four interpolation functions are enough to fit the deformation of the cantilever plate accurately. WX1, WX2, WX3, and WX4 are taken to calculate local vortex. First of all, attack angle caused by four interpolation functions dWX1/dY, dWX2/dY, dWX3/dY, and dWX4/dY are calculated, respectively. Substitute equation (2) into dWX1/dY, dWX2/dY, dWX3/dY, and dWX4/dY; then, functions of the attack angle become related to θ. These four functions are set up as K1, K2, K3, and K4, which are substituted into equations (4a) and (4b), respectively. A0 and A calculated by equations (4a) and (4b) are substituted into equation (3) to solve the corresponding local vortex of top four interpolation functions: γWX1, γWX2, γWX3, and γWX4.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. Thus, top four interpolation functionsVX1, VX2, VX3, and VX4 are selected to calculate the local vortex. Firstly, angle of attack functions VX1/U∞, VX2/U∞, VX3/U∞, and VX4/U∞of four interpolation functions are calculated, respectively. And substitute equation (2) into VX1/U∞, VX2/U∞, VX3/U∞, and VX4/U∞. Then, angle of attack functions are translated into functions related to θ. These four functions are set up as Q1, Q2, Q3, and Q4, which are substituted into equations (6a) and (6b), respectively. A0′ and An′ calculated by equations (6a) and (6b) are substituted into equation (5) to solve the corresponding local vortex: γVX1, γVX2, γVX3, and γVX4 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 γz, which can be expressed as linear superposition of the local vortex strength calculated by abovementioned eight interpolation functions. γz is written as follows:(14)γz=γWX1+γWX2+γWX3+γWX4+γVX1+γVX2+γVX3+γVX4=pw1⋅wx,b,t+pw2⋅wx,0,t+pw3⋅wx,b2,t+pw4⋅wx,b2−3b6,t+pv1⋅∂wx,b,t∂t+pv2⋅∂wx,0,t∂t+pv3⋅∂wx,b/2,t∂t+pv4⋅∂wx,b/2−3b/6,t∂t,where(15a)pw1=123−36Yb+10−63U∞b×Yb−Y2b2+3−72U∞bbY−1,(15b)pw2=−26−63+123+36YbU∞bYb−Y2b2+7+32U∞bbY−1,(15c)pw3=483Yb16−243U∞bYb−Y2b2+23U∞bbY−1,(15d)pw4=363−723YbU∞bYb−Y2b2−33U∞bbY−1,(15e)pv1=12−43Yb−Y2b23/2+3−32−4Yb×Yb−Y2b2−12bY−1,(15f)pv2=−43−12Yb−Y2b23/2−12bY−1+11+32−4YbYb−Y2b2,(15g)pv3=−163Yb−Y2b23/2+8Yb+23−4×Yb−Y2b2−bY−1,(15h)pv4=24Yb−Y2b23/23−3Yb−Y2b23.According to equations (1) and (14), the aerodynamic force Δp can be expressed as(16)Δp=ρU∞γz.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.


## 2.4. Aerodynamic Correction and Error Analysis

Value of itemb/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. Thus, it is necessary to add a correlation coefficient inb/Y−1. The 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 itemb/Y−1 in equation (15a) is corrected asb−Y/Y+0.05b. At the moment, the value of b−Y/Y+0.05b at Y=0 is equal to the value of b/Y−1 at Y=0.95b. The value of these two functions at the trailing edge portion of the plate changes little, as shown in Figure 11. The corrected aerodynamic expression is denoted as Δp′.Figure 11
Correction of the coefficient of the aerodynamic force.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. Thus, 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 Δpp 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Δpp is as follows:(17)Δpp=Δp′1−Ma∞2+1/2Δp′⋅1−1−Ma∞2.Aerodynamic forceΔpp is the linear superposition of aerodynamic forces calculated by several interpolation functions. Moreover, inflow air must satisfy hypotheses of irrotational and nonviscous. Thus, the aerodynamic force Δpp calculated by equation (17) is an approximate result. There 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. The essential reason why the lift force can be generated is to change the attack angle of the wing, which indirectly affects the aerodynamic force. Thus, 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 Figure12. A thin shell model of 0.01 meter thickness is constructed by stretching the spline curve we drew to 10 meters along the spanwise direction. Assuming the elasticity modulus of the shell structure in subsonic, compressible, and nonsticky turbulence flow field is infinite. The velocity of the airflow is U∞=200m/s. After numerical simulation, pressure cloud pictures of the upper and lower surface can be obtained, as shown in Figures 13 and 14.Figure 12
The coordinate graph of the spline curve.Figure 13
Pressure of the lower surface of the airfoil.Figure 14
Pressure of the upper surface of the airfoil.The lift force can be calculated by the pressure difference between the upper and lower surface of the airfoil, as shown in Figure15.Figure 15
Comparison between the theoretical result and the ANSYS result.If we take the spline curve in Figure12 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:(18)wx,0,t=0.376,wx,b2,t=0.083,wx,b2−3b6,t=0.581,wx,b,t=−0.3,wx,b2+3b6,t=−0.45.In calculation of the lift force generated only by deformation of the mean camber line, velocities in the aerodynamic force expressionΔppare taken as zero, namely,(19)∂wx,0,t∂t=0,∂wx,b/2,t∂t=0,∂wx,b/2−3b/6,t∂t=0,∂wx,b,t∂t=0,∂wx,b/2+3b/6,t∂t=0.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.The 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.


## 3. Establishment and Dispersion of Aerodynamic Equations for the Cantilever Plate

The schematic diagram of the wing is shown in Figure1, which is simplified as a cantilever plate in subsonic flow. The 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:(20a)ux,y,z,t=u0x,y,t+zφxx,y,t,(20b)vx,y,z,t=v0x,y,t+zφyx,y,t,(20c)wx,y,z,t=w0x,y,t,where u0, v0, and w0 are displacements along X, Y, and Z directions of points on the neutral plane of the blade, and φx and φy are rotation angles around y-axis and x-axis.Strains of the von Karman plate, which are computed by using the nonlinear strain-displacement relation, are obtained as follows:(21a)εxx=∂u0∂x+12∂w∂x2+z∂φx∂x=εxx0+zεxx1,(21b)εyy=∂v0∂y+12∂w∂y2+z∂φy∂y=εyy0+zεyy1,(21c)γxy=∂u0∂y+∂v0∂x+∂w∂x∂w∂y+z∂φx∂y+∂φy∂x=γxy0+zγxy1,(21d)γxz=φx+∂w∂xγyz=φy+∂w∂y.The constitutive relation of the isotropic material is expressed as follows:(22a)σxx=E1−υ2εxx+υεyy,(22b)σyy=E1−υ2εyy+υεxx,(22c)τxy=E21+υγxy,(22d)τxz=E21+υγxz,(22e)τyz=E21+υγyz,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:(23a)QxQy=∫−h/2h/2τxzτyzdz,(23b)NxNyNxy=∫−h/2h/2σxxσyyτxydz,(23c)MxMyMxy=∫−h/2h/2σxxσyyτxyzdz.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:(24)Ud=12∫ΩNxxεxx0+Nyyεyy0+Nxyγxy0+Mxεxx1+Myεyy1+Mxyεxy1+Qxzγxz+Qyzγyzdxdy.The kinetic energy and the virtual work done by external forces are as follows:(25)T=12∫Ω∫−h/2h/2ρcu˙2+v˙2+w˙2dxdydz,(26)WP=∫ΩΔP⋅w−c⋅dWdtdxdy,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:(27)∫t1t2δT−δUd+δWpdt=0.The lateral displacement of the thin plate is more obvious compared with others. Thus, displacementsu0 and v0 on the middle surface and rotation angle φx and φyare neglected. Substituting equations (24)–(26) into equation (27), nonlinear dynamic equation of the plate in the lateral direction is obtained as follows:(28)112π2E∂/∂xϕxx,y,t+∂2/∂x2wx,y,th2+2υ+112π2E∂/∂yϕyx,y,t+∂2/∂y2wx,y,th2+2υ+Eυh∂/∂xwx,y,t∂/∂ywx,y,t∂2/∂y∂xwx,y,t−υ2+1+Eh∂/∂xwx,y,t∂/∂xwx,y,t∂2/∂x2wx,y,t−υ2+1+1/2υ∂/∂ywx,y,t2hE∂2/∂x2wx,y,t−υ2+1+1/2∂/∂xwx,y,t2hE∂2/∂x2wx,y,t−υ2+1+π2E∂/∂xwx,y,t2∂2/∂y2wx,y,th12×2+2υ+υ∂/∂xwx,y,t∂2/∂y∂xwx,y,tEh∂/∂ywx,y,t−υ2+1+∂/∂ywx,y,t∂2/∂y2wx,y,tEh∂/∂ywx,y,t−υ2+1+1/2υ∂/∂xwx,y,t2Eh∂2/∂y2wx,y,t−υ2+1+1/2∂/∂ywx,y,t2Eh∂2/∂y2wx,y,t−c∂/∂twx,y,t+Δp=ρch∂2/∂t2wx,y,t.Ii is calculated by(29)I0I1I2=∫−h/2h/2ρc1zz2dz.The boundary conditions of the fix end of the plate are(30a)x=0w=0,φx=0,φy=0,(30b)y=0w=0,φx=0,φy=0.The boundary conditions of the free end of the plate are(31a)x=aw=0,Mx=0,Mxy=0,Qx=0,(31b)y=bw=0,My=0,Mxy=0,Qy=0.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(32a)XM1=coshk1⋅x−cosk1⋅x−β1sinhk1⋅x−sink1⋅x,(32b)XM2=coshk2⋅x−cosk2⋅x−β2sinhk2⋅x−sink2⋅x,(32c)YM1=1,(32d)YM2=31−2Yb,where k1, k2, β1, and β2 are calculated by(33)k1=1.875a,k2=4.694a,β1=1.875a,β2=sinhk1⋅a−sink1⋅acoshk1⋅a+cosk1⋅a.The lateral displacement of the plate can be expressed as follows:(34)w=x1⋅XM1⋅YM1+x3⋅XM1YM2.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:(35a)x¨1=−1470740.224x13+575174.9082x32−113819521.7x32x1+0.01384683053x33−74.29541441x1+0.3027136145U∞2x3+0.001143287344x12x3+5383.891352x12+0.04079⋅U∞⋅x˙3−0.02857⋅2⋅c⋅x˙1,(35b)x¨3=0.0005710592x13+0.0004853467x32−0.08089898x32x1−110064645x33−11668.03765x3+0.2032087116U∞2x3−113739904.7x12x3+1149545.143x3x1−0.02855⋅2⋅c⋅x˙3+0.04076⋅U∞⋅x˙1.Table 1
Value of parameters.
Physical quantityNumerical valuea5 mb1 mh0.02 mρ1.28 kg/m3E102 GPaρc1750kg/m3c10N⋅s/mυ0.3In equations (35a) and (35b), x1 and x3 are generalized coordinates of the first-order and the second-order modals. The 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. Thus, amplitude of the vibrating structure increases continuously by absorbing energy from the airflow, which may lead to the occurrence of flutter.


## 4. 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 x1=0.001, x3=0.001, x˙1=0, and x˙3=0.When the inflow velocityU∞ is chosen as 72m/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 the plane w,v, where v is the transverse vibration velocity. The adjacent trajectory gradually becomes a closed loop, as shown in Figure 16(b).Figure 16
Stable state of the system when the inflow velocity U∞ is chosen as 72 m/s: (a) time history diagram of vibration amplitude on plane t,x1 and (b) phase portrait on plane w,v.
(a)(b)When the inflow speedU∞ is set as 80.43m/s, the system keeps in a critical state between stability and instability, as shown in Figure 17(a). The amplitude of the plate at X=5 and Y=0.8converges to a constant. The 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). The phase portrait is an isolated closed loop, which do not have closed rails nearby. This is called a limit cycle. At the moment, the critical flutter velocity is reached.Figure 17
The system is kept in a critical state of stable and unstable when the inflow velocity U∞ is chosen as 80.43 m/s: (a) time history diagram of vibration amplitude on plane t,x1 and (b) phase portrait on plane w,v.
(a)(b)When the inflow velocityU∞ is increased from 80.43m/s to 89m/s gradually, the system appears to diverge and becomes instable, as shown in Figure 18(a). The 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). The phase portrait is far from a closed loop.Figure 18
A divergent and unstable state of the system is kept in a critical state of stable and unstable when the inflow velocity U∞ is chosen as 89 m/s: (a) time history diagram of vibration amplitude on plane t,x1 and (b) phase portrait on plane w,v.
(a)(b)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. Thus, when the inflow velocity reaches a certain value, there will be periodic vibration of equal amplitudes, as shown in Figure 17. This 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 as80.43m/s, the system vibrates periodically with equal amplitudes. When the inflow velocity increases to a value smaller than 89m/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 velocityU∞ 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). The critical flutter velocity increases as the thickness of the plate increases.Figure 19
Influences of different factors on critical flutter velocities: (a) plate thickness; (b) aspect ratio; (c) damping.
(a)(b)(c)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 Table1. The 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 Table1, 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.


## 5. 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. The 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. The results showed that when inflow velocity reached the critical value, the limit cycle occurred. The 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.



---

*Source: 8841590-2020-06-19.xml*