AIAA 2007-317

45th AIAA Aerospace Sciences Meeting and Exhibit 8 - 11 January 2007, Reno, Nevada

Probabilistic Collocation: An Efficient Non-Intrusive Approach For Arbitrarily Distributed Parametric Uncertainties G.J.A. Loeven,∗†J.A.S. Witteveen,∗ and H. Bijl‡ Delft University of Technology, The Netherlands

Complex flow and fluid-structure interaction simulations require efficient uncertainty quantification methods. In addition, a good property of an uncertainty quantification method is non-intrusiveness, meaning that existing deterministic solvers can be used for uncertainty quantification. In this paper the Probabilistic Collocation method is introduced, which combines the non-intrusiveness of the Stochastic Collocation method with the exponential convergence for arbitrary probability distributions of the Galerkin Polynomial Chaos method. Due to the non-intrusiveness and exponential convergence, the Probabilistic Collocation method requires only a few collocation points for a high accuracy. For the one dimensional piston problem the efficiency of the Probabilistic Collocation method is compared with the Galerkin Polynomial Chaos method, the Non-Intrusive Polynomial Chaos method and the Stochastic Collocation method. The strength of the Probabilistic Collocation method is demonstrated by solving steady flow around a NACA0012 airfoil with an uncertain free stream velocity using a commercial flow solver. Different possibilities of representing the stochastic solution are demonstrated to show the potential use of uncertainty quantification.

I.

Introduction

Increasing computer facilities and more advanced algorithms lead to more accurate solutions. The solutions appear physically correct, however, only little is known with respect to their reliability. Physical uncertainties, which are always present, can lead to significant variations in the results. Neglecting these inherent uncertainties is not acceptable with the current state of technology. A trade-off has to be made between accuracy of the deterministic solution and the inclusion of uncertainties. This can be achieved by dividing the computational budget in a deterministic part and a stochastic part. When a system is known to incorporate uncertainties, it seems logical to do several less accurate deterministic runs to obtain information about the uncertainty propagation and therefore the reliability of the solution. However, this not common practice in engineering yet. Real life computations are often computationally intensive. This holds especially for computational fluid dynamics problems due to the large numbers of degrees of freedom. Uncertainty quantification requires an additional method to obtain the stochastic properties, this leads to an increase in the computational work compared to the deterministic case. Efficiency of the uncertainty quantification method is, therefore, of great importance. Several methods exist to quantify the uncertainties. A straightforward method is Monte Carlo simulation, which is often too expensive, since a lot of samples are required for reasonable accuracy. A commonly used more efficient method is the Galerkin Polynomial Chaos method, which shows exponential convergence with respect to the polynomial order. The method is based on the Homogeneous Chaos theory of Wiener 1 and was proposed by Ghanem and Spanos.2 To find the polynomial coefficients a Galerkin projection is used, ∗ Ph.

D. researcher, Faculty of Aerospace Engineering, P.O. Box 5058, 2600 GB Delft, The Netherlands, Member AIAA [email protected] ‡ Full Professor, Faculty of Aerospace Engineering, P.O. Box 5058, 2600 GB Delft, The Netherlands, Member AIAA. † Email:

1 of 14 American Institute of Aeronautics and Astronautics

Copyright © 2007 by G.J.A. Loeven, J.A.S. Witteveen, and H. Bijl. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

which results in a coupled system of deterministic equations. The method is therefore intrusive, meaning the existing code has to be edited to be able to solve the coupled system. For practical applications a non-intrusive method is preferred, since it enables engineers to use their existing deterministic solver as a black box for uncertainty quantification. That is why in practice Monte Carlo simulation is still extensively applied. To benefit from the efficiency of the Galerkin Polynomial Chaos method, Walters, Hosder, and Perez3, 4 developed the Non-Intrusive Polynomial Chaos method and applied it to computational fluid dynamics problems. A similar approach is followed by Reagan, Najm, Ghanem, and Knio5 called Non-Intrusive Spectral Projection. Both methods use samples to estimate the coefficients of the polynomial chaos expansion. Another non-intrusive method is the Stochastic Collocation method, which was proposed by Mathelin and Hussaini.6, 7 The probability distribution of the uncertain parameter forms the basis of a transformation to an artificial space. Collocation points are chosen in this artificial space and with Lagrange interpolation the probability distribution of the solution is constructed. Mathelin and Hussaini7 showed a significant decrease in computational time compared to the Galerkin Polynomial Chaos method for a quasi-1D nozzle flow. A previous study8 showed that the Stochastic Collocation method only shows exponential convergence for uniformly distributed uncertain parameters. The Stochastic Collocation method results in a set of decoupled deterministic equations, making the method non-intrusive. In this paper the Probabilistic Collocation method is presented. This method combines the idea of a chaos transformation like the Polynomial Chaos method and the collocational approach of the Stochastic Collocation method. The Probabilistic Collocation method presented in this paper is a generalization of the Stochastic Collocation method. Since the optimal collocation points with corresponding weights are computed based on the input distribution, exponential convergence is obtained with respect to the polynomial order for arbitrary probability distributions. Babuˇska, Nobile, and Tempone9 proved that for elliptic partial differential equations the Probabilistic Collocation method is equivalent to the Galerkin Polynomial Chaos method. The basis of the Probabilistic Collocation method is formed by Lagrange interpolation and Gaussian quadrature, with the same chaos transformation as the Polynomial Chaos method. The Probabilistic Collocation method computes the Galerkin projection numerically using Gaussian quadrature, resulting in a decoupled system of equations. The method is therefore non-intrusive. A major advantage is that the Probabilistic Collocation method can be used for all kinds of non-linearities whereas the Galerkin based Polynomial Chaos methods show difficulties.7, 10, 11 Numerical results are shown for the one dimensional piston problem. A comparison is made on efficiency between the mentioned methods. A problem of the collocation methods is the fact that for multiple uncertain parameters the collocation points are obtained using the tensor products of univariate polynomials. This leads to a rapid increase of the number of collocation points. Xiu and Hesthaven12 have investigated sparse collocation grid approaches to decrease the number of points for multivariate approximations. To demonstrate how uncertainties can efficiently be quantified using the Probabilistic Collocation method combined with a commercial flow solver, the method is applied to steady flow around a NACA0012 airfoil with uncertain free stream velocity. The flow is solved using the FineTM/Hexa flow solver of Numeca Int. Three possibilities of representing the stochastic solution are demonstrated to show the potential use of uncertainty quantification. First the probability density functions of quantities of interest like drag or lift are shown. Secondly, the mean and variance pressure field are computed to see in what part of the flow domain the pressure is most sensitive to an uncertain free stream velocity. And thirdly, the mean skin friction and pressure on the airfoil surface are plotted with 99% uncertaintybars, which are obtained from the distribution function of the solutions. The paper is organized as follows: section II reviews some efficient non-intrusive methods for uncertainty quantification. Section III introduces the Probabilistic Collocation method, which is compared on efficiency using the linear piston problem with existing methods in section IV. In section V the strength of the Probabilistic Collocation method is demonstrated for steady flow around a NACA0012 airfoil, using the FineTM/Hexa flow solver as a black box. The conclusions are drawn in section VI.

