Analysis of a Nonlinear Aeroelastic System with Parametric Uncertainties Using Polynomial Chaos Expansion

Abstract

Aeroelastic stability remains an important concern for the design of modern structures such as wind turbine rotors, more so with the use of increasingly flexible blades. A nonlinear aeroelastic system has been considered in the present study with parametric uncertainties. Uncertainties can occur due to any inherent randomness in the system or modeling limitations, and so forth. Uncertainties can play a significant role in the aeroelastic stability predictions in a nonlinear system. The analysis has been put in a stochastic framework, and the propagation of system uncertainties has been quantified in the aeroelastic response. A spectral uncertainty quantification tool called Polynomial Chaos Expansion has been used. A projection-based nonintrusive Polynomial Chaos approach is shown to be much faster than its classical Galerkin method based counterpart. Traditional Monte Carlo Simulation is used as a reference solution. Effect of system randomness on the bifurcation behavior and the flutter boundary has been presented. Stochastic bifurcation results and bifurcation of probability density functions are also discussed.

1. Introduction

Fluid-structure interaction can result in dynamic instabilities like flutter. Nonlinear parameters present in the system can stabilize the diverging growth of flutter oscillations to a limit cycle oscillation (LCO). Sustained LCO can lead to fatigue failure of rotating structures such as wind turbine rotors. Hence, it is an important design concern in aeroelastic analysis. Moreover, there is a growing interest in understanding how system uncertainties in structural and aerodynamic parameters and initial conditions affect the characteristics of such dynamical response.

Uncertainty quantification in a stochastic framework with stochastic inputs has traditionally been analyzed with Monte Carlo simulations (MCSs). To apply this procedure one should use the distribution of the input parameters to generate a large number of realizations of the response. Probability density function (PDF) and other required statistics are then approximated from these realizations; however, it is computationally expensive, especially for large complex problems. Hence, there is a need to develop alternate approaches which are computationally cheaper than direct MCS procedure. Perturbation method is a fast tool for obtaining the response statistics in terms of its first and second moments [1]. The statistical response is determined by expanding the stochastic parameters around their mean via a Taylor series [2]. The application of this method is, however, limited to small perturbations and does not readily provide information on high-order statistics [3, 4]. The resulting system of equations becomes extremely complicated beyond second-order expansions as shown in the literature. Sensitivity method is a more economical approach, based on the moments of samples, but it is less robust and depends strongly on the modeling assumptions [5]. Another approach based on expanding the inverse of the stochastic operator in a Neumann series is also limited to small fluctuations only; even combining with the Monte Carlo method also seems to result in a computationally prohibitive algorithm for complex systems [4].

Polynomial chaos expansion (PCE) is a more effective approach, pioneered by Ghanem and Spanos [4], proposed first in the structural mechanics finite elements area. It is a spectral representation of the uncertainty in terms of orthogonal polynomials. The stochastic input is represented spectrally by employing orthogonal polynomial functionals from the Askey scheme as basis in the random space. The original homogeneous PCE was based on Hermite polynomials from the Askey family [6]. It can give optimal exponential convergence for Gaussian inputs [7]. A standard Galerkin projection is applied along the random dimensions to obtain the weak form of the equations. The resulting deterministic systems are solved using standard techniques to solve for each random mode [8]. Galerkin polynomial chaos expansion (Galerkin PCE) based approaches have been examined extensively with different basis functions to model several uncertain flow and flow-induced instability problems [9, 10].

Galerkin PCE (also called intrusive approach) modifies the governing equations to a coupled form in terms of the chaos coefficients. These equations are usually more complex and arriving at them is quite often a tedious task for some choices of the uncertain parameters. In order to avoid these, several uncoupled alternatives have been proposed. These are collectively called nonintrusive approaches. In a nonintrusive polynomial chaos method a deterministic solver is used repeatedly as in Monte Carlo simulation. The Probabilistic Collocation (PC) method is such a nonintrusive polynomial chaos method in which the problem is collocated at Gauss quadrature points in the probability space [11, 12]. The deterministic solutions are performed at these collocation points. The nonintrusive polynomial chaos method proposed by Walters and coworkers [13–15] is based on approximating the polynomial chaos coefficients. A similar approach called nonintrusive Spectral Projection has been used by Reagan et al. [16]. Pettit and Beran [17, 18] have also used a stochastic projection technique to compute the chaos expansion coefficients in an aeroelastic system. When multiple uncertain parameters are involved the collocation grids are constructed using tensor products of one-dimensional grids. Thus, the number of collocation points and therefore the number of required deterministic solutions increases rapidly. As an alternative, sparse grid collocation approaches can be implemented [19–21].

The intrusive and nonintrusive PCE approaches and their implementation to an aeroelastic model with structural nonlinearity are discussed in detail in the subsequent sections.

2. Nonlinear Aeroelastic Model

