Numerical investigation into discontinuity-induced bifurcations in an aeroelastic system with coupled non-smooth nonlinearities

The present study focuses on investigating the bifurcation characteristics of a pitch–plunge aeroelastic system possessing coupled non-smooth nonlinearities, both in structural and aerodynamic fronts. To this end, a freeplay nonlinearity is considered in the stiffness of the pitch degree-of-freedom. The effects of dynamic stall arising due to large instantaneous angles-of-attackare incorporated using the semi-empirical Leishman–Beddoes aerodynamic model. A systematic response analysis is carried out to discern the bifurcation characteristics of the aeroelastic system considering the airspeed as the system parameter. At low airspeeds, a series of dynamical transitions, including aperiodic responses, occur predominantly due to the structural freeplay nonlinearity while the flow remains attached to the surface of the wing. However, beyond a critical value of airspeed, the system response is dominated by high amplitude pitch-dominated limit-cycle oscillations, which can be attributed to stall flutter. It is demonstrated that the freeplay gap plays a key role in combining the effects of structural and aerodynamic nonlinearities. At higher values of the freeplay gap, interesting discontinuity-induced bifurcation scenarios, such as grazing and boundary equilibrium bifurcations arise due to coupled nonlinear interactions, which can significantly impact the safety of the aeroelastic system.


Introduction
Aeroelastic systems, such as flexible wings, helicopter rotor blades, and wind turbines among many others, are subjected to coupled interactions between inertial, elastic, and aerodynamic forces often with non-negligible nonlinearity, both in the structure or/and aerodynamics [8]. These nonlinear aeroelastic systems exhibit phenomenologically rich dynamical signatures, such as Hopf bifurcations leading to self-sustained limit cycle oscillations (LCOs) [18,34], transition to chaos [1,3], grazing bifurcations [31]. Identifying the onset of such dynamical transitions is important from the standpoint of structural safety as self-sustained oscillations can induce fatigue damage over time, eventually leading to structural failure. Therefore, a bifurcation analysis can directly benefit the design of aeroelastic systems through an appropriate understanding of the underlying fluid-structure interaction dynamics.
The dynamical signatures of the system responses are impacted significantly by the type of nonlinearity, emerging from either the structural [18], or the aero-dynamic counterpart [11] of the aeroelastic system, or even from the combination of both [3,17]. Structural nonlinearities can be broadly categorized into distributed or concentrated types [1]. For large amplitude oscillations, the effects of distributed structural nonlinearities become pronounced and are usually incorporated using a polynomial expression for the stiffness term. On the other hand, concentrated structural nonlinearities, often mathematically modeled using piecewise linear functions, can have a significant effect even in the case of low amplitude oscillations. Freeplay nonlinearity is one such concentrated nonlinearity that has received wide attention in the literature [1,23,32] due to its prevalence and impact on aeroelastic systems. Several studies have typically considered an isolated case of non-smooth structural nonlinearity to show its impact on the system and observed phenomenological transitions in the structural response. Recent studies, such as Vasconcellos et al. [31] and Monfared et al. [26] have considered Theodorsen's unsteady and quasi-steady aerodynamic models, respectively, to estimate the loads on a pitch-plunge airfoil with freeplay nonlinearity and observed grazing bifurcations owing to the non-smooth structural nonlinearity. Verstraelen et al. [35] mathematically and experimentally demonstrated the co-existence of two-domain and three-domain LCOs through a grazing bifurcation scenario in aeroelastic systems with freeplay in pitch.
On the other hand, the aerodynamic nonlinearity can significantly alter the aeroelastic responses when the structure is subjected to flow at high effective AoA under dynamic stall condition [25]. An aeroelastic instability in the dynamic stall regime, known as stall flutter, can give rise to high amplitude selfsustained oscillations. Dimitriadis and Li [9] experimentally investigated an aeroelastic system under dynamic stall and observed the route to stall flutter to occur via a subcritical Hopf bifurcation. Although highfidelity Navier-Stokes simulations [3] can accurately capture the nonlinear aerodynamic effects and the flow separation characteristics, this approach is associated with a prohibitive computational cost. Hence, several semi-empirical aerodynamic models [16] were developed to estimate the aerodynamic loads under dynamic stall conditions at a considerably low computational cost. In this context, the semi-empirical LB model is proven to effectively capture the different dynamic stall events except at very high AoA. The indicial form of the model was modified into a state-space represen-tation by Leishman [20] and Leishman and Nguyen [21], making it convenient for the nonlinear dynamic analysis of an aeroelastic system under the influence of dynamic stall without the loss of generality. LB model is non-smooth in nature; the effect of the associated non-smoothness on the response characteristics of an aeroelastic system subjected to dynamic stall has been investigated by Galvanetto et al. [13]. The authors observed a regime of LCOs with increasing amplitude that transition to aperiodic responses at the threshold of static stall event.
A large section of the existing literature has only focused on investigating the aeroelastic systems with isolated nonlinearity either in the structure or in the flow. However, the presence of coupled nonlinearities can give rise to radically different dynamics that are not observed otherwise [7,19,27]. The presence of discontinuity in both structural and aerodynamic models adds more complexities to analyzing such systems and characterizing the underlying dynamics. Vasconcellos et al. [33] experimentally investigated stall-induced aeroelastic responses in the presence of cubic hardening type of distributed nonlinearity in pitch DoF and have shown that the preset angles (angle set in wind-off conditions) alter the Hopf bifurcation onset. Candon et al. [5] characterized the interactions between the freeplay and aerodynamic nonlinearities at transonic regimes by investigating their higher-order spectra. A recent study by Bethi et al. [2] investigated a pitch-plunge airfoil with various forms of distributed structural nonlinearities under dynamic stall conditions and reported a variety of transition scenarios, such as period-doubling cascade beyond a critical airspeed. Sai Vishal et al. [36] established the route to stall flutter in a system possessing freeplay nonlinearity in the structural component and subjected to nonlinear aerodynamic loads by invoking the theory of synchronization. However, the authors did not explore the response signatures arising from coupling two non-smooth nonlinearities (freeplay and dynamic stall) from a bifurcation perspective and only a synchronization route was suggested.
Discontinuity-induced-bifurcations (DIBs) can be detrimental to the structure as they involve abrupt jumps in the response dynamics between coexisting attractors and have a significant impact even at low amplitudes of oscillations. Kalmár-Nagy et al. [17] showed that aeroelastic systems with non-smooth nonlinearities in both structure and flow exhibited border collision and rapid bifurcations. The authors used a simplified bilinear model fitted to the coefficient of lift data from static stall experiments to calculate the aerodynamic loads. Lelkes et al. [22], using the same aerodynamic model, investigated the impact of a vibration absorber on the stability of equilibria and observed a shift in the onset of rapid bifurcation to higher airspeeds. However, this bilinear model does not account for the delay in the flow separation and added lift generation from the formation of dynamic stall vortex and hence, provides limited insight into the system. A more realistic approach would be to consider a dynamic stall model in order to incorporate the effects of flow separation and vortex shedding. To the best of the authors' knowledge, the bifurcation behavior of such an aeroelastic system consisting of freeplay nonlinearity in the structure and a non-smooth dynamic stall model has largely remained unexplored in the literature. The present study is devoted to addressing this aspect. Another important parameter in aeroelastic systems involving freeplay nonlinearity is the freeplay gap, which can significantly alter the DIB behavior. Studies by Verstraelen et al. [35] have underscored the pivotal role played by the freeplay gap in dictating the DIBs in aeroelastic systems with linear aerodynamics. Extending this understanding to coupled non-smooth nonlinearities, therefore, becomes a necessary step to understand the associated DIB scenario better. The present study aims to take up this analysis as well.
To that end, this study focuses on numerically investigating the dynamical signatures of an aeroelastic system subjected to coupled non-smooth nonlinearities in the flow, as well as in the structure. A pitch-plunge airfoil, possessing freeplay nonlinearity in the pitch DoF and subjected to dynamic stall in the aerodynamics, is considered. The nonlinearity in aerodynamics is captured by using the LB model. A systematic response analysis is carried out with airspeed as the bifurcation parameter. Furthermore, the effect of the freeplay gap on the discontinuity-induced bifurcation scenario is investigated and the effects of mass-ratio and uncoupled frequency-ratio on the bifurcation behavior are also explored.
The rest of the paper is organized as follows. Section 2 presents the mathematical model of the aeroelastic system and validation studies of the present LB model. The detailed bifurcation analysis of the pitch-plunge aeroelastic system in the presence of isolated nonlinearities in the structure and the flow, as well as combined nonlinearities, are presented in Sects. 3.1, 3.2 and 3.3, respectively. The effect of the freeplay gap on the bifurcation behavior and the impact of structural parameters on the onset of dynamic stall are presented next in Sects. 3.4 and 3.5, respectively. Finally, the salient findings from this study are summarized in Sect. 4.