II.

Several Efficient Uncertainty Quantification Methods

Different methods exist to quantify uncertainties. Monte Carlo simulation is a natural way to quantify uncertainties, however, a lot of samples are required to obtain reasonable accuracy. In this paper the focus is on methods that increase the efficiency by approximating the stochastic response of the solution based on

2 of 14 American Institute of Aeronautics and Astronautics

the input distribution of the uncertain parameters. Next to efficiency, non-intrusiveness is considered to be the second important property, since advanced existing solvers can be used as a black box for uncertainty quantification. Several efficient non-intrusive methods exist, among others the Non-Intrusive Polynomial Chaos method3–5 and the Stochastic Collocation method.6, 7 In this section the methods are explained using the following general stochastic differential equation L (a (ω)) u(x, t, ω) = S(x, t),

(1)

where u(x, t, ω) is the solution and L (a (ω)) is a differential operator which contains space and time differentiation and can be nonlinear and depends on a random parameter a (ω). The solution u(x, t, ω) is a function of space x ∈ D ⊂ Rd , time t and the random event ω ∈ Ω. The complete probability space is then given by (Ω, F, P ), with Ω the set of outcomes, F ⊂ 2Ω the σ-algebra of events and P : F → [0, 1] a probability measure. S(x, t) is a space and time dependent source term. The random event ω is introduced by the presence of an uncertain parameter a(ω). A.

Galerkin and Non-Intrusive Polynomial Chaos

The Polynomial Chaos methods results in a spectral representation of the stochastic response of the solution and high order approximations of the mean and variance. Based on the Homogeneous Chaos theory of Wiener1 the original Galerkin Polynomial Chaos method was developed by Ghanem and Spanos. 2 The Polynomial Chaos framework has been extended to obtain exponential convergence for arbitrary distributions 13–15 using numerical techniques to construct an optimal set of orthogonal polynomials. It has successfully been applied to fluid mechanics by Walters and Huyse.16 An advantage of the Polynomial Chaos method is the exponential convergence with respect to the polynomial order, a disadvantage is the intrusiveness due to the coupled system of equations that has to be solved. A lot of research in Non-Intrusive Polynomial Chaos 3–5 methods is done to be able to use existing solvers for uncertainty quantification. The Polynomial Chaos methods starts with the polynomial chaos expansion of the solution u(x, t, ω) and all uncertain parameters, with the infinite sum truncated at a finite number of M + 1 terms u(x, t, ω) ≈

M X

ui (x, t)Ψi (ξ(ω)).

(2)

i=0

This expansion is a spectral expansion in the vector of random variables ξ(ω) with the random trial basis {Ψi }. Equation (2) divides the random variable u(x, t, ω) into a deterministic part, the coefficients u i (x, t) and a stochastic part, the polynomial chaoses Ψi (ξ(ω)). The total number of expansion terms is M + 1, which is determined by the dimension n of the vector of random variables ξ(ω) (i.e. the number of uncertain parameters) and the highest order p of the polynomials {Ψi } M +1=

(n + p)! . n!p!

(3)

The basis {Ψi } is a set of polynomials orthogonal with respect to the probability density function of the input uncertainty. For some standard distributions like the normal or uniform distribution, Hermite or Legendre polynomials are used.13 For other commonly used distributions like the lognormal distribution, the orthogonal polynomials are constructed numerically using the Gram-Schmidt algorithm. 15 Substituting the polynomial chaos expansion (2) into the differential equation (1) results in L (a (ω))

M X i=0

ui (x, t)Ψi (ξ(ω)) ≈ S(x, t).

(4)

Here, the parameter a has to be expanded as well, however for convenience in this notation it is omitted. A Galerkin projection on each basis {Ψk } is applied, to ensure that the truncation error is orthogonal to the functional space spanned by {Ψi }: + * M X k = 0, 1, . . . , M. (5) L (a (ω)) ui (x, t)Ψi , Ψk = hS, Ψk i , i=0

3 of 14 American Institute of Aeronautics and Astronautics

