Analysis of a piecewise linear aeroelastic system with and without tuned vibration absorber

The dynamics of a two-degrees-of-freedom (pitch–plunge) aeroelastic system is investigated. The aerodynamic force is modeled as a piecewise linear function of the effective angle of attack. Conditions for admissible (existing) and virtual equilibria are determined. The stability and bifurcations of equilibria are analyzed. We find saddle-node, border collision and rapid bifurcations. The analysis shows that the pitch–plunge model with a simple piecewise linear approximation of the aerodynamic force can reproduce the transition from divergence to the complex aeroelastic phenomenon of stall flutter. A linear tuned vibration absorber is applied to increase stall flutter wind speed and eliminate limit cycle oscillations. The effect of the absorber parameters on the stability of equilibria is investigated using the Liénard–Chipart criterion. We find that with the vibration absorber the onset of the rapid bifurcation can be shifted to higher wind speed or the oscillations can be eliminated altogether.


Introduction
Nonlinear aeroelastic phenomena affect several types of aeroelastic systems such as flexible wings, helicopter rotor blades and wind turbines. Nonlinear aeroelasticity studies the interactions between inertial, elastic and aerodynamic forces on flexible structures that are exposed to airflow and feature non-negligible nonlinearity [1]. The theory of aeroelasticity is extensively covered in the literature [2][3][4].
The sources of nonlinearity in aeroelastic systems include geometric nonlinearity, structural nonlinearity, flow separation, friction, free-play in actuators, backlash in gears, nonlinear control laws, oscillating shock waves and other nonlinear phenomena [5,6]. A comprehensive study for such nonlinearities was presented by Lee et al. in [7] together with the derivation of the equations of motion of a 2D airfoil oscillating in pitch and plunge. Dowell et al. [8] summarized the physical basis and the effect of nonlinear aeroelasticity on the flight and its association with limit cycle oscillations (LCO). The effect of structural nonlinearities on the dynamical behavior of the system was investigated in [9,10]. Aerodynamic nonlinearities were studied by Dowell et al. in [11] using the describing function method. A combination of structural nonlinearity and the nonlinear ONERA stall aerodynamic model was investigated in the aeroelastic response of a nonrotating helicopter blade in the work of Tang et al. [12]. The airfoil model consisted of a NACA 0012 profile (the same used in this paper) with three types of nonlinearities: nonlin-ear structure linear aerodynamics, linear structure nonlinear aerodynamics and nonlinear structure with nonlinear aerodynamics. CFD-based aeroelastic investigations are described in [13][14][15].
Gilliatt et al. [16] studied structural and aerodynamic nonlinearities arising from stall conditions. Experimental investigation of a NACA0012 airfoil undergoing stall flutter oscillations in a low-speed wind tunnel was presented by Dimitriadis et al. in [17]. Santos et al. [18] carried out an experimental and numerical study of a NACA 0012 airfoil under the influence of structural and aerodynamic nonlinearities due to dynamic stall effects at high angles of attack. Jian et al. [19] developed a first-order, state-space model by combining a geometrically exact, nonlinear beam model with nonlinear ONERA-EDLIN dynamic stall model to investigate high-aspect-ratio flexible wings.
Nonlinear aeroelastic forces can produce complex bifurcation scenarios leading to chaotic oscillations. Sarkar et al. [20] observed a period-doubling route to chaos in the case of a 2-DOF aeroelastic model using the ONERA nonlinear aerodynamic dynamic stall model. In a recent study, Bose et al. [21] confirmed the Ruelle-Takens-Newhouse quasi-periodic route to chaos in a nonlinear aeroelastic system.
Flexible aeroelastic systems can lose stability by flutter or divergence depending on the system parameters. When flexible airfoils experience wind excitation, dynamic stability loss (called flutter instability), triggered by a Hopf bifurcation, may occur. The flow velocity for which the instability starts is called flutter velocity.
The starting point of the paper by Kalmár-Nagy et al. [22] was the observation that the coefficient of lift function can be well approximated by a piecewise linear function. Piecewise linear functions are often used to model free-play, nonlinear stiffness, nonlinear aerodynamic forces, and hysteresis nonlinearity in aeroelastic systems [10,[23][24][25]. Free-play nonlinearity is often studied in aeroelasticity, as it is a common occurrence within the actuated control surface [26]. Sales et al. [27] studied an aeroviscoelastic system with a nonsmooth, free-play-type nonlinearities in their control surface. The dynamic aeroelastic response of multi-segmented hinged wings with bilinear stiffness characteristics was studied theoretically and experimentally in [28]. Dimitriadis et al. [29] investigated a complex mathematical wing model in unsteady flow with a piecewise linear nonlinearity.
Piecewise nonlinear functions are also used to analyze aeroelastic systems with complex piecewise nonlinear structural stiffness [30,31]. Sun et al. [32] approximated the static lift coefficient by a fourth-order piecewise polynomial function. A piecewise aerodynamic flutter equation was established by Goodman in [33], which uses a piecewise aerodynamic interpolation function. Replacing the nonlinearities of a dynamical system with piecewise linear makes the problem more analytically tractable. The phase space of piecewise linear systems consists of regions, each of which has linear dynamic equations of motion [34]. Magri et al. [35] proved that a typical airfoil section applying nonsmooth definition of the dynamic stall model can generate a nonsmooth Hopf bifurcation, similar to a rapid bifurcation.
Once the dynamics of the system is understood (at least on a rudimentary level), engineering questions can be addressed. Important questions are how to reduce the oscillation amplitude for a given range of parameters or how to control the flutter instability in an aeroelastic system. The linear tuned vibration absorber (LTVA), developed by Frahm [36], is a classical passive device for controlling flutter instabilities. This device consists of a small lumped mass attached to the primary structure through a linear spring and a damper. The first detailed study of properties and optimization of LTVAs was presented by Den Hartog [37]. If its parameters are correctly tuned, the LTVA can significantly shift the flutter speed [38]. The parameter optimization of nonlinear tuned vibration absorber for flutter control was carried out by Mahler et al. in [39]. Kassem et al. [40] presented an analytical and experimental study of an active dynamic vibration absorber for flutter control.
This work investigates the dynamics of a 2-DOF piecewise linear aeroelastic system for a broad range of the angle of attack. In Sect. 2, dimensional governing equations of the 2-DOF aeroelastic system are presented. The piecewise linear modeling of quasi-steady aerodynamic forces is summarized in Sect. 3. The piecewise linear models of lift coefficient as the function of the effective angle of attack were determined for NACA 0009, 0012, 23012 profiles. Section 4 is devoted to the nondimensionalization of the governing equations. In Sect. 5, the equilibrium points of the aeroelastic system are determined. The stability of the equilibria is investigated in Sect. 6. In Sect. 7, the classical and discontinuity-induced equilibrium bifurcations are studied. In the case of NACA 0012 and 23012 stall

The aeroelastic model
We analyze the two-degrees-of-freedom (2-DOF) aeroelastic system shown in Fig. 1. Coordinates y and α describe the vertical (plunge) displacement (positive downward) and angular (pitch) displacement (positive in the clockwise direction), respectively. The semichord of the airfoil is denoted by b. In this model, we assume that the center of gravity (denoted by G) is located at three quarters of the chord length. The elastic axis passes through the center of gravity. The mass of the wing is m, and I cg is the moment of inertia about the center of gravity. The linear spring constant is k y , and the damping is c y for the plunge DOF. The linear spring constant is k α , and the damping is c α for the pitch DOF. The aerodynamic lift force (positive upward) and moment (positive in the clockwise direction) applied at the airfoils aerodynamic center (denoted by A) are L and M. The system equations are given by mÿ + c yẏ + k y y = −L (C l (α eff )) , The aerodynamic lift force L and moment M are functions of the lift coefficient C l . The lift coefficient is The system parameters used in this study were taken from the paper by Gilliatt et al. [41] and are summarized in Table 1.

Aerodynamic forces
The aerodynamic lift force and moment on the righthand side of Eq. (1) are assumed to be proportional to the lift coefficient, i.e., where ρ is the air density, S is the wing span, and b is the semichord. The measured static lift coefficient versus effective angle of attack for different airfoil sections [42,43] is illustrated in Fig. 2 Parameter α + stall characterizes the stall condition for positive α eff at which lift starts to decrease as α eff is increased. Parameter α + switch corresponds to the switching point at which the slope of C l starts to increase again. The negative angles of attack α − stall and α − switch are defined analogously. We also impose the following continuity constraints to hold Symmetric profiles give rise to an odd coefficient of lift function with parameters The piecewise linear models of the lift coefficient for the NACA 0012, 0009 and 23012 profiles are shown in Fig. 2 were determined by least square method using piecewise linear model (4) and continuity constraints (5). The parameters are given in Table 2. We point out that the NACA 0012, 0009 profiles are symmetric, while NACA 23012 is asymmetric.
In the next section, the governing nondimensional equations are derived for symmetric airfoil profiles (for the nonsymmetric case the derivation can be done similarly).

Nondimensional equations for symmetric profiles
Nondimensionalization results in a better understanding of the physical phenomenon (via the nondimensional groups obtained) and in a reduction of system parameters [44]. Following [22] we choose the length scale L, the timescale T , and the nondimensional freestream velocity μ > 0 as These scales yield the nondimensional plungeŷ = y/L, the nondimensional time τ = t/T . The derivative with respect to the nondimensional time is denoted by () = d()/dτ . The nondimensional effective angle of attack can now be expressed as The nondimensional form of Eq. (1) iŝ where the nondimensional parameters p 1 , p 2 , p 3 and p 4 are Using trilinear approximation (4) of the lift coefficient for symmetric profiles (see Eq. (6)), we obtain the following equations: for α stall ≤ |α eff | ≤ α switcĥ and for α switch ≤ |α eff | We introduce the state vector the system matrices and translation vectors Equations (11)-(13) can now be concisely written as a collection of affine systems, together with their domains of validitẏ where the domains of subsystems (17)-(21) are given by These domains are separated by the switching planes defined by With these pairwise disjoint sets, the full state space R 4 can be decomposed as

Equilibrium points
The equilibrium points of system (17)-(21) are formally given byẋ = 0. Clearly, an equilibrium point only exists if it is inside the corresponding domain of validity. Following Di Bernardo et al. [45] we refer to such an equilibrium point as admissible. We call an equilibrium point virtual if it is not admissible (i.e., it satisfiesẋ = 0 but falls outside the domain of validity). Virtual equilibria play an important role in organizing the dynamics of the corresponding region [46]. The admissible equilibrium points of system (17)-(21) are We express the admissibility conditions embedded in Eqs. (30 )- (32) in terms of the system parameters.
E 0 is always admissible for symmetric airfoil profiles, From admissibility condition (34) we obtain that equi- One can utilize continuity condition (5) and symmetry condition (6) to express μ 1 as Let us define (provided c 2 > 0) From condition (35) we find that equilibrium E + 2 is The admissible equilibria of the symmetric aeroelastic models

Classical and discontinuity-induced bifurcations
The results in this section pertain to structural and aerodynamic parameters of the model with a NACA  (Tables 1 and 2). This profile is symmetric; therefore, we only state the results for equilibria with the + superscript. In piecewise linear aeroelastic system (17)-(21) discontinuity-induced bifurcations [45,46,48,49] are also observed in addition to their classical counterparts: border collision and rapid bifurcation. Border collision is a boundary equilibrium bifurcation where virtual equilibrium points become admissible. The rapid bifurcation is the abrupt appearance of stable periodic oscillations [50,51]. The rapid bifurcation of a piecewise linear system is also called a focuscenter-limit cycle bifurcation, describing its bifurcation scenario [52]. As shown in Sect. 6, loss of stability of equilibria corresponds to the values of dimensionless freestream velocity where μ 1 = 0.215, μ RB = 0.304 and μ A = 0.321. The first classical bifurcation occurs when equilibrium E 0 loses stability at μ 1 = 0.215. At this bifurcation point the roots of the characteristic equation R 0 (λ, μ 1 ) = 0 (Eq. (46)) are λ 1 = 0, λ 2 = −0.059, λ 3 = −0.081 − 0.996i, From these eigenvalues we conclude that equilibrium E 0 loses its stability by undergoing a saddle-node bifurcation.  Fig. 7.
The third bifurcation is a rapid bifurcation [53] of equilibrium E + 1 at μ RB = 0.304. At the rapid bifurcation point the roots of the characteristic equation At this point the pair of complex conjugate roots λ 1 , λ 2 of the characteristic equation R 1 (λ, μ) = 0 crosses the imaginary axis with positive speed, similarly to the Hopf bifurcation for smooth systems. Figure 8 demonstrates that at the rapid bifurcation the stable spiral equilibrium E + 1 becomes unstable, and a finite-amplitude limit cycle is born. We note that the limit cycle amplitude decreases by increasing the bifurcation parameter μ.
The fourth bifurcation of system (17)-(21) is a border collision (BC) of equilibria E + 1 and E + 2 at μ 2 = 0.391. At this border collision the admissible equilibria E + 1 and E + 2 become virtual by simultaneously crossing the switching plane + 2 . This type of border collision (when two admissible equilibria simultaneously cross the same switching plane and become virtual) is also called as a nonsmooth fold bifurcation of piecewise systems [45].Typical trajectories of system (17 )- (21) before and after the second-border collision bifurcation (μ 2 = 0.391) are illustrated in Fig. 9.
The bifurcation diagram of system (17)-(21) is shown in Fig. 10, where the solid line indicates the stable, and the dashed line the unstable equilibrium points. The bifurcation points are denoted by dots.
Even though the derivation here does not include formulae for the NACA 0009 and NACA 23012 profiles, we nevertheless show their bifurcation diagrams in Fig. 11.
The above results were confirmed numerically using Mathematica. The built-in event location method with a maximum step size of 0.001 was used to solve the piecewise linear system. For a given μ, a onedimensional Poincaré section (the set of points for which α = 0) was computed, and these sections were spliced together to obtain the diagrams as a function of μ (see Fig. 12). The numerical results were overlaid on top of the bifurcation diagrams of Figs. 10 and 11.

Stall flutter reduction with a linear tuned vibration absorber
Stall flutter is a dynamic aeroelastic instability resulting in unwanted oscillatory loads on the structure [54]. Stall flutter of aircraft wings is caused by nonlinear aerodynamic characteristics at very high angles of attack, when a partial or complete separation of the fluid flow from the airfoil occurs periodically during the oscillation [3,55]. Static divergence can cause asymmetric stall flutter oscillations around nonzero pitch angles [17,56]. Piecewise linear model (17) (Table 2) Increasing the nondimensional freestream velocity above the rapid bifurcation point (for NACA 0012 profile aerodynamic parameters) the system starts to oscillate. To eliminate these stall flutter oscillations we applied a linear tuned vibration absorber (LTVA) to the primary system, see Fig. 13. The absorber is attached along the mid-chord of the airfoil at distance z from the center of gravity (G). It is composed of a mass m ltva , a spring of linear stiffness k ltva and a linear damper c ltva . The equations of motion of the aeroelastic system with the absorber are mÿ +c yẏ +k y y − f (y,ẏ, α,α, h,ḣ) =−L (C l (α eff )) , f (y,ẏ, α,α, h,ḣ) = M (C l (α eff )) ,  y,ẏ, α,α, h,ḣ Equations (63) and (64) are nondimensionalized by a length scale L, a timescale T and a nondimensional freestream velocity μ > 0, given in Eq. (7). We define the damping ratio ξ = c ltva /c y , the stiffness ratio η = k ltva /k y and the mass ratio ε = m ltva /m. The nondimensional distance of the absorber from the center of gravity is defined as ζ = z/L. The nondi-mensional displacement of the absorber is denoted bŷ h = h/L. The nondimensional equation will be the followinĝ wherê f (ŷ,ŷ , α, α ,ĥ,ĥ ) = ξ p 1 ĥ − (ŷ − ζ α ) and the nondimensional coupling ratio is We extend the state vector as We define the new system matrices The equations of the extended system arė where the domains of subsystems (71)-(75) are given by These domains are separated by the switching planes defined by With these pairwise disjoint sets, the full state space R 6 can be decomposed as

Equilibrium points
The admissible equilibrium points of subsystems (71)-(75) are The expression of equilibria F 0 , F ± 1 , F ± 2 and their range of validity (Eqs. (84)-(86)) determine the admissible ranges depending on the bifurcation parameter μ as F 0 is always admissible for symmetric airfoil profiles, (87)

Effects of the absorber parameters
In this section we investigate the effects of the absorber parameters on the stability of equilibria F 0 , F ± 1 , F ± 2 . To reduce the number of parameters we fixed the mass ratio ε = 0.1. The other absorber parameters were varied in the following intervals: First, we investigate the asymptotic stability of these equilibria of the aeroelastic system with the LTVA. To determine the asymptotic stability of the equilibrium points, we study the characteristic polynomial of the coefficient matrices Q k defined in Eq. (69) where the parameters ε, ζ, η, ξ , but this dependence is not explicitly shown to keep the notation short. Applying the Liénard-Chipart stability criterion (see Eq. (45)), the aeroelastic system with the linear tuned vibration absorber is asymptotically stable if the following conditions are fulfilled q 2 (k, μ) > 0, q 4 (k, μ) > 0, q 6 (k, μ) > 0, (97) 5 (k, μ) > 0, (98) where the Hurwitz determinants 3 (k, μ) and 5 (k, μ) are calculated according to Eq. (45).
The range of asymptotic stability μ ∈ [μ 1 , μ RB ) of equilibria F ± 1 can, however, be extended by tuning the parameters of the absorber. In other words, we can shift μ RB to larger μ values by appropriate absorber parameters. The stable regions of equilibria F ± 1 in the 2-dimensional parameter space (η, μ) for different ξ and ζ values, determined by conditions (97)-(98), are illustrated in Fig. 14. The stable regions of equilibria F ± 1 can also be visualized in the 3-dimensional parameter space (ξ, η, μ) for ζ ∈ {0, 0.01, 0.05, −0.01}, see Fig. 15. The range of admissibility and instability μ ∈ (μ A , μ 2 ] of equilibrium F ± 2 is the same as that of E ± 2 . This is the consequence of which means that the zero of q 6 (2, μ) coincides with that of a 4 (2, μ) at μ = μ A . Even though the derivation here does not include formulae for NACA 23012 profile, we nevertheless show the stable regions of equilibria F + 1 in Figs. 16 and 17.

Quenching oscillations with the absorber
In this section we demonstrate that oscillations can be eliminated in certain regions of the parameter space of the absorber. Once again we use the structural and aerodynamic parameters of the model with NACA 0012 profile (Tables 1 and 2). Using absorber parameters we draw the bifurcation values μ 1 , μ RB , μ 2 as the functions of η see Fig. 18a. For parameter values (101) and η = 0.05 the rapid bifurcation is eliminated and the limit cycle oscillations (see Fig. 18b). For parameter values (101) and η = 0.12 the value of μ RB is increased and the limit cycle oscillations still exist if μ ∈ [μ RB , μ 2 ) (see Fig. 18c). Typical trajectories with and without absorber are shown in Fig. 19.

Conclusions
A piecewise linear aeroelastic system with and without a tuned vibration absorber was studied. A piecewise linear model was utilized using experimental data for the lift coefficient versus the angle of attack for NACA 0012, NACA 0009 and NACA 23012 airfoil profiles. The equations of motion for the system were nondimensionalized for symmetric airfoil profiles. The nondimensional freestream velocity was introduced as the system bifurcation parameter. The equilibria of the piecewise linear aeroelastic system were determined analytically as a function of the bifurcation parameter. The admissibility conditions of the equilibria were analyzed. The stability of those equilibrium points was studied by applying the Liénard-Chipart criterion. The bifurcation analysis found classical and discontinuity-induced bifurcations as saddle-node, border collision and rapid bifurcation. The saddle-node bifurcation corresponds to the static divergence of the aeroelastic system. The rapid bifurcation is the abrupt appearance of stable periodic stall flutter oscillations of the aeroelastic system. The bifurcation analysis showed that the pitch-plunge model with a simple piecewise linear approximation of the aerodynamic force can reproduce the complex aeroelastic phenomenon of stall flutter. Finally, the effect of absorber parameters was investigated. The stall flutter oscillations were eliminated using a well-tuned absorber. This parametric study of the linear tuned vibration absorber predicts that not only active [32,54] but also passive vibration reduction methods can be applied to suppress stall flutter.
Acknowledgements Open access funding provided by Budapest University of Technology and Economics (BME). The publication of the work reported herein has been supported by the ÚNKP-19-3 New National Excellence Program of the Ministry for Innovation and Technology of Hungary. The research reported in this paper was supported by the Higher Education Excellence Program of the Ministry of Human Capacities in the frame of the Water Sciences & Disaster Prevention research area of BME (BME FIKP-VÍZ). The research reported in this paper has been supported by the National Research, Development and Innovation Fund (TUDFO/51757/2019-ITM, Thematic Excellence Program).
Funding Open access funding provided by Budapest University of Technology and Economics (BME).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.