Mathematical model
The structure pertaining to the aeroelastic system is considered to be a 2-D NACA 0012 airfoil with pitch and plunge DoFs constrained using torsional and translational springs attached at the elastic axis (see Fig. 1). The pitch angle is denoted by α, and the plunge displacement, non-dimensionalised using the semi-chord b = c/2 (where c is chord length), is denoted by . The non-dimensional velocity and acceleration terms for the pitch and plunge motions are obtained by differentiating the corresponding responses with respect to the non-dimensional time τ = Vt/b, where V is the dimensional airspeed and t is the dimensional time. The elastic axis is located at a h b from the mid-chord point and the center of gravity is located at x α b from the elastic axis. The stiffness coefficients along the pitchplunge DoFs is denoted by k and k α , respectively. The governing equations of motion for the given aeroelastic system in non-dimensional form is given by Here, U (= V/bω α ) is the non-dimensional airspeed; ω α is the uncoupled pitch natural frequency; is the ratio of uncoupled pitch (ω α ) and plunge (ω ) natural frequencies. μ represents the mass ratio given by  m/πρb 2 , where m is the mass of the airfoil and ρ is the density of air. r α is the radius of gyration about the elastic axis given by I α /mb 2 , where I α is the moment of inertia about pitch DoF. The non-dimensional stiffness along the plunge and pitch DoFs-G( ) and R(α) -are represented as functions of and α, respectively. In the current investigation, the plunge stiffness is considered to be a linear function of , given by G( ) = . The freeplay nonlinearity in the non-dimensional pitch stiffness R(α) is defined using a piece-wise linear function as given by Here, δ represents the region of zero stiffness, known as freeplay gap. The variation of R(α) with the change in α is shown schematically in Fig. 2. The aerodynamic load coefficients C L and C m are estimated using the state-space representation of the LB formulation [20,21], consisting of a set of first order ordinary differential equations (ODEs) that model the aerodynamic loads associated with the different flow modules under dynamic stall; see Eq. (4).
Here, f is a 12 × 1 vector of nonlinear functions, . , x 12 ] T are the twelve aerodynamic states (see Appendix A), U is the non-dimensional airspeed introduced earlier, M is the Mach number,α and q denote the effective pitch angle and effective nondimensional pitch rate, given by [8] α = tan −1 sin α + cos α cos α − sin α , since a h = −1/2. Note that q is a nonlinear function of the acceleration so that Eq. (4) features mass nonlinearity. Although the resulting equations of motion can be integrated numerically using an iterative scheme for , this approach is unwieldy for bifurcation analysis. An alternative is to linearize Eq. (6) by assuming that 2 is small such that This alternative is chosen for the present work, the difference between Eqs. (6) and (7) being of the order of 10 −8 for the frequencies and amplitudes considered here.
The aerodynamic load coefficients are given by where c, N , m represent the non-dimensional coefficients of the chordwise force, normal force, and pitching moment, respectively. The LB model uses the frame of reference where the loads perpendicular to and along the airfoil chord are calculated. The lift coefficient (C L ) is computed as The states x 1x 8 are associated with the unsteady attached flow regime, modified from the Wagner's function formulation by taking the Mach number (M) and flow compressibility into account (refer to Appendix A for the non-dimensional ODEs of the respective states). The normal force and moment coefficients associated with the unsteady attached flow are given by where M is the Mach number, T I = c/a, where a is the speed of sound and A 1 , 4 and b 5 are LB model constants (corresponding values are provided in the Appendix A). C N α is the slope of the normal force coefficient curve against pitch angle at attached flow conditions. K α , K q , K α M and K q M are empirical constants; the corresponding expressions are provided in the Appendix A. The superscript I indicates the impulsive load coefficient and C indicates the circulatory load coefficient terms. A chord value of 0.637 m is chosen in this study following McAlister et al. [24].
The states x 9 , x 10 and x 12 model the flow separation regime and represent the delayed normal force component, trailing edge separation point location, and the delayed version of the trailing edge separation point location. The delay in x 12 is incorporated in order to improve the representation of the pitching moment due to trailing edge separation during flow reattachment. The aerodynamic load coefficients for the corresponding flow regime are given by where the superscript f denotes separated circulatory loads.
The state x 11 accounts for the extra lift generated due to the formation of leading edge vortex, when the value of x 9 crosses an experimentally obtained critical normal force (C N 1 ) value. K 0 , K 1 and K 2 are coefficients related to the position of the aerodynamic center and the shape of the moment break at stall that are estimated from static experiments. The normal force and moment coefficients generated by the leading edge vortex are given by Here, T vl is the experimentally obtained value of time taken for the vortex to travel one chord, while τ v is the vortex time that starts when |x 9 | = C N 1 and progresses with non-dimensional time τ as |x 9 | is increasing. The value of τ v is reset to 0 when |x 9 | < C N 1 .
It is worthwhile to mention that the boundary at |x 9 | = C N 1 introduces non-smoothness in the aerodynamic model. The values of empirical constants C N α , K 0 , K 1 , K 2 and T vl are dependent on M and are provided in Appendix A. The expressions for the non-dimensional derivatives of the states x 9 − x 12 are also provided in Appendix A. The total values of the aerodynamic load coefficients are obtained from summing the separated circulatory, impulsive and vortex contributions, such that The present LB model is validated with existing experimental results on a NACA 0012 airfoil at dynamic stall conditions. The coefficient of moment (C m ) and the coefficient of normal force (C N ) calculated from the LB model, used in this study [13], is validated against the experimental results of McAlister et al. [24] at M = 0.3 for an airfoil, sinusoidally pitching with the kinematics-α(τ ) = 12 • + 10 • sin(κτ )with the value of the reduced frequency κ = 0.0976 (see Fig. 3a, b). It is observed that the present results show a close agreement with the reference experimental results in this M regime. Evidently, the LB model is seen to be capable of capturing the aerodynamic loads under dynamic stall conditions with an agreeable accuracy.
The full list of the discontinuity boundaries of the model is It should be noted that most of these discontinuities are modeling artifacts. Although dynamic stall is an abrupt phenomenon comprising a series of rapid transitions among different stages, it has not yet been established with confidence that it is discontinuous in the strict mathematical sense. The shedding of the Leading Edge Vortex, corresponding to |x 9 | = C N 1 in the LB model, is probably the most abrupt phenomenon occurring during dynamic stall but the time scale of this event is yet to be established for general cases. The LB model is quite successful in modeling the aerodynamic load response during dynamic stall, even in the challenging post-shedding phase, and is therefore used here. Some aspects of the model can be smoothed, such as Kirchoff's load computation for low Mach numbers [4]. As the Mach number is significant in the present study, the original LB formulation is used, with the discontinuous Kirchoff function of Eq. (28) (although one aspect of this function is smoothed over, see later).