This deterministic set of M + 1 equations is coupled and can be solved using standard numerical techniques. Here a Block-Gauss-Seidel algorithm is employed to obtain the coefficients of the expansion. From Eq. (5) the coefficients ui (x, t) are computed. The idea of the non-intrusive approaches3–5 is to estimate the coefficients of the polynomial chaos expansion based on few deterministic solves. The Non-Intrusive Polynomial Chaos3, 4 takes for a polynomial chaos of order M a set of M + 1 vectors ξ i = ξ1 , ξ2 , . . . , ξn i for i = 0, 1, 2, . . . , M in the random space. For one normally distributed parameter (n = 1) the vector of a 4th order approximation is given by ξ = 0, 1, −1, 2, −2. The Non-Intrusive Spectral Projection5 uses Latin Hypercube sampling to estimate the coefficients. For each of these samples the deterministic code is run. The polynomial coefficients of expansion (2) are obtained by solving the following linear system:      Ψ0 (ξ 0 ) Ψ1 (ξ 0 ) . . . ΨM (ξ 0 ) u0 (x, t) u(x, t, ξ 0 )       Ψ0 (ξ 1 ) Ψ0 (ξ 1 ) . . . ΨM (ξ 1 )   u1 (x, t)   u(x, t, ξ 1 )       . (6) .. .. .. .. .. ..   =  .      . . . . . Ψ0 (ξ M ) Ψ0 (ξ M ) . . . ΨM (ξ M ) uM (x, t) u(x, t, ξM )

The intrusive Galerkin approach results in the exact polynomial coefficients ui (x, t), while the nonintrusive approach yields approximations of the polynomial coefficients ui (x, t). The stochastic solution u(x, t, ω) is reconstructed using Eq. (2). The mean µu and the variance σu2 of the solution are determined using: µu = u 0 , σu2 =

M X i=1

(7)

 ui (x, t)2 Ψ2i .

(8)

These expressions follow from the definition of the mean and variance. B.

Stochastic Collocation

A different spectral approach is the Stochastic Collocation method. It was developed by Mathelin and Hussaini6 to reduce the costs of the Galerkin Polynomial Chaos method in case of nonlinear equations. The method was successfully applied to a quasi-1D nozzle with uncertain inlet conditions and nozzle shape. Mathelin and Hussaini showed a substantial decrease of computing time using the Stochastic Collocation method compared to the Polynomial Chaos method. For each Legendre collocation point the problem is solved deterministically. The distribution function of the solution is reconstructed using polynomial interpolation. The mean and variance are obtained using Gauss-Legendre quadrature. The Stochastic collocation method is explained using the general differential equation (1). In case of the Stochastic Collocation method the distribution function of the random variable is projected from [0, 1] on the domain [−1, 1], which is called the α domain, by the linear transformation Fu (u) = 2Fu (u) − 1

,

(9)

where Fu (u) is the projected distribution function on the α domain [−1, 1] and Fu (u) the distribution function on [0, 1]. The solution in the α domain is u(x, t, α). The α domain is a stochastic space which is defined according to a standard domain of orthogonal polynomials [−1, 1]. From the α domain N p collocation points αi are taken. The method proposed by Mathelin and Hussaini6 uses Np Gauss-Legendre quadrature points and Lagrange interpolating polynomials of order Np − 1 for the function approximation. The solution u(x, t, α) is approximated by the following expansion u(x, t, α) ≈

Np X

ui (x, t)hi (α) ,

(10)

i=1

with ui (x, t) the values of u(x, t, α) at the collocation points αi and hi (α) interpolating polynomials of degree Np − 1, with hi (αj ) = δij . Transformation (9) is applied to the general differential equation (1), after

4 of 14 American Institute of Aeronautics and Astronautics

which expansion (10) is substituted. A Galerkin projection on each basis {hl } is applied to make the error orthogonal to the functional space spanned by {hi }: *

L (a (ω))

Np X

ui (x, t)hi , hl

i=1

+

= hS, hl i ,

l = 1, . . . , Np .

(11)

The Galerkin projection (11) is approximated using Gaussian quadrature. For a general inner product hf (α), g(α)i of two functions f (α) and g(α) Gaussian quadrature results in: hf (α), g(α)i = =

Np Np Np XX X

fi gj hi (αk )hj (αk )wk ,

i=1 j=1 l=1

Np Np Np XX X

fi gj δik δjk wk ,

i=1 j=1 k=1

=

Np X

fk gk wk ,

(12)

k=1

where wk are the quadrature weights corresponding to the collocation points αk . The resulting set of equations is fully decoupled. The mean µu and the variance σu2 of the stochastic solution can be determined using: µu =

Np X

1 2 ui (x, t)wi ,

(13)

i=1

σu2 =

Np X i=1

2  Np X 2 1 1 (ui (x, t)) wi −  ui (x, t)wi  , 2

2

(14)

i=1

where wi are the weights corresponding to the collocation points αi . These relations are derived from the definition of the mean and variance. The most important difference of the Stochastic Collocation method with the Polynomial Chaos method is the transformation to the artificial space α. The polynomial expansion exists in the α space and is, therefore, not a polynomial chaos expansion.

III.

The Probabilistic Collocation Method

