# Higher Period Stochastic Bifurcation of Nonlinear Airfoil Fluid-Structure Interaction

Authors: Jeroen A. S. Witteveen; Hester Bijl
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

The higher period stochastic bifurcation of a nonlinear airfoil fluid-structure interaction system is analyzed
using an efficient and robust uncertainty quantification method for unsteady problems. The computationally efficient numerical approach achieves a constant error with a constant number of samples in time. The robustness of the method is assured by the extrema diminishing concept in probability space. The numerical results demonstrate that the system is even more sensitive to randomness at the higher period bifurcation than in the first bifurcation point. In this isolated point in parameter space
the clear hierarchy of increasing importance of the random nonlinearity parameter, initial condition, and natural frequency ratio, respectively, even suddenly reverses. Disregarding seemingly less important random parameters based on a preliminary analysis can, therefore, be an unreliable approach for reducing the number of relevant random input parameters.


---

## Body

## 1. Introduction

It is widely know that the behavior of nonlinear dynamical systems is highly sensitive to small variations. Examples of significant effects of varying initial conditions and model parameters in time-dependent problems can be found in many branches of science and engineering. In turbulence modeling and nonlinear stability theory of transition it is recognized that uncertainty in the initial conditions has a substantial effect on the long-term solution [1–3]. The inherent sensitivity of meteorological and atmospheric models for weather prediction results in a rapid loss of simulation accuracy over time [4, 5]. Stochastic parameters also affect the voltage oscillations in the electric circuit of a nonlinear transistor amplifier [6]. In this paper, the aeronautical application of the effect of randomness on the bifurcation of a nonlinear aeroelastic wing structure is analyzed. Physical uncertainties are encountered in this kind of fluid-structure systems due to varying atmospheric conditions, wear and tear, and production tolerances affecting material properties and the geometry. Compared to the deterministic case the stochastic bifurcation can lead to an earlier onset of unstable flutter behavior, which can cause fatigue damage and structural failure.Fluid-structure interaction systems can be modeled deterministically using detailed finite-element method (FEM) structural discretizations and high-fidelity unsteady computational fluid dynamics (CFD) simulations. This computationally highly intensive approach is usually too expensive for performing many deterministic simulations required in a flutter analysis study. In flutter analysis the structure is, therefore, usually modeled by rigid airfoil mass-spring systems for wing structures and by plate equations for plate-like designs [7, 8]. Structural nonlinearity is then modeled by a cubic nonlinear spring, since structures commenly behave as a cubic stiffness hardening spring [9]. In this framework the flow forces are taken into account in the governing structural equations by source terms prescribed by aerodynamic models. These simplifications are even more frequently used in the stochastic analysis of aeroelastic systems, since each additional random input parameter contributes to the dimensionality of the parameter domain under consideration.The stochastic bifurcation behavior of aeroelastic systems has previously been studied using perturbation techniques, Monte Carlo simulation, Polynomial Chaos formulations, and a range of other numerical and analytical methods. The perturbation approach [10] has been used by Poirion to obtain a first-order approximation of the flutter probability of a bending-torsion structural model, see [11, 12]. A second moment perturbation-based stochastic finite-element method has been applied by Liaw and Yang [13] to determine the effect of uncertainties on panel flutter. The vibration of a hydrofoil in random flow has been considered using a stochastic pertubation approach by Carcaterra et al. [14]. Monte Carlo simulations [15] have, for example, been used by Lindsley et al. [16, 17] to study the periodic response of nonlinear plates under supersonic flow subject to randomness. Poirel and Price [18] have studied random bending-torsion flutter equations with turbulent flow conditions and a linear structural model using also a Monte Carlo-type approach. The stochastic postflutter behavior of limit cycle oscillations has been studied by Beran et al. [19] using Monte Carlo sampling.Other uncertainty quantification methodologies have, for example, been employed in an investigation of nonlinear random oscillations of aeroservoelastic systems by Poirion [20] using random delay modeling of control systems. Choi and Namachchivaya [21] have used nonstandard reduction through stochastic averaging in nonlinear panel flutter under supersonic flow subject to random fluctuations in the turbulent boundary layer to find response density functions. De Rosa and Franco [22] have predicted the stochastic response of a plate subject to a turbulent boundary layer using numerical and analytical approaches. Frequency domain methods have been considered for solving linear stochastic operator equations by Sarkar and Ghanem [23].Polynomial Chaos methods [24–28] are, in general, computationally efficient alternatives for the detailed and quantitative probabilistic modeling of physical uncertainties by Monte Carlo simulation. However, in dynamic simulations the Polynomial Chaos method usually requires a fast increasing expansion order to maintain a constant accuracy in time. This leads to a fast increasing sample size in the more practical nonintrusive Polynomial Chaos formulations [29–33], which are based on the polynomial interpolation of an in general small number of samples. Resolving the asymptotic stochastic effect in a postflutter analysis using nonintrusive Polynomial Chaos can, however, lead to a very high number of required deterministic simulations. This effect is especially profound in problems with an oscillatory solution in which the frequency of the response is affected by the random parameters [3, 34]. Pettit and Beran [35] have demonstrated that the Polynomial Chaos expansion is subject to energy loss in representing the periodic response of a bending-torsion flutter model for long integration times. They have also found that the wavelet based Wiener-Haar expansion of Le Maître et al. [36] loses its accuracy less rapidly. Also other multielement Polynomial Chaos formulations have been shown to postpone the resolution problems [37]. Millman et al. [38] have proposed a Fourier-Chaos expansion for oscillatory responses in application to a bending-torsion flutter problem subjected to Gaussian distributions. The effectivity of intrusive and nonintrusive Polynomial Chaos methods has been compared for a single-degree-of-freedom pitching airfoil stall flutter system in [39, 40].Another special Polynomial Chaos formulation for oscillatory problems was recently also developed to maintain a constant accuracy in time with a constant polynomial order [41, 42]. The nonintrusive approach is based on normalizing the oscillatory samples in terms of their phase. The uncertainty quantification interpolation of the samples is then performed at constant phase, which eliminates the effect of frequency differences on the increase of the required sample size [43]. The method is proven to result in a bounded error as function of the phase with a constant number of samples for periodic responses and under certain conditions also in a bounded error in time [44]. The formulation was also extended to multifrequency responses of continuous structures by using a wavelet decomposition preprocessing step [45]. Application of the method to an elastically mounted airfoil showed that this fluid-structure interaction system is sensitive to small variations at the bifurcation from a stable solution to a period-1 limit cycle oscillation [43]. A period-1 motion refers to an oscillation that repeats itself after a 2π-orbit around a fixed point in phase space.In this paper the latter uncertainty quantification approach is employed to analyze the stochastic higher period bifurcation of an aeroelastic airfoil with nonlinear structural stiffness. It is demonstrated that the fluid-structure interaction system is even more sensitive to randomness at the higher period bifurcation than in the previously considered first bifurcation point. The resulting general mathematical formulation of the uncertainty quantification problem is given in Section2. The efficient and robust uncertainty quantification method for unsteady problems based on extrema diminishing interpolation of oscillatory samples at constant phase used to resolve the stochastic bifurcation behavior numerically is introduced in Section 3. As is common in stochastic flutter analysis, the aeroelastic system is modeled by a two-dimensional rigid airfoil with two degrees of freedom in pitch and plunge, and cubic nonlinear spring stiffness. The aerodynamic loads are computed using an aerodynamic model as described in Section 4. Randomness is introduced in terms of three random parameters in the system and its initial conditions. The effect of uncertainty in the ratio of natural pitch and plunge frequencies is resolved in Section 5. A random nonlinear spring parameter is considered in Section 6. The effect of randomness in the initial condition of the pitch angle is investigated in Section 7. The main findings are summarized in Section 8.