Results and discussion
Three test cases are investigated in this section: • A system with freeplay in pitch and linear aerodynamics, obtained from the low-amplitude, linearized version of the Leishman-Beddoes model. • A system with linear structure and fully nonlinear Leishman-Beddoes aerodynamic loads. • A system with both freeplay in pitch and nonlinear aerodynamic loads, which is the combination of the previous two systems.
Studying the effects of the structural and aerodynamic nonlinearities separately sheds light on the dynamic behavior of the combined nonlinear system since characteristics of the responses of the first two systems are found in the third. The tools used for the dynamic investigation are stability analysis of fixed points, numerical

Non-smoothness only in the structure
In this section, we consider a pitch-plunge aeroelastic system with structural freeplay nonlinearity, subjected to linear aerodynamics. The structural parameters are chosen from Lee et al. [18] and are shown in Table 1.
The aerodynamic loads come from the linearized version of the Leishman-Beddoes model of the previous section. The flow is assumed to be attached at all times such that only the first eight aerodynamic states are required. The Eq. (5) for the effective pitch angle is linearized tô while the non-dimensional pitch rate is given by Eq. (7). The aeroelastic equations of motion are reduced to Substituting these latest expressions in the equations of motion (1) and (2) and combining with Eq. (27) finally leads to Q(U, M) and q(U, M). Of all the Machdependent parameters in Table 4, only the lift-curve slope C N α and K 0 affect the linearized aerodynamic model. The value of C N α for Mach number less than 0.3 is obtained from the Prandtl-Glauert rule while K 0 , which determines the moment of the lift around the aerodynamic center, is set to zero at very low Mach numbers in this work. 1 Therefore, K 0 is set to decrease exponentially from K 0 (0.3) = 0.0125 at M = 0.3 to 0 at M = 0.15 and to remain zero down to This treatment is slightly arbitrary but it ensures the continuity of the system eigenvalues over the entire airspeed range of interest. The values of C N α and K 0 for all Mach numbers equal to or higher than 0.3 are obtained by interpolating the data in Table 4. The presence of freeplay in a linear system splits the phase space into three piecewise linear subdomains, as shown in Fig. 2. Inside the freeplay region, |α| ≤ δ/2, the response dynamics are determined by the linear system, which has zero stiffness in pitch and is referred to as the inner linear system. The pitch axis lies on the c/4; as mentioned earlier, the aerodynamic pitching moment around this point has been set to 0 for M ≤ 0.15. At higher Mach numbers, the aerodynamic pitching moment is non-zero and determined by Eq. (24) or Table 4. For K 0 = 0, the matrix Q(U, M) has a zero eigenvalue, which means that the fixed point lying at x F = 0 is neutrally stable. This situation occurs at M ≤ 0.15, which corresponds to U ≤ 2.5 given the choice of system parameters. For K 0 = 0, Q(U, M) has a positive real eigenvalue since it has no structural stiffness in pitch to counterbalance the aerodynamic pitching moment. This means that for all U > 2.5, the fixed point x F = 0 is an unstable spiral. This condi-tion is known as static divergence in linear aeroelasticity and the divergence airspeed is denoted by U D i . All response trajectories are pushed away and cross the freeplay boundary. As a consequence, a boundaryequilibrium bifurcation occurs at U D i = 2.5. Furthermore, linear flutter occurs at another critical value of the airspeed, U F i = 1.04, which means that the single pair of complex conjugate eigenvalues of matrix Q(U, M) becomes purely imaginary and the fixed point x F = 0 becomes a neutrally stable focus. Circles whose pitch amplitude is equal to δ/2 graze the two discontinuity boundaries. Hence, a grazing bifurcation occurs at the flutter speed of the inner linear system. The flutter speed, U F i , flutter frequency ω F i and divergence speed U D i of the inner linear system are given in Table 2.
Outside the freeplay region, |α| ≥ δ/2, the system dynamics are described by which is essentially a linear system with a nonlinear fixed point, given by The α component of this fixed point is plotted in Fig. 4a; it lies exactly on the freeplay boundary for airspeeds up to U D i when it starts moving away from this boundary. The stability of this fixed point is determined by the stability of the linear system x = Q o (U, M)x, which will be referred to as the outer linear system. The matrix Q o (U, M) has no positive real eigenvalues throughout the speed range of interest. The system undergoes flutter at U F o = 5.53, as shown in Table 2, when one pair of complex conjugate eigenvalues becomes purely real, as shown in Fig. 4b. Therefore, beyond U F o , both the  In order to calculate numerically the complete dynamics of the system, Eq. (22) are solved using the pseudo-arclength continuation scheme in [8]. Given initial guesses for a point x 0 on a limit cycle and for the period of the cycle, T 0 , at a chosen airspeed U 0 , a Newton system is set up such that ⎛ where Δx, ΔT , ΔU are improvements to x 0 , T and U, f is the system's nonlinear function, x(T ) is the value of the system states obtained by means of a time integration of Eq. (22) over T seconds and with initial conditions x 0 , s is the pseudo-arclength and x −1 , x −1 are the states and their derivatives of a point on a nearby limit cycle that has already been determined accurately. The time integration is carried out by means of a 4-5 order Runge-Kutta scheme with event detection in order to determine accurately the crossings of the freeplay boundaries. The derivatives are determined numerically using a forward finite difference scheme. The Newton system is set up and solved repeatedly until convergence is achieved and then new direction vectors ∂x/∂s, ∂ T /∂s, ∂U/∂s are calculated and a new limit cycle is pinpointed. Limit cycles branches are continued until a chosen amplitude is exceeded; the stability of each limit cycle is determined using Floquet theory and folds and branch points are detected. Figure 5a plots the variation with airspeed of the maximum value of α over a complete cycle. This maximum corresponds to the amplitude for symmetric limit cycles but not for asymmetric ones. There is a primary branch of symmetric limit cycles that emanates from the grazing point that occurs at the flutter speed of the inner linear system, U F i , and with period 2π/ω F i . The branch is initially unstable and propagates in the decreasing airspeed direction but it then undergoes a fold and becomes stable. At U = 1.74 there is a branch point that causes the primary branch to become unstable and gives rise to a secondary asymmetric branch. The latter is again unstable and propagates in the decreasing U direction until it folds, becomes stable, and changes direction. It folds again at a higher airspeed and rejoins the primary branch at the second branch point, lying at U = 3.2. The amplitude of the primary branch continues to increase with U and tends to infinity as the flutter speed of the outer linear system is approached. The figure also plots the α coordinate of the stable fixed point that appears as a result of the boundary-equilibrium bifurcation mentioned earlier. Figure 5b plots the variation of the non-dimensional period, T = 2π V/ωb, of the two limit cycle branches with airspeed.
Both branches are related to grazing bifurcations. The primary branch is generated when the circles of the fluttering inner linear system graze the freeplay boundary, as shown in Fig. 6a which plots a phase plane section of the unstable limit cycle at the grazing point. The secondary branch also grazes the freeplay boundary but in a different way, as seen in Fig. 6b. The cycle features two loops, a large one that crosses both freeplay boundaries and a small one that grazes the boundary at −δ/2; this limit cycle lies near the first fold of the secondary  Fig. 6b, which grazes +δ/2, is also a valid solution of the system. Finally, it must be mentioned that Floquet analysis did not detect any period doubling or torus bifurcation.