The Probabilistic Collocation method is based on the idea of a chaos transformation which is also used in the Polynomial Chaos methods. Exponential convergence of the Probabilistic Collocation is shown numerically in section IV. In the Probabilistic Collocation method first a chaos version of Lagrange interpolation is used for the distribution function approximation. With this special Lagrange interpolation it is possible to describe the input distribution with two collocation points only. Just as for the Polynomial Chaos method, where two polynomial coefficients describe the input distribution. Secondly, Gaussian quadrature is used to compute the Galerkin projection and the integration of the approximation of the distribution function. Like ordinary Gaussian quadrature using Np points is able to integrate polynomials up to degree 2Np − 1 exactly, the Gaussian quadrature presented here integrates polynomial chaoses exactly. By using Gaussian quadrature a decoupled set of equations and a high order approximation of the mean and variance are obtained. Probabilistic Collocation expansion with chaos Lagrange interpolation Each variable depending on the uncertain input parameter is expanded as follows: u(x, t, ω) =

Np X

ui (x, t)hi (ξ(ω)) ,

i=1

5 of 14 American Institute of Aeronautics and Astronautics

(15)

where ui (x, t) is the solution u(x, t, ω) at the collocation point ωi and hi is the Lagrange interpolating polynomial chaos corresponding to the collocation points ωi . Np is the number of collocation points. Np − 1 is the order of the interpolating polynomial chaos. The difference with ordinary Lagrange interpolation is that here a polynomial chaos is used instead of ordinary polynomials. The Lagrange interpolating polynomial is a function in terms of the random variable ξ(ω). The random variable ξ(ω) is chosen on the standard domains [−1, 1], [0, ∞) or (−∞, ∞), such that the uncertain input parameter a(ω) is a linear transformation of ξ(ω): a(ω1 )ξ(ω0 ) − a(ω0 )ξ(ω1 ) a(ω0 ) − a(ω1 ) a(ω) = + ξ(ω) = a ˜0 + a ˜1 ξ(ω), (16) ξ(ω0 ) − ξ(ω1 ) ξ(ω0 ) − ξ(ω1 ) which follows from the expansion of a(ω) similar to equation (15). The Lagrange interpolating polynomial chaos is the polynomial chaos hi (ξ(ω)) of order Np − 1 that passes through the Np collocation points. It is given by: Np Y ξ(ω) − ξ(ωk ) hi (ξ(ω)) = , (17) ξ(ωi ) − ξ(ωk ) k=1 k6=j

with hi (ξ(ωj )) = δij . The collocation points are chosen such that they correspond to the Gaussian quadrature points used to integrate the function u(x, t, ω) in the ω domain. The solution has to be integrated in order to obtain for instance the mean or variance. To find the Gaussian quadrature points and weights the procedure below is followed. Computing Gaussian quadrature points with corresponding weights A powerful method to compute Gaussian quadrature rules is by means of the Golub-Welsch algorithm. 17 This algorithm requires the recurrence coefficients18 of polynomials which are orthogonal with respect to the weighting function of the integration. Exponential convergence for arbitrary probability distributions is obtained when the polynomials are orthogonal with respect to the probability density function of ξ, so w(ξ) = fξ (ξ). The recurrence coefficients are computed using the discretized Stieltjes procedure, 19 which is a stable method for arbitrary distribution functions. For convenience of notation the argument ω is omitted from here on. First the recurrence coefficients have to be computed. Orthogonal polynomials satisfy the three-term recurrence relation: πi+1 (ξ) = (ξ − αi )πi (ξ) − βi πi−1

i = 2, 3, . . . , Np ,

π0 (ξ) = 0, π1 (ξ) = 1,

(18) N

p where αi and βi are the recurrence coefficients determined by the weighting function w(ξ) and {πi (ξ)}i=1 is a i i−1 set of (monic) orthogonal polynomials with πi (ξ) = ξ + O(ξ ), i = 1, 2, . . . , Np . The recurrence coefficients are given by the Darboux’s formulae:18

(ξπi , πi ) (πi , πi ) (πi , πi ) βi = (πi−1 , πi−1 )

i = 1, 2, . . . , Np ,

αi =

i = 2, 3, . . . , Np ,

(19)

where (·, ·) denotes an inner product. The first coefficient β1 is given by (π1 , π1 ). Gander and Karp19 showed that discretizing the weighting function leads to a stable algorithm. Therefore the discretized Stieltjes procedure is used to obtain the recurrence coefficients. Hereto, first the weighting w(ξ) is discretized by wN (ξ) =

N X j=1

wj δ(ξ − ξj ),

wj > 0,

(20)

where δ is the Dirac delta function. Stieltjes’ procedure starts with i = 1. With (19) the first coefficient P α1 is computed, β1 = N j=1 wj . Now π2 (ξ) is computed by (18) using α1 and β1 . This is repeated for 6 of 14 American Institute of Aeronautics and Astronautics

i = 2, 3, . . . , Np . When continuous weighting functions are considered Np  N , for discrete measures Np ≤ N . The inner product is defined as (p(ξ), q(ξ)) =

Z

p(ξ)q(ξ)wN (ξ)dξ = S

N X

wj p(ξj )q(ξj ),

(21)

j=1

for two functions p(ξ) and q(ξ). From the recurrence coefficients αi and βi , i = 1, 2, . . . , Np , the collocation points ξi and corresponding weights wi are computed using the Golub-Welsch algorithm.17 With the recurrence coefficients the following matrix is constructed:   √ α1 β2 √   √   β2 α 2 β3 ∅   √ √   β α β 3 3 4   (22) J = . . . . .. .. ..     p p   ∅ βNp −1 αNp −1 βNp   p βNp αN p

The eigenvalues of J are the collocation points ξi , i = 1, . . . , Np , which are the roots of the polynomial of 2 order Np from the set of the constructed orthogonal polynomials. The weights are found by w i = β1 v1,i , i = 1, . . . , Np , where v1,i is the first component of the normalized eigenvector corresponding to eigenvalue ξ i . Now the collocation points ξi in the ξ-domain are known. They are mapped to the ω-domain using the distribution function of ξ. The optimal collocation points ωi are found by ωi = Fξ (ξi ),

