The Effect of Structural Parameters on Response Characteristics of Aeroelastic System in the Presence of Dynamic Stall

The eﬀect of structural paramters on the response and aerodynamic stiﬀness characteristics of the free aeroelastic system under the inﬂuence of dynamic stall is investigated adopting CFD (Computational Fluid Dynamics) method. The equilibrium angle of the spring and the structural stiﬀness are taken as parameters of interest. Systems with small equilibrium angles enter the symmetric limit-cycle state more quickly after a Hopf bifurcation and experience dynamic stall in both directions, rather than slowly decreasing in minimum angle of attack and remaining in the asymmetric limit-cycle state before dynamic stall in the opposite direction, as is the case with systems with large spring equilibrium angles. Thus, aerodynamic stiﬀness of system with large equilibrium angles can be more signiﬁcantly inﬂuenced by the change in aerodynamic moment characteristics at the minimum angle of attack. Furthermore, by increasing the initial angular velocity, we ﬁnd that the system response all becomes symmetric limit cycle and therefore the aerodynamic stiﬀness appears to have a monotonically increasing characteristic. As to the eﬀect of structural stiﬀness, it is found that the limit cycle amplitude ﬁrst increases with structural stiﬀness after bifurcation, then the amplitude is unchanged with varying structural stiﬀness at higher Mach number. Energy maps show that the parametric distribution of the energy transfer contributes to this phenomenon. Moreover, when entering the symmetric limit cycle state, the structural stiﬀness no longer has a signiﬁcant eﬀect on the aerodynamic stiﬀness of the system, as the increase The eﬀect of structural parameters on aeroelastic system in the aerodynamic stiﬀness is determined solely by the increase in dynamic pressure without the eﬀect of changes in moment characteristics.


Abstract
The effect of structural paramters on the response and aerodynamic stiffness characteristics of the free aeroelastic system under the influence of dynamic stall is investigated adopting CFD (Computational Fluid Dynamics) method. The equilibrium angle of the spring and the structural stiffness are taken as parameters of interest. Systems with small equilibrium angles enter the symmetric limit-cycle state more quickly after a Hopf bifurcation and experience dynamic stall in both directions, rather than slowly decreasing in minimum angle of attack and remaining in the asymmetric limit-cycle state before dynamic stall in the opposite direction, as is the case with systems with large spring equilibrium angles. Thus, aerodynamic stiffness of system with large equilibrium angles can be more significantly influenced by the change in aerodynamic moment characteristics at the minimum angle of attack. Furthermore, by increasing the initial angular velocity, we find that the system response all becomes symmetric limit cycle and therefore the aerodynamic stiffness appears to have a monotonically increasing characteristic. As to the effect of structural stiffness, it is found that the limit cycle amplitude first increases with structural stiffness after bifurcation, then the amplitude is unchanged with varying structural stiffness at higher Mach number. Energy maps show that the parametric distribution of the energy transfer contributes to this phenomenon. Moreover, when entering the symmetric limit cycle state, the structural stiffness no longer has a significant effect on the aerodynamic stiffness of the system, as the increase

