Identifying route to stall flutter through stochastic bifurcation analysis

The interaction of an elastic structure such as an airfoil and fluid flow can give rise to nonlinear phenomenon such as limit cycle oscillations, period doubling or chaos. These phenomena are indicated by a change in the stability behaviour of the dynamical known as bifurcations. Presence of viscous effects in the fluid flow can give rise to flow separation which causes a stability change in the system that is identified to happen via a Hopf bifurcation. In such cases, the airfoil exhibits limit cycle oscillations which are torsionally dominant, known as stall flutter. Despite identifying the route to stall flutter under uniform flow conditions, investigating a stall problem under stochastic wind has received minimal attention. The ability of fluctuating flows to change the stability boundaries and disrupt the route to flutter, compels the need to carry out a stochastic analysis of stalling airfoils. Study of stall flutter in such systems under the influence of a time varying sinusoidal gust is undertaken and the route to flutter is identified by carrying out a stochastic bifurcation analysis.


Introduction
Aeroelasticity is the study of interaction of aerodynamic and elastic forces that arise when an elastic structure such as airfoil is subjected to a fluid flow [1,2].The interaction of these forces can give rise to complex nonlinear phenomenon such as stall, limit cycle oscillations or chaos [2,3] and are indicated by a change in the qualitative features in the dynamical system known as bifurcations.Bifurcation is a change in the stability behaviour of a nonlinear dynamical system that could result in qualitatively different dynamical responses, leading to topological changes in the phase space, as one or more parameter(s) which the system is dependent on is/are varied [4,5].
A nonlinear study of aeroelastic systems accounts for presence of viscous effects, freeplay, hardening and softening springs, etc [1,2].Viscous effects in the fluid flow can give rise to flow separation [6] which causes large negative moment associated with dynamic stall followed by flow reattachment that could induce torsionally dominant limit cycle oscillations.This phenomenon is known as stall flutter and the route to stall flutter happens via a Hopf bifurcation [7].Stall flutter can be observed in systems such as wind turbines, helicopter blades, turbomachinery, etc.The onset of a limit cycle oscillations is undesirable in such systems as these oscillations can induce fatigue damage to the structure that can affect the structural integrity or the system may even abruptly fail [8][9][10].Consequently, investigations into the bifurcation and stability characteristics of aeroelastic systems under deterministic flows have received substantial attention in the literature.In reality, the fluid MATEC Web of Conferences 211, 02011 (2018) https://doi.org/10.1051/matecconf/201821102011VETOMAC XIV flow is random and a mean flow assumption would not be sufficient enough for bifurcation analysis.A stochastic bifurcation analysis is undertaken to address this issue.This paper is devoted towards investigating the stochastic bifurcation scenarios in a pitching airfoil undergoing stall flutter.A one degree of freedom pitching airfoil is considered.The pitching stiffness is incorporated by attaching a torsional spring at the elastic axis.The input flow is modelled to be a randomly time varying quantity.With the mean speed as the bifurcation parameter a stochastic bifurcation analysis is systematically carried out.The study concludes by providing insights into the pre-flutter oscillations such as intermittency in light of the stochastic attractors.
The organization of this paper is as follows.Section 2 describes the mathematical model of the stochastic stall flutter problem.The dynamical responses obtained and the bifurcation scenarios in a stalling airfoil along with the associated discussions are provided in Section 3. The salient features of this study are summarized and concluded in Section 4. A one degree of freedom model of an airfoil as shown in Fig. 1 is considered, since the predominant mode of vibration in case of stall flutter is pitch [7].Above shown is a NACA 0012 airfoil having chord length 'c'.'E' is the elastic axis, U is oncoming fluid flow and 'α' is the angle of attack.Kα and Cα are the stiffness and damping coefficient of the spring which enables pitching.The elastic axis is assumed to be at 1/4 th of the chord i.e. at a distance of c/4 from the leading edge.

Mathematical model
The equation of motion of a pitch airfoil in presence of a linear structural damping coefficient ζα and nonlinear cubic spring stiffness βα, can be expressed in non-dimensional form as where, (′) indicates derivative with respect to non-dimensional time τ which is defined as τ = (t ̅ )/b, where b is the airfoil semi-chord.U is the non-dimensional wind velocity defined as U =  ̅ /  .The structural parameters λ and rα are the mass ratio and the radius of gyration respectively which are defined as λ= / 2 and rα = Iα/ ( 2 ) where 'm' is the mass per unit span, Iα is the mass moment of inertia of the airfoil, ρ is the free-stream density and Cm is the aerodynamic moment coefficient.
During the dynamic stall of an airfoil, various nonlinear interactions take place and the flow field is highly unsteady due to the growth, evolution and subsequent shedding of leading edge vortex [6].Accurate modelling of the flow requires solving Navier-Stokes equations, which is a tedious task, so an available semi empirical model, Beddoes Leishman model [11] is employed for analysis.Beddoes () =    () +    () where,    () is the unsteady lift coefficient due to trailing edge flow separation,    () is the forcing term accounting for increased lift due to leading edge vortex and    ()    () are the moment coefficients due to the trailing edge separation and dynamic stall phase.The speed of the oncoming flow as discussed is fluctuating with time.So, it is modelled using the function, where, h(τ) accounts for the randomness in the airflow and is represented as where,  is the intensity of the fluctuations,   is the mean speed of airflow,   () is a time varying frequency of the sinusoid, such that   () =   +  1 ().Here,   is the mean frequency of oscillations, () is a random number and  1 is the scaling frequency.