i = 1, . . . , Np .

(23)

For a uniformly distributed parameter the Probabilistic Collocation method results in the Stochastic Collocation method,6 if ξ(ω) is chosen to be U (−1, 1) distributed. Application to a general differential equation Expansion (15) is substituted into the general differential equation (1). A Galerkin projection on each basis {hk (ξ(ω))} is applied: * + Np X L (a (ω)) ui (x, t)hi , hk = hS, hk i , k = 1, . . . , Np . (24) i=1

This projection is approximated using Gaussian quadrature, with optimal collocation points and corresponding weights based on the input distribution. The result is a deterministic system of equations which is fully decoupled. The mean and variance of the solution are found by µu =

Np X

ui (x, t)wi ,

(25)

i=1

σu2 =

Np X i=1

2



(ui (x, t)) wi − 

Np X i=1

2

ui (x, t)wi  ,

(26)

where wi are the weights corresponding to the collocation points ωi . These relations are derived from the definition of the mean and variance. For example if an input parameter is normally distributed and the solution is approximated using a 4th order expansion. Table 1 shows the five collocation points and weights corresponding to a 4 th order Probabilistic Collocation expansion for a normally distributed random variable. They are computed using the method described above where ξ is a N (0, 1) distributed random variable. For an uncertain parameter a with a N (µ, σ 2 ) distribution the deterministic solver has to be run for the following input values ai = Fa−1 (ωi ),

i = 1, . . . , 5.

7 of 14 American Institute of Aeronautics and Astronautics

(27)

PSfrag replacements where ωi the collocation points in Table 1. Evaluating the problem deterministically at the collocation points results in five solutions ui . From these five solutions the distribution function can be reconstructed using expansion (15). The mean and variance are obtained using equations (25) and (26). Table 1. Five collocation points and weights corresponding to a 4th order Probabilistic Collocation expansion for a normally distributed random variable.

Collocation point ω1 0.00213853121130 ω2 0.08760906885846 ω3 0.50000000000000 ω4 0.91239093114154 ω5 0.99786146878870

Weight 0.01125741132772 0.22207592200562 0.53333333333333 0.22207592200562 0.01125741132772

w1 w2 w3 w4 w5

A disadvantage of the presented Probabilistic Collocation method compared to the Galerkin Polynomial Chaos method is that the number of required collocation points increases rapidly with the number of uncertain parameters. When the number of uncertain parameters increases the collocation points are generally found using tensor products of one-dimensional polynomials. Improvements on this area have been performed by Xiu and Hesthaven12 for the Stochastic Collocation method. An extension to the Probabilistic Collocation method has to be investigated.

IV.

Efficiency of the Probabilistic Collocation Method

rag replacements L

q(t),x























































 

























 





















































 





































































































p



 Fluid

m

k

amb

Piston position q

The test problem is the linear piston problem8, 20 indicated in Figure 1. The linear piston is chosen since it is a test problem which is easily evaluated, but with a nonlinear dependence between the solution and the uncertain parameter. It consists of a one-dimensional fluid domain on one side bounded by a mass 1

k − 10% k = 1.429 k + 10%

0.5 0

-0.5 -1 0

5

10

15

20

25

time t Figure 2. Realizations using different values for the spring stiffness k of the linear piston problem for k = µk = 1.429 and k = µk ± 10%.

Figure 1. Linear piston problem.

which is attached to a spring. The fluid is modeled using the linearized isentropic Euler equations. The piston position q(t) and velocity q(t) ˙ are determined by the mass of the piston m, the spring stiffness k, the length of the fluid domain L, the fluid state U = (ρ, ρu)T , the ambient pressure pamb , and the initial conditions. At t = 0 the fluid is at rest and the piston state is q(0) = 1 and q(0) ˙ = 0. A second-order central finite volume discretization is used for the flow. The time integration is performed using an ESDIRK-4 scheme.20 The complete equations can be found in Ref. 8. The mass of the piston is taken m = 2, which together with the mean value for k leads to a interesting interaction between fluid and structure. The mean value for k is µk = 1.429 and corresponding to U (µk ± 10%) for a uniform distribution, the variance is set to σk2 = 6.8068 · 10−3 . Figure 2 shows that for a 10% variation of k a significant change of the solution is present. The solution is considered the piston position q at t = 10. The Probabilistic Collocation method is compared with the Galerkin and Non-Intrusive Polynomial Chaos methods and the Stochastic Collocation method. The error convergence with respect to the amount of computational work is used to compare the efficiency. The work is expressed as the number of times a deterministic problem needs to be solved. For the Galerkin Polynomial Chaos method the amount of 8 of 14 American Institute of Aeronautics and Astronautics

work is I ∗ (M + 1) where M is the number of polynomial chaos expansion terms of the approximation and I the number of Block-Gauss-Seidel iterations required to obtain the polynomial coefficients. The BlockGauss-Seidel algorithm makes use of the structure of the resulting matrix. The stopping criterion for the iterations is a residual of 10−8 . For the computations in this paper about 3–5 iterations were required. The Non-Intrusive Polynomial Chaos method requires M + 1 deterministic solves to approximate the coefficients. The Stochastic and Probabilistic Collocation methods both use one deterministic computation for every collocation point. The errors are defined by σ2 − σ2 µu − µu,ref u u,ref (28) ε µu = and εσu2 = , 2 σu,ref µu,ref