Figure 1 shows a schematic plot of the two degree-of-freedom pitch-plunge aeroelastic system and also the notations used in the analysis. The aeroelastic equations of motion for the linear system have been derived by Fung [22]. For nonlinear restoring forces such as with cubic springs in both pitch and plunge, the mathematical formulation is given by Lee et al. [23] as follows:𝜀^′′+𝑥_𝛼𝛼^′′+2𝜁_𝜀𝜔𝑈𝜀^′+(𝜔𝑈)²(𝜀+𝛽_𝜀𝜀³1)=−𝐶𝜋𝜇^𝐿𝑥(𝜏),_𝑎𝑟²_𝛼𝜀^′′+𝛼^′′𝜁+2^𝛼𝑈𝛼^′+1𝑈²(𝛼+𝛽_𝛼𝛼³2)=𝜋𝜇𝑟²_𝛼𝐶_𝑀(𝜏). (2.1) The above equations are shown in the nondimensional form. The nondimensional parameters are given below. The plunge deflection is considered positive in the downward direction and the pitch angle about the elastic axis is denoted positive nose up. Elastic axis is located at a distance 𝑎_ℎ𝑏 from the midchord where 𝑏 is the half chord. Let us also use 𝑣 as the wind velocity ℎ as the plunge deflection. Among the nondimensional quantities, 𝜀=ℎ/𝑏 = nondimensional displacement of the elastic axis point; 𝜏=𝑣𝑡/𝑏 = nondimensional time; 𝑈=𝑣/(𝑏𝜔_𝛼)= the nondimensional velocity (also called reduced velocity); 𝜔=𝜔_𝜀/𝜔_𝛼, 𝜔_𝜀 and 𝜔_𝛼 are the natural frequencies of the uncoupled plunging and pitching modes respectively. In the structural part, 𝜁_𝜀 and 𝜁_𝛼 are the damping ratios in plunge and pitch respectively, 𝑟_𝛼 is the radius of gyration about the elastic axis, and 𝜇 is the airfoil mass ratio defined as 𝑚/(𝜋𝜌𝑣). 𝛽_𝛼 and 𝛽_𝜀 denote coefficients of cubic spring in pitch and plunge respectively. For incompressible, inviscid flow, Fung [22] gives the expressions for unsteady lift and pitching moment coefficients, 𝐶_𝐿(𝜏) and 𝐶_𝑀(𝜏):𝐶_𝐿(𝜏)=𝜋(𝜀^′′−𝑎_ℎ𝛼^′′+𝛼^′)+2𝜋{𝛼(0)+𝜀^′1(0)+[2−𝑎_ℎ]𝛼^′(0)}𝜙(𝜏)+2𝜋^𝜏₀𝜙(𝜏−𝜎)[𝛼^′(𝜎)𝜀^′′1(𝜎)+[2−𝑎_ℎ]𝛼^′′𝐶(𝜎)]𝑑𝜎,_𝑀1(𝜏)=𝜋[2+𝑎_ℎ]×{𝛼(0)+𝜀^′1(0)+[2−𝑎_ℎ]𝛼^′1(0)}𝜙(𝜏)+𝜋[2+𝑎_ℎ]×^𝜏₀𝜙(𝜏−𝜎){𝛼^′(𝜎)+𝜀^′′1(𝜎)+[2−𝑎_ℎ]𝛼^′′+𝜋(𝜎)}𝑑𝜎2𝑎_ℎ(𝜀^′′−𝑎_ℎ𝛼^′′1)−[2−𝑎_ℎ]𝜋2𝛼^′−𝜋𝛼16^′′ (2.2) The Wagner function 𝜙(𝜏) is given by:𝜙(𝜏)=1−𝜓₁e^−𝜖1𝜏−𝜓₂e^−𝜖2𝜏. (2.3) Values for the constants are, 𝜓₁= 0.165, 𝜓₂= 0.335, 𝜖₁= 0.0455 and 𝜖₂ = 0.3 [24]. Introducing the following new variables 𝑤₁,𝑤₂,𝑤₃,and𝑤₄ [23], the original integrodifferential equations for aeroelastic system given by (2.1) are reformulated: 𝑤₁=^𝜏₀𝑒^{−𝜖1(𝜏−𝜎)}𝑤𝛼(𝜎)𝑑𝜎,₂=^𝜏₀𝑒^{−𝜖2(𝜏−𝜎)}𝑤𝛼(𝜎)𝑑𝜎,₃=^𝜏₀𝑒^{−𝜖1(𝜏−𝜎)}𝑤𝜀(𝜎)𝑑𝜎,₄=^𝜏₀𝑒^{−𝜖2(𝜏−𝜎)}𝜀(𝜎)𝑑𝜎. (2.4) Now a set of autonomous differential equations of the form 𝐗^′=𝑓(𝐗) are obtained as, 𝐗={𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅,𝑥₆,𝑥₇,𝑥₈} = {𝛼,𝛼^′,𝜀,𝜀^′,𝑤₁,𝑤₂,𝑤₃,𝑤₄}.

Figure 1 

The schematic of a symmetric airfoil with pitch and plunge degrees-of-freedom.