Non-smoothness only in the flow
In this case, the structure is linear so that the effects of the non-smooth dynamic stall model can be studied in isolation. Galvanetto et al. [13] analyzed the bifurcation behavior of a similar system, choosing the Mach number as the bifurcation parameter; they observed a period-doubling cascade, culminating in aperiodic oscillations. However, in the present study, the nondimensional airspeed U is chosen as the bifurcation parameter. The density is chosen to be constant so that the Mach number varies with U and the parameters of the LB model at each U value are calculated by interpolation from the data provided in Table 4 in Appendix A.
Adding the equations of motion to the aerodynamic equations of the LB model given in Eq. (4) results in an augmented system of the form    symmetric. Therefore, the system is in equilibrium when all the displacements are zero; the flow is fully attached in this condition, which means that the separation point lies at the trailing edge. Consequently, the fixed point of the system is a 16 × 1 vector of zeros, except for the 10th and 12th states that are equal to 1; these states represent the position of the separation point. The next step is to calculate the stability in the Lyapunov sense of this fixed point by linearizing the system around x F in the form where Δx is a small departure from x F and ∂f/∂x| x F is the Jacobian evaluated at x F . Figure 7 plots the variation of the real and imaginary parts of the eigenvalues of the Jacobian, showing clearly that one pair of complex conjugate eigenvalues becomes purely imaginary at U H = 5.53. This is the well-known condition of the Hopf bifurcation and signifies that a branch of limit cycles will emanate from the Hopf point U H , with frequency ω H = 35.4 rad/s. Note that the U H is identical to the flutter speed of the outer linear system in Table 2. This is logical since the structural system is now linear and has stiffness K α in the entire phase space. Its stability around the fixed point is determined by very low-amplitude aerodynamic loads, which are linearized. Therefore, close to the fixed point, the present system and the outer linear system of Sect. 3.1 are identical.
Strictly speaking, the nonlinear function f is not differentiable around x F because of the discontinuity in slope of Eq. (28) atα = 0. For the present case, the slope of f (α, α 1 ) jumps from 2.6×10 −3 to −2.6×10 −3 asα crosses from 0 − to 0 + . In order to calculate the Jacobian around x F the slope d f /dα was set to 0 at α = 0. This is a negligible modification of the LB model, considering that the same slope can take values over ±30 at higher values of the effective pitch angle. It also renders the model more physical because there is no reason why the derivative of the lift curve slope of an airfoil should be discontinuous across α = 0. As a consequence, the system does indeed undergo a Hopf bifurcation at x F and the resulting limit cycle branch can be determined by means of numerical continuation, up to a point.
The shooting-based pseudo-arc length continuation scheme mentioned previously is used again in order to determine the dynamics of the system with linear structure and nonlinear aerodynamics. At small amplitudes of oscillation, there are no discontinuities in the system, other than the one occurring atα = 0 and dealt with already. The system becomes non-smooth first when |x 9 | = C N 1 , which is grazed before the other discontinuity boundaries as the oscillation amplitude increases. Figure 8 plots the variation of limit cycle amplitude and period with U close to the Hopf point. It can be seen that the bifurcation is supercritical, as the resulting limit cycle branch is stable and propagates in the increasing airspeed direction. Both the amplitude and frequency increase smoothly with airspeed, without undergoing any folds, until x 9 grazes C N 1 , as seen in Fig. 8c. The grazing limit cycle is plotted in a phase plane section in Fig. 8d. At that point, the numerical continuation procedure detects a branch point; it then breaks down because the system develops memory. The time response from initial conditions over a cycle depends on the time  instance at which |x 9 | exceeds C N 1 , that is the time the vortex is shed. The vortex shedding time now becomes an additional system state, which is not associated with an equation of motion. If the vortex time is retained as an initial condition, the Newton system of the continuation process becomes overdetermined; if it is ignored, the Newton system fails to converge. It may be possible to write an additional equation of motion for the vortex state but such a development is not the focus of the present paper. As a consequence, the dynamics of the system after the grazing bifurcation are determined by means of long time integrations only. Figure 9 plots a bifurcation diagram, obtained by calculating the maxima of α(τ ) over a long time duration at each airspeed, by means of a 4-5 order Runge-Kutta scheme with very small tolerance of 10 −12 and no event detection. In Appendix B it is shown that this way of solving the equations of motion gives identical results to using event detection. The same conclusion was reached by Liu et al. [23] and Galvanetto et al. [13]. The initial behavior of the system is identical to the one seen in Fig. 8. After the branch point, the response undergoes a transition to aperiodic oscillations. At around U = 6.7, the system transitions again to periodic high-amplitude oscillations typical of a fully developed stall flutter. This second transition occurs when all the discontinuity boundaries are exceeded by a significant amount during an oscillation cycle so that no grazing or near-grazing events are occurring. The time integrations are stopped at U = 8 when the amplitude in pitch reaches 30 • .
It is not clear whether the aperiodic oscillations observed in Fig. 9 are physical. Experimentally observed aperiodic stall flutter oscillations and coexisting limit cycles have been reported in the literature [10] but at very low Mach numbers and the nature of the coexisting limit cycles were different. As the numerical continuation scheme could not proceed past the branch point, the cause of the aperiodic oscillations is unclear. It may be that the multitude of discontinuous boundaries results in numerous grazing events so that the limit cycle becomes a torus but this explanation could not be confirmed. The present section deals with the bifurcation behavior of the pitch-plunge aeroelastic system featuring both freeplay nonlinearity in the pitch stiffness and aerodynamic nonlinearity due to dynamic stall during large amplitude oscillations. The freeplay gap is set to δ = 1 • , as in Sect. 3.1. At low amplitudes of oscillation and airspeeds, the Leishman-Beddoes model behaves nearly linearly, so that the dynamics of the system are still determined by the inner and outer linear systems of Sect. 3.1. The flutter and divergence characteristics of Table 2 and the fixed point positions and dynamics of Fig. 4 are still relevant. The high-amplitude dynamics are dictated by the outer linear system and nonlinear aerodynamic loads, just like the system in Sect. 3.2. Figure 10 plots the limit cycle amplitude and period of this system, calculated using numerical continuation. The plot stops when the amplitude in x 9 reaches C N 1 at U = 5.47, a branch point occurs and the numerical continuation scheme breaks down, as mentioned in the previous section. Figure 10a, b are nearly identical to those of Fig. 5a, b for the system with linear aerodynamics; the only significant difference being the occurrence of the branch point at |x 9 | = C N 1 . The stable fixed points are also nearly identical, while the boundaryequilibrium bifurcation occurs at the same airspeed i.e. at U = 2.5. As the numerical continuation scheme breaks down at the third branch point, the subsequent dynamics of the system are analyzed by means of long time integrations. Figure 11 plots the bifurcation diagram of the system, obtained using long time integrations at each airspeed. The numerical scheme was again the Runge-Kutta 4-5 method with the same low tolerance. The first limit cycle oscillations occur at around U = 0.9 and they have small amplitude, as shown in Fig. 12a, b. Furthermore, the flow is fully attached; Fig. 12c shows that x 10 = 1 at all times, i.e. the separation point lies at the trailing edge, and that x 9 does not increase beyond ±0.04, i.e. |x 9 | << C N 1 so that no vortex shedding can occur. The limit cycle occurs due to the grazing bifurcation occurring at U F o = 1.04 as discussed in Sect. 3.1. The corresponding frequency spectrum (Fig. 12c) shows that the response has a single dominant frequency along with its super-harmonic, representative of a period-1 limit cycle [23,31].
As U is gradually increased up to U = 3.4, a secondary limit cycle branch is created at the first branch point, resulting in asymmetric limit cycles that feature a small loop that grazes the freeplay boundary, as seen in Fig. 13a, b. At higher airspeeds, this secondary branch merges with the primary branch, and only symmetric limit cycles can occur, see Fig. 13c, d. Again, the amplitudes of these limit cycles are low enough to ensure that the aerodynamic loads are nearly linear and no stall or dynamic stall occurs. Figure 14 plots x 9 versus x 10 phase portraits up to U = 5.4, showing that x 10 = 1 and |x 9 | < C N 1 at all times.
At U = 5.47, x 9 starts to graze the C N 1 boundary, as shown in Fig. 15. At U = 5.45, x 9 approaches this boundary but does not graze it yet, see Fig. 15a, b. Furthermore, x 10 does not drop below 0.96, which means that the flow is nearly fully attached. At U = 5.6, x 9 crosses the C N 1 boundary and α grazes the α 1 boundary, see Fig. 15c, d. Now, the trailing edge separation point position, x 10 oscillates between 0.4 and 1, so that a significant portion of the wing experiences separated flow. Furthermore, the crossing of the C N 1 boundary means that the vortex-induced aerodynamic loads are no longer equal to zero at all times. As a consequence, the flow oscillates between attached conditions, light stall, and deep stall and the oscillations start to become quasi-periodic. Light stall signifies that only trailing edge separation occurs. Deep stall means that a leadingedge vortex is shed [9,13]. At U = U F o = 5.53 the stable fixed points shown in Fig. 10a become unstable. Response trajectories starting near the fixed point spiral outwards until they cross the freeplay boundary. Then, they increase abruptly in amplitude and graze or cross both the |x 9 | = C N 1 and |α| = α 1 boundaries. The steady-state response is an aperiodic limit cycle since no periodic limit cycles can exist between U = 5.47 and U = 6.6. Even though the bifurcation occurring at U F o involves a pair of complex eigenvalues becoming purely imaginary and the sys- tem is nonlinear, it cannot really be called a Hopf. No limit cycle branch is generated at the fixed point and the response jumps directly to a high-amplitude aperiodic oscillation. It is worth summarizing that, between U = 5.47 and 5.53, the system undergoes multiple grazing bifurcations and a Hopf-like bifurcation. As the airspeed is increased further, the response becomes completely aperiodic until U = 6.6; see Fig. 16a, b. The aperiodic attractors have been categorized as quasiperiodic and chaotic states using the topological invariants, such as Largest Lyapunov Exponent and correlation dimension; see Appendix C. In Fig. 16a, the range of α values is split into five sub-domains; S 1 lies inside the freeplay region, S 2,3 inside the δ/2 ≤ |α| ≤ α 1 region and S 4,5 in the |α| > α 1 region. It can be seen that the time response of α oscillates aperiodically between all five of these regions. The boundary crossing picture of Fig. 16a is not complete because there are other discontinuous boundaries that do not depend on α, see Eq. (21). The x 9 versus x 10 phase plot (Fig. 16c) shows multiple loops, the smallest of which stay within x 10 > 0.9 and |x 9 | < C N 1 , which means that the flow is nearly attached. The vast majority of the loops cross the C N 1 boundary, so that vortex shedding is occurring over most of the cycles. Similar in (b) represent the boundaries of sub-domains S4 and S5 located at static stall angle (α 1 = ± 11.96 • ) and the red dashed lines represent the boundaries of sub-domains S1, S2 and S3, respectively aperiodic responses were observed and attributed to an aperiodic attractor in [13]. The response transitions from aperiodic to periodic as the airspeed is increased to U = 6.6 (Fig. 17a), marking the onset of high-amplitude periodic and symmetric stall-induced LCOs; see Fig. 17b. This event can be attributed to deep stall as the amplitudes of the pitch response reach a value considerably higher than the static stall angle; see Fig. 17c. The value of x 10 decreases to nearly zero at every half cycle and |x 9 | exceeds ±C N 1 by a significant amount. This means that not only does the trailing edge separation point move up towards the leading edge, but a strong vortex is also shed twice every cycle, causing additional vortex-induced lift and pitching moment, as seen in Fig. 17d, e. Flow reattachment begins when the vortex clears the trailing edge, after the sharp drop in the pitching moment. Separation, vortex shedding, and reattachment occur alternatively on the upper and lower sides of the airfoil.
The effect of freeplay nonlinearity on the system dynamics is predominant at lower amplitudes of oscillations. As shown in Sect. 3.1, the oscillation amplitude due to freeplay remains small and only starts to increase significantly as the flutter speed of the outer linear system is approached. At this condition, the freeplay-induced oscillations start to cause light stall earlier than the nonlinear aerodynamics alone would have caused it. Deep stall, on the other hand, appears virtually unaffected by freeplay. However, if a system has only freeplay nonlinearity, the response amplitude is directly proportional to the freeplay gap [6]. This means that increasing the freeplay gap will increase the amplitude range in which freeplay has a significant effect and, in the case of the system with both freeplay and nonlinear aerodynamics, will increase the interaction between the two nonlinearities [35]. This serves as an impetus to investigate the effect of increasing the freeplay gap and the following section presents the same.