where the subscript u,ref indicates the reference solution. The reference solution is a 20 th order polynomial chaos computation. The method was validated for lower orders using a Monte Carlo simulation with 10 6 samples. The convergence of the methods is much faster than Monte Carlo simulation, therefore very quickly the accuracy is below the accuracy of the Monte Carlo simulation. Probabilistic Collocation versus Galerkin Polynomial Chaos

First the Probabilistic Collocation method is compared with the Galerkin Polynomial Chaos method to demonstrate the exponential convergence with respect to the polynomial chaos order and the increase in efficiency due to the decoupling of the equations. Figure 3 shows the convergence of both methods, first with respect to polynomial chaos order (cf. Figure 3(a)). The convergence is the same for both methods with respect to polynomial chaos order, this holds for arbitrarily distributed uncertain parameters. Figure 3(b) shows the convergence with respect to the amount of computational work, expressed in the number of deterministic solves. When the computational work is considered, the strength of the Probabilistic Collocation approach comes forward. The amount of work is equal to the number of collocation points, while for the rag replacementsGalerkin Polynomial Chaos method a coupled system of equations has to be solved. This is, however, a very efficient case for Galerkin Polynomial Chaos, since the test problem does for example not contain any non-polynomial nonlinearities. For problems nonlinearities (especially non-polynomial) involved in the PSfragwith replacements equations, the Galerkin approach can become very inefficient.11 The Probabilistic Collocation method has no problems dealing with nonlinear equations, since the problem is solved deterministically for each collocation point. The deterministic problems can be implemented as a black box.

PCM: PCM: GPC: GPC:

100

mean var mean var

10−5

Error

Error

100

10−10

PCM: PCM: GPC: GPC:

mean var mean var

10−5

10−10

0

1 2 3 4 Polynomial chaos order p

5

6

5

10

15

20

25

Work

(a)

(b)

Figure 3. Convergence for the Probabilistic Collocation method (PCM) and the Galerkin Polynomial Chaos method (GPC), with respect to the polynomial chaos order(a) and the computational work(b).

To show that the Probabilistic Collocation method results in exactly the same polynomial chaos approx-

9 of 14 American Institute of Aeronautics and Astronautics

imation as the Galerkin Polynomial Chaos, the monomial coefficients of the expansions are compared. Np X

u(ωi )hi (ξ(ω)) =

M X

ui Φi (ξ(ω)) = u ˆ0 + u ˆ1 ξ(ω) + u ˆ2 ξ(ω)2 + . . . + u ˆM ξ(ω)M

(29)

i=0

i=1

Table 2 shows the monomial coefficients of the expansion for the piston position at t = 10. To compute the coefficients of the polynomial chaos method exactly the system of equations is solved directly. For general use of the Galerkin Polynomial Chaos the system is not solved directly, then the Block-Gauss-Seidel algorithm is used. For elliptic partial differential equations Babuˇska, Nobile, and Tempone9 have proved that the Galerkin Polynomial Chaos method and the Probabilistic Collocation method are equivalent. Table 2. The monomial coefficients of the 4th order expansion (29) of the piston position at t = 10 obtained from the Galerkin Polynomial Chaos method and the Probabilistic Collocation method.

Coefficient u ˆ0 u ˆ1 u ˆ2 u ˆ3 u ˆ4

Galerkin Polynomial Chaos -0.73204167739370 0.11007272407613 0.00969183126072 -0.00069309123364 -0.00000788645601

Probabilistic Collocation -0.73204167739370 0.11007272407613 0.00969183126072 -0.00069309123364 -0.00000788645601

Non-Intrusive Polynomial Chaos -0.73204167739370 0.11008631168539 0.00969155246323 -0.00069926645535 -0.00000775982742

Probabilistic Collocation versus Non-Intrusive Polynomial Chaos The Non-Intrusive Polynomial Chaos method results in approximations of the polynomial chaos coefficients. The accuracy with respect to the amount of computational work cannot be better than the Probabilistic Collocation method , since this method yields the exact polynomial chaos coefficients. The approximated monomial coefficients for the piston position q at t = 10 are included in Table 2. The first coefficient is accurately approximated, the other coefficients show some discrepancy. The influence of these small differences between the exact and approximated coefficients on the convergence is clearly shown in Figure 4. PSfrag replacements It must be stated that the convergence depends heavily on the choice of the sampling vectors ξ i . For this comparison the choice of Walters, Hosder, and Perez3, 4 is followed. For the case of a normally distributed spring stiffness k the sampling vector is ξ = 0, 1, −1, 2, −2. Reagan, Najm, Ghanem, and Knio 5 address the importance of the choice of the sampling vector as well. They are not convinced that M + 1 samples are enough to accurately estimate the expansion coefficients. They, therefore, use many more samples computed using Latin Hypercube sampling.

Error

100

PCM: mean PCM: var NIPC: mean NIPC: var

10−5

10−10

1

2

3

4 Work

5

6

7

Figure 4. Convergence for the Probabilistic Collocation method (PCM) and the Non-Intrusive Polynomial Chaos method (NIPC) with respect to the computational work.

10 of 14 American Institute of Aeronautics and Astronautics

Probabilistic Collocation versus Stochastic Collocation

rag replacementsFigure 5 shows the convergence with respect to the amount of computational work, which is for both methods equal to the number of collocation points. Figure 5(a) shows that when ξ(ω) is taken U (−1, 1) the PSfrag replacements Probabilistic Collocation is equivalent to the Stochastic Collocation method. The Probabilistic Collocation method is a generalization of the Stochastic Collocation method, since exponential convergence is maintained for arbitrarily distributed parameters. This is due to the fact that the optimal collocation points and corresponding weights are computed based on the input distribution. Figure 5(b) shows the convergence for a normally distributed parameter. This is the general picture for non-uniform distribution.