Explicitly, the system looks like,𝑥^′₁=𝑥₂,𝑥^′₂=(𝑐₀𝑁−𝑑₀𝑀)(𝑐₁𝑑₀−𝑐₀𝑑₁),𝑥^′₃=𝑥₄,𝑥^′₄=(−𝑐₁𝑁+𝑑₁𝑀)(𝑐₁𝑑₀−𝑐₀𝑑₁),𝑥^′₅=𝑥₁−𝜖₁𝑥₅,𝑥^′₆=𝑥₁−𝜖₂𝑥₆,𝑥^′₇=𝑥₃−𝜖₁𝑥₇,𝑥^′₈=𝑥₃−𝜖₂𝑥₈, (2.5) where 𝑀=𝑐₂𝑥₄+𝑐₃𝑥₂+𝑐₄𝑥₃+𝑐₅𝑥³₃+𝑐₆𝑥₁+𝑐₇𝑥₅+𝑐₈𝑥₆+𝑐₉𝑥₇+𝑐₁₀𝑥₈−𝑓(𝜏),𝑁=𝑑₂𝑥₂+𝑑₃𝑥₁+𝑑₄𝑥³₁+𝑑₅𝑥₄+𝑑₆𝑥₃+𝑑₇𝑥₅+𝑑₈𝑥₆+𝑑₉𝑥₇+𝑑₁₀𝑥₈−𝑔(𝜏). (2.6) The coefficients 𝑐₀⋯𝑐₁₀ and 𝑑₀⋯𝑑₁₀ depend on the system parameters, and their expressions along with 𝑓(𝜏) and 𝑔(𝜏) are given in the appendix.

3. Uncertainty Quantification and Polynomial Chaos Expansion

It is increasingly being felt among the aeroelastic community that aeroelastic analysis should include the effect of parametric uncertainties. This can potentially revolutionize the present design concepts with higher rated performance and can also reshape the certification criteria. Nonlinear dynamical systems are known to be sensitive to physical uncertainties, since they often amplify the random variability with time. Hence, quantifying the effect of system uncertainties on the aeroelastic stability boundary is crucial. Flutter, a dynamic aeroelastic instability involves a Hopf bifurcation where a damped (stable response) oscillation changes to a periodic oscillatory response at a critical wind velocity. In a linear system the post flutter response can grow in an unbounded fashion [22]. System parametric uncertainties can significantly affect the onset and properties of bifurcation points. The importance of stochastic modeling of these uncertainties is that they quantify the effect of the uncertainties on flutter and bifurcation in a probabilistic sense and gives the response statistics in a systematic manner.

The original homogeneous polynomial chaos expansion [4] is based on the homogeneous chaos theory of Wiener [6, 25]. This is based on a spectral representation of the uncertainty in terms of orthogonal polynomials. In its original form, it employs Hermite polynomials as basis from the generalized Askey scheme and Gaussian random variables. Spectral polynomial chaos-based approaches with other basis functions have also been used in the recent past in various unsteady flow and flow-structure interaction problems of practical interest [8, 26, 27].

3.1. Classical Galerkin Polynomial Chaos Approach

In the classical Galerkin-PCE approach, the polynomial chaos expansion of the system response is substituted into the governing equation and a Galerkin error minimization in the probability space is followed. This results in a set of coupled equations in terms of the polynomial chaos coefficients. The resulting system is deterministic, but it is significantly modified to a higher order and complexity depending on the order of chaos expansion and system nonlinearity. After solving this set of coefficient equations, they are substituted back to get the system response.

As per the Cameron-Martin theorem [28], a random process 𝑋(𝑡,𝜃) (as function of random event 𝜃) which is second-order stationary can be written as𝑋(𝑡,𝜃)=𝑎₀𝜓₀+^∞_𝑖1=1𝑎_i1𝜓₁(𝜉_𝑖1(𝜃))+^∞_𝑖1^𝑖₌₁¹_𝑖2=1𝑎_𝑖1𝑖2𝜓₂(𝜉_i1(𝜃),𝜉_i2+(𝜃))_{∞𝑖1𝑖=11𝑖2𝑖=12𝑖3=1}(𝑎_i1𝑖2𝑖3)𝜓₃(𝜉_𝑖1(𝜃),𝜉_𝑖2(𝜃),𝜉_𝑖3(𝜃))+⋯, (3.1) where 𝜓_𝑛(𝜉_𝑖1,𝜉_𝑖2,…,𝜉_𝑖𝑛) denotes the Hermite polynomial of order 𝑛 in terms of 𝑛-dimensional independent standard Gaussian random variables 𝝃=(𝜉_𝑖1,𝜉_𝑖2,…,𝜉_𝑖𝑛) with zero mean and unit variance. The above equation is the discrete version of the original Wiener polynomial chaos expansion, and the continuous integrals are replaced by summations. For notational convenience equation (3.1) can be rewritten as𝑋(𝑡,𝜃)=^∞_𝑗=0𝑎_𝑗(𝑡)Φ_𝑗(𝝃(𝜃)). (3.2) There is an one-to-one relationship between the 𝜓's and Φ's and also 𝑎_𝑗's and 𝑎_𝑗's in (3.1) and (3.2). In the original form, chaos expansion uses Hermite polynomials (𝜓_𝑛s). The form of the one-dimensional Hermite polynomials is given as follows.𝜓₀𝜓=1,₁𝜓=𝜉,₂=𝜉²𝜓−1,₃=𝜉³𝜓−3𝜉,₄=𝜉⁴−6𝜉²𝜓+3,_𝑛=𝜉𝜓_𝑛−1−(𝑛−1)𝜓_𝑛−2. (3.3) One can also use orthogonal polynomials from the generalized Askey scheme for some standard nonGaussian input uncertainty distributions such as gamma and beta [8]. For any arbitrary input distribution, a Gram-Schmidt orthogonalization can be employed to generate the orthogonal family of polynomials [29].

