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

Authors: Ajit Desai; Sunetra Sarkar
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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.


---

## Body

## 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

Figure1 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)ε′′+xαα′′+2ζɛω¯Uε′+(ω¯U)2(ɛ+βɛɛ3)=-1πμCL(τ),xarα2ɛ′′+α′′+2ζαUα′+1U2(α+βαα3)=2πμrα2CM(τ).
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 ahb from the midchord where b is the half chord. Let us also use v as the wind velocity h as the plunge deflection. Among the nondimensional quantities, ɛ=h/b = nondimensional displacement of the elastic axis point; τ=vt/b = nondimensional time; U=v/(bωα)= 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, rα is the radius of gyration about the elastic axis, and μ is the airfoil mass ratio defined as m/(πρv). βα 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, CL(τ) and CM(τ):(2.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α′′
The Wagner function ϕ(τ) is given by:(2.3)ϕ(τ)=1-ψ1e-ϵ1τ-ψ2e-ϵ2τ.
Values for the constants are, ψ1= 0.165, ψ2= 0.335, ϵ1= 0.0455 and ϵ2 = 0.3 [24]. Introducing the following new variables w1,w2,w3,andw4 [23], the original integrodifferential equations for aeroelastic system given by (2.1) are reformulated: (2.4)w1=∫0τe-ϵ1(τ-σ)α(σ)dσ,w2=∫0τe-ϵ2(τ-σ)α(σ)dσ,w3=∫0τe-ϵ1(τ-σ)ɛ(σ)dσ,w4=∫0τe-ϵ2(τ-σ)ɛ(σ)dσ.
Now a set of autonomous differential equations of the form X′=f(X) are obtained as, X={x1,x2,x3,x4,x5,x6,x7,x8} = {α,α',ɛ,ɛ',w1,w2,w3,w4}.Figure 1
The schematic of a symmetric airfoil with pitch and plunge degrees-of-freedom.Explicitly, the system looks like,(2.5)x1′=x2,x2′=(c0N-d0M)(c1d0-c0d1),x3′=x4,x4′=(-c1N+d1M)(c1d0-c0d1),x5′=x1-ϵ1x5,x6′=x1-ϵ2x6,x7′=x3-ϵ1x7,x8′=x3-ϵ2x8,
where (2.6)M=c2x4+c3x2+c4x3+c5x33+c6x1+c7x5+c8x6+c9x7+c10x8-f(τ),N=d2x2+d3x1+d4x13+d5x4+d6x3+d7x5+d8x6+d9x7+d10x8-g(τ).
The coefficients c0⋯c10 and d0⋯d10 depend on the system parameters, and their expressions along with f(τ) and g(τ) 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 X(t,θ) (as function of random event θ) which is second-order stationary can be written as(3.1)X(t,θ)=â0ψ0+∑i1=1∞âi1ψ1(ξi1(θ))+∑i1=1∞∑i2=1i1âi1i2ψ2(ξi1(θ),ξi2(θ))+∑i1=1∞∑i2=1i1∑i3=1i2(âi1i2i3)ψ3(ξi1(θ),ξi2(θ),ξi3(θ))+⋯,
where ψn(ξi1,ξi2,…,ξin) denotes the Hermite polynomial of order n in terms of n-dimensional independent standard Gaussian random variables ξ=(ξi1,ξi2,…,ξin) 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(3.2)X(t,θ)=∑j=0∞aj(t)Φj(ξ(θ)).
There is an one-to-one relationship between the ψ's and Φ's and also âj's and aj's in (3.1) and (3.2). In the original form, chaos expansion uses Hermite polynomials (ψns). The form of the one-dimensional Hermite polynomials is given as follows.(3.3)ψ0=1,ψ1=ξ,ψ2=ξ2-1,ψ3=ξ3-3ξ,ψ4=ξ4-6ξ2+3,ψn=ξψn-1-(n-1)ψn-2.
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α(t,ξ(θ)), governed by Gaussian random variables ξ (ξ can always be normalized as a standard Gaussian one), can be approximated by the following truncated series:(3.4)α(t,ξ(θ))=∑j=0pα̂j(t)Φj(ξ(θ))
Note that, here the infinite upper limit of (3.2) is replaced by p, called the order of the expansion. For n number of random variable and polynomial order np, p is given by the following [26]:(3.5)p=(n+np)!n!np!-1.
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]:(3.6)£[α(t,θ)]+βα[α(t,θ)]3=0,
here £ is a linear differential operator.Equation (3.4) is now rewritten for a single random variable ξ as(3.7)α(t,ξ(θ))=∑j=0pα̂j(t)Φj(ξ(θ)).
Here Φj’s are now Hermite polynomials ψj 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(3.8)βα=β¯α+ξβ̃α=∑l=01βαlΦl,
with, βα0=β¯α and βα1=β̃α.Substituting the chaos expansion terms, (3.7) and (3.8) in (3.6),(3.9)£[∑j=0pα̂j(t)Φj]+∑l=01βαlΦl[∑j=0pα̂j(t)Φj]3=0.
The cubic nonlinear function can be expressed in the following form:(3.10)[∑j=0pα̂j(t)Φj]3=∑i=0p∑m=0p∑n=0pα̂iα̂mα̂nΦiΦmΦn.
Substituting (3.10) into (3.9) and simplifying, we get,(3.11)£[∑j=0pα̂j(t)Φj]+∑l=01βαlΦl[∑i=0p∑m=0p∑n=0pα̂iα̂mα̂nΦiΦmΦn]=0.
Using Galerkin projection on (3.11) by taking 〈·,Φk〉, for k=0,1,…,p,(3.12)£[α̂k(t)]+1〈Φk2〉[∑l=01∑i=0p‍∑m=0p‍∑n=0pβαlα̂iα̂mα̂n〈ΦlΦiΦmΦnΦk〉]=0.
The expected value operator 〈·〉, called the inner product, is defined as,(3.13)〈Φl⋯Φk〉=∫-∞∞Φl⋯Φkω(ξ)dξ.
For Hermite polynomials the weighting function ω(ξ) is the Gaussian probability density function. For single random variable case it is given as(3.14)ω(ξ)=(12π)e-((1/2)ξ2).
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 L2 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)〈ΦlΦk〉=〈Φl2〉δlk,δlk is the Kronecker delta function, given as:(3.16)δlk={1ifl=k,0otherwise.
The inner product terms in (3.12) 〈ΦlΦiΦmΦnΦk〉 and 〈Φk2〉 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〈Φ1Φ1〉=1〈Φ0Φ0Φ0Φ0Φ0〉=1〈Φ2Φ2〉=2〈Φ0Φ0Φ1Φ1Φ2〉=2〈Φ3Φ3〉=6〈Φ0Φ1Φ1Φ2Φ2〉=10〈Φ4Φ4〉=24〈Φ0Φ1Φ2Φ2Φ3〉=48The 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〈ΦiΦj〉=0 for i≠j, hence the expansion coefficients can be directly evaluated as(3.17)α̂j(t)=〈α(t,ξ(θ))Φj〉〈Φj2〉.
The denominator in (3.17) can be shown to satisfy 〈Φj2〉=j! for nonnormalized Hermite polynomials [31]. So the key step in projecting α(t,ξ(θ)) along the polynomial chaos basis is the evaluation of 〈α(t,ξ(θ))Φj〉. The inner product is given by the following integral:(3.18)〈α(t,ξ(θ))Φj〉=∫-∞∞α(t,ξ(θ))Φjω(ξ)dξ.
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 α(t,θ) are then used to estimate the deterministic coefficients, α̂j(t)’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.


## 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 X(t,θ) (as function of random event θ) which is second-order stationary can be written as(3.1)X(t,θ)=â0ψ0+∑i1=1∞âi1ψ1(ξi1(θ))+∑i1=1∞∑i2=1i1âi1i2ψ2(ξi1(θ),ξi2(θ))+∑i1=1∞∑i2=1i1∑i3=1i2(âi1i2i3)ψ3(ξi1(θ),ξi2(θ),ξi3(θ))+⋯,
where ψn(ξi1,ξi2,…,ξin) denotes the Hermite polynomial of order n in terms of n-dimensional independent standard Gaussian random variables ξ=(ξi1,ξi2,…,ξin) 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(3.2)X(t,θ)=∑j=0∞aj(t)Φj(ξ(θ)).
There is an one-to-one relationship between the ψ's and Φ's and also âj's and aj's in (3.1) and (3.2). In the original form, chaos expansion uses Hermite polynomials (ψns). The form of the one-dimensional Hermite polynomials is given as follows.(3.3)ψ0=1,ψ1=ξ,ψ2=ξ2-1,ψ3=ξ3-3ξ,ψ4=ξ4-6ξ2+3,ψn=ξψn-1-(n-1)ψn-2.
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α(t,ξ(θ)), governed by Gaussian random variables ξ (ξ can always be normalized as a standard Gaussian one), can be approximated by the following truncated series:(3.4)α(t,ξ(θ))=∑j=0pα̂j(t)Φj(ξ(θ))
Note that, here the infinite upper limit of (3.2) is replaced by p, called the order of the expansion. For n number of random variable and polynomial order np, p is given by the following [26]:(3.5)p=(n+np)!n!np!-1.
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]:(3.6)£[α(t,θ)]+βα[α(t,θ)]3=0,
here £ is a linear differential operator.Equation (3.4) is now rewritten for a single random variable ξ as(3.7)α(t,ξ(θ))=∑j=0pα̂j(t)Φj(ξ(θ)).
Here Φj’s are now Hermite polynomials ψj 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(3.8)βα=β¯α+ξβ̃α=∑l=01βαlΦl,
with, βα0=β¯α and βα1=β̃α.Substituting the chaos expansion terms, (3.7) and (3.8) in (3.6),(3.9)£[∑j=0pα̂j(t)Φj]+∑l=01βαlΦl[∑j=0pα̂j(t)Φj]3=0.
The cubic nonlinear function can be expressed in the following form:(3.10)[∑j=0pα̂j(t)Φj]3=∑i=0p∑m=0p∑n=0pα̂iα̂mα̂nΦiΦmΦn.
Substituting (3.10) into (3.9) and simplifying, we get,(3.11)£[∑j=0pα̂j(t)Φj]+∑l=01βαlΦl[∑i=0p∑m=0p∑n=0pα̂iα̂mα̂nΦiΦmΦn]=0.
Using Galerkin projection on (3.11) by taking 〈·,Φk〉, for k=0,1,…,p,(3.12)£[α̂k(t)]+1〈Φk2〉[∑l=01∑i=0p‍∑m=0p‍∑n=0pβαlα̂iα̂mα̂n〈ΦlΦiΦmΦnΦk〉]=0.
The expected value operator 〈·〉, called the inner product, is defined as,(3.13)〈Φl⋯Φk〉=∫-∞∞Φl⋯Φkω(ξ)dξ.
For Hermite polynomials the weighting function ω(ξ) is the Gaussian probability density function. For single random variable case it is given as(3.14)ω(ξ)=(12π)e-((1/2)ξ2).
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 L2 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)〈ΦlΦk〉=〈Φl2〉δlk,δlk is the Kronecker delta function, given as:(3.16)δlk={1ifl=k,0otherwise.
The inner product terms in (3.12) 〈ΦlΦiΦmΦnΦk〉 and 〈Φk2〉 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〈Φ1Φ1〉=1〈Φ0Φ0Φ0Φ0Φ0〉=1〈Φ2Φ2〉=2〈Φ0Φ0Φ1Φ1Φ2〉=2〈Φ3Φ3〉=6〈Φ0Φ1Φ1Φ2Φ2〉=10〈Φ4Φ4〉=24〈Φ0Φ1Φ2Φ2Φ3〉=48The 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〈ΦiΦj〉=0 for i≠j, hence the expansion coefficients can be directly evaluated as(3.17)α̂j(t)=〈α(t,ξ(θ))Φj〉〈Φj2〉.
The denominator in (3.17) can be shown to satisfy 〈Φj2〉=j! for nonnormalized Hermite polynomials [31]. So the key step in projecting α(t,ξ(θ)) along the polynomial chaos basis is the evaluation of 〈α(t,ξ(θ))Φj〉. The inner product is given by the following integral:(3.18)〈α(t,ξ(θ))Φj〉=∫-∞∞α(t,ξ(θ))Φjω(ξ)dξ.
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 α(t,θ) are then used to estimate the deterministic coefficients, α̂j(t)’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,ah=-0.5,xα=0.25,ζα=0,ζɛ=0,rα=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 the8th-order flutter system to an 8×(p+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 t=7000 at which the solutions are well past their transients and stationary. The reduced speed considered here is U=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, att=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 Figure6. 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, att=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 U=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 atU=6.42.Figure 8
Uncertain nonlinear stiffness: comparison of a typical time history with12th order PCE and MCS for ξ=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 Figure11 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 U=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 U=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 atU=6.52.Figure 12
Uncertain viscous damping: comparison of the PDFs with increasing order of PCE at nondimensional time 1400 atU=6.52.Uncertain viscous damping: comparison of the PDFs with increasing order of PCE atU=6.52 with nondimensional time (a) t=5000, (b) t=7800.
(a)(b)Uncertain viscous damping: a typical time history with (a) The15th order PCE and MCS at ξ=2.3 and reduced speed U=6.52. (b) A close-up of (a).
(a)(b)A typical realization time history with PCE along with its deterministic counterpart are presented in Figure14. 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 Figure15 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 U=6.5, the double-peak behavior of the PDF indicates the two different LCO amplitudes around which most of the realizations are concentrated. Towards U=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 (Figure10) as U=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.Uncertain viscous damping: (a) CDF, (b) PDF of the critical flutter point.
(a)(b)


## 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.



---

*Source: 379472-2010-07-18.xml*