Effect of freeplay gap on the onset of dynamic stall
The bifurcation plot of Fig. 11 may be seen as an approximate superposition of the bifurcation plots of Figs. 5 and 9. This indicates a minimal interaction between the nonlinearities arising from the structure and the aerodynamics due to the fact that the freeplay gap is very small. In fact, it is so small that, when the system is undergoing dynamic stall oscillations, the freeplay is nearly linearized. In this section, it is shown that increasing the freeplay gap causes higher degrees of interaction between the two nonlinearities. The bifurcation plots for the cases δ = 2 • , 3 • and 4 • , focusing on the speed range U = 5 − 7, are shown in Fig. 18a-c. These particular airspeeds are chosen because the system response transitions into stall flutter in this range. It is observed that the qualitative nature of the dynamics is similar in all cases. However, the amplitude of the oscillations reaches the static stall angle (indicated using black lines) at progressively lower airspeeds as δ is increased. For example, the amplitude of oscillations reaches α 1 at U = 5.6 when δ = 1 • (refer to Fig. 11 in Sect. 3.3) while for δ = 4 • (Fig. 18c), it reaches it at U = 5.3. The onset of large-amplitude stall flutter is also shifted to lower airspeeds as the freeplay gap is increased. It is worth noting that the width of the region of aperiodic oscillations reduces as δ is increased. The qualitative nature of the response is similar for all δ in the chosen airspeed range: small amplitude oscillations increasing to light stall, followed by higher amplitude aperiodic oscillations involving vortex shedding and finally periodic deep stall oscillations. However, the amplitudes of each of these phenomena and the airspeed ranges in which they occur are quantitatively affected by the increase in δ.
The effect of the freeplay gap can be seen more clearly when the bifurcation plot is calculated using δ 3042 S. Vishal et al. and not U as a bifurcation parameter. Figure 19 depicts the bifurcation plots for distinct airspeeds U = 5, 5.5, 6, 6.5 as function of δ. The chosen airspeeds include all the dynamical transitions that occur in the range U = 5 -7. It can be seen that, at U = 5 and 5.5, increasing the freeplay gap from 0 • to 5 • changes the nature of the steady-state response from static to aperiodic or quasi-periodic stall flutter oscillations. For U = 6 all responses are quasi-periodic, only the amplitude is affected by the freeplay gap. For U = 6.5, the response evolves from aperiodic to periodic stall flutter. The most important conclusion from Fig. 19 is that increasing the freeplay gap can cause stall flutter to occur at an airspeed significantly lower than the Hopf speed of the system without freeplay, which lies at U = 5.53.
The phase plots of x 9 and x 10 corresponding to δ = 0 • , 1 • , 2 • and 3 • are given in Fig. 20a-d, respectively, for U = 5.7, in order to demonstrate further the impact of freeplay on dynamic stall onset. Figure 20a represents the case δ = 0 • ; flow separation is seen to be negligible, and no vortex formation can be observed as the x 9 values still lie within the C N 1 boundary (marked by dashed lines in black). However, at δ = 1 • , the phase portrait of x 9 and x 10 is of quasi-periodic nature, such  that the response switches between light stall and deep stall; see Fig. 20b. As δ is increased to values greater than 1 • , the phase portraits are indicative of aperiodic dynamics with the response reaching further into the deep stall regime (see Fig. 20c, d). The extent of trailing edge flow separation is also observed to increase with freeplay gap as x 10 reaches values close to zero, representing an almost completely separated flow.

