Higher Period Stochastic Bifurcation of Nonlinear Airfoil Fluid-Structure Interaction

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.

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 Section 2. 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 to 𝑛_𝐚 uncorrelated second-order random input parameters 𝐚(𝜔)={𝑎₁(𝜔),…,𝑎_𝑛a(𝜔)}∈𝐴 with parameter space 𝐴∈ℝ^𝑛𝐚, which governs an oscillatory response 𝑢(𝐱,𝑡,𝐚)

with operator ℒ and source term 𝑆 defined on domain 𝐷×𝑇×𝐴, and appropriate initial and boundary conditions. The spatial and temporal dimensions are defined as 𝐱∈𝐷 and 𝑡∈𝑇, respectively, with 𝐷⊂ℝ^𝑑, 𝑑={1,2,3}, and 𝑇=[0,𝑡_max]. A realization of the set of outcomes Ω of the probability space (Ω,ℱ,𝑃) is denoted by 𝜔∈Ω=[0,1]^𝑛𝐚, with ℱ⊂2^Ω the 𝜎-algebra of events and 𝑃 a probability measure.

Here we consider a nonintrusive uncertainty quantification method 𝑙 which constructs a weighted approximation 𝑤(𝐱,𝑡,𝐚) of response surface 𝑢(𝐱,𝑡,𝐚) based on 𝑛_𝐬 deterministic solutions 𝑣_𝑘(𝐱,𝑡)≡𝑢(𝐱,𝑡,𝐚_𝑘) of (2.1) for different parameter values 𝐚_𝑘≡𝐚(𝜔_𝑘) for 𝑘=1,…,𝑛_𝐬. The samples 𝑣_𝑘(𝐱,𝑡) can be obtained by solving the deterministic problem

for 𝑘=1,…,𝑛_s, using standard spatial discretization methods and time marching schemes. A nonintrusive uncertainty quantification method 𝑙 is then a combination of a sampling method 𝑔 and an interpolation method ℎ. Sampling method 𝑔 defines the 𝑛_s sampling points {𝐚_𝑘}^𝑛s_𝑘=1 and returns the deterministic samples 𝐯(𝐱,𝑡)={𝑣₁(𝐱,𝑡),…,𝑣_𝑛s(𝐱,𝑡)}. Interpolation method ℎ constructs an interpolation surface 𝑤(𝐱,𝑡,𝐚) through the 𝑛_𝐬 samples 𝐯(𝐱,𝑡) as an approximation of 𝑢(𝐱,𝑡,𝐚). We are eventually interested in an approximation of the probability distribution and statistical moments 𝜇_u𝑖(𝐱,𝑡) of the output 𝑢(𝐱,𝑡,𝐚), which can be obtained by sorting and weighted integration of 𝑤(𝐱,𝑡,𝐚):

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 Section 3.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 method 𝑙 evaluates integral (2.3) by dividing parameter space 𝐴 into 𝑛_e non-overlapping simplex elements 𝐴_𝑗⊂𝐴:

Here we consider a multielement Polynomial Chaos method based on Newton-Cotes quadrature points and simplex elements [47]. A piecewise polynomial approximation 𝑤(𝐱,𝑡,𝐚) is then constructed based on 𝑛_s deterministic solutions 𝑣_𝑗,𝑘(𝐱,𝑡)=𝑢(𝐱,𝑡,𝐚_𝑗,𝑘) for the values of the random parameters 𝐚_𝑗,𝑘 that correspond to the 𝑛_s̃ Newton-Cotes quadrature points of degree 𝑑 in the elements 𝐴_𝑗:

where 𝑐_𝑗,𝑘 is the weighted integral of the Lagrange interpolation polynomial 𝐿_𝑗,𝑘(𝐚) through Newton-Cotes quadrature point 𝑘 in element 𝐴_𝑗:

for 𝑗=1,…,𝑛_e and 𝑘=1,…,𝑛_s̃. Here, second-degree Newton-Cotes quadrature with 𝑑=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 𝐴 for the adaptive scheme consists of the minimum of 𝑛_eini=𝑛_a! simplex elements and 𝑛_sini=3^𝑛a 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 𝐴_𝑗 are adaptively refined using a refinement measure 𝜌_𝑗 based on the largest absolute eigenvalue of the Hessian 𝐻_𝑗, as a measure of the curvature of the response surface approximation in the elements, weighted by the probability 𝑓_𝑗 contained by the elements

(a)

(b)

(c)

Figure 1 

Discretization of two-dimensional parameter space 𝐴 using 2-simplex elements and second-degree Newton-Cotes quadrature points given by the dots.

with ∑^𝑛e_𝑗=10𝑥0200𝑑𝑓_𝑗=1. The stochastic grid refinement is terminated when convergence measure 𝛿_𝑛e is smaller than a threshold value 𝛿_𝑛e<𝛿 where

with 𝜇_w(𝐱,𝑡) and 𝜎_w(𝐱,𝑡) the mean and standard deviation of 𝑤(𝐱,𝑡,𝜔), or when a maximum number of samples 𝑛_s is reached. Convergence measure 𝛿_𝑛e 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 into 𝑛_ẽ=2^𝑛a 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

where the arguments 𝐱 and 𝑡 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 𝐚_𝑘 results in oscillatory samples 𝑣_𝑘(𝑡)=𝑢(𝐚_𝑘), of which the phase 𝑣_𝜙𝑘(𝑡)=𝜙(𝑡,𝐚_𝑘) is a well-defined monotonically increasing function of time 𝑡 for 𝑘=1,…,𝑛_s.

In order to interpolate the samples 𝐯(𝑡)={𝑣₁(𝑡),…,𝑣_𝑛a(𝑡)} at constant phase, first, their phase as function of time 𝐯_𝜙(𝑡)={𝑣_𝜙1(𝑡),…,𝑣_𝜙𝑛a(𝑡)} is extracted from the deterministic solves 𝐯(𝑡). Second, the time series for the phase 𝐯_𝜙(𝑡) are used to transform the samples 𝐯(𝑡) into functions of their phase 𝐯(̂𝐯_𝜙(𝑡)) according to

for 𝑘=1,…,𝑛_s, see Figure 2. Third, the sampled phases 𝐯_𝜙(𝑡) are interpolated to the function 𝑤_𝜙(𝑡,𝐚):

(a)

(b)

Figure 2 

Oscillatory samples as function of time and phase.

as approximation of 𝜙(𝑡,𝐚). Finally, the transformed samples 𝐯(̂𝐯_𝜙(𝑡)) are interpolated at a constant phase 𝜑∈𝑤_𝜙(𝑡,𝐚) to

Repeating the latter interpolation for all phases 𝜑∈𝑤_𝜙(𝑡,𝐚) results in the function 𝑤(𝑤_𝜙(𝑡,𝐚),𝐚). The interpolation 𝑤(𝑤_𝜙(𝑡,𝐚),𝐚) is then transformed back to an approximation in the time domain 𝑤(𝑡,𝐚) as follows:

The resulting function 𝑤(𝑡,𝐚) is an approximation of the unknown response surface 𝑢(𝑡,𝐚) as function of time 𝑡 and the random parameters 𝐚(𝜔). The actual sampling 𝑔 and interpolation ℎ is performed using the extrema diminishing uncertainty quantification method 𝑙 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 𝜀̂(𝜑,𝐚)=|𝑤(𝜑,𝐚)−𝑢(̂𝜑,𝐚)| as function of phase 𝜑 for periodic responses according to

where 𝛿 is defined by

The error 𝜀(𝑡,𝐚)=|𝑤(𝑡,𝐚)−𝑢(𝑡,𝐚)| is also bounded in time under certain conditions, see [44].