PCM: mean PCM: var SC: mean SC: var

100

Error

Error

100

PCM: mean PCM: var SC: mean SC: var

10−5

10−5

10−10 10−10 1

2

3

4 Work

5

6

7

1

(a) Uniform distribution

2

3

4 Work

5

6

7

(b) Normal distribution

Figure 5. Convergence for the Probabilistic Collocation method (PCM) and the Stochastic Collocation method (SC) with respect to computational work for uniformly(a) and normally(b) distributed uncertain parameters. The result for the normal distribution is general for all non-uniform distributions.

V.

Application to steady flow around a NACA 0012 airfoil

A flow around a NACA0012 airfoil at an angle of attack of 5 degrees is considered. The Reynolds number is set to 3 · 106 . The deterministic computations are performed using the FineTM/Hexa solver by Numeca Int. on a grid of 50.000 cells. The grid layout is shown in Figure 6.

(a)

(b)

Figure 6. The computational mesh layout(a) and a detailed view of the airfoil(b).

11 of 14 American Institute of Aeronautics and Astronautics

The flow is modeled by the Reynolds-Averaged Navier-Stokes equations using the Spalart-Allmaras turbulence model. The air properties are at 0m ISA. The free stream velocity is assumed to be normally distributed, with mean µM = 0.3 and a coefficient of variation of 5%. The coefficient of variation is defined 2 as COV = σ/µ, this results in a variance of σM = 2.25e−4 . The uncertainty is propagated using the Probabilistic Collocation method, with five collocation points. The five collocation points are given in Table 3 with the corresponding solutions for the lift and the drag. The flow solver is run deterministically for every collocation point. After the runs were completed the deterministic solutions were used to compute the mean, variance and distribution function of properties of interest according equations (25), (26), and (15). Table 3. The input Mach number for the deterministic solves and the solutions for the drag and lift for the five collocation points.

Collocation point ω1 ω2 ω3 ω4 ω5

Mach number 0.25714 0.27966 0.30000 0.32033 0.34285

Drag [N] 25.08 29.59 34.04 38.85 44.60

Lift [N] 1099 1305 1510 1733 2001

0.12

3

0.1

2.5

0.08

2 fL (L)

fD (D)

The computed mean drag is µD = 34.14N with a variance of σD = 3.42N, this results in a coefficient of variation of 10%. The problem is therefore behaving exactly as expected, with the drag force quadratically dependent on the free stream velocity.PSfrag The lift shows the same behavior, with a mean of µ L = 1515N, a replacements variance of σL = 158N and a coefficient of variation of 10%. The probability density functions of the drag rag replacementsand lift are shown in Figure 7.

0.06

1.5

0.04

1

0.02

0.5

0

25

30

35 40 Drag, N

45

0

50

×10−3

1000

1200 1400 1600 1800 2000 2200 Lift, N (b)

(a)

Figure 7. The probability density functions of the drag(a) and the lift(b) obtained from a 4 th order Probabilistic Collocation computation, with FineTM/Hexa as deterministic solver.

Another way to visualize the resulting stochastic solution is by observing the mean and variance of the pressure field around the airfoil. Figure 8 shows the pressure fields obtained from a 4th order Probabilistic Collocation computation. The mean pressure around the airfoil is shown in Figure 8(a) and the variance of the pressure field in (b). From the variance field it is clear that the input uncertainty results in a higher uncertainty near the leading edge of the airfoil. If one is interested in the pressure distribution or skin friction on the airfoil surface, the uncertainty can be represented by uncertaintybar plots. Figure 9 shows the pressure distribution on the upper and lower surface with the uncertaintybars indicating the interval which contains 99% of all possible values. This 12 of 14 American Institute of Aeronautics and Astronautics

(a)

(b)

PSfrag replacements

Figure 8. The mean(a) and the variance(b) pressure field obtained from a 4th order Probabilistic Collocation computation, with FineTM/Hexa as deterministic solver.

interval is obtained from the distribution function of the pressure at each position on the airfoil. It can be seen in the figure that an uncertain free stream velocity leads to the highest variation in pressure on the upper surface near the leading edge. The uncertaintybars for the skin friction along the airfoil surface are shown in Figure 10. The skin friction shows large variations along the entire airfoil, for an uncertain free stream velocity. rag replacements 0.85

80

Skin friction, N/m2

Pressure P/P0

0.9 0.95 1

upper lower

60

40

20

1.05 0

0.2

0.4 0.6 x/c coordinate

0.8

0

1

0.2

0.4 0.6 x/c coordinate

0.8

1

Figure 10. Skin friction on the upper and the lower surface of the NACA0012 airfoil. The mean is indicated by the solid line, the uncertaintybars show the interval which contains 99% of all possible values.

Figure 9. Pressure distribution on the upper and the lower surface of the NACA0012 airfoil. The mean is indicated by the solid line, the uncertaintybars show the interval which contains 99% of all possible values.

VI.

0

Conclusions

The Probabilistic Collocation is shown to be a very efficient method for quantification of parametric uncertainties with arbitrary distributions. It is a generalization of the Stochastic Collocation method, which converges only exponentially for uniformly distributed parameters. The generalization is achieved by employing a chaos approach similar to the Polynomial Chaos methods. The convergence with respect to the polynomial chaos order of the approximation of the Probabilistic Collocation method and the Galerkin Polynomial Chaos method are the same. It is shown numerically that the resulting Polynomial Chaos ap13 of 14 American Institute of Aeronautics and Astronautics