## 2. Mathematical Formulation of the Uncertainty Quantification Problem

Consider a dynamical system subject tona uncorrelated second-order random input parameters a(ω)={a1(ω),…,ana(ω)}∈A with parameter space A∈ℝna, which governs an oscillatory response u(x,t,a)(2.1)ℒ(x,t,a;u(x,t,a))=S(x,t,a),
with operator ℒ and source term S defined on domain D×T×A, and appropriate initial and boundary conditions. The spatial and temporal dimensions are defined as x∈D and t∈T, respectively, with D⊂ℝd, d={1,2,3}, and T=[0,tmax]. A realization of the set of outcomes Ω of the probability space (Ω,ℱ,P) is denoted by ω∈Ω=[0,1]na, with ℱ⊂2Ω the σ-algebra of events and P a probability measure.Here we consider a nonintrusive uncertainty quantification methodl which constructs a weighted approximation w(x,t,a) of response surface u(x,t,a) based on ns deterministic solutions vk(x,t)≡u(x,t,ak) of (2.1) for different parameter values ak≡a(ωk) for k=1,…,ns. The samples vk(x,t) can be obtained by solving the deterministic problem(2.2)ℒ(x,t,ak;vk(x,t))=𝒮(x,t,ak),
for k=1,…,ns, using standard spatial discretization methods and time marching schemes. A nonintrusive uncertainty quantification method l is then a combination of a sampling method g and an interpolation method h. Sampling method g defines the ns sampling points {ak}k=1ns and returns the deterministic samples v(x,t)={v1(x,t),…,vns(x,t)}. Interpolation method h constructs an interpolation surface w(x,t,a) through the ns samples v(x,t) as an approximation of u(x,t,a). We are eventually interested in an approximation of the probability distribution and statistical moments μui(x,t) of the output u(x,t,a), which can be obtained by sorting and weighted integration of w(x,t,a):(2.3)μui(x,t)≈μwi(x,t)=∫A‍w(x,t,a)ifa(a)da.
This information can be used for reducing design safety factors and robust design optimization, in contrast to reliability analysis in which the probability of failure is determined [46].


## 3. An Efficient Uncertainty Quantification Method for Unsteady Problems

The efficient uncertainty quantification formulation for oscillatory responses based on interpolation of scaled samples at constant phase is developed in Section3.2. The robust extrema diminishing uncertainty quantification method based on Newton-Cotes quadrature in simplex elements employed in the unsteady approach is first presented in the next section.


### 3.1. Robust Extrema Diminishing Uncertainty Quantification

A multielement uncertainty quantification methodl evaluates integral (2.3) by dividing parameter space A into ne non-overlapping simplex elements Aj⊂A:(3.1)μwi(x,t)=∑j=1ne∫Aj‍w(x,t,a)ifa(a)da.
Here we consider a multielement Polynomial Chaos method based on Newton-Cotes quadrature points and simplex elements [47]. A piecewise polynomial approximation w(x,t,a) is then constructed based on ns deterministic solutions vj,k(x,t)=u(x,t,aj,k) for the values of the random parameters aj,k that correspond to the ñs Newton-Cotes quadrature points of degree d in the elements Aj:(3.2)μwi(x,t)=∑j=1ne‍∑k=1ñs‍cj,kvj,k(x,t)i,
where cj,k is the weighted integral of the Lagrange interpolation polynomial Lj,k(a) through Newton-Cotes quadrature point k in element Aj:(3.3)cj,k=∫Aj‍Lj,k(a)fa(a)da,
for j=1,…,ne and k=1,…,ñs. Here, second-degree Newton-Cotes quadrature with d=2 is considered in combination with adaptive mesh refinement in probability space, since low-order approximations are more effective for approximating response response surfaces with singularities. The initial discretization of parameter space A for the adaptive scheme consists of the minimum of neini=na! simplex elements and nsini=3na samples, see Figure 1. The example of Figure 1 for two random input parameters can geometrically be extended to higher dimensional probability spaces. The elements Aj are adaptively refined using a refinement measure ρj based on the largest absolute eigenvalue of the Hessian Hj, as a measure of the curvature of the response surface approximation in the elements, weighted by the probability fj contained by the elementsDiscretization of two-dimensional parameter spaceA using 2-simplex elements and second-degree Newton-Cotes quadrature points given by the dots.
(a)(b)(c)(3.4)fj=∫Aj‍fa(a)da,
with ∑j=1ne‍fj=1. The stochastic grid refinement is terminated when convergence measure δne is smaller than a threshold value δne<δ̅ where(3.5)δne=max(∥μw⌊ne/2⌋(x,t)-μwne(x,t)∥∞∥μwne(x,t)∥∞,∥σw⌊ne/2⌋(x,t)-σwne(x,t)∥∞∥σwne(x,t)∥∞),
with μw(x,t) and σw(x,t) the mean and standard deviation of w(x,t,ω), or when a maximum number of samples n̅s is reached. Convergence measure δne can be extended to include also higher statistical moments of the output.In elements where the quadratic second-degree interpolation results in an extremum other than in a quadrature point, the element is subdivided intoñe=2na subelements with a linear first-degree Newton-Cotes approximation of the response without performing additional deterministic solves. It is proven in [44] that the resulting approach satisfies the extrema diminishing (ED) robustness concept in probability space(3.6)minA(w(a))≥minA(u(a))⋀maxA(w(a))≤maxA(u(a)),∀u(a),
where the arguments x and t are omitted for simplicity of the notation. The ED property leads to the advantage that no non-zero probabilities of unphysical realizations can be predicted due to overshoots or undershoots at discontinuities in the response. Due to the location of the Newton-Cotes quadrature points the deterministic samples are also reused in successive refinements and the samples are used in approximating the response in multiple elements.