Any stochastic process 𝛼(𝑡,𝝃(𝜃)), governed by Gaussian random variables 𝝃 (𝜉 can always be normalized as a standard Gaussian one), can be approximated by the following truncated series:𝛼(𝑡,𝝃(𝜃))=^𝑝_𝑗=0𝛼_𝑗(𝑡)Φ_𝑗(𝝃(𝜃)) (3.4) Note that, here the infinite upper limit of (3.2) is replaced by 𝑝, called the order of the expansion. For 𝑛 number of random variable and polynomial order 𝑛_𝑝, 𝑝 is given by the following [26]:𝑝=(𝑛+𝑛_𝑝)!𝑛!𝑛_𝑝!−1. (3.5) We demonstrate the Galerkin-PCE approach for a generalized dynamical system for a single random variable case, that is, with a random cubic stiffness. Let us write the governing equation with cubic nonlinearity in the following form [27]:£[𝛼(𝑡,𝜃)]+𝛽_𝛼[𝛼(𝑡,𝜃)]³=0, (3.6) here £ is a linear differential operator.

Equation (3.4) is now rewritten for a single random variable 𝜉 as𝛼(𝑡,𝜉(𝜃))=^𝑝_𝑗=0𝛼_𝑗(𝑡)Φ_𝑗(𝜉(𝜃)). (3.7) Here Φ_𝑗’s are now Hermite polynomials 𝜓_𝑗 as shown in (3.3).

If the cubic spring constant 𝛽_𝛼 is assumed to be a Gaussian random variable with mean 𝛽_𝛼 and standard deviation 𝛽_𝛼, it can be characterized by𝛽_𝛼=𝛽_𝛼+𝜉𝛽_𝛼=¹_𝑙=0𝛽_𝛼𝑙Φ_𝑙, (3.8) with, 𝛽_𝛼0=𝛽_𝛼 and 𝛽_𝛼1=𝛽_𝛼.

Substituting the chaos expansion terms, (3.7) and (3.8) in (3.6),£[^𝑝_𝑗=0𝛼_𝑗(𝑡)Φ_𝑗]+¹_𝑙=0𝛽_𝛼𝑙Φ_𝑙[^𝑝_𝑗=0𝛼_𝑗(𝑡)Φ_𝑗]³=0. (3.9) The cubic nonlinear function can be expressed in the following form:[^𝑝_𝑗=0𝛼_𝑗(𝑡)Φ_𝑗]³=^𝑝𝑝_𝑖=0^𝑝_𝑚=0𝑛=0𝛼_𝑖𝛼_𝑚𝛼_𝑛Φ_𝑖Φ_𝑚Φ_𝑛. (3.10) Substituting (3.10) into (3.9) and simplifying, we get,£[^𝑝_𝑗=0𝛼_𝑗(𝑡)Φ_𝑗]+¹_𝑙=0𝛽_𝛼𝑙Φ_𝑙[^𝑝𝑝_𝑖=0^𝑝_𝑚=0𝑛=0𝛼_𝑖𝛼_𝑚𝛼_𝑛Φ_𝑖Φ_𝑚Φ_𝑛]=0. (3.11) Using Galerkin projection on (3.11) by taking ⟨⋅,Φ_𝑘⟩, for 𝑘=0,1,…,𝑝,£[𝛼_𝑘1(𝑡)]+⟨Φ²_𝑘⟩[^1𝑝_𝑙=0^𝑝_𝑖=0^𝑝_𝑚=0𝑛=0𝛽_𝛼𝑙𝛼_𝑖𝛼_𝑚𝛼_𝑛⟨Φ_𝑙Φ_𝑖Φ_𝑚Φ_𝑛Φ_𝑘⟩]=0. (3.12) The expected value operator ⟨⋅⟩, called the inner product, is defined as,⟨Φ_𝑙⋯Φ_𝑘⟩=^∞_−∞Φ_𝑙⋯Φ_𝑘𝜔(𝜉)d𝜉. (3.13) For Hermite polynomials the weighting function 𝜔(𝜉) is the Gaussian probability density function. For single random variable case it is given as1𝜔(𝜉)=(√2𝜋)𝑒^{−((1/2)𝜉2)}. (3.14) The Hermite polynomials are orthogonal with respect to this weighting function in the Hilbert space. The polynomial chaos forms a complete orthogonal basis in the 𝐿₂ space of real-valued functions depending on the Gaussian random variables; hence the inner product of two orthogonal polynomial can be replaced by the identity⟨Φ_𝑙Φ_𝑘⟩=⟨Φ²_𝑙⟩𝛿_𝑙𝑘, (3.15)𝛿_𝑙𝑘 is the Kronecker delta function, given as:𝛿_𝑙𝑘=1if𝑙=𝑘,0otherwise. (3.16) The inner product terms in (3.12) ⟨Φ_𝑙Φ_𝑖Φ_𝑚Φ_𝑛Φ_𝑘⟩ and ⟨Φ²_𝑘⟩ can be evaluated analytically before-hand and substituted in the equation. The resulting system becomes a deterministic differential equation in terms of the chaos coefficients. Depending on the type of nonlinearity, the number of random variables and the number of expansion terms, evaluating the inner products could be tedious. In the present study, they are computed by numerical integration by using Gauss-Hermite quadrature rule and verified analytically by using symbolic mathematical solver Mathematica. Some typical nonzero inner-products are given in Table 1.