The phases 𝐯_𝜙(𝑡) are extracted from the samples based on the local extrema of the time series 𝐯(𝑡). 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 𝜀_𝑘 for at least one additional cycle length. The functions for the phases 𝐯_𝜙(𝑡) in the whole time domain 𝑇 are constructed by identifying all successive cycles of 𝐯(𝑡) and linear extrapolation to 𝑡=0 and 𝑡=𝑡_max 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 𝑡=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 𝑣_𝜙𝑘(𝑡) cannot be extracted from one of the samples 𝑣_𝑘(𝑡) for 𝑘=1,…,𝑛_s, then uncertainty quantification interpolation ℎ is directly applied to the time-dependent samples 𝐯(𝑡).

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 Section 4.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 Figure 3, 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] by

Figure 3 

The two-degree-of-freedom airfoil flutter model.

where 𝛼(𝜏) is the pitch angle and 𝜉(𝜏)=ℎ/𝑏 is the nondimensional version of the plunge displacement ℎ of the elastic axis, with 𝑏=𝑐/2 the half-chord, and initial conditions 𝛼(0)=𝛼₀ and 𝜉(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 𝑚/𝜋𝜌𝑏², with 𝑚 the airfoil mass, and 𝜌 the air density. The radius of gyration about the elastic axis is 𝑟_𝛼, where elastic axis is located at a distance 𝑎_h𝑏 from the mid-chord point, and the mass center is located at a distance 𝑥_𝛼𝑏 from the elastic axis. The bifurcation parameter is the ratio of time scales of the structure and the flow defined as 𝑈^∗=𝑈/(𝑏𝜔_𝛼), with 𝑈 the free stream velocity. The primes denote differentiation with respect to nondimensionalized time 𝜏=𝑈𝑡/𝑏. The expressions for the aerodynamic force and moment coefficients, 𝐶_L(𝜏) and 𝐶_M(𝜏) are given by Fung [49] as

where 𝜙(𝜏) is the Wagner function

with the constants 𝜓₁=0.165, 𝜓₂=0.335, 𝜀₁=0.0455, and 𝜀₂=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]

with

where a vector {𝑥_𝑖}^8𝑖=1 of new variables is defined as

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, 𝑎_h=−0.5, 𝑥_𝛼=0.25, 𝑟_𝛼=0.5, 𝜉₀=0, 𝛽_𝜉=0, and 𝜁_𝛼=𝜁_𝜉=0 as in [48].

