Stochastic Bifurcation Analysis of an Elastically Mounted Flapping Airfoil

The present paper investigates the effects of noisy flow fluctuations on the fluid-structure interaction (FSI) behaviour of a span-wise flexible wing modelled as a two degree-of-freedom elastically mounted flapping airfoil. In the sterile flow conditions, the system undergoes a Hopf bifurcation as the free-stream velocity exceeds a critical limit resulting in a stable limit-cycle oscillation (LCO) from a fixed point response. On the other hand, the qualitative dynamics changes from a stochastic fixed point to a random LCO through an intermittent state in the presence of irregular flow fluctuations. The probability density function depicts the most probable system state in the phase space. A phenomenological bifurcation (P-bifurcation) analysis based on the transition in the topology associated with the structure of the joint probability density function (pdf) of the response variables has been carried out. The joint pdf corresponding to the stochastic fixed point possesses a Dirac delta function like structure with a sharp single peak around zero. As the mean flow speed crosses the critical value, the joint pdf bifurcates to a crater-like structure indicating the occurrence of a P-bifurcation. The intermittent state is characterized by the co-existence of the unimodal as well as the crater like structure.


Introduction
Fluid-structure interaction (FSI) is known to lead to undesirable phenomena like aero-elastic flutter which is usually detrimental for engineering systems.On the contrary, biological flyers take the advantage of FSI to augment their propulsive efficiency by exploiting the coupling between the flexible wings and the surrounding unsteady flow-field.A proper understanding of the role of FSI in natural flights may directly benefit biologically-inspired design of Micro Aerial Vehicles (MAVs) [1].Moreover, the performance of very light weight flapping wing MAVs is significantly altered in the presence of gusty flow-field.
Flapping wing MAVs have received significant research attention in recent years due to their potential application in surveillance and environmental monitoring [2,3].Extensive experimental and computational studies have been carried out to understand the unsteady aerodynamics of flapping wings [4][5][6].However, in most of these studies, the effects of flow fluctuations on the flapping flight are ignored by considering a uniform flow condition.However, MAVs often need to operate in severe weather conditions.Therefore, the consideration of realistic flow fluctuations is crucial to investigate the performance of flapping wing MAVs.
There are substantial literature available regarding the effect of wind gust on the dynamical behaviour and structural safety of aircrafts at high Reynolds number regime [7][8][9].However, MAVs are more susceptible to gusty flows due to their low inertia and low flight speed.The studies investigating the effect of realistic flow fluctuations in flapping wing aerodynamics are limited.Lian and Shyy [10] observed that a flapping airfoil can significantly alleviate the impact of wind gust as compared to fixed airfoil thus potentially benefiting the MAV design.Shyy et al. [11] compared the aerodynamic performance of rigid and flexible airfoils subjected to gusty flows by carrying out numerical simulations as well as experiments and observed that a flexible airfoil is able to maintain high lift to drag ratio as compared to a rigid airfoil.Thus, flexibility can make the flapping flight gust tolerant.Watkins et al. [12] reported that the wind gust made up of small scale eddies produce uneven lift distribution over a flapping wing, causing a rolling motion.On the other hand, gusts having large scale eddy motion cause pitch motion.Lian and Shyy [13] observed that both the lift and drag coefficients are altered significantly under gust environment.Wind tunnel tests are also carried out to study the impact of gusty flows on MAVs [14,15].However, it is to be noted that most of these studies have been carried out based on rigid wing assumption, where wing flexibility is ignored.The inclusion of flow fluctuations in the flapping wing aerodynamic models poses several computational challenges as the difference in the gust frequency and the flapping frequency leads to a multi-scale problem [10].The different time scales present in the problem increases the computational cost.The fluid-elastic coupling adds to the difficulty level.Although, wind gusts encompass a broad range of frequencies and amplitudes, majority of the available studies con-MATEC Web of Conferences 148, 08001 (2018) https://doi.org/10.1051/matecconf/201814808001ICoEV 2017 sider a simplified model of harmonic gust with a single frequency [10,16].The aerodynamic performance of flexible airfoils under irregular flow-fluctuations is still not well understood.Therefore, a stochastic bifurcation analysis is essential to carry out to understand the effects of irregular fluctuations in the flow on the dynamical stability characteristics of the coupled system.The present study is focussed to investigate the behaviour of flexible wing section under irregular gust.A two degrees of freedom model with pitch-plunge flexibility has been considered to investigate the dynamical behaviour.The irregular gust is modelled as stochastic input and the wing output behaviour is examined qualitatively through a Phenomenological bifurcation (P-bifurcation) analysis [17].P-bifurcation behaviour is indicative of stochastic bifurcation.
In this present study, the large deformation of a spanwise flexible wing has been modelled by incorporating cubic nonlinearity in the rotational spring along the pitch (torsion) degree of freedom.The nonlinear structural model has been coupled with an inviscid flow solver using a weak coupling strategy to build the FSI framework.The equations related to flow are solved by unsteady vortex lattice method (UVLM).A stochastic bifurcation analysis has been carried out for the two degree-of-freedom pitch-plunge airfoil subjected to noisy flow fluctuations in inviscid condition, with the mean speed of the flow taken to be the bifurcation parameter.
The remainder of the paper consists of the following sections: Section 2 describes the FSI framework and the gust model; P-bifurcation results are presented in Section 3. The paper ends with the concluding remarks at Section 4.