Table 1 

The Galerkin approach is also called the intrusive approach as it modifies the system governing equations in terms of the chaos coefficients. The modification results into a higher order and much more complex form. As a result, this approach may become computationally quite expensive.

3.2. Nonintrusive Projection Method

A number of nonintrusive variants of PCE have been developed to counter the disadvantages of the classical Galerkin method. Stochastic projection is one of them [4, 30]. In the present study, a stochastic projection-based approach is used to evaluate the chaos coefficients. Here, the chaos expansions are not substituted in the governing equations; instead samples of the solutions are used (using low-order deterministic simulations) to evaluate the coefficients directly using a projection formula. As a result, this approach can utilize the existing deterministic code and hence the name nonintrusive. The random process is approximated by a truncated series, as shown in (3.7).

The Hermite polynomials are statistically orthogonal, that is, they satisfy ⟨Φ_𝑖Φ_𝑗⟩=0 for 𝑖≠𝑗, hence the expansion coefficients can be directly evaluated as𝛼_𝑗(𝑡)=⟨𝛼(𝑡,𝜉(𝜃))Φ_𝑗⟩⟨Φ²_𝑗⟩. (3.17) The denominator in (3.17) can be shown to satisfy ⟨Φ²_𝑗⟩=𝑗! for nonnormalized Hermite polynomials [31]. So the key step in projecting 𝛼(𝑡,𝜉(𝜃)) along the polynomial chaos basis is the evaluation of ⟨𝛼(𝑡,𝜉(𝜃))Φ_𝑗⟩. The inner product is given by the following integral:⟨𝛼(𝑡,𝜉(𝜃))Φ_𝑗⟩=^∞_−∞𝛼(𝑡,𝜉(𝜃))Φ_𝑗𝜔(𝜉)d𝜉. (3.18) A Gauss-Hermite quadrature will be suitable for evaluating the above as the domain is (−∞,∞) and the weight function is Gaussian PDF. The quadrature points are the zeros of the Hermite polynomials of chosen order. A number of deterministic runs are performed at the quadrature points which is much lower than the full Monte Carlo simulations. We refer to this step as a pseudo-Monte Carlo simulation approach. In the pseudo-MCS approach the samples of 𝛽_𝛼 are generated from the corresponding 𝜉 values which are the Gauss-Hermite quadrature points. The realizations of the system response 𝛼(𝑡,𝜃) are then used to estimate the deterministic coefficients, 𝛼_𝑗(𝑡)’s, in (3.17) using the Gauss-Hermite quadrature rule. It should also be noted that for each evaluation of the inner product integral a convergence study is done by gradually increasing the number of quadrature points.

4. Results and Discussions

The main focus of the present study is quantifying the effect of system uncertainties on the bifurcation behavior and the flutter boundary of the nonlinear aeroelastic system. A fourth order variable step Runge-Kutta method is employed for the time integration. The main bifurcation parameter in a flutter system is the nondimensional wind velocity, also called the reduced velocity. In a linear aeroelastic system, the response changes to an exponentially growing solution from a stable damped oscillation at some critical wind velocity, known as the linear flutter speed. Nonlinear aeroelastic system can stabilize the response at the post-flutter regime to limit cycle oscillations [23] and the critical point becomes a Hopf bifurcation point. With a cubic nonlinearity, both supercritical and subcritical Hopf bifurcations are possible [32]. The latter case is observed for a softening cubic spring. Here in the stochastic analysis, we focus on the supercritical case. A deterministic bifurcation diagram with the following parameter values [23] is shown in Figure 2:𝜇=100,𝜔=0.2,𝑎_ℎ=−0.5,𝑥_𝛼=0.25,𝜁_𝛼=0,𝜁_𝜀=0,𝑟_𝛼=0.5,𝛽_𝛼=3,𝛽_𝜀=0. The variation of the limit cycle oscillation (LCO) amplitude is plotted with reduced velocities. Bifurcation occurs at the corresponding linear flutter speed of 6.285, and the observation match well with the earlier results [33]. At the post flutter velocities, limit cycle oscillations are observed and the amplitude of the LCOs increase as the reduced velocity increases.