The system parameters that are assumed to be uncertain are 𝜔, 𝛽_𝛼, and 𝛼₀. The randomness of these three parameters is described by a symmetric unimodal beta distribution with 𝛽₁=𝛽₂=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 𝛼₀ has an interval of 𝛼₀∈[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 𝛼₀ 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 𝑈^∗∈[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 𝑈^∗=6.25 the system response is a decaying oscillation to 𝛼=0. At 𝑈^∗=6.25 the system exhibits a supercritical Hopf bifurcation to a stable period-1 limit cycle oscillation with an increasing amplitude for increasing 𝑈^∗. At the second bifurcation point 𝑈^∗=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 𝑈^∗={5,10,15}. At 𝑈^∗=5 both 𝛼 and 𝜉 are decaying oscillations to the stable fixed point (𝛼,𝜉)=(0,0). A period-1 oscillation can be identified for 𝛼 and 𝜉 at 𝑈^∗=10. For 𝑈^∗=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.

(a)

(b)

(c)

Figure 6 

Time histories of pitch 𝛼(𝜏) and plunge 𝜉(𝜏) in the three regimes of the deterministic airfoil flutter system.

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 𝑈^∗={6.25;10;13.42;15} are considered. This corresponds to the first (𝑈^∗=6.25) and second (𝑈=13.42) deterministic bifurcation point, and the period-1 (𝑈^∗=10) and higher period (𝑈^∗=15) regime. The case of 𝑈^∗=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 𝑛_s in the stochastic simulations for the different values of 𝑈^∗ 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 𝑛_s=10³ samples.

Table 1 

Relation between convergence measure 𝛿_𝑛e and 𝐿_∞ error 𝜀_𝐿∞ for the mean 𝜇_𝛼 and standard deviation 𝜎_𝛼 at 𝑈^∗=6.25.

Table 2 

Relation between convergence measure 𝛿_𝑛e and 𝐿_∞ error 𝜀_𝐿∞ for the mean 𝜇_𝛼 and standard deviation 𝜎_𝛼 at 𝑈^∗=10.

Table 3 

Relation between convergence measure 𝛿_𝑛e and 𝐿_∞ error 𝜺_𝐿∞ for the mean 𝜇_𝛼 and standard deviation 𝜎_𝛼 at 𝑈^∗=13.42.

Table 4 

Relation between convergence measure 𝛿_𝑛e and 𝐿_∞ error 𝜺_𝐿∞ for the mean 𝜇_𝛼 and standard deviation 𝜎_𝛼 at 𝑈^∗=15.

(a)

(b)

(c)

(d)

Figure 7 

Time histories of the mean 𝜇𝛼(𝜏) and standard deviation 𝜎𝛼(𝜏) of the pitch angle 𝛼 due to the random frequency ratio 𝜔.

(a)

(b)

(c)

(d)

Figure 8 

Response surface of the pitch angle 𝛼 at 𝜏=2000 as function of the random frequency ratio 𝜔.

(a)

(b)

(c)

(d)

Figure 9 

Bifurcation of the probability density function of the pitch angle 𝛼 due to the random frequency ratio 𝜔.

For 𝑈^∗=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 𝑈^∗=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 𝑈^∗=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 𝑈^∗. 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 point 𝑈^∗=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 𝑈^∗=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 𝑈^∗=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 samples 𝑛_s used in the stochastic simulations depends significantly on the value of bifurcation parameter 𝑈^∗. In Tables 1–4 the convergence 𝛿_ne and the error 𝜀_𝐿∞ 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 𝐿_∞ error𝜀_𝐿∞ is defined for 𝜇_𝛼 as

and equivalently for 𝜎_𝛼. The method is highly efficient in the periodic regimes, in which 𝑛_s=3 samples is already sufficient to match the Monte Carlo results based on 𝑛_s=10³ samples. This holds even for the oscillatory response surface in the higher period case of 𝑈^∗=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 𝑈^∗=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 𝑈^∗=10. The multimodal PDF evolves further into a distribution with 6 peaks at approximately 𝛼={±5,±11,±16} due to the higher period motion at 𝑈^∗=15. At the second deterministic bifurcation point 𝑈^∗=13.42 the PDF is in an intermediate state between the approximately symmetric multimodal distributions of 𝑈^∗=10 and 𝑈^∗=15.

The stochastic behavior of the system is also shown in Figure 10 for the three random parameters in terms of the bifurcation of the maximum standard deviation 𝜎_𝛼max in the asymptotic range defined as

(a)

(b)

(c)

Figure 10 

Bifurcation of the maximum of the standard deviation 𝜎𝛼max for the three random parameters.

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 𝑈^∗=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 𝑛_s in the simulations for random 𝛽_𝛼 is again determined based on convergence studies.

(a)

(b)

(c)

(d)

Figure 11 

Time histories of the mean 𝜇𝛼(𝜏) and standard deviation 𝜎𝛼(𝜏) of the pitch angle 𝛼 due to the random nonlinearity parameter 𝛽𝛼.

(a)

(b)

(c)

(d)

Figure 12 

Response surface of the pitch angle 𝛼 at 𝜏=2000 as function of the random nonlinearity parameter 𝛽𝛼.

(a)

(b)

(c)

(d)

Figure 13 

Bifurcation of the probability density function of the pitch angle 𝛼 due to the random nonlinearity parameter 𝛽𝛼.

For 𝑈^∗=6.25, 𝑈^∗=10, and 𝑈^∗=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 𝑈^∗=6.25, which suggests that randomness in 𝛽_𝛼 does not lead to an earlier bifurcation. For both 𝑈^∗=10 and 𝑈^∗=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 𝑈^∗ 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 𝑈^∗=13.42. The stochastic system shows at 𝑈^∗=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 𝑈^∗=13.42. The adaptive method resolves also this discontinuous response accurately and the other response surface approximations require again only 𝑛_s=3 deterministic simulations to match the Monte Carlo results.

The PDF in Figure 13 is also significantly distorted from the unimodal input distribution at 𝑈^∗=13.42 only. For 𝑈^∗=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 𝑈^∗ 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 𝑈^∗=6.25 and 𝑈^∗=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 𝑈^∗=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 𝛼₀(𝜔)

The results for randomness in the pitch initial condition 𝛼₀ are given in Figures 14–16. The mean 𝜇_𝛼 shows a decaying oscillation to zero for 𝑈^∗=6.25 and periodic oscillations for 𝑈^∗=10 and 𝑈^∗=15, which closely resemble the deterministic results of Figure 6. The standard deviation 𝜎_𝛼 also decays to zero for 𝑈^∗=6.25, and oscillates around only 𝜎_𝛼=1.4 and 𝜎_𝛼=0.5 at 𝑈^∗=10 and 𝑈^∗=15, respectively. For 𝑈^∗=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 𝑈^∗=13.42 only.

(a)

(b)

(c)

(d)

Figure 14 

Time histories of the mean 𝜇𝛼(𝜏) and standard deviation 𝜎𝛼(𝜏) of the pitch angle 𝛼 due to the random initial condition 𝛼0.

(a)

(b)

(c)

(d)

Figure 15 

Response surface of the pitch angle 𝛼 at 𝜏=2000 as function of the random initial condition 𝛼0.

(a)

(b)

(c)

(d)

Figure 16 

Bifurcation of the probability density function of the pitch angle 𝛼 due to the random initial condition 𝛼0.

The bifurcation of the maximum standard deviation in Figure 10(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 𝛼₀ 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 𝛼₀ 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 𝛼₀. 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 point 𝑈^∗=13.42 the maximum standard deviation peaks for random 𝛼₀ at a higher value of 𝜎_𝛼max=10.0 compared to random 𝜔. At 𝑈^∗=13.42 the effect of randomness in the parameters 𝛽_𝛼 and 𝛼₀ on 𝜎_𝛼max is, therefore, larger than that of random 𝜔, while in the rest of the bifurcation domain the parameters 𝛽_𝛼, 𝛼₀, 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 𝛼₀.

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 𝛼₀. 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 𝑈^∗=13.42. The output variability also starts to increase from the trivial solution at an earlier position compared to the first deterministic bifurcation point 𝑈^∗=6.25.

The effect of uncertainty in the nonlinear stiffness parameter 𝛽_𝛼 is approximately constant beyond 𝑈^∗=6.25 at 𝜎_𝛼max=1 due to its effect on the limit cycle oscillation amplitude. At the second deterministic bifurcation point 𝑈^∗=13.42 random 𝛽_𝛼 results in a sudden peak in the output standard deviation of 𝜎_𝛼max=10.3. The random initial condition 𝛼₀ reaches its maximum output standard deviation of 𝜎_𝛼max=10.0 also in the second bifurcation point. In the rest of the bifurcation domain 𝛼₀ 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 𝛼₀ 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 𝛼₀ due to the singularity at 𝑈^∗=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.

Acknowledgments

The presented work is supported by the NODESIM-CFD project (Non-Deterministic Simulation for CFD based design methodologies); a collaborative project funded by the European Commission, Research Directorate-General in the 6th Framework Programme, under contract AST5-CT-2006-030959.