Introduction
In an aeroelastic system such as a helicopter rotor system, the unsteady flow field and the flexible structure are coupled much tightly and strongly. The complex response of this system will aggravate the damage to the structure and greatly affect the performance of the aircraft. Therefore, aeroelasticity is an important nonlinear dynamics problem in aircraft design. Among the tremendous flow phenomena in unsteady flow fields for aircraft, dynamic stall is one of the most complex phenomena and highly non-linear. In the event of dynamic stall, the stall angle will be delayed compared to the static situation. In addition to this, the formation, convection and shedding of large scale leading edge vortex can occur during the dynamic stall [1], which results in lift overshoot and a sharp drop after vortex shedding. In addition, as the airfoil is pitching down, the flow will recover from the massive separation state to the unstalled state through a flow reattachment process. This also creates a hysteresis in aerodynamic moment and lift. The non-linearity of the aerodynamic loads due to these complex flow phenomena as well as the structural non-linearity [2] contribute to the complexity of the aeroelastic response.
Some researchers focus on the conditions under which stall-induced oscillations can occur. As the general structural damping only plays a dissipative role, the onset of limit cycle oscillation is linked to the work done by aerodynamic forces in one cycle [3]. Carta et al. [4] and Oates et al. [5] derived the aerodynamic damping coefficient to determine the stability of single-degreeof-freedom system. If negative aerodynamic damping is present in the system, self-sustained limit cycle oscillations may occur. Bowles et al. [6] proposed a "time-resolved" aerodynamic damping coefficient concept based on the Hilbert transform to analyze the link between instantaneous flow changes and aerodynamic damping coefficients during the pitching cycle. Bhat et al. [7] returned to the concept of work done and energy transfer and use data on energy transfer during forced oscillations of the airfoil for different parameters to explore which conditions are possible for self-sustained aeroelastic oscillations. On this basis, Menon et al [8] used energy transfer maps between flow and airfoil systems to understand and predict the dynamics of free-oscillating aeroelastic systems. On the other hand, the aerodynamic nonlinearity can have a great effect on the characteristics of aeroelastic responses. Gilliatt et al. [9] used a quasi-steady aerodynamic model containing a cubic nonlinear term to analyze the two-degree-of-freedom aeroelastic response and explored the effect of aerodynamic non-linearity on internal resonance using a multi-scale method.
Vasconcellos et al. [10] considered both the effects of structural and aerodynamic non-linearities. They found that these two non-linearities may dominate in different free-stream velocity intervals. Dimitriadis et al. [11] investigated the bifurcation behavior of aeroelastic responses experimentally in the presence of non-linearity due to dynamic stall. Two bifurcations occurred in their study, the first of which resulted in symmetric limit-cycle oscillations (accompanied by dynamic stalls in both directions), while the second bifurcation allowed separation to occur only on one side of the wing, resulting in an asymmetric limit cycle oscillation. Similar limit-cycle phenomena dominated by dynamic stall in different directions are also present in this paper.
In addition to affecting the amplitude characteristics of the aeroelastic system, aerodynamic non-linearity also affects the frequency characteristics of the system through aerodynamic stiffness effects. The effect of aerodynamic stiffness was investigated in free-oscillation aeroelastic systems at low-to-moderate Reynolds numbers in Ref. [12]. In their study, the appearance of laminar trailing edge separation or laminar separation bubbles has changed the aerodynamic stiffness characteristics and thus the frequency characteristics of the system. While in the present paper, the dynamic stall accompanying the large scale turbulent separation is then the dominant factor in the variation of the aerodynamic stiffness effect. Actually, the aerodynamic stiffness effect and the subharmonic resonances phenomenon was investigated in the previous studies ([13] and [14]) of the authors. In these studies, it was found that there is a non-monotonic variation in the aerodynamic stiffness of the system as the free-stream flow velocity increases. The different trends in aerodynamic stiffness characteristics are thought to be related to the fact that the minimum AOA (Angle Of Attack) of the system is in different regions, as the aerodynamic moments can vary considerably in different oscillation ranges. These characteristics are revisited in Section 4.1 in this paper. However, these studies are only based on a few specific parameters and the characteristics of the change in aerodynamic stiffness due to the influence of dynamic stall need to be investigated in more depth under a wider range of different conditions.
In present paper, a parametric study of the effect of the structural parameters on the aeroelastic response and aerodynamic stiffness is conducted to fully explore the mechanism by which dynamic stall affects the change in aerodynamic stiffness under different conditions. A CFD solver is coupled with the aeroelastic equation to obtain the aeroelastic responses, or simulating the forced oscillations with given amplitudes and frequency in order to construct the energy map to understand the amplitude response of the aeroelastic system.
The paper is organized as follows: Section 2 introduces the numerical method, including the aeroelastic model and CFD solver. Section 3 describe two analysis methods, which are the equivalent linearization analysis and method of energy map. Section 4 first introduce the aerodynamic stiffness characteristics of system with dynamic stall, then discusses the effect of equilibrium angle and structural stiffness coefficients of the spring. The effect of initial condition is also investigated and merged into Section 4.2 because it has a similar mechanism to the equilibrium angle which affects the aerodynamic stiffness. Section 5 concludes the paper with the final discussions 2 Numerical method 2.1 Aeroelastic model This paper investigates the aeroelastic response of an airfoil pitching around a fixed spring equilibrium angle θ 0 . As shown in Figure 1, the angle of attack is α, which is the angle between chord line and horizontal line, while the torsional angle is φ , which is the difference angle between the angle of attack α and the spring equilibrium angle θ 0 . Therefore, the moment generated by the torsional spring is K φ φ. The aerodynamic moment exerted by the surrounding fluid on the airfoil is M . The free-stream velocity is U ∞ .
According to Lagrangian equation, the motion equation of the system is given by Where C φ and K φ are damping coefficient and stiffness coefficient of the spring, respectively. I α is the moment of inertia of the system. M is the aerodynamic moment, which can be obtained from a CFD solver of ANSYS FLUENT 17.2.
Crank-Nicolson method is employed here as the time advance method, which is unconditionally stable and second-order in time. The method is not limited in the selection of time step and relatively simple for solving compared to Runge-Kutta method. According to the method, the expressions of torsional angle and its derivatives are given by Rearranging the equation 2, the first and second derivatives are given bẏ In this case, the first and second derivatives of the responses can be expressed by the value at old time step and the responses at current time which is φ n+1 . At time n+1, substituting equations 3 into 1, the following equation are obtained The specific terms of the effective stiffness are given by And the effective loads are given by Finally, we can solve the responses φ by equation 4 coupled with the CFD solver. This is achieved through the user define function (UDF) files of ANSYS FLUENT, where the CFD solver provides aerodynamic lift and moment for calculating Eq. 4 while the obtained aeroelastic responses data are sent to the CFD solver to update the motion of the airfoil.