Computational methodology
In the present FSI framework, the 2-DOF pitch-plunge flexible flapping system has been coupled (using a partitioned approach based 'weak coupling' strategy) with a potential flow solver where both flow and structure solvers exchange information explicitly at each time step.

Pitch-plunge flexible flapping model
The nonlinear structural model consists of a rigid NACA 0012 airfoil elastically supported by a cubic nonlinear rotational spring and a translational spring along pitch and plunge degrees of freedom (see Fig. 1).The structural model takes limited mode spanwise flexibility into account along the translational plunge (y) and rotational pitch (α) direction.The elastic center where the axis system is located is at a distance 0.25b towards the foil leading edge from the mid-chord with the chord length being c = 2b.The center of mass (c.m.) is located at a distance x α b from the elastic axis toward the foil trailing edge.The direction of positive lift (L) and moment (M) is indicated in Fig. 1.The nondimensionalized structural equations of motion are given as [18], x α r α 2 h + α + 2 ( (2) The non-dimensional variables are defined as follows: Here, b is the semichord, r α b is the radius of gyration, x α b is the distance between the center of mass and elastic center; ω ξ and ω α are the uncoupled natural frequencies of plunge and pitch respectively; cubic nonlinear co-efficient β α determines the extent of non-linearity in the spring stiffness along the pitch degree of freedom; v ∞ is the free-stream; C l and C m are the unsteady lift and moment coefficients about the aerodynamic center of the airfoil and are computed at each time step using an unsteady vortex lattice method (UVLM) based potential flow solver.The unsteady potential solver details are given in the next subsection.The structural responses are obtained by numerically integrating Eq. 1-2 using an explicit fourth-order Runge-Kutta method.The time step for integrating structural equations of motion is taken to be equal to that of the flow solver.

Flow solver
An 2D unsteady vortex lattice method (UVLM) based potential flow solver is implemented following the Hess and Smith approach to estimate the aerodynamic loads in the coupled system governing equations [19,20].The flow solver assumes the flow to be inviscid, incompressible and irrotational and is governed by the Euler equation, where u is the flow velocity, p is the pressure and ρ is the fluid density.A schematic of the UVLM approach has been presented in Fig. 2. Unlike the classical analytical approach (Theodorsen's Method [21]), in this method, the unsteady wake is not assumed to be rigid and is free to evolve with its own local velocities.Hence, the actual shape of the wake is captured.The wake consists of discretized wake elements.Furthermore, the actual shape of the airfoil is considered to compute the flow-field variables.The airfoil surface is divided into 'N' number of small segments called panels.In the present study, the number of panels (N) is chosen to be 400 based on a panel convergence study.Each panel is represented using two types of singularity elements, sources and vortices.The source singularity strength is considered to be constant over a particular panel (q j , j = 1, 2, ..., N) and the vorticity strength is considered to be constant over all the panels (τ).For unsteady flows, the time dependent source strengths and the time dependent vorticity strength are denoted as (q j ) k and τ k respectively for the time step index k.The velocity at any point (x,y) in the flow-field is the vector sum of velocities of the undisturbed flow (freestream) and the disturbance field due to the presence of the oscillating body and also the wake behind the body.The boundary condition that the surface of the body is a streamline is satisfied by taking the summation of velocities induced by body bound singularities, free-stream and wake vortices to be equal to zero in the direction normal to the surface at each discretized body element or panel as, (4) Here, A n i j and B n i j are the influence co-efficients at the i th panel control point by a unit strength source distribution and unit strength vorticity distribution on the j th panel respectively.v str denotes the unsteady stream velocity at body fixed co-ordinate system including airfoil motion.(τ w ) k is the uniform vorticity distribution over the wake panel attached to the trailing edge at t k .(B n i,n+1 ) k and (C n im ) k are the normal velocity induced at the i th panel control point by unit strength vorticity distribution on the shed wake panel and unit strength m th wake vortex at time-step t k respectively.Γ m and Γ m−1 are the overall circulation on the airfoil surface at time-steps t m and t m−1 respectively.The control point is situated at the mid point of each panel where the zero normal flow boundary condition is satisfied.The unknown body bound vortex strength, unknown source strengths and the immediately shed wake vortex at the current time step are computed using the no-normal boundary condition along with the additional condition of Kelvin's circulation theorem given by, where ∆ k is the length of the wake panel attached to the trailing edge at t k .Kelvin's theorem states that the total circulation in the flow-field must be preserved and that any changes in the circulation about the body are balanced by an equal and opposite circulation added in the wake in the form of new vortices.These vortices influence the local velocity field significantly, and as a result the forces on the airfoil at any instant become a function of the past motion of the airfoil.All the wake vortices, shed from the trailing edge of the airfoil, are carried downstream by their local velocities induced by the free-stream, body bound vortices and the other shed vortices.In unsteady potential flow, the calculation of pressure at any point on the body is done by using unsteady Bernoulli's equation [20].The pressure coefficients at the i th panel control point at time-step t k can be written as, where v t is the tangential velocity induced at the i th panel; (φ i ) k and (φ i ) k−1 are the velocity potential of the i th panel at time-steps t k and t k−1 respectively.The small effective angle of attack of the system justifies the use of the potential flow solver in the present study.The present UVLM based potential flow solver has been quantitatively validated by comparing the peak lift coefficients with the Navier-Stokes (N-S) results reported by Young [22] as well as with lumped vortex method (LVM) results computed by an inhouse code (see Fig. 3).LVM is a similar panel method based potential flow solver [23].In LVM, each panel has a lumped vortex at the control point situated at one-fourth of its length.The zero normal flow boundary condition is satisfied at the collocation point situated at three-fourth of its length.It can be seen that the present results show an excellent match with the LVM results and are also in close agreement with the N-S results especially at low κh values.

Noisy flow fluctuations
The fluctuating flow field is defined in the present study as, where u(τ) is the irregular flow fluctuation modelled as Ornstein-Uhlenbeck (O-U) process [24,25].The O-U is a stationary Gauss-Markov process and is defined by the following stochastic differential equation: du(τ) = −∆ωu(τ) dτ + ∆ωq dW(τ); (8) where W(τ) is a Wiener process, ∆ω and q are noise parameters.The correlation function of the O-U process is given by, The O-U process is simulated by numerically integrating the stochastic differential equation using the Euler-Maruyama method [26].It is to be noted that ∆ω is independently varied keeping q fixed to vary the correlation of the O-U noise.The variation of the correlation function of the O-U noise with different values of ∆ω for q = 1 is presented in Fig. 4. It can be seen that as the ∆ω is increased to 5 from 1, the correlation length of the process is decreased significantly.Thus, this noise model allows us to control the time-scale present in the fluctuations by changing ∆ω of the process.The time histories of the noisy inflow with different noise correlation (∆ω = 1 and 5) have been presented in Fig. 5.It is observed that the time-scales present in the noisy inflow drastically change as the correlation length of the noise is varied.A long time-scale gusty inflow is obtained at ∆ω =1; it becomes a short time-scale gusty inflow at ∆ω = 5 as the correlation of the O-U noise is decreased.

Bifurcation in uniform flow
The bifurcation behaviour of a pitch-plunge flexible flapping system in the presence of uniform flow is well investigated in the literature.The bifurcation diagram of the pitch response of the coupled system is presented in Fig. 6 considering the non-dimensional free-stream velocity (U * ) to be the control parameter.It is seen that the coupled system response damps down to a fixed equilibrium point (h = 0) below a critical value of the control parameter U * .The system undergoes a supercritical Hopf bifurcation at U * = U * cr and stable limit-cycle oscillations (LCO) emerge at U * > U * cr .The Hopf bifurcation point (U * cr ) is observed to be 3.75 from Fig. 6.For U * > U * cr LCOs persists and the amplitude of LCOs increase with U.
The pitch response at U * = 5 has been presented in Fig. 7 which shows a typical LCO in the post Hopf-bifurcation regime.

Bifurcation in fluctuating flow
This section presents the results of the bifurcation analysis of the flexible flapping system in the presence of irregular flow fluctuations.Typical atmospheric wind-spectra modelled using Dryden or von Kármán spectra are known to capture the long time scales.Hence, a long time-scale gust (O-U noise with ∆ω = 1, q = 1) has been imposed on the present system.The mean of the fluctuating flow speed (U * mean ) has been taken to be the bifurcation parameter.The pitch time histories are presented at different U * mean values in Fig. 8.It can be seen that the coupled system response does not transition from a damped oscillatory state to well developed limit cycle oscillations abruptly beyond a particular critical mean flow velocity.Instead, the system response gradually transforms from a stochastic fixed point (at U * mean =2) to a randomly modulated LCO (at U * mean =6) through irregular alternations between two qualitatively different dynamics as the bifurcation parameter is increased.In the dynamical systems theories, this kind of transitional dynamics is known as intermittency.In the presence of long time-scale fluctuations, an 'on-off' type intermittency phenomenon [27,28] has been observed where the system response is seen to irregularly alter between an 'off' state (very low amplitude aperiodic fluctuations) and an 'on' state (high amplitude periodic bursts) at U * mean = 4.0 and 4.5.As the mean flow speed is gradually increased, the segments of periodic bursts are seen to be more pronounced.Eventually the periodic oscillations completely dominate the time history as random LCOs on further increasing the mean flow speed to U * mean = 6.
To obtain a better insight into the system dynamics, a P-bifurcation analysis is carried out considering U * mean as the control parameter.Conventionally, P-bifurcation is defined by the qualitative changes in the structure of the stationary joint-pdf of the response variables and their instantaneous time derivatives, as the bifurcation parameter is varied [17,29,30].Figure 9 shows the typical joint pdfs p(α, α) for the pitch response and its instantaneous time derivatives for some representative values of U * mean .It can be clearly seen that at U * mean = 2.0, the joint p(α, α) appears to have a unimodal structure centered about the origin; see Fig. 9(a).Figure 9(b) shows the j-pdf structure for U * mean = 4.5.It is observed that while a peak at the origin still exists, a ring like structure appears to form around the peak at the origin.It is observed that the amplitude of the peak at the origin has decreased from what is observed in Fig. 9(a).It is quite obvious that the strength of the attractor associated with the annular ring is significantly weaker than the one at the origin.As U * mean is increased to 5, it can be observed from Fig. 9(c) that the ring like structure has become more prominent and the peak at the center has become smaller.Finally, at U * mean = 6.0, it can be seen from Fig. 9(d) that the joint pdf resembles a crater like structure with the peak at the origin having completely disappeared.Clearly, one can see that there have been topolog- The underlying physics associated with these joint pdfs can be better appreciated from an inspection of the contour plots of the joint pdfs and shown in Fig. 10.Note that these pdf contours provide a measure of the time spent by the trajectories in a volume element of the state space which is indicative of different stochastic attractors.The contour plots of p(α, α) for U * mean = 2.0 reveal a circular structure about the origin with very small radius (Figure 10(a)) as the trajectories spend most of the time around the zero state.At U * mean = 4.5, the contours shown in Fig. 10(b) reveal that the radius of the attractor at the origin has increased due to the presence of periodic bursts denoting the wider base of the corresponding jpdf.As U * mean is increased to 5, two attractors-one at the origin and another the annular ring are seen to co-exist characterizing the 'off' and the 'on' states respectively in the intermittent phase; see Fig. 10(c).However, as the annular ring like attractor emerges reflecting prominent periodic bursts, the contour level values of the existing attractor at the origin gets reduced.It signifies the fact that the probability of the trajectories spending time in the zero state becomes lower.Finally, Fig. 10(d) shows that at U * mean = 6.0, the attractor at the origin has disappeared and only the annular ring like attractor exists characterizing the random limit-cycle oscillations.

Concluding remarks
A stochastic bifurcation analysis has been carried out for a two degree-of-freedom pitch-plunge airfoil subjected to long time-scale noisy flow fluctuations, with the mean speed of the flow taken to be the bifurcation parameter.A P-bifurcation analysis based on studying the changes in the topology associated with the structure of the joint pdf of the response variables has been presented.The presence of long time-scale fluctuations in the flow result in an 'onoff' type intermittency phenomenon in flow regimes prior to Hopf bifurcation.The system exhibits intermittency at flow regimes where the phase space is characterized by two co-existing attractors.This is evidently reflected in the joint pdf of response variables.The effect of correlation on the noise induced intermittency is currently being investigated by the authors.

Figure 1 .
Figure 1.Schematic of airfoil in pitch and plunge degrees of freedom.

Figure 6 .
Figure 6.Bifurcation diagram of pitch response for uniform flow condition.