### 3.2. Efficient Uncertainty Quantification Interpolation at Constant Phase

Polynomial Chaos methods usually require a fast increasing number of samples with time to maintain a constant accuracy. Performing the uncertainty quantification interpolation of oscillatory samples at constant phase instead of at constant time results, however, in a constant accuracy with a constant number of samples. Assume, therefore, that solving (2.2) for realizations of the random parameters ak results in oscillatory samples vk(t)=u(ak), of which the phase vϕk(t)=ϕ(t,ak) is a well-defined monotonically increasing function of time t for k=1,…,ns.In order to interpolate the samplesv(t)={v1(t),…,vna(t)} at constant phase, first, their phase as function of time vϕ(t)={vϕ1(t),…,vϕna(t)} is extracted from the deterministic solves v(t). Second, the time series for the phase vϕ(t) are used to transform the samples v(t) into functions of their phase v̂(vϕ(t)) according to(3.7)v̂k(vϕk(t))=vk(t),
for k=1,…,ns, see Figure 2. Third, the sampled phases vϕ(t) are interpolated to the function wϕ(t,a):Oscillatory samples as function of time and phase.
(a)(b)(3.8)wϕ(t,a)=h(vϕ(t)),
as approximation of ϕ(t,a). Finally, the transformed samples v̂(vϕ(t)) are interpolated at a constant phase φ∈wϕ(t,a) to(3.9)ŵ(φ,a)=h(v̂(φ)).
Repeating the latter interpolation for all phases φ∈wϕ(t,a) results in the function ŵ(wϕ(t,a),a). The interpolation ŵ(wϕ(t,a),a) is then transformed back to an approximation in the time domain w(t,a) as follows:(3.10)w(t,a)=ŵ(wϕ(t,a),a).
The resulting function w(t,a) is an approximation of the unknown response surface u(t,a) as function of time t and the random parameters a(ω). The actual sampling g and interpolation h is performed using the extrema diminishing uncertainty quantification method l based on Newton-Cotes quadrature in simplex elements described in the previous section.This uncertainty quantification formulation for oscillatory responses is proven to achieve a bounded errorε̂(φ,a)=|ŵ(φ,a)-û(φ,a)| as function of phase φ for periodic responses according to(3.11)ε̂(φ,a)<δ,∀φ∈ℝ,a∈A,
where δ is defined by(3.12)ε̂(φ,a)<δ,∀φ∈[0,1],a∈A.
The error ε(t,a)=|w(t,a)-u(t,a)| is also bounded in time under certain conditions, see [44].The phasesvϕ(t) are extracted from the samples based on the local extrema of the time series v(t). A trial and error procedure identifies a cycle of oscillation based on two or more successive local maxima. The selected cycle is accepted if the maximal error of its extrapolation in time with respect to the actual sample is smaller than a threshold value ε̅k for at least one additional cycle length. The functions for the phases vϕ(t) in the whole time domain T are constructed by identifying all successive cycles of v(t) and linear extrapolation to t=0 and t=tmax before and after the first and last complete cycle, respectively. The phase is normalized to zero at the start of the first cycle and a user-defined parameter determines whether the sample is assumed to attain a local extremum at t=0. The interpolation at constant phase is restricted to the time domain that corresponds to the range of phases that is reached by all samples in each of the elements. If the phase vϕk(t) cannot be extracted from one of the samples vk(t) for k=1,…,ns, then uncertainty quantification interpolation h is directly applied to the time-dependent samples v(t).


## 3.1. Robust Extrema Diminishing Uncertainty Quantification

A multielement uncertainty quantification methodl evaluates integral (2.3) by dividing parameter space A into ne non-overlapping simplex elements Aj⊂A:(3.1)μwi(x,t)=∑j=1ne∫Aj‍w(x,t,a)ifa(a)da.
Here we consider a multielement Polynomial Chaos method based on Newton-Cotes quadrature points and simplex elements [47]. A piecewise polynomial approximation w(x,t,a) is then constructed based on ns deterministic solutions vj,k(x,t)=u(x,t,aj,k) for the values of the random parameters aj,k that correspond to the ñs Newton-Cotes quadrature points of degree d in the elements Aj:(3.2)μwi(x,t)=∑j=1ne‍∑k=1ñs‍cj,kvj,k(x,t)i,
where cj,k is the weighted integral of the Lagrange interpolation polynomial Lj,k(a) through Newton-Cotes quadrature point k in element Aj:(3.3)cj,k=∫Aj‍Lj,k(a)fa(a)da,
for j=1,…,ne and k=1,…,ñs. Here, second-degree Newton-Cotes quadrature with d=2 is considered in combination with adaptive mesh refinement in probability space, since low-order approximations are more effective for approximating response response surfaces with singularities. The initial discretization of parameter space A for the adaptive scheme consists of the minimum of neini=na! simplex elements and nsini=3na samples, see Figure 1. The example of Figure 1 for two random input parameters can geometrically be extended to higher dimensional probability spaces. The elements Aj are adaptively refined using a refinement measure ρj based on the largest absolute eigenvalue of the Hessian Hj, as a measure of the curvature of the response surface approximation in the elements, weighted by the probability fj contained by the elementsDiscretization of two-dimensional parameter spaceA using 2-simplex elements and second-degree Newton-Cotes quadrature points given by the dots.
(a)(b)(c)(3.4)fj=∫Aj‍fa(a)da,
with ∑j=1ne‍fj=1. The stochastic grid refinement is terminated when convergence measure δne is smaller than a threshold value δne<δ̅ where(3.5)δne=max(∥μw⌊ne/2⌋(x,t)-μwne(x,t)∥∞∥μwne(x,t)∥∞,∥σw⌊ne/2⌋(x,t)-σwne(x,t)∥∞∥σwne(x,t)∥∞),
with μw(x,t) and σw(x,t) the mean and standard deviation of w(x,t,ω), or when a maximum number of samples n̅s is reached. Convergence measure δne can be extended to include also higher statistical moments of the output.In elements where the quadratic second-degree interpolation results in an extremum other than in a quadrature point, the element is subdivided intoñe=2na subelements with a linear first-degree Newton-Cotes approximation of the response without performing additional deterministic solves. It is proven in [44] that the resulting approach satisfies the extrema diminishing (ED) robustness concept in probability space(3.6)minA(w(a))≥minA(u(a))⋀maxA(w(a))≤maxA(u(a)),∀u(a),
where the arguments x and t are omitted for simplicity of the notation. The ED property leads to the advantage that no non-zero probabilities of unphysical realizations can be predicted due to overshoots or undershoots at discontinuities in the response. Due to the location of the Newton-Cotes quadrature points the deterministic samples are also reused in successive refinements and the samples are used in approximating the response in multiple elements.