Impact of system parameters on the onset of dynamic stall
In the previous section, the freeplay gap was varied but the structural system parameters remained constant. Naturally, these parameters play a primordial role in dictating the aeroelastic response of the system as they affect the flutter speeds and frequencies of the inner and outer linear systems and, hence, the airspeeds at which freeplay-induced and stall-induced oscillations occur as well as the frequency of these oscillations. Many studies have demonstrated the dependence of various nonlinear aeroelastic systems on the structural system parameters, see [3,14] for recent examples. The effects of two important system parameters, namely, the mass ratio (μ) and frequency ratio ( ) on the bifurcation characteristics and flutter onset, are investigated by keeping the other structural parameters constant, as already given in Table 1. First, μ is incremented in the range 50-250, and the resulting changes in the system dynamics and onset of stall flutter are studied. It is important to study the effect of the value of μ on the dynamic characteristics of the inner and outer linear systems. Figure 21a plots the variation of the flutter speeds of the two linear systems with μ, while Fig. 21b presents the variation of the flutter frequency. It can be seen that, as μ increases, both flutter speeds increase while both flutter frequencies decrease. As the wing becomes heavier with respect to the air and the springs, higher airspeeds are required for the aerodynamic loads to be important enough to destabilize the system, while the wind-off natural frequencies decrease.
A combined bifurcation plot for μ = 50, 100, 150, 200, 250 is presented in Fig. 22a and the corresponding zoomed sections are provided in Fig. 22b-e. It is observed that the onset of stall flutter is postponed as the value of μ is increased. This can be attributed to the fact that the flutter speeds of both the inner and outer linear systems increase with μ, as mentioned in the previous paragraph. No significant change in the qualitative nature of the system was observed but all oscillations are moved to higher airspeeds, all stall flutter amplitudes are decreased (the amplitudes of the freeplayinduced oscillations depend mostly on the size of the freeplay gap), while the airspeed range of the aperiodic oscillations increases between μ = 50 and 100 and decreases steadily for μ > 100. Table 3 details the airspeed ranges of the three main types of response encountered.
Next, the frequency ratio ( ) is incremented in the range of 0.3-0.99 while the other parameters are kept constant. Figure 23 plots the effect on the flutter speeds and frequencies of the inner and outer linear systems. The figure shows that both the flutter condition of the inner linear system increases linearly with , since this parameter determines the stiffness of the only spring in the system. On the other hand, the flutter speed of the outer linear system first decreases and then increases with . At around = 0.8, the outer linear system  starts to flutter at a lower airspeed than the inner. This means that, for ≥ 0.8, only stall can cause selfexcited oscillations. The inner linear system may flutter at a higher airspeed than the outer but it is still unstable due to the positive real eigenvalue at U > 2.5. Consequently, in the speed range between U = 2.5 and U F o the system has a stable fixed point just outside the freeplay region, whose position is given by Eq. (25). At all speeds above U F o , the system response is stall flutter, as a Hopf bifurcation occurs around the fixed point at U F o .
Similar to the previous case, a combined bifurcation plot for = 0.3, 0.5, 0.7, 0.8 and 0.99 is presented in Fig. 24a and the zoomed sections of the marked inset (U = 0 − 3) are provided in Fig. 24b-f. In all cases, the initial conditions used in the time integrations lay outside the freeplay range. It is observed that the onset of freeplay-induced LCOs is delayed while the onset of stall flutter is shifted to lower airspeeds as the value of increases. The response dynamics of the system are significantly affected by the change in the value.