CFD Solver Setting and Validation
The computational domain is a circle zone as shown in Fig. 2, and the NACA0012 airfoil lies in the center of the circle. The O-type grid is adopted to discrete the flow domain. The independence of the results from the grid resolution was examined with different grid sizes. The final grid utilized 20328 cells, with 169 nodes around the airfoil while the location of the first row of cells keep y + ≤ 1.
The whole flow domain undergoes the same rigid motion as the airfoil which is specified by the UDF file according to Eq. 4. The incompressible Navier-Stokes equations are numerically solved using the pressure-based solver of FLUENT 17.2. For the boundary conditions, the no-slip wall and the Pressurefar-field condition are applied on the airfoil and the far field, respectively. The Pressure-Implicit with Splitting of Operators (PISO) algorithm is selected as the pressure-velocity coupled scheme. The second order upwind method is adopted for spatial discretization and first order implicit method is adopted for transient formulation. The Transition Shear Stress Transport (SST) model is selected as the turbulent model. The Transition SST model is based on the coupling of the SST k −ω transport equations with two other transport equations, one for the intermittency and one for the transition. It is shown in Ref. [15] that this transition model overtakes the Spalart-Allmaras model in predicting the complex variation of the boundary layer in reattachment process.
The convergence criteria for the residuals were O(10 −6 ) in magnitude. There are 236 time-steps per forcing period, and 50 inner iterations are used for convergence of each time step. In the cases of grid independent test, four oscillation periods are calculated in each case and the last period is selected for comparison. In aeroelastic computation, at least 35 forcing periods are calculated for each case. The calculation time will be extended to ensure the aeroelastic response enters the steady state and repeats several cycles if necessary.
Grid independent test was performed, and three computational domains with different grid resolutions were selected (coarse, fine, and very fine grids). These domains have G1=1.11 × 10 4 , G2=2.03 × 10 4 , G3=4.15 × 10 4 cells corresponding to 123, 169, and 259 nodes around the airfoil, respectively. Fig.  3 is provided to examine the grid sensitivity of the computational domain (Re=1 × 10 6 , Ma=0.072 and k=0.1). As shown in Fig. 3, except for the result of coarse grid G1, the lift coefficient of fine grid G2 is very close to that of very fine grid G3. Moreover, the experiment result [16] is also provided for comparison. The results of G2 and G3 show good agreement with the experiment even in the reattachment process of dynamic stall. 3 Analysis methods

Equivalent linearization analysis on aerodynamic stiffness
In Ref. [14], an analysis method is used to explain qualitatively the reasons for the different aerodynamic stiffness characteristics caused by the different aerodynamic moment characteristics. This method is based on equivalent linearization method in Ref. [17]. First, we assume that the aeroelastic response is a pure sinusoid of the form the nonlinear aerodynamic moment is given by Following the Harmonic Balance methodology, this function can be expanded as a first order Fourier series so that where Since the moment is on the right hand size of Eq.2, move it to the left of the equation, then the equivalent linear aerodynamic stiffness is given by in practical computation process of b 1 , the term sin(ωt) is replaced with α−A0 A , and the integration is performed in the last period of the simulation result.
Since the aerodynamic stiffness based on the equivalent linearization is the −M (α,α) sin(ωt)dt, we could focus on the variation of integrand −M (α,α) sin(ωt) in a cycle to elucidate the link between the trend in aerodynamics stiffness and the variation in moment characteristics. In Ref. [14], by analyzing the variation of this integrated term over a period with the Mach number at different stages, the causes of the trend in the different stages of aerodynamic stiffness are qualitatively explained. In this paper, this method is also used to analyze the differences in aerodynamic stiffness due to differences in moment characteristics for different spring structural parameters.