Figure 2 

Deterministic flutter and bifurcation diagram with cubic nonlinearity (supercritical Hopf bifurcation).

We now consider random variations in the system parameters and investigate the influence on the overall dynamics. We consider only single uncertain parametric variation in this paper, that is, a single random variable model. First, the hardening cubic spring constant is considered to be a Gaussian random variable with mean 𝛽_𝛼 = 3 and standard deviation 𝛽_𝛼= 0.3. All other parameters are assumed to be deterministic. Figure 3 shows bifurcation behavior with the cubic stiffness as random, it now has a range of possible LCO amplitudes for each reduced velocity and the onset of flutter is unaffected. The standard deviation, that is, the amplitude variation range increases as reduced speed increases.

Figure 3 

Uncertain nonlinear stiffness: stochastic bifurcation diagram.

A Galerkin PCE approach is used to quantify the propagation of this uncertainty on the response. The Galerkin approach modifies the 8th-order flutter system to an 8×(𝑝+1) order system. It also involves calculating the complex fifth-order inner product terms as shown earlier. As a result, the solution process is computationally intensive for the nonlinear system in question. After solving for the chaos coefficients, in the post processing stage, the coefficients are substituted back to the expansion form to get the stochastic response. Probability density functions (PDFs) and other required statistics can then be readily obtained. The time histories of the first few random modes in pitch are plotted in Figure 4. The zeroth-order mode is the mean; one can also see that the contribution of higher-order random modes is gradually diminishing. A representative PDF is shown in Figure 5 for increasing order of chaos expansion terms. PDFs are calculated at time 𝑡=7000 at which the solutions are well past their transients and stationary. The reduced speed considered here is 𝑈=6.42, close to the deterministic bifurcation point. The figure also presents results from a standard MCS with 12000 samples as a comparison reference. One can see how increasing the order of expansion the CPU time for the solution is getting magnified. Results are presented up to the 12th order of expansion at which the PCE results match well with that of MCS. However, the simulation time also approaches to that of the reference MCS. While calculating the CPU time for the Galerkin-PCE approach the inner products computation and post processing of results are not taken into account.

Figure 4 

Galerkin-PCE: behavior of the first few random modes.

Figure 5 

Galerkin-PCE: PDF comparisons for increasing order of chaos expansions, at $t=7000$ with $U=6.42$.

Now the nonintrusive projection approach is followed using a Gauss-Hermite quadrature. Galerkin-PCE and nonintrusive results are compared in terms of their accuracy and simulation time in Figure 6. A good match with MCS is seen for the 12th order of expansion. Once again a standard MCS with 12000 samples is used as a reference solution (in other cases too we have used a standard 12000 samples MCS as reference). However, nonintrusive approach is seen to be much faster than Galerkin-PCE for the same level of accuracy as indicated in Figure 6. With this, the computational disadvantage of the conventional Galerkin based PCE for nonlinear systems is demonstrated. Henceforth, this approach will not be used for further simulations in this paper.

Figure 6 

Galerkin-PCE: PDF comparisons for intrusive and nonintrusive PCE, at $t=7000$ with $U=6.42$.

The response realization time histories for a few samples of random variable 𝜉 are plotted in Figure 7. The response time histories show difference in amplitude but not in phase. A typical realization time history obtained with the 12th order PCE along with its deterministic counterpart is compared in Figure 8. The match is perfect even at long time. Amplitude response PDFs as a function of reduced velocities (bifurcation parameter) are shown in Figure 9. They represent single peak monotonic behavior as all the realizations give finite amplitude LCOs. Effectively, the PDFs are not undergoing any qualitative change or bifurcations. Close to 𝑈=6.4 the PDF looks sharper and narrower as most realizations are going towards the same limit cycle amplitude. As the speed increases, the PDF is broader and less sharp, indicating that the realization amplitudes are spread over a wider band of amplitudes.

Figure 7 

Uncertain nonlinear stiffness: five different realizations time histories at $U=6.42$.

Figure 8 

Uncertain nonlinear stiffness: comparison of a typical time history with $12$th order PCE and MCS for $\xi =1.5$ and $U=6.42$.

Figure 9 

Uncertain nonlinear stiffness: amplitude response PDF as a function of reduced speed.