Results and discussion
The model is validated against results obtained in [11] where the change in Cm with change in angle of attack of a one dof pitch airfoil undergoing dynamic stall is shown in Fig. 2.

Fig. 2. Variation of Cm with respect to 𝛼𝛼
The equations ( 2) and (3) are calculated using the coefficients modelled in different modules and the equation of motion for pitch (1) is solved using Runge-Kutta method in MATLAB using the ODE45 function with a sample time step of 0.1s with the empirical constants and equations ( 2 As observed from Fig. 3, when the incoming flow speed U is lower than the critical speed Ucr, the pitch response decays down to a stable attractor at the origin (see Fig. 3a).However, when U ≥ Ucr the response transforms itself into sustained LCOs (see Fig. 3b) via a supercritical Hopf bifurcation.Similar observations have been reported in the literature [6,12].It is worth noting that the Hopf bifurcation route to stall flutter well explains the onset of instability only when the flows are uniform.However, realistic flows often possess small temporal fluctuations.Consequently, input flows to aeroelastic systems are modelled as a second order random process [13].These fluctuations can disrupt the Hopf bifurcation route to flutter [3] and also alter the stability regimes [13].Noting the significance of modelling the flow to be a random process, the governing aeroelastic equation (Eq.1) is solved using the sinusoidal gust (modelled using Eq.4) as input.With mean flow speed as bifurcation, sample time histories are obtained and shown in Fig. 4. A visual inspection of the time responses in Fig. 4 reveal that they are qualitatively distinct from the single step Hopf bifurcation encountered in the deterministic counterpart.Fig. 4a.
shows an intermittent switching between varying amplitude periodic oscillations.When Um approaches Ucr the system exhibits random LCOs where in the noise drives the system between two attractors (as shown in Fig. 4b).Indeed, there is a regime of pre-flutter oscillations that possess sporadic bursts of periodicity amidst aperiodic fluctuations.These oscillations are called intermittency [5] and are encountered in dynamical systems in the presence of a noisy input parameter [13,14].To the best of author's knowledge, the present study is the first to show an intermittency regime presaging stall flutter.On one hand, it can be argued that the noisy system parameter (here the input wind) drives the system from one MATEC Web of Conferences 211, 02011 (2018) https://doi.org/10.1051/matecconf/201821102011VETOMAC XIV basin of attractor to another in a random fashion resulting in the formation of alternating dynamics in the same time response [15].On the other hand, recent studies by Venkatramani et al. [14] have demonstrated that the cumulative unsteady wake effects play a significant role in the appearance of intermittency.In the case of our present study, wherein the flow is unsteady and separated, further investigations are required to pinpoint on the mechanism behind the appearance of intermittency.Increasing Um above the critical speed, large amplitude LCOs are encountered indicating the onset of stall flutter.
In the light of stochastic responses under consideration, it is imperative that any comments on stability or bifurcation characteristics would demand a stochastic bifurcation analysis.It is noteworthy that for noisy dynamical systems, demarcating regimes of decaying oscillations against those of LCOs are not well defined.As seen earlier, the presence of noisy parameter can drive the dynamics back and forth about different attractor regimes in an unpredictable fashion.Consequently, typical characterization of stochastic responses are in terms of its joint-probability density function (j-pdf) of its state variables.As a bifurcation parameter is varied, the structure of the j-pdf undergoes a visual topological change, marking a change in the stability characteristics of the system -and is called Phenomenological or Ptype bifurcations.By estimating the j-pdf of the state variable, one can statistically discern the average time spent by trajectories of a dynamical system in a volume element of the state space.Therefore, using the obtained time responses in Fig. 4, the j-pdfs of the pitch and pitch rate are computed and are shown in Fig. 5.
A visual inspection of Fig. 5a indicates that the j-pdf has a centred peak about the mean amidst a circular ring like structure with relatively lower strength.This is reflective of an alternating dynamic between a weakly stable limit cycle attractor amidst a relatively strong fixed point attractor.Indeed, the corresponding time history shown in Fig. 4a represents the airfoil undergoing an intermittent response.As the mean flow speed increases, the occurrence of periodic bursts increases, marked by an increase in the strength of the limit cycle attractor along with a drop in the strength of the fixed point attractor; see Fig. 5b.Further increase in Um results in the complete destruction of the fixed point attractor and the limit cycle attractor alone prevails as shown in Fig. 5c.J-pdfs are only a qualitative measure for intermittency and a quantitative representation is required to ascertain the regimes of intermittency and random LCOs.Quantitative measures such as Shannon entropy as shown in [15] will provide deeper insights into stochastic stability regimes and can be explored in future studies.
Fig. 3a.Time history for U < Ucr Fig. 3b.Time history for U ≥ Ucr

Fig. 5 .
Fig. 5. Joint probability density function of pitch responses for (a) Um<<Ucr (b) Um<Ucr (c) Um>Ucr Leishman model divides the stalling of an airfoil into different flow modules which are, attached flow module, vortex formation phase, vortex shedding phase and the flow reattachment module.The lift and moment coefficients MATEC Web of Conferences 211, 02011 (2018) https://doi.org/10.1051/matecconf/201821102011VETOMAC XIV are modelled for each module and the total lift and moment coefficients are given by the relation,