proximations of both methods are equivalent. The amount of work of the Galerkin Polynomial Chaos method is more, since a coupled systems of equations has to be solved. For non-linear equations the efficiency gain of the Probabilistic Collocation method with respect to Galerkin Polynomial Chaos becomes significantly. The Non-intrusive Polynomial Chaos method results in approximations of the polynomial coefficients. For the same amount of work, the Non-Intrusive Polynomial Chaos method can, therefore, not obtain the same accuracy of the Probabilistic Collocation method. Furthermore, the non-intrusiveness of the Probabilistic Collocation method enables the use of sophisticated commercial codes. The method is demonstrated using the FineTM/Hexa solver of Numeca Int. A steady flow is considered around a NACA0012 airfoil with uncertain free stream velocity. A normal distribution is assumed with a coefficient of variation of 5%. The stochastic solution is represented in three different ways. The probability density functions of quantities of interest, like the drag and lift, are shown. Furthermore, the mean and variance pressure fields are given to see places in the computational domain which are sensitive to the uncertain input parameter. Finally, uncertaintybar plots are shown for the pressure and skin friction distribution on the airfoil surface. The stochastic results demonstrate how uncertainty quantification produces results which provide more information about the solution with respect to uncertain input parameters.

References 1 Wiener,

N., “The Homogeneous Chaos,” Amer. J. Math., Vol. 60, 1938, pp. 897–936. R. G. and Spanos, P. D., Stochastic Finite Elements: A Spectral Approach, Dover, 1991. 3 Walters, R. W., “Towards Stochastic Fluid Mechanics via Polynomial Chaos,” Proceedings of the 41 st AIAA Aerospace Sciences Meeting and Exhibit, AIAA paper 2003–0413, Reno, January 2003. 4 Hosder, S., Walters, R. W., and Perez, R., “A Non-Intrusive Polynomial Chaos Method for Uncertainty Propagation in CFD Simulations,” Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, AIAA paper 2006–891, Reno, January 2006. 5 Reagan, M. T., Najm, H. N., Ghanem, R. G., and Knio, O. M., “Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection,” Combustion and Flame, Vol. 132, 2003, pp. 545–555. 6 Mathelin, L. and Hussaini, M. Y., “A Stochastic Collocation algorithm for uncertainty analysis,” Tech. Report NASA/CR2003-212153, NASA Langley Research Center, 2003. 7 Mathelin, L., Hussaini, M. Y., and Zang, T. A., “Stochastic approaches to uncertainty quantification in CFD simulations,” Numerical Algorithms, Vol. 38, 2005, pp. 209–236. 8 Loeven, G. J. A., Witteveen, J. A. S., and Bijl, H., “(Student paper) Efficient Uncertainty Quantification in Computational Fluid-Structure Interaction,” Proceedings of the 8th AIAA Non-Deterministic Approaches Conference, AIAA paper 2006–1634, Newport, May 2006. 9 Babuˇ ska, I., Nobile, F., and Tempone, R., “A Stochastic Collocation Method for Elliptic Partial Differential Equations with Random Input Data,” ICES report 05-48, 2005. 10 Debusschere, B. J., Najm, H., P´ ebay, P. P., Knio, O. M., Ghanem, R. G., and Maˆıtre, O. P. L., “Numerical challenges in the use of Polynomial Chaos representations for stochastic processes,” SIAM, J. Sci. Comput., Vol. 26, No. 2, 2001, pp. 698–719. 11 Witteveen, J. A. S. and Bijl, H., “Using Polynomial Chaos for Uncertainty Quantification in Problems with Nonlinearities,” Proceedings of the 8th AIAA Non-Deterministic Approaches Conference, AIAA paper 2006–2066, Newport, May 2006. 12 Xiu, D. and Hesthaven, J., “High order collocation methods for differential equations with random inputs,” SIAM J. Sci. Comput., Vol. 27, No. 3, 2005, pp. 1118–1139. 13 Xiu, D. and Karniadakis, G. E., “The Wiener-Askey Polynomial Chaos For Stochastic Differential Equations,” SIAM J. Sci. Comput., Vol. 24, No. 2, 2002, pp. 619–644. 14 Wan, X. and Karniadakis, G. E., “Beyond Wiener-Askey expansions: Handling arbitrary PDFs,” J. Sci. Comput., Published online: December 23 2005. 15 Witteveen, J. A. S. and Bijl, H., “Modeling Arbitrary Uncertainties Using Gram-Schmidt Polynomial Chaos,” Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, AIAA paper 2006–896, Reno, January 2006. 16 Walters, R. W. and Huyse, L., “Uncertainty analysis for fluid mechanics with applications,” ICASE Report No. 2002-1, 2002. 17 Golub, G. H. and Welsch, J. H., “Calculation of Gauss quadrature rules,” Mathematics of Computation, Vol. 23, No. 106, 1969, pp. 221–230. 18 Gautschi, W., “Orthogonal polynomials (in Matlab),” Journal of Computational and Applied Mathematics, Vol. 178, 2005, pp. 215–234. 19 Gander, M. J. and Karp, A. H., “Stable computation of high order Gauss quadrature rules using discretization for measures in radiation transfer,” Journal of Quantitative Spectroscopy Radiative Transfer , Vol. 68, 2001, pp. 213–223. 20 van Zuijlen, A. H. and Bijl, H., “Implicit and explicit higher order time integration schemes for structural dynamics and fluid-structure interaction computations,” Computers and Structures, Vol. 83, 2004, pp. 93–105. 2 Ghanem,

14 of 14 American Institute of Aeronautics and Astronautics