Conclusions
The present study investigates the response dynamics of a pitch-plunge aeroelastic system possessing discontinuous nonlinearities in both the structure and the aerodynamics. First, the bifurcation behavior of the system with freeplay only is analyzed, followed by that of the system with nonlinear aerodynamics only. Then, the system is subjected to the combined effects of dynamic stall and freeplay nonlinearity in the pitch degree of freedom. Finally, the impact of three important structural system parameters-freeplay gap (δ), mass ratio (μ), and frequency ratio ( ), on the response characteristics of the system are investigated. The important findings of this study are the following: • The bifurcation plot of the system with both freeplay and aerodynamic nonlinear is approximately a superposition of the bifurcation plots of the systems with isolated nonlinearities. Freeplayinduced and stall-induced effects interact only within a narrow airspeed region around the flutter speed of the outer linear system, which is also the   Hopf condition of the system without freeplay. Past this airspeed, the freeplay is effectively linearized, since the oscillation amplitude is much higher than the freeplay gap. Freeplay only affects the dynamics of the system at low airspeeds, causing lowamplitude oscillations. • Increasing the size of the freeplay gap also increases the interaction between the two nonlinearities. The amplitude of the freeplay-induced oscillations is proportional to the freeplay gap, so that, if this gap is large enough, dynamic stall can occur at airspeeds significantly lower than the flutter speed of the outer linear system. • Increasing the value of the mass ratio μ does not alter qualitatively the system's dynamics, even though all limit cycles occur at higher airspeeds and have slightly lower amplitudes. The range of airspeeds in which aperiodic LCOs can occur shrinks drastically as μ increases. • Varying the ratio of plunge-to-pitch natural frequency modifies significantly the behaviour of the system. At particular values of , the flutter speeds of the inner and outer linear systems are inversed so that no freeplay-induced oscillations can occur.
The system investigated in this work features a very high number of discontinuous boundaries: two boundaries due to the freeplay and 11 boundaries due to the Leishman-Beddoes model. The number and nature of the discontinuities have made it impossible to carry out numerical continuation beyond the grazing of the first Leishman-Beddoes boundary. Future work will be aimed at developing a numerical continuation strategy that can handle the sudden appearance of an additional vortex-related system state when deep stall occurs. from Wallonie-Bruxelles International, Belgium towards this research.