## 3.2. Efficient Uncertainty Quantification Interpolation at Constant Phase

Polynomial Chaos methods usually require a fast increasing number of samples with time to maintain a constant accuracy. Performing the uncertainty quantification interpolation of oscillatory samples at constant phase instead of at constant time results, however, in a constant accuracy with a constant number of samples. Assume, therefore, that solving (2.2) for realizations of the random parameters ak results in oscillatory samples vk(t)=u(ak), of which the phase vϕk(t)=ϕ(t,ak) is a well-defined monotonically increasing function of time t for k=1,…,ns.In order to interpolate the samplesv(t)={v1(t),…,vna(t)} at constant phase, first, their phase as function of time vϕ(t)={vϕ1(t),…,vϕna(t)} is extracted from the deterministic solves v(t). Second, the time series for the phase vϕ(t) are used to transform the samples v(t) into functions of their phase v̂(vϕ(t)) according to(3.7)v̂k(vϕk(t))=vk(t),
for k=1,…,ns, see Figure 2. Third, the sampled phases vϕ(t) are interpolated to the function wϕ(t,a):Oscillatory samples as function of time and phase.
(a)(b)(3.8)wϕ(t,a)=h(vϕ(t)),
as approximation of ϕ(t,a). Finally, the transformed samples v̂(vϕ(t)) are interpolated at a constant phase φ∈wϕ(t,a) to(3.9)ŵ(φ,a)=h(v̂(φ)).
Repeating the latter interpolation for all phases φ∈wϕ(t,a) results in the function ŵ(wϕ(t,a),a). The interpolation ŵ(wϕ(t,a),a) is then transformed back to an approximation in the time domain w(t,a) as follows:(3.10)w(t,a)=ŵ(wϕ(t,a),a).
The resulting function w(t,a) is an approximation of the unknown response surface u(t,a) as function of time t and the random parameters a(ω). The actual sampling g and interpolation h is performed using the extrema diminishing uncertainty quantification method l based on Newton-Cotes quadrature in simplex elements described in the previous section.This uncertainty quantification formulation for oscillatory responses is proven to achieve a bounded errorε̂(φ,a)=|ŵ(φ,a)-û(φ,a)| as function of phase φ for periodic responses according to(3.11)ε̂(φ,a)<δ,∀φ∈ℝ,a∈A,
where δ is defined by(3.12)ε̂(φ,a)<δ,∀φ∈[0,1],a∈A.
The error ε(t,a)=|w(t,a)-u(t,a)| is also bounded in time under certain conditions, see [44].The phasesvϕ(t) are extracted from the samples based on the local extrema of the time series v(t). A trial and error procedure identifies a cycle of oscillation based on two or more successive local maxima. The selected cycle is accepted if the maximal error of its extrapolation in time with respect to the actual sample is smaller than a threshold value ε̅k for at least one additional cycle length. The functions for the phases vϕ(t) in the whole time domain T are constructed by identifying all successive cycles of v(t) and linear extrapolation to t=0 and t=tmax before and after the first and last complete cycle, respectively. The phase is normalized to zero at the start of the first cycle and a user-defined parameter determines whether the sample is assumed to attain a local extremum at t=0. The interpolation at constant phase is restricted to the time domain that corresponds to the range of phases that is reached by all samples in each of the elements. If the phase vϕk(t) cannot be extracted from one of the samples vk(t) for k=1,…,ns, then uncertainty quantification interpolation h is directly applied to the time-dependent samples v(t).


## 4. A Nonlinear Airfoil Fluid-Structure Interaction System

The nonlinear airfoil flutter model used to simulate the airfoil fluid-structure interaction system is given in Section4.1. The deterministic bifurcation behavior is briefly considered in Section 4.2.


### 4.1. A Two-Degree-of-Freedom Airfoil Flutter Model