Energy maps
This method was proposed by Ref. [8]. Multiplication ofα to equation 1 as follows: results in a simple energy balance describing the growth of energy in the oscillator ⇒ d dt where e f is the energy extracted by the airfoil from the flow and e d is the energy lost to structural damping. Integrating this equation over an oscillation cycle with time period T , where the beginning of the cycle is at time t n , results in an equation for the change of amplitude over the oscillation cycle where E is the net energy gain of the structure over one cycle. Assuming a purely sinusoidal form for α andα , i.e., α = A sin ωt , results in the following equation for the change in amplitude over an oscillation cycle Equilibrium is attained when A(t n + T ) = A(t n ), and this corresponds to When structural damping coefficients is zero, the condition is changed to The nondimensional energy transfer coefficient is is the moment coefficients. ρ is the fluid density, S is the reference area of the wing and c is the chord length of the airfoil. In order to understand the amplitude response of the dynamical system as a function of the various parameters, a map of energy transfer as a function of the parameters of interest can be established for analysis. First, the process starts with conducting simulations of airfoils that are forced to undergo sinusoidal pitching oscillations over a range of amplitudes and frequencies, while holding the other parameters such as θ 0 and Ma fixed. The AOA and moment data obtained in the simulations are used to calculated the energy transfer coefficients. Then, the condition E * f = 0 can be used to determine the presence or absence of limit-cycle oscillations with a particular combination of amplitude and frequency. The validity of these energy maps rests primarily on the assumption that the free oscillations are sinusoidal. Although the free oscillations of the highly nonlinear system in this paper are not always strictly sinusoidal under all parameters, the energy maps could still provide a general trend with some accuracy.

Results and discussions 4.1 The aerodynamic stiffness revisited
In Ref. [14], aerodynamic stiffness in free aeroelastic system with low structural stiffness coefficient was investigated. It was found that the aerodynamic stiffness increases first, then decreases and finally increases with the free-stream velocity. In this paper, we first reproduce this aerodynamic stiffness characteristic by performing calculations using the following parameters shown in Table Table 1 Parameters for reproducing the aerodynamic stiffness characteristics of low structural stiffness systems.
1. In the actual calculation, we set the parameter as the standard sea level Mach number corresponding to these velocities. Here the structural stiffness coefficient was set to be slightly lower than that in Ref. [14] to account for the generality of this aerodynamic stiffness characteristics.
Since the natural frequency of the system ω n = (K φ + K a )/I α (K a is the aerodynamic stiffness coefficients), it can well represent the change of aerodynamic stiffness. Firstly, in the upper part of Figure 4, it can be noticed that three intervals are divided according to the trend of the system frequency. The minimum angle of attack for the response is also drawn in the lower part of Figure 4 to allow a clearer comparison of where the minimum angle of attack lies for different intervals. It can be observed that the two critical points between the three intervals are the zero degree angle and the stall angle in the opposite direction. Crossing both of these critical points implies significant changes in aerodynamic moment characteristics around the minimum AOA, which could lead to different characteristics of aerodynamic stiffness. The typical aerodynamic moment coefficients in these three intervals are shown in Figure 4 as black curves. The effect of the variation of moment characteristics on the aerodynamic stiffness is qualitatively analyzed in Ref. [14] using equivalent linearization analysis and is not repeated here.
It is worth noting that due to the presence or absence of dynamic stall in the opposite direction, two kinds of limit-cycle responses with different characteristics appear in the first two stages and in the third stages. We compare the phase plane plots, lift characteristics and moment characteristics of the two cases with different limit-cycle properties in Fig. 5, and we can find that the limit cycle response changes from an asymmetric type in the first two stages to a symmetric type in the third stage due to the appearance of the dynamic stall in the opposite direction.