Next, we consider the viscous damping ratio in pitch (𝜁_𝛼) to be uncertain (in the earlier part damping was put to be zero) and all the other parameters deterministic. This case is potentially more interesting than the earlier one. The damping ratio is assumed to be a Gaussian random variable with mean 𝜁_𝛼 =0.1 and standard deviation 𝜁̃_𝛼 = 0.01. Figure 10 shows the bifurcation behavior with random damping ratio. The firm line gives the deterministic bifurcation behavior. The other bifurcation branches are for the two different extreme realizations of the random damping. Thus they represent the boundaries of the possible random variations of the bifurcation behavior (stochastic).

Figure 10 

Uncertain viscous damping: stochastic bifurcation diagram.

The major difference between the uncertain damping and the earlier considered uncertain stiffness is that, variation in damping can show phase shifting behavior in the response realizations. This is presented in Figure 11 where five different realizations are shifted in phase from each other. This behavior becomes more pronounced as time increases. As a result, response PDFs can now show bimodal behavior especially at large times. A few representative PDFs are plotted now at different reduced velocities. Figure 12 shows the response PDF at 𝑈=6.52 and time = 1400. Though the time level is past the transients it is not large and the phase shift is not yet very pronounced. The response PDF shows a single peak pattern. A reasonably good match with MCS is obtained within the 12th order of chaos expansion. However, at higher time levels as the phase shifting becomes stronger, the PDFs start to look distorted from their single peak behavior and goes towards a double peak pattern. Figure 13 shows the PDF at 𝑈=6.52 and time 5000 and 7800. In the first case, a double-peak bimodal PDF is just emerging as shown in Figure 13(a). In this case a 12th order expansion is not sufficient to capture the response accurately; a 15th order expansion gives better accuracy. At higher time 7800, the response PDF is more towards a two-peak bimodal shape as is seen in Figure 13(b). However, even a 15th order chaos expansion does not give the required accuracy. If one considers a different reduced velocity, the bimodal behavior can appear at some different time levels. The important observations from these figures are two-fold. First is the gradual double-peak behavior with increasing time. The second one is the apparent mismatch between the MCS and PCE results which seems to be increasing again with time. The reason for the first is nothing but the increasing phase shifting between the realizations time histories. However, for the second, the reason for the mismatch is the long time degeneracy which is shown in Figure 14. This mismatch can be improved by using higher order chaos expansions.

Figure 11 

Uncertain viscous damping: five different realizations time history at $U=6.52$.

Figure 12 

Uncertain viscous damping: comparison of the PDFs with increasing order of PCE at nondimensional time 1400 at $U=6.52$.

(a)

(b)

Figure 13 

Uncertain viscous damping: comparison of the PDFs with increasing order of PCE at $U=6.52$ with nondimensional time (a) $t=5000$, (b) $t=7800$.

(a)

(b)

Figure 14 

Uncertain viscous damping: a typical time history with (a) The $15$th order PCE and MCS at $\xi =2.3$ and reduced speed $U=6.52$. (b) A close-up of (a).

A typical realization time history with PCE along with its deterministic counterpart are presented in Figure 14. One can clearly see a degeneracy in the time history which starts around time levels close to 6000. PCE can show such type of degenerate behavior in capturing LCO response [17], especially at large times. As a counter measure, one can increase the order of the chaos expansion. However, this can only push the degeneracy to a later time but can not solve it entirely. Nonpolynomial based chaos approaches have been attempted in the recent past towards this end [17]. An unsteady adaptive stochastic finite elements method, developed by Witteveen and Bijl [34–36] has also been used successfully. This approach is based on time-independent parametrization. This achieves a constant accuracy in time with a constant number of samples. In this method interpolation of oscillatory samples is based on constant phase instead of a constant time.

The amplitude response PDFs for the uncertain damping case is shown in Figure 15 for different reduced velocities. Here the LCO amplitudes are captured after the initial transients have died down but before the time degeneracy has started. A nonmonotonic behavior is clearly indicated; some realizations are going to damped oscillation and others give LCO amplitudes scattered within the domain boundary. At 𝑈=6.5, the double-peak behavior of the PDF indicates the two different LCO amplitudes around which most of the realizations are concentrated. Towards 𝑈=6.6, all realizations give finite amplitude LCO, thus essentially they are of the same type. The PDF shows a single-peak monotonous behavior. Therefore, the PDFs of the response amplitude have clearly gone through a qualitative change here, in other words, a bifurcation.

Figure 15 

Uncertain viscous damping: amplitude response PDF as a function of reduced velocities.

For the uncertain damping case, we also see that the critical reduced velocity at which flutter can occur, has come down from its corresponding deterministic value. This value can be read off the bifurcation plot (Figure 10) as 𝑈=6.35. This is the lowest extrema of the critical points. The cumulative distribution function (CDF) and the PDF of the critical points are shown in Figure 16. The CDF can directly give the probability of flutter (in other words, probability of failure) at any given reduced velocity.

(a)

(b)

Figure 16 

Uncertain viscous damping: (a) CDF, (b) PDF of the critical flutter point.

5. Conclusions