The two-degree-of-freedom model for the pitch and plunge motion of an airfoil used here, see Figure3, has also been studied deterministically, for example, by Lee et al. [48] and stochastically using Fourier Chaos by Millman et al. [38]. The aeroelastic equations of motion with cubic restoring springs in both pitch and plunge are given in [48] byFigure 3
The two-degree-of-freedom airfoil flutter model.(4.1)ξ''+xαα''+2ζξω̅U*ξ'+(ω̅U*)2(ξ+βξξ3)=-1πμCL(τ),xαrα2ξ''+α''+2ζαU*α'+1U*2(α+βαα3)=2πμrα2CM(τ),
where α(τ) is the pitch angle and ξ(τ)=h/b is the nondimensional version of the plunge displacement h of the elastic axis, with b=c/2 the half-chord, and initial conditions α(0)=α0 and ξ(0)=ξ0. The nonlinear spring constants in plunge and pitch are, respectively, βξ and βα. EEquivalently, the viscous damping coefficients are ζξ and ζα. The ratio of natural frequencies is given by ω̅=ωξ/ωα, where ωξ and ωα are the natural frequencies of the uncoupled plunging and pitching modes, respectively. The mass ratio μ is defined as m/πρb2, with m the airfoil mass, and ρ the air density. The radius of gyration about the elastic axis is rα, where elastic axis is located at a distance ahb from the mid-chord point, and the mass center is located at a distance xαb from the elastic axis. The bifurcation parameter is the ratio of time scales of the structure and the flow defined as U*=U/(bωα), with U the free stream velocity. The primes denote differentiation with respect to nondimensionalized time τ=Ut/b. The expressions for the aerodynamic force and moment coefficients, CL(τ) and CM(τ) are given by Fung [49] as(4.2)CL(τ)=π(ξ''-ahα''+α')+2π{α(0)+ξ'(0)+[12-ah]α'(0)}ϕ(τ)+2π∫0τϕ(τ-σ)[α'(σ)+ξ''(σ)+(12-ah)α''(σ)]dσ,CM(τ)=π(12+ah){α(0)+ξ'(0)+(12-ah)α'(0)}ϕ(τ)+π(12+ah)∫0τϕ(τ-σ){α'(σ)+ξ''(σ)+(12-ah)α''(σ)}dσ+π2ah(ξ''-ahα'')-(12-ah)π2α'-π16α'',
where ϕ(τ) is the Wagner function(4.3)ϕ(τ)=1-ψ1e-ε1τ-ψ2e-ε2τ,
with the constants ψ1=0.165, ψ2=0.335, ε1=0.0455, and ε2=0.3 given by Jones [50]. Based on (4.1) to (4.3), the following set of first-order ordinary differential equations for the motion of the airfoil is derived in [48](4.4)x1'=x2,x2'=c0H-d0Pd0c1-c0d1,x3'=x4,x4'=-c1H+d1Pd0c1-c0d,x5'=x1-ε1x5,x6'=x1-ε2x6,x7'=x3-ε1x7,x8'=x3-ε2x8,
with(4.5)P=c2x4+c3x2+c4x3+c5x33+c6x1+c7x5+c8x6+c9x7+c10x8-f(τ),H=d2x2+d3x1+d4x13+d5x4+d6x3+d7x5+d8x6+d9x7+d10x8-g(τ),
where a vector {xi}i=18 of new variables is defined as(4.6)x1=α,x2=α',x3=ξ,x4=ξ',x5=w1,x6=w2,x7=w3,x8=w4,w1=∫0τe-ϵ1(τ-σ)α(σ)dσ,w2=∫0τe-ϵ2(τ-σ)α(σ)dσ,w3=∫0τe-ϵ1(τ-σ)ξ(σ)dσ,w4=∫0τe-ϵ2(τ-σ)ξ(σ)dσ.
Following [48], the solution is determined numerically until τ=2000 using the explicit fourth order Runge-Kutta method with a time step of Δτ=0.1, which is approximately 1/256 of the smallest period. The other parameter values are chosen to be μ=100, ah=-0.5, xα=0.25, rα=0.5, ξ0=0, βξ=0, and ζα=ζξ=0 as in [48].The system parameters that are assumed to be uncertain areω̅, βα, and α0. The randomness of these three parameters is described by a symmetric unimodal beta distribution with β1=β2=2 to limit their parameter range to a finite domain with vanishing probability at the interval boundaries. The ratio of natural frequencies ω̅ has a mean of μω̅=0.2 and an interval of μω̅∈[0.15;0.25]. The input probability density function for ω̅ is shown as an example in Figure 4. For a hard spring model with βα>0 the system exhibits a stable limit cycle oscillation beyond the first bifurcation point [51]. The mean of the nonlinear stiffness parameter βα(ω) is chosen to be μβα=100 to limit the pitch angle α to the domain in which the aerodynamic model is valid, and the interval is set to βα(ω)∈[90,110]. The initial condition α0 has an interval of α0∈[9,11] degrees around mean μα0=10°. The resulting coefficients of variation for the random parameters are cvω̅=11.2%, cvβα=4.48%, and cvα0=4.48%. The effect of the random parameters ω̅, βα, and α0 is analyzed in Sections 5 to 7. First the deterministic bifurcation behavior is explored in the next section.Figure 4
Input probability density function for the ratio of natural frequenciesω̅.


### 4.2. Deterministic Bifurcation Behavior

The deterministic bifurcation plot for the aeroelastic system given by (4.1) to (4.3) is shown in Figure 5 as function of bifurcation parameter U*∈[5,15] in terms of the angles of attack α for which α'=0 in the asymptotic range. In what follows the first deterministic bifurcation point of U*=6.25 the system response is a decaying oscillation to α=0. At U*=6.25 the system exhibits a supercritical Hopf bifurcation to a stable period-1 limit cycle oscillation with an increasing amplitude for increasing U*. At the second bifurcation point U*=13.42 the response shows an abrupt bifurcation to a higher period limit cycle oscillation. The oscillation amplitude continues to increase beyond the second bifurcation point.Figure 5
Deterministic bifurcation plot of the nonlinear airfoil flutter system.Three typical time histories of pitchα and plunge ξ in the three different regimes are given in Figure 6 as function of nondimensional time τ for U*={5,10,15}. At U*=5 both α and ξ are decaying oscillations to the stable fixed point (α,ξ)=(0,0). A period-1 oscillation can be identified for α and ξ at U*=10. For U*=15 the angle of attack α exhibits a higher period oscillation with a higher amplitude, while the plunge deflection ξ maintains a period-1 oscillation. The mean, standard deviation, and probability distribution of the more interesting pitch degree of freedom α is, therefore, considered in the following stochastic flutter analysis. The plunge ξ is used to extract the phase of the oscillation.Time histories of pitchα(τ) and plunge ξ(τ) in the three regimes of the deterministic airfoil flutter system.
(a)(b)(c)


## 4.1. A Two-Degree-of-Freedom Airfoil Flutter Model

The two-degree-of-freedom model for the pitch and plunge motion of an airfoil used here, see Figure3, has also been studied deterministically, for example, by Lee et al. [48] and stochastically using Fourier Chaos by Millman et al. [38]. The aeroelastic equations of motion with cubic restoring springs in both pitch and plunge are given in [48] byFigure 3
The two-degree-of-freedom airfoil flutter model.(4.1)ξ''+xαα''+2ζξω̅U*ξ'+(ω̅U*)2(ξ+βξξ3)=-1πμCL(τ),xαrα2ξ''+α''+2ζαU*α'+1U*2(α+βαα3)=2πμrα2CM(τ),
where α(τ) is the pitch angle and ξ(τ)=h/b is the nondimensional version of the plunge displacement h of the elastic axis, with b=c/2 the half-chord, and initial conditions α(0)=α0 and ξ(0)=ξ0. The nonlinear spring constants in plunge and pitch are, respectively, βξ and βα. EEquivalently, the viscous damping coefficients are ζξ and ζα. The ratio of natural frequencies is given by ω̅=ωξ/ωα, where ωξ and ωα are the natural frequencies of the uncoupled plunging and pitching modes, respectively. The mass ratio μ is defined as m/πρb2, with m the airfoil mass, and ρ the air density. The radius of gyration about the elastic axis is rα, where elastic axis is located at a distance ahb from the mid-chord point, and the mass center is located at a distance xαb from the elastic axis. The bifurcation parameter is the ratio of time scales of the structure and the flow defined as U*=U/(bωα), with U the free stream velocity. The primes denote differentiation with respect to nondimensionalized time τ=Ut/b. The expressions for the aerodynamic force and moment coefficients, CL(τ) and CM(τ) are given by Fung [49] as(4.2)CL(τ)=π(ξ''-ahα''+α')+2π{α(0)+ξ'(0)+[12-ah]α'(0)}ϕ(τ)+2π∫0τϕ(τ-σ)[α'(σ)+ξ''(σ)+(12-ah)α''(σ)]dσ,CM(τ)=π(12+ah){α(0)+ξ'(0)+(12-ah)α'(0)}ϕ(τ)+π(12+ah)∫0τϕ(τ-σ){α'(σ)+ξ''(σ)+(12-ah)α''(σ)}dσ+π2ah(ξ''-ahα'')-(12-ah)π2α'-π16α'',
where ϕ(τ) is the Wagner function(4.3)ϕ(τ)=1-ψ1e-ε1τ-ψ2e-ε2τ,
with the constants ψ1=0.165, ψ2=0.335, ε1=0.0455, and ε2=0.3 given by Jones [50]. Based on (4.1) to (4.3), the following set of first-order ordinary differential equations for the motion of the airfoil is derived in [48](4.4)x1'=x2,x2'=c0H-d0Pd0c1-c0d1,x3'=x4,x4'=-c1H+d1Pd0c1-c0d,x5'=x1-ε1x5,x6'=x1-ε2x6,x7'=x3-ε1x7,x8'=x3-ε2x8,
with(4.5)P=c2x4+c3x2+c4x3+c5x33+c6x1+c7x5+c8x6+c9x7+c10x8-f(τ),H=d2x2+d3x1+d4x13+d5x4+d6x3+d7x5+d8x6+d9x7+d10x8-g(τ),
where a vector {xi}i=18 of new variables is defined as(4.6)x1=α,x2=α',x3=ξ,x4=ξ',x5=w1,x6=w2,x7=w3,x8=w4,w1=∫0τe-ϵ1(τ-σ)α(σ)dσ,w2=∫0τe-ϵ2(τ-σ)α(σ)dσ,w3=∫0τe-ϵ1(τ-σ)ξ(σ)dσ,w4=∫0τe-ϵ2(τ-σ)ξ(σ)dσ.
Following [48], the solution is determined numerically until τ=2000 using the explicit fourth order Runge-Kutta method with a time step of Δτ=0.1, which is approximately 1/256 of the smallest period. The other parameter values are chosen to be μ=100, ah=-0.5, xα=0.25, rα=0.5, ξ0=0, βξ=0, and ζα=ζξ=0 as in [48].The system parameters that are assumed to be uncertain areω̅, βα, and α0. The randomness of these three parameters is described by a symmetric unimodal beta distribution with β1=β2=2 to limit their parameter range to a finite domain with vanishing probability at the interval boundaries. The ratio of natural frequencies ω̅ has a mean of μω̅=0.2 and an interval of μω̅∈[0.15;0.25]. The input probability density function for ω̅ is shown as an example in Figure 4. For a hard spring model with βα>0 the system exhibits a stable limit cycle oscillation beyond the first bifurcation point [51]. The mean of the nonlinear stiffness parameter βα(ω) is chosen to be μβα=100 to limit the pitch angle α to the domain in which the aerodynamic model is valid, and the interval is set to βα(ω)∈[90,110]. The initial condition α0 has an interval of α0∈[9,11] degrees around mean μα0=10°. The resulting coefficients of variation for the random parameters are cvω̅=11.2%, cvβα=4.48%, and cvα0=4.48%. The effect of the random parameters ω̅, βα, and α0 is analyzed in Sections 5 to 7. First the deterministic bifurcation behavior is explored in the next section.Figure 4
Input probability density function for the ratio of natural frequenciesω̅.


## 4.2. Deterministic Bifurcation Behavior

The deterministic bifurcation plot for the aeroelastic system given by (4.1) to (4.3) is shown in Figure 5 as function of bifurcation parameter U*∈[5,15] in terms of the angles of attack α for which α'=0 in the asymptotic range. In what follows the first deterministic bifurcation point of U*=6.25 the system response is a decaying oscillation to α=0. At U*=6.25 the system exhibits a supercritical Hopf bifurcation to a stable period-1 limit cycle oscillation with an increasing amplitude for increasing U*. At the second bifurcation point U*=13.42 the response shows an abrupt bifurcation to a higher period limit cycle oscillation. The oscillation amplitude continues to increase beyond the second bifurcation point.Figure 5
Deterministic bifurcation plot of the nonlinear airfoil flutter system.Three typical time histories of pitchα and plunge ξ in the three different regimes are given in Figure 6 as function of nondimensional time τ for U*={5,10,15}. At U*=5 both α and ξ are decaying oscillations to the stable fixed point (α,ξ)=(0,0). A period-1 oscillation can be identified for α and ξ at U*=10. For U*=15 the angle of attack α exhibits a higher period oscillation with a higher amplitude, while the plunge deflection ξ maintains a period-1 oscillation. The mean, standard deviation, and probability distribution of the more interesting pitch degree of freedom α is, therefore, considered in the following stochastic flutter analysis. The plunge ξ is used to extract the phase of the oscillation.Time histories of pitchα(τ) and plunge ξ(τ) in the three regimes of the deterministic airfoil flutter system.
(a)(b)(c)


## 5. Random Natural Frequency Ratioω̅(ω)

First the effect of randomness in the ratio of natural frequenciesω̅ is resolved. The results are presented in terms of the time histories of the mean μα(τ) and standard deviation σα(τ) of the pitch angle α in Figure 7. The bifurcation of the system is illustrated in Figure 8 by the response surface of α as function of the random parameter ω̅ at τ=2000. In Figure 9 the P-bifurcation behavior of the probability density function (PDF) of α also at τ=2000 is given. In all three figures the bifurcation parameter values U*={6.25;10;13.42;15} are considered. This corresponds to the first (U*=6.25) and second (U=13.42) deterministic bifurcation point, and the period-1 (U*=10) and higher period (U*=15) regime. The case of U*=5 also considered in the previous section is not shown here, since the system response in the prebifurcation domain is equal to the trivial solution. The required number of sampling points ns in the stochastic simulations for the different values of U* is established after performing an convergence study which is summarized as an example in Tables 1–4. The results are compared to converged Monte Carlo reference solutions based on ns=103 samples.Table 1
Relation between convergence measureδne and L∞ error εL∞ for the mean μα and standard deviation σα at U*=6.25.
nensMeanμαStandard deviationσαconv.δneerrorεL∞conv.δneerrorεL∞13—1.515·10-1—1.533·100251.151·10-14.115·10-27.117·10-14.293·10-1493.334·10-28.088·10-33.190·10-11.097·10-18176.332·10-31.768·10-36.915·10-25.138·10-2Table 2
Relation between convergence measureδne and L∞ error εL∞ for the mean μα and standard deviation σα at U*=10.
nensMeanμαStandard deviationσαconv.δneerrorεL∞conv.δneerrorεL∞13—1.677·10-3—1.743·10-3251.814·10-32.141·10-41.791·10-34.635·10-4492.533·10-41.092·10-44.536·10-42.335·10-48171.318·10-41.237·10-42.291·10-41.983·10-4Table 3
Relation between convergence measureδne and L∞ error εL∞ for the mean μα and standard deviation σα at U*=13.42.
nensMeanμαStandard deviationσαconv.δneerrorεL∞conv.δneerrorεL∞13—3.176·10-1—3.224·10-1252.666·10-11.905·10-13.389·10-12.117·10-1491.679·10-11.024·10-11.685·10-11.629·10-18175.335·10-21.136·10-17.391·10-21.272·10-116331.103·10-16.984·10-31.253·10-11.060·10-2Table 4
Relation between convergence measureδne and L∞ error εL∞ for the mean μα and standard deviation σα at U*=15.
nensMeanμαStandard deviationσαconv.δneerrorεL∞conv.δneerrorεL∞13—7.315·10-3—7.503·10-3256.319·10-39.960·10-47.444·10-31.619·10-3499.979·10-42.598·10-41.615·10-33.905·10-48172.636·10-41.675·10-43.608·10-43.812·10-4Time histories of the meanμα(τ) and standard deviation σα(τ) of the pitch angle α due to the random frequency ratio ω̅.
(a)(b)(c)(d)Response surface of the pitch angleα at τ=2000 as function of the random frequency ratio ω̅.
(a)(b)(c)(d)Bifurcation of the probability density function of the pitch angleα due to the random frequency ratio ω̅.
(a)(b)(c)(d)ForU*=6.25 the mean μα of the pitch angle shows a decaying oscillation to zero and the standard deviation approaches the steady asymptotic value of σα=0.423 after an initial increase from the deterministic initial condition in Figure 7(a). The decaying mean is caused by a combination of decaying and periodic realizations as can be concluded from the response surface of Figure 8(a). The non-zero asymptotic value of the standard deviation also indicates that due to the randomness in ω̅ the system is already stochastically bifurcated in the deterministic bifurcation point U*=6.25. The onset of the stochastic bifurcation occurs, therefore, at a lower value of the bifurcation parameter than in the deterministic case. As a consequence a deterministic flutter analysis predicts a later start of unstable behavior by neglecting the variability in system parameters, which can lead to disastrous effects by defining the flight envelope based on a too optimistic deterministic flutter boundary.In the period-1 regime at U*=10 the mean μα exhibits a decaying oscillation due the fully periodic response. The resulting frequency differences lead to increasing phase differences in time and increasingly to realizations of opposite sign, which cancel each other resulting in a decaying mean pitch. The standard deviation reaches a significantly higher steady asymptotic value of σα=4.8 due to the increased amplitudes of the limit cycle oscillation at higher values of U*. The effect of ω̅ on the frequency of the response can be derived from the oscillatory response surface of Figure 8(b). The deterministic oscillation period shape of α shown in Figure 6(b) can also be recognized in the shape of the response surface.At the second deterministic bifurcation pointU*=13.42 the mean μα and standard deviation σα show an irregular behavior with only a slowly decaying mean and a large asymptotic standard deviation of approximately σα=7. This is a result of the discontinuity in the response of α in Figure 8(c) caused by the deterministic bifurcation present at μω̅=0.2. On the left and the right of the discontinuity at ω̅=0.2 the higher period and period-1 shape function can be recognized in the response, respectively, which suggests a subcritical Hopf bifurcation as function of ω̅. For U*=15 in the higher period regime μα and σα give again a decaying oscillation and a steady asymptotic value of σα=8.8, respectively. The time histories of μα and σα are initially more complex than in the period-1 regime of U*=10 due to the higher period behavior of the realizations. In the response surface of Figure 8(d) the deterministic higher period shape of α shown Figure 6(c) can again be identified.The required number of samplesns used in the stochastic simulations depends significantly on the value of bifurcation parameter U*. In Tables 1–4 the convergence δne and the error εL∞ with respect to the Monte Carlo reference solutions μαMC(τ) and σαMC(τ) are given. The convergence measure δne used for μα and σα separately is defined by (3.5) and the L∞ errorεL∞ is defined for μα as(5.1)εL∞=∥μα(τ)-μαMC(τ)∥∞∥μαMC(τ)∥∞
and equivalently for σα. The method is highly efficient in the periodic regimes, in which ns=3 samples is already sufficient to match the Monte Carlo results based on ns=103 samples. This holds even for the oscillatory response surface in the higher period case of U*=15. At the deterministic bifurcation points the adaptive method robustly captures the singularity in the response surface by automatically refining near the bifurcation in probability space.The resulting bifurcation behavior of the PDF ofα at τ=2000 is shown in Figure 9. At U*=6.25 the PDF is already bifurcated from a delta function in the stochastic pre-bifurcation domain to a unimodal PDF with the highest probability at α=0°. The PDF develops into a multimodal distribution with peaks at α=±8 due to the oscillatory behavior of the response at U*=10. The multimodal PDF evolves further into a distribution with 6 peaks at approximately α={±5,±11,±16} due to the higher period motion at U*=15. At the second deterministic bifurcation point U*=13.42 the PDF is in an intermediate state between the approximately symmetric multimodal distributions of U*=10 and U*=15.The stochastic behavior of the system is also shown in Figure10 for the three random parameters in terms of the bifurcation of the maximum standard deviation σαmax in the asymptotic range defined asBifurcation of the maximum of the standard deviationσαmax for the three random parameters.
(a)(b)(c)(5.2)σαmax=maxτ∈[1500,2000](σα(τ)).
It can be seen in Figure 10(a) that the first bifurcation of the maximum standard deviation σαmax starts at an earlier location than the first deterministic bifurcation. In the period-1 regime in between the two deterministic bifurcations σαmax gradually increases due to the increasing limit cycle oscillation amplitude in combination with the random frequency. At the second deterministic bifurcation point the standard deviation reaches a local maximum of σαmax=8.0 and it continues to increase at a higher rate beyond U*=13.42.


## 6. Random Nonlinearity Parameterβα(ω)

Next the effect of a random nonlinearity parameter for the pitch degree of freedomβα on the stochastic behavior of the system is considered. The results for the mean and standard deviation, the response surface, and the PDF are shown in Figures 11–13. The required number of samples ns in the simulations for random βα is again determined based on convergence studies.Time histories of the meanμα(τ) and standard deviation σα(τ) of the pitch angle α due to the random nonlinearity parameter βα.
(a)(b)(c)(d)Response surface of the pitch angleα at τ=2000 as function of the random nonlinearity parameter βα.
(a)(b)(c)(d)Bifurcation of the probability density function of the pitch angleα due to the random nonlinearity parameter βα.
(a)(b)(c)(d)ForU*=6.25, U*=10, and U*=15 the random parameter βα has a qualitatively different effect on the system than ω̅. Both the mean μα and standard deviation σα decay in this case to zero for U*=6.25, which suggests that randomness in βα does not lead to an earlier bifurcation. For both U*=10 and U*=15 the mean shows an oscillatory behavior, which closely resembles the deterministic time histories of Figures 6(b) and 6(c). The standard deviation has for these two cases a low constant value of approximately σα=0.7 and σα=0.3, respectively. The response surfaces of Figures 12(b) and 12(d) are also nonoscillatory, which indicates that βα has little effect on the oscillation frequency. It can be concluded that randomness in βα has for these values of U* a small effect of the system behavior. This can be understood from the fact that the nonlinearity parameter has only a significant effect on the limit cycle oscillation amplitude. It can, therefore, be expected that a random βα has a small effect in the prebifurcation domain, and that the effect in the periodic regimes is constant in time.However, the effect of randomβα is significant in the second deterministic bifurcation point U*=13.42. The stochastic system shows at U*=13.42 an irregular behavior with a non-decaying mean μα and a large standard deviation oscillating around approximately σα=7, which is comparable to the results for random ω̅. The sudden large effect of βα is caused by the discontinuity in the response at μβα=100. So, even the randomness in parameter βα, which has in general a small effect on the response, becomes important at the conditions of the second deterministic bifurcation point U*=13.42. The adaptive method resolves also this discontinuous response accurately and the other response surface approximations require again only ns=3 deterministic simulations to match the Monte Carlo results.The PDF in Figure13 is also significantly distorted from the unimodal input distribution at U*=13.42 only. For U*=6.25 the histogram shows a delta function PDF, which indicates that the stochastic bifurcation for random βα has not yet started in the first deterministic bifurcation point. This observation is confirmed by the bifurcation of the maximum standard deviation σαmax in Figure 10(b). The stochastic bifurcation of σαmax coincides with the location of the first deterministic bifurcation, which suggests that βα has no effect on the value of U* at which the unstable behavior starts, since bifurcation is initially a linear phenomenon. On the other hand, the nonlinearity parameter does have an effect on the limit cycle oscillation amplitude in the period-1 regime between U*=6.25 and U*=13.42, where the standard deviation is approximately constant at a value of σαmax=1. In accordance with the previous results the standard deviation reaches a maximum in the deterministic bifurcation point of σαmax=10.3. Beyond U*=13.42 the standard deviation drops to the value σαmax=1 of the period-1 domain. Whether the response is period-1 or higher period does, therefore, not affect the influence of βα on the oscillation amplitude.


## 7. Random Initial Conditionα0(ω)

The results for randomness in the pitch initial conditionα0 are given in Figures 14–16. The mean μα shows a decaying oscillation to zero for U*=6.25 and periodic oscillations for U*=10 and U*=15, which closely resemble the deterministic results of Figure 6. The standard deviation σα also decays to zero for U*=6.25, and oscillates around only σα=1.4 and σα=0.5 at U*=10 and U*=15, respectively. For U*=13.42 the mean and standard deviation shows again a sudden irregular behavior with the standard deviation oscillating around σα=7 due to the discontinuity in the response surface of Figure 15(c). The PDF of α of Figure 16 shows also a multimodal character for U*=13.42 only.Time histories of the meanμα(τ) and standard deviation σα(τ) of the pitch angle α due to the random initial condition α0.
(a)(b)(c)(d)Response surface of the pitch angleα at τ=2000 as function of the random initial condition α0.
(a)(b)(c)(d)Bifurcation of the probability density function of the pitch angleα due to the random initial condition α0.
(a)(b)(c)(d)The bifurcation of the maximum standard deviation in Figure10(c) gives in the largest part of the bifurcation parameter domain a two times higher value of σαmax than for random βα with the identical input coefficient of variation of cvβα=4.48%. The system is, therefore, twice as sensitive to randomness in α0 than to random βα. The random initial condition results actually in a variation of the initial phase of the realizations. Phase differences in the response have a large effect on the stochastic behavior as we have observed for random ω̅. However, the phase differences do not increase in time for random α0 in contrast to the case with randomness in the ratio of natural frequencies. The gradually increasing σαmax due to random ω̅ is, therefore, in the majority of the bifurcation parameter range larger than the effect of random α0. This effect is not only caused by the larger input coefficient of variation for ω̅ but also mainly by the increasing phase differences in time.However, in the second deterministic bifurcation pointU*=13.42 the maximum standard deviation peaks for random α0 at a higher value of σαmax=10.0 compared to random ω̅. At U*=13.42 the effect of randomness in the parameters βα and α0 on σαmax is, therefore, larger than that of random ω̅, while in the rest of the bifurcation domain the parameters βα, α0, and ω̅, respectively, have a clear hierarchy of increasing importance. Even while ω̅ has twice the coefficient of variation compared to the other parameters, at the second deterministic bifurcation point the singularity results in larger variation in the response surface for βα and α0.In deterministically already highly computationally intensive problems subject to a large number of random input parameters, the actual uncertainty analysis is usually performed for a subset of the most important random input parameters only, which is selected based on preliminary results for a limited number of parameter settings. The current results should warn the reader that this can be a dangerously unreliable approach, since the importance of the random input parameters can highly depend on the chosen bifurcation parameter value. In isolated points in parameter space a clear relative importance of the random parameters can even suddenly reverse, such that none of the parameters can be disregarded in advance. A multidimensional treatment of the combined effect of multiple random parameters in stochastic aeroelastic applications will, therefore, be considered in future work.


## 8. Conclusions

The higher period stochastic bifurcation of a nonlinear airfoil flutter model is studied numerically. The fluid-structure interaction model consists of a two-degree-of-freedom rigid airfoil with cubic nonlinear springs and an aerodynamic model to determine the fluid loads in pitch and plunge. The employed uncertainty quantification method for unsteady problems is robust and efficient due to the extrema diminishing interpolation of oscillatory samples at constant phase.The effect on the time history of the pitch angleα is considered for randomness in the ratio of natural pitch and plunge frequencies ω̅, a nonlinear spring parameter βα, and the initial condition of the pitch angle α0. The random natural frequency ratio ω̅ affects the frequency of the response, which results in a gradual increase of the maximum standard deviation of the pitch angle in the asymptotic range to σαmax=8.0° in the second deterministic bifurcation point of U*=13.42. The output variability also starts to increase from the trivial solution at an earlier position compared to the first deterministic bifurcation point U*=6.25.The effect of uncertainty in the nonlinear stiffness parameterβα is approximately constant beyond U*=6.25 at σαmax=1 due to its effect on the limit cycle oscillation amplitude. At the second deterministic bifurcation point U*=13.42 random βα results in a sudden peak in the output standard deviation of σαmax=10.3. The random initial condition α0 reaches its maximum output standard deviation of σαmax=10.0 also in the second bifurcation point. In the rest of the bifurcation domain α0 results approximately in a two times higher output randomness than βα due to its effect on the phase of the response.Despite the largest effect ofω̅ in the majority of the bifurcation domain, βα and α0 are the most important sources of randomness at the second deterministic bifurcation point. This is caused by the larger variance in the response surface for βα and α0 due to the singularity at U*=13.42. Reducing the number of random input parameters based on preliminary results for a limited number of parameter settings can, therefore, give unreliable results, since the order of relative parameter importance can reverse in isolated singular points in parameter space.



---

*Source: 394387-2009-06-14.xml*