The effect of equilibrium angle of the spring
This paper first discusses the effect of the equilibrium angle of the spring on the aerodynamic stiffness characteristics of the aeroelastic system. In fact, the spring equilibrium position somehow determines the oscillation range of the angle of attack of the airfoil. If the influence of aerodynamic forces on the aeroelastic system is relatively small, the angle of attack is usually limited to near the equilibrium position due to the spring's restoring force. And based on the analysis in the previous section, we know that the oscillation range has a significant effect on the aerodynamic moment characteristics of the airfoil, which also changes the aerodynamic stiffness characteristics of the system. In  Table 1. Figure 6 shows different response characteristics of system with different equilibrium angles and free-stream Mach number. In Fig. 6, if the system response eventually decays to a fixed value, it is represented by a circle, and if the system exhibits limit-cycle oscillations, it is represented by a square symbol. Based on the classification of limit cycles with and without symmetry presented in the previous section, we also use two colors to distinguish between the two kinds of LCOs.
Usually, when the free-stream velocity is low, the system has not yet undergone a Hopf bifurcation, and the presence of fixed-point response is normal. However, we can find that, when the springs equilibrium angle of attack is 10 and 25 degrees, the fixed-point response also occurs at higher free-stream velocities (red circles in Fig.6), even after the limit cycle oscillation has already occurred. The trend of the response with free-stream Mach number for one of the systems that contain this special fixed point is shown in Fig.7. In Fig. 7, when the springs equilibrium AOA is 10 degrees, a Hopf bifurcation occurs at Ma = 0.03 as the free-stream Mach number increases, followed by limit cycle oscillation responses. However, at Ma=0.15, the system reaches a new equilibrium of about -13 degrees different from the previous one of about 10 degrees. That means α = −13.17 • is a stable fixed point when Ma=0.15.
Let us compare the two special fixed points of the system with different equilibrium AOAs. In Fig. 8, the time series, phase plane plot and moment coefficients of the two cases are shown. It can be observed that although the transient of the second case(θ 0 = 25 • , Ma=0.19) had approached the symmetric limit cycle, both eventually converge to a similar fixed point. These fixed points are located at the minimum angle of attack and the rate of change of the moment coefficient relative to the angle of attack is small, which facilitates the system to reach equilibrium.
Here we assume that this particular fixed point at different spring equilibrium angles is constant as α 0 and the moment coefficient at equilibrium is C m0 . Based on the assumptions above, it is possible to find the critical Mach number Ma c at which this fixed-point response can be obtained for other spring equilibrium angles. To achieve equilibrium, it is necessary to satisfy We could get Then the value of k and α 0 can be obtained by substituting the two sets of parameters obtained earlier above that make stable fixed points exist into equation 20. The purple curve satisfying equation 20 is plotted in Figure 9 on the basis of Figure 6. Around the curve, we found several parameters (purple circle in Fig.9) with which the simulations were able to find stable fixed points. Surprisingly, this curve is roughly the dividing line between symmetric and asymmetric limit-cycle response areas and also between the two areas with different aerodynamic stiffness characteristics.
Next, we analyze the effect of spring equilibrium angle on the aeroelastic response and aerodynamic stiffness characteristics of the system by comparing several sets of calculations. First, two more spring equilibrium angles that do not differ much from the 15 degrees in section 4.1 are selected for comparison, which are 10 degrees and 20 degrees.
The response frequencies of the system for different free-stream Mach numbers at three spring equilibrium angles are given in Figure 10 Fig. 9 The parameter curve that allows the existence of a negative stable fixed point.
characterizes the variation of the aerodynamic stiffness. Also Figure 10(b) gives the amplitude variation of the system, where the three data points near the xaxis correspond to the appearance of the special fixed points mentioned above. Figure 10(a) shows that the variation of the aerodynamic stiffness characteristics is similar for the three spring equilibrium angles of attack, all of which exhibit three phases, ie., growth, reduction and then growth again. At lower Mach numbers, the natural frequencies of the three are the same, because the aerodynamic stiffness has less influence at this stage and it is mainly the structural stiffness that dominates the frequency. It is also observed that the larger the spring balance angle, the later the system enters the third phase. As analyzed above, the entry into the third stage is mainly due to the appearance of the dynamic stall in the opposite direction. Therefore, a larger spring equilibrium angle requires a higher aerodynamic load for the system to reach the dynamic stall angle in the opposite direction. A comparison of the results at the three spring equilibrium angles shows that the so-called third phase (the phase where the aerodynamic stiffness grows again) is also divided into two parts. In the second part, the trend of the aerodynamic stiffness with the free-stream Mach number becomes consistent for the different equilibrium angles. This issue will be explored in more depth in the next section. Next, we will compare the aeroelastic response and aerodynamic stiffness characteristics for two spring equilibrium angles that differ significantly. Figure  11 gives the AOA, amplitude and frequency variation characteristics of the system response for spring equilibrium angles of 3 degrees and 20 degrees respectively. The special fixed-point responses mentioned above are not plotted in the diagram as they exist only in very narrow intervals and are not the focus here.
As can be observed from Figure 11, the amplitude and frequency variation characteristics of the system appear very different for large differences in equilibrium angles. In Figure 11(a), with a small spring equilibrium angle like 3 degrees, the fixed-point response continues until the higher Mach number of 0.11 before the bifurcation occurs, compared to the previous example. Moreover, the system quickly enters the symmetric limit cycle state after the bifurcation. Conversely, at a spring equilibrium angle of 20 degrees, as the Mach number increases, there is a period after the bifurcation when the system has an asymmetric limit-cycle response because the amplitude of the system has not yet grown sufficiently large to allow the minimum angle of attack to reach the dynamic stall angle in the opposite direction.
The reason why the case with spring equilibrium angle of 3 degrees enters the symmetrical limit cycle state earlier is that their equilibrium angle is closer to the negative dynamic stall angle. This makes it easier to have dynamic stall in both directions at the same time as the amplitude increases. Therefore, it can be seen that the occurrence of dynamic stall phenomena can have a very significant impact on the variation of the system response. To further observe the effect of dynamic stall on the aeroelastic response, the aeroelastic responses at a spring equilibrium angle of 3 degrees are simulated with different initial conditions. By changing the initial conditions, the system response is allowed to enter any region of the phase plane artificially and thus be subject to different attractors. Figure 12(a) gives a comparison of the change in amplitude with and without initial angular velocity. As shown in Figure 12(b), there is a stable fixed point and a stable symmetric limit cycle for the system at low Mach numbers. As the Mach number increases, a stable asymmetric limit cycle and a stable symmetric limit cycle co-exist, as in Fig. 12(c). Finally, when the Mach number becomes high enough, there is only one symmetrical limit cycle as in Figure 12(d). Here, we note that the initial angular velocity allows the system to be attracted to the symmetric limit cycle attractor at low Mach numbers. As shown in Figure 13, by the comparison for the lift and moment characteristics of the two limit cycle responses at Ma = 0.13, it can be seen that the presence of a symmetric limit cycle is always accompanied by the presence of a dynamic stall in the negative direction.
The comparison of response frequencies of the cases with or without initial angular velocity is also shown in Figure 14. It is found that, when the initial angular velocity is set, the response frequency of the system increases monotonically with Mach number. Because the response type of the system in this case is always a symmetric limit cycle, the change in aerodynamic stiffness is dominated by the increase in overall aerodynamic load due to the increase in Mach number. Note that the trend in frequency variation with and without initial angular velocity is the same at low Mach numbers, but there is a significant difference in the aeroelastic response between the two. At low Mach numbers, the response of an aeroelastic system without initial angular velocity is the fixed point. Actually, the frequencies obtained here are the result of a spectral analysis of the transient oscillation response of the system.