The bifurcation behavior of a nonlinear pitch-plunge flutter problem with uncertain system parameters has been studied. The problem is a simple model problem to understand the mechanism of nonlinear flutter in a stochastic framework. The parameters which have been assumed to be random could attribute their uncertainties to laboratory testing conditions. Moreover, a cubic nonlinear stiffness is used for various sources of analytic nonlinearities; they often represent different control mechanisms and could face modeling uncertainties.

The classical Galerkin Polynomial Chaos method and the nonintrusive Projection method are applied to capture the propagation of uncertainty through the nonlinear aeroelastic system. The focus of this work is to investigate the performance of these techniques and to see how the aeroelastic stability characteristics are altered due to the random effects. The Monte Carlo solution is used as reference solution. The computational cost of the Galerkin Polynomial Chaos method is seen to be very high and subsequently only the Projection method based on Gauss-Hermite quadrature is used for the analysis. The effect of uncertain cubic structural nonlinearity and viscous damping parameter are investigated. Uncertainty in the cubic stiffness does not alter the bifurcation (flutter) point, it only affects the amplitudes of the periodic response in the post flutter stage. The PDF behavior also does not show any qualitative changes. On the other hand, uncertainty in damping affects the bifurcation point. It can lower the onset of flutter; the PDF of the response amplitude also undergoes a qualitative change. In other words, a bifurcation of the response PDF takes place. The results highlight the risk induced by parametric uncertainty and importance of uncertainty quantification in nonlinear aeroelastic systems.

The uncertain damping case by polynomial chaos suffer from long time degeneracy, as is also discussed in the literature. The degeneracy can be controlled by using higher order chaos expansions, though this cannot be a permanent solution. For the uncertain nonlinear stiffness, the problem of time degeneracy is not encountered.

Appendix

The coefficients introduced in Section 2 are used from [23] and are reproduced here for the sake of completeness: 𝑐₀1=1+𝜇,𝑐₁=𝑥_𝛼−𝑎_ℎ𝜇,𝑐₂=2𝜁_𝜀𝜔𝑈+2𝜇(1−𝜓₁−𝜓₂𝑐),₃=1+2(1/2−𝑎_ℎ)(1−𝜓₁−𝜓₂)𝜇,𝑐₄=(𝜔𝑈)²+2𝜇(𝜓₁𝜖₁+𝜓₂𝜖₂𝑐),₅=𝛽_𝜀(𝜔𝑈)²,𝑐₆=2𝜇[(1−𝜓₁−𝜓₂1)+(2−𝑎_ℎ)(𝜓₁𝜖₁+𝜓₂𝜖₂𝑐)],₇=2𝜇𝜓₁𝜖₁1[1−(2−𝑎_ℎ)𝜖₁],𝑐₈=2𝜇𝜓₂𝜖₂1[1−(2−𝑎_ℎ)𝜖₂𝑐],₉2=−𝜇𝜓₁𝜖²₁,𝑐₁₀2=−𝜇𝜓₂𝜖²₂,𝑑₀=𝑥_𝛼𝑟²_𝛼−𝑎_ℎ𝜇𝑟²_𝛼,𝑑₁=1+1+8𝑎²_ℎ8𝜇𝑟²_𝛼,𝑑₂𝜁=2^𝛼𝑈+1−2𝑎_ℎ2𝜇𝑟²_𝛼−(1+2𝑎_ℎ)(1−2𝑎_ℎ)(1−𝜓₁−𝜓₂)2𝜇𝑟²_𝛼𝑑₃=1𝑈²−(1+2𝑎_ℎ)(1−𝜓₁−𝜓₂)𝜇𝑟²_𝛼−(1+2𝑎_ℎ)(1−2𝑎_ℎ)(𝜓₁𝜖₁+𝜓₂𝜖₂)2𝜇𝑟²_𝛼,𝑑₄=𝛽_𝛼𝑈²,𝑑₅=−(1+2𝑎_ℎ)(1−𝜓₁−𝜓₂)𝜇𝑟²_𝛼,𝑑₆=−(1+2𝑎_ℎ)(𝜓₁𝜖₁+𝜓₂𝜖₂)𝜇𝑟²_𝛼,𝑑₇=−(1+2𝑎_ℎ)(𝜓₁𝜖₁[1−(1/2−𝑎_ℎ)𝜖₁])𝜇𝑟²_𝛼,𝑑₈=−(1+2𝑎_ℎ)(𝜓₂𝜖₂[1−(1/2−𝑎_ℎ)𝜖₂])𝜇𝑟²_𝛼,𝑑₉=(1+2𝑎_ℎ)𝜓₁𝜖²₁𝜇𝑟²_𝛼,𝑑₁₀=(1+2𝑎_ℎ)𝜓₂𝜖²₂𝜇𝑟²_𝛼,2𝑓(𝜏)=𝜇1((2−𝑎_ℎ)𝛼(0)+𝜀(0))(𝜓₁𝜖₁𝑒^−𝜖1𝜏+𝜓₂𝜖₂𝑒^−𝜖2𝜏),𝑔(𝜏)=−(1+2𝑎_ℎ)𝑓(𝜏)2𝑟²_𝛼. (A.1)