Conflict of interest
The authors declare that they have no conflict of interest.

A Aerodynamic state-space ODEs
The set of ODEs used to calculate the states pertaining to the unsteady attached flow regime are given by The LB model uses the Kirchhoff theory to calculate the load coefficients corresponding to the flow separation regime. The point of trailing edge separation is used to determine the loss in the normal force coefficient with respect to the ideal flow scenario. The position of the trailing edge separation point is given by Here, α 1 is the point where f (α 1 , α 1 ) = 0.7, which is approximately equal to the static stall angle. It should be noted that α 1 serves as a discontinuous boundary that is responsible for splitting the phase space into two additional domains apart from the existing domains due to the presence of freeplay in the structure. S 1 and S 2 are constants that are obtained from experiments for each airfoil type. The ODEs corresponding to the flow separation regime are given by T P , where T P , T f and T f 0 are time delay constants obtained from dynamic stall experiments. The parameters T f and α 1 vary as the flow detaches and re-attaches. The vortex shedding phase begins when the value of |x 9 | ≥ C N 1 marking the onset of flow separation. In this regime, the parameters vary such that Here α 1 0 is the experimentally obtained static stall angle of incidence, and δ α 1 is a parameter dependent on airfoil shape that is used to capture the point of static stall angle with better accuracy during a dynamic event. The flow reattachment process begins when |x 9 | < C N 1 and the parameters T f and α 1 in this regime are defined as where T vl is the experimentally obtained value of time taken for a vortex to travel one chord. The ODE that provides the solution to the state x 11 is given by c v is the derivative of the vortex feed c v that determines the strength of vortex induced normal force given by The parameter (T v ) that controls the change in x 11 also changes according to the flow condition and is given by In the flow separation phase and during flow reattachment (|x 9 | < C N 1 ), T v = T v 0 . Finally, the 12 aerodynamic states x 1x 12 defined earlier and the structural states depicting the pitch and plunge, velocity and acceleration terms, [α , α , , ] together form a state-space system of total 16 ODEs i.e h = [x 1 , x 2 , ..., x 16 ]. The parameters α 1 0 , δ α 1 , S 1 , S 2 , T P , T f 0 , T vl and T v 0 are also dependent on M and the parameter values for each M concerned with this study are provided in Table 4. Note that the parameter values are obtained from are obtained from Galvanetto et al. [13].
The intermediate values of the empirical parameters between any two Mach numbers are estimated using the polynomial cubic hermite interpolation technique for the bifurcation study. The expressions of the empirical constants K α , K q , K α M and K q M are given by

B Validity of the chosen numerical integration scheme
The integration method adopted can affect the solutions of the system. It was verified by Galvanetto et al. [13] that such stringent values of tolerance can visibly negate the errors arising from integration for discontinuous systems (albeit the computational time may marginally increase), by comparing results obtained through event detection and without event detection.
To ensure an accurate solution with the fourth-order Runge-Kutta integration scheme, we have considered a stringent tolerance of 10 −12 (both in absolute and relative levels of tolerance) for the numerical integration using an adaptive time-stepping approach. To validate the accuracy of the chosen numerical integration scheme, we compare the aperiodic responses presented in our manuscript for two different tolerance values with the solutions obtained with event detection (see Fig. 25). It is evident from Fig. 25 that a less stringent tolerance of 10 −6 can influence the errors in capturing the discontinuous boundaries. However, the stringent tolerance of 10 −12 considered in this study seems to be as effective as the event detection scheme in accurately capturing the dynamics. Therefore, we adopt this approach as it reduces the computational cost significantly. Figure 26 presents the reconstructed phase-portraits corresponding to the pitching responses of the aeroelastic system with combined structural and aerodynamic nonlinearities as a representative case for three different flow velocities, where the system exhibits quasi-periodic (U = 5.8), weakly chaotic (U = 6.5), and periodic (U = 7.0) response, respectively. The phase-space reconstruction [29] is carried out based on a optimum time delayτ , estimated using the method of average mutual information [12]. Although the phasespace attractors can qualitatively distinguish between the periodic and the aperiodic states, quantitative measures are required to precisely identify the nature of aperiodicity. To that end, we characterize the aperiodic responses by estimating the quantitative topological measures of the corresponding reconstructed phaseportraits, namely the largest Lyapunov exponent (LLE) and the correlation dimension; see Fig. 27. Lyapunov exponent is the quantitative measure of the exponential rate at which an infinitesimal perturbation to a trajectory of a system grows or decays in the state space and is a measure of the sensitivity of the system to the initial conditions. LLE is calculated in this study using the Rosenstein algorithm [28]. It can be observed from the first column of Fig. 27 that a positive LLE is estimated for U = 6.5, representing chaotic dynamics as the trajectories diverge exponentially within a bounded volume of the phase space. It is to be noted that the very small positive value of LLE (= 0.00033) represents weak chaos. On the other hand, the quasi-periodic and periodic dynamics are categorized by zero (for U = 5.8) and negative LLE (for U = 7.0), respectively.

C Characterization of the aperiodic dynamics
Correlation dimension based on Grassberger-Procaccia algorithm [15] has been determined next to confirm the dynamical signatures of the attractors observed in the reconstructed phase-portraits. It helps to identify the chaotic oscillations, characterized by the presence of a strange attractor in the phase space with a non-integer correlation dimension. As presented in the second column of Fig. 27, chaos is characterized by a correlation dimension of 2.51 for U = 6.5, depicting the fractal nature of the chaotic signal. However, the quasi-periodic and periodic dynamics are characterized by the integer correlation dimensions of 2 and 1, respectively, for U = 5.8 and U = 7.0. The correlation dimension of 2 for a quasi-periodic attractor corresponds to the two-dimensional toroidal attractor in the phase space. On the other hand, a one-dimensional periodic attractor is denoted by a correlation dimension of 1.