The effect of the spring stiffness
Similar to last section, several structural stiffness coefficients close to that in section 4.1 were selected for simulation to compare their aeroelastic response in Figure 15, the aeroelastic response amplitude and frequency as a function of Mach number for three sets of simulations with structural stiffness coefficients of 7.26, 10.6 and 15.6 N · m/rad are presented. All of the spring equilibrium angles of the systems in this section are 15 degrees. Similar to Figure 10(a), all three cases exhibit three phases, i.e. growing, then reducing, and then growing again. In this section, however, there are still some differences in the natural frequencies of the system at low Mach numbers due to differences in the structural stiffness. In contrast, in the corresponding region of Fig. 10(a), the frequencies of the three cases overlap, as they share the same structural parameters. Furthermore, in Fig. 15(b), we find that the amplitude of the limit cycle increases with structural stiffness when the Mach number is 0.03, i.e., at the onset of the limit-cycle responses in the three sets of simulations for different structural stiffnesses. Here we use the method of energy map to illustrate the emergence of this phenomenon. Firstly, flow field simulations are carried out for oscillating airfoils with a Mach number of 0.03, a mean angle of attack of 13.4 degrees and 130 sets of different oscillation frequencies ω(5rad/s∼30rad/s) and oscillation amplitudes (6deg∼24deg) as oscillation parameters. 13.4 degrees is the value obtained by averaging the mean AOA of the aeroelastic response of each set of simulation at Ma = 0.03. Then the energy map is obtained using the simulation data and the concept mentioned in section 3.2. E * f =0 means that the energy extracted by the airfoil from the free stream over one cycle is zero. This is a condition for the existence of a stable limit cycle response. Looking at the zone where E * f = 0 in Figure 16(a), one can see that the amplitude of the oscillation increases monotonically as the frequency of the oscillation increases. The frequencies and amplitudes of the three limit-cycle responses for different structural stiffnesses at Ma = 0.03 are also plotted on Figure 16(a) (the shape and color of the symbols are the same as in Fig. 15), which indicates that they lie exactly around the E * f = 0 line. Therefore, as the natural frequency of the system is higher due to larger structural stiffness coefficient, larger amplitude limit-cycle responses appear, depending on the characteristics of the energy map of the system at low Mach numbers (0.03). The energy map at Mach 0.07 in Figure 16 (b) are also constructed by performing another set of simulations (the mean angle of attack is about 11 degrees). The region where the energy transfer is about 0 is parallel to the x-axis, which just provides the reason why the limit cycle amplitudes of the three systems are the same when the Mach number is 0.07 in Figure 15 (b).
The similarity between the Figure 10 magnitude of the Mach number, and grows at a certain slope. In fact, when we also look at Figure 15(b), we can see that the amplitude of the system response has also reached a state of saturation. In other words, as the Mach number increases, the response amplitude does not increase as it did in the previous period, but stops increasing or increases very little. The same phenomenon was also found in Ref. [8], which describes this behavior as a loss of memory of the system for the structural parameters of the spring. Figure 17 (a-b) gives the moment coefficients and moments at three Mach numbers for a structural stiffness coefficients of 10.6 N · m/rad as the system enters amplitude saturation. As can be seen from Figure 17(a), the magnitude of the moment coefficient of the system remains almost constant as the Mach number increases. This is due to the fact that the range of oscillation is almost constant as the system enters amplitude saturation state, and the compressibility effect of the Mach number change is not yet apparent. However, the increase in moment value due to an increase in Mach number is more pronounced because of the increase in dynamic pressure 1 2 ρU 2 ∞ on the system, as shown in Figure 17(b). As mentioned in references [13] and [14], the integral value of M sin(ωt) can in some way indicate the magnitude of the aerodynamic stiffness of the system. Figure 17(c-d) shows that the peak of M sin(ωt) does increase with Mach number and therefore the overall integral value, while actually the integral value of C M sin(ωt) does not increase much. The change in moments brings about a change in the aerodynamic stiffness coefficients. In this phase, however, the change in moment is almost exclusively related to the increase in dynamic pressure, so that the increase in aerodynamic stiffness is monotonically increasing, and the slope of the increase is almost identical for different structural stiffnesses and for different spring equilibrium angles. Assume that K a0 = K a /( 1 2 ρU 2 ∞ Sc) = f (C M , ), i.e. K a0 is only related to C M . Therefore, the natural frequency is Where b is the constant term. since K a0 remains constant in the amplitude saturation state, the natural frequency is roughly linear with respect to the Mach number. This also explains the linear increase in natural frequency with Mach number in the simulation results for increasing the initial angular velocity in the previous section. The studies above and Ref. [14] are only based on simulations with small structural stiffnesses, so next we will add three sets of the larger structural stiffnesses for simulation, which are 30, 50 and 100 N · m/rad. The results are shown in Figure 18. Firstly, from Figure 18(a) we can see that as the structural stiffness coefficient increases, it is clear that the natural frequency of the system increases from a base of higher structural frequencies. In addition to this, the variation in the natural frequency of the system becomes smaller and tends to increase more monotonically for larger structural stiffness coefficients. The slope of the increase in frequency with Mach number is also the same for each structural stiffness coefficient during the amplitude saturation phase, which is in agreement with the previous analysis.
One phenomenon that becomes more apparent in Figure 18(b) as higher structural stiffness coefficient data is added is that the amplitude at the onset of the limit cycle increases as the structural stiffness coefficient increases. In addition, for high structural stiffness coefficients, a constant amplitude appears at a certain interval of the Mach number. It appears that the oscillation range of the system no longer varies in this interval, but a detailed analysis of the oscillation angle of attack of the system is still required.  Fig. 19 Maximum, minimum and mean AOA of case with structural stiffness coefficients of (a)10.6 N · m/rad and (b)100 Nm/rad.
From Figure 19(b) it can be concluded that the range of oscillations of the system is not unchanged during this interval when the structural stiffness factor is 100 N · m/rad, but that both the maximum and minimum angles of attack are decreasing simultaneously and at exactly the same rate. In Fig.  19(a), correspondingly, at a structural stiffness factor of 10.6 N · m/rad, the maximum angle of attack of the system remains constant, while the minimum angle of attack is continuously decreasing and therefore its overall amplitude is increasing.
Next, we analyze the reasons for the smaller variation in aerodynamic stiffness at high structural stiffness coefficients. It is more evident that the natural frequency of the system decreases to a lesser extent after reaching a local peak (compare the portion of the curve for structural stiffness of 100 at 0.31-0.35 Mach and the portion of the curve for structural stiffness of 10.6 N · m/rad at The effect of structural parameters on aeroelastic system  Figure 18(a)). The moment and M sin(ωt) of cases with two structural stiffness coefficients around the local peak of natural frequency are compared in Figure 20. From the Figure 20(a-b), we can see that as the Mach number increases, the moment near the angle of attack minimum extends to the left and downwards, which is the key reason why the aerodynamic stiffness becomes less. As shown in Figure 20(c-d), this extension corresponds to an increase in the area enclosed by the part of the M sin(ωt) curve less than 0 and the x axis. Comparing Fig. 20(a) and Fig. 20(b), it can be seen that when the structure is more rigid, the moment curve at the minimum of the angle of attack is more blunt, thus making less of the moment become negative. It is believed to account for the smaller drop in the aerodynamic stiffness drop phase when the structural stiffness coefficients are greater. The difference in moment characteristics around minimum AOA is actually related to the oscillating frequency of the airfoil.

Conclusion
Based on the CFD method, this paper presents a parametric study of the response characteristics and aerodynamic stiffness variation characteristics of a single-degree-of-freedom free aeroelastic system under the influence of dynamic stall.
Firstly, different types of aeroelastic response of the system with increasing free-stream Mach number were found to exist for different spring equilibrium angles. With large equilibrium angles, the response of the system first produces an asymmetric limit-cycle response through a Hopf bifurcation that persists for a Mach interval, followed by the presence of a symmetric limit cycle as the minimum angle of attack continues to decrease to a negative dynamic stall angle. Due to the change in the form of the limit cycle, the system appears to have different moment characteristics and consequently the aerodynamic stiffness first increases, then decreases and then increases. In contrast, with small equilibrium angles the bifurcation velocity is greater and the response rapidly enters the symmetric limit-cycle state after the bifurcation due to the closer proximity to the negative dynamic stall angle. Due to the specific influence of the moment characteristics of the asymmetric limit cycle on the aerodynamic stiffness variation, large differences in aerodynamic stiffness characteristics can occur with different equilibrium AOAs. Furthermore, by increasing the initial angular velocity, we find that the aerodynamic stiffness appears to increase monotonically because the system response all becomes symmetric limit cycle. In addition, a special fixed-point response stabilized around the minimum angle of attack is found at the transition threshold between the asymmetric and symmetric limit cycles.
The influence of the structural stiffness of the system on the aeroelastic response and aerodynamic stiffness characteristics is then investigated. It is first observed that the limit cycle amplitude increases with structural stiffness when the system first bifurcates and then is unaffected by structural stiffness as the Mach number increases. By the analysis on the energy maps of the system at different Mach numbers, it is found that differences in the parametric distribution of the energy transfer coefficients at different Mach numbers contribute to this phenomenon. The reasons why systems with different structural stiffnesses converge in aerodynamic stiffness as they enter the symmetric limit cycle interval are also discussed. In this interval, the increase in the aerodynamic stiffness of the system is determined solely by the increase in dynamic pressure without the effect of changes in moment characteristics. Finally, the reasons for the smaller variation in aerodynamic stiffness at high structural stiffness coefficients are also analyzed. This is because the aerodynamic moments coefficient curve is blunter at the minimum angle of attack for high structural stiffness.