Performance impact of flow and geometric variations for a turbine blade using an adaptive NIPC method

Accepted Manuscript Performance impact of flow and geometric variations for a turbine blade using an adaptive NIPC method

Zhiheng Xia, Jiaqi Luo, Feng Liu

PII: DOI: Reference:

S1270-9638(18)31940-0 https://doi.org/10.1016/j.ast.2019.04.025 AESCTE 5117

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

3 September 2018 22 January 2019 16 April 2019

Please cite this article in press as: Z. Xia et al., Performance impact of flow and geometric variations for a turbine blade using an adaptive NIPC method, Aerosp. Sci. Technol. (2019), https://doi.org/10.1016/j.ast.2019.04.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Performance impact of flow and geometric variations for a turbine blade using an adaptive NIPC method Zhiheng Xiaa , Jiaqi Luob,1,∗, Feng Liuc a Department

of Aeronautics and Astronautics, Peking University, Beijing 100871, China of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China c Department of Mechanical and Aerospace Engineering, University of California, Irvine, CA 92697-3975, United States b School

Abstract Geometric variations of the realistic blades due to the limited machining precision and flow variations at the outlet boundary due to the operational condition changes for turbomachinery blades cannot be totally eliminated, the significant effects of which require to be taken into account in the aerodynamic shape design and optimization. The present paper investigates the performance impact of flow and geometric variations by using the polynomial chaos. The adaptive non-intrusive polynomial chaos (NIPC) model based on an adaptive sparse grid technique is firstly introduced, the response performance of which is validated through function experiments and uncertainty quantification of performance change for a three-dimensional turbine blade. Then the adaptive NIPC is used to statistically evaluate the adiabatic efficiency change of the turbine blade due to the geometric and back pressure variations separately. The performance change exhibits nonlinear dependence on the geometric and back pressure variations. Moreover, the geometric variations induce more performance deteriorations to the turbine blade in the present study. Finally, the simultaneous performance impact of geometric and back pressure variations are investigated. The results are presented in detail and compared with those obtained from the uncoupled studies, demonstrating that these two kinds of uncertain effects to the performance change are nonlinearly dependent. Keywords: uncertainty quantification, turbomachinery, geometric variation, flow variation, adaptive sparse grid, non-intrusive polynomial chaos Preprint submitted to Aerospace Science and Technology

April 19, 2019

1. Introduction Turbomachines are desired to work at the specific operation conditions in industry and the real profiles of turbomachinery blades are expected to be the same as design. However, in the real world, the flow at the inlet and outlet of turbomachinery cannot always maintaine as the design conditions. For example, flow distortion is quite common in transonic compressors, especially at the operation conditions near stall [1,2]. The inlet flow of the turbine deviates from the design conditions because of the high turbulence and spanwise non-uniformity at the outlet of the upstream combustor [3,4]. Furthermore, the real blade profiles are usually different from the nominal ones due to the limited machining precision [5-7]. The flow and geometric variations inevitably influence the aerodynamic performance of turbomachinery blades. Since 1950s, the studies on performance impact of flow and geometric variations have become more and more active. At the beginning of this century, Garzon et al. [8] successfully evaluated the performance change of turbomachinery blades due to the manufacturing variations by using computational fluid dynamics (CFD), releasing the researchers from expensive and time-consuming experiments. In the past decade, more and more attentions have been paid on the aerodynamic performance uncertainty quantification (UQ) for turbomachinery blades and airfoils using CFD [9-20]. Loeven et al. [9] introduced a non-intrusive approach to study the effects of uncertain free stream velocity to the aerodynamic performance of NACA0012 airfoil. Bhattrai et al. [10] employed the gradient-enhanced kriging method to evaluate the impact of the variations of inlet flow velocity, density, temperature, and wall temperature on the aerodynamic force acting on a hypersonic trailing-edge flap model. Zheng et al.[11] employed the novel Bernstein-polynomial-based method to propagate geometric uncertainty for hypersonic vehicles. Loeven et al. [12] used the probabilistic collocation method to evaluate the performance impact of the uncertain inlet total pressure profile for a transonic compressor rotor, NASA Rotor 37. The results ∗ Corresponding 1 Email

author. address: [email protected] (J. Luo)

2

showed that the mass flow rate was dramatically sensitive to the change of total pressure profile. Pecnik et al. [13] studied the performance impact of the variations of boundary layer transition and turbulence inlet conditions for a transonic turbine guide vane. Currently, some UQ methods have been widely applied: the direct Monte Carlo simulation (MCS), sensitivity-based method such as the method of moment (MM) and polynomial chaos (PC) model [14-20]. Grazon et al. [8] determined the converged statistics by the direct MCS from two thousand samples. If the direct MCS method is employed in robust design and optimization, it would be time-consuming to obtain the optimum. The sensitivity-based method can be applied to small-scale UQ studies, for which the performance change exhibits almost linear or quadratic dependence on the perturbation. The sensitivity-based method can support quite accurate results with high efficiency. By the MM method, Putko et al. [16] determined the uncertainty propagations of geometric and flow variations separately in the quasi-one-dimensional flow. Recently, Luo et al. [17,18] used the adjoint method to calculate the second order sensitivities, by which the evaluations of performance impact of manufacturing variations for a turbine blade were successfully achieved with significantly reduced computational cost. However, it should be noticed that the sensitivity-based method is only useful for the small-scale UQ problems. When the perturbation increases to be considerable, the evaluations of performance change by using the sensitivity-based method lack of sufficient accuracy. In such situations, a substitutional method is necessary to propagate performance uncertainty for the large-scale UQ problems. Due to the easy implementation and significantly reduced computational cost compared with the direct MCS, PC methods (PCs) have been widely used. Essentially, PCs are surrogate model methods. Once the model is constructed based on a series of samples, function response can then be achieved and thus the MCS can be utilized to determine the statistics. PCs can be classified by intrusive polynomial chaos (IPC) and non-intrusive polynomial chaos (NIPC). Although one can directly get the PC representation of model outputs by a one-time shot solution, one has to reformulate the governing equations. On 3

the contrary, NIPC regards the original code as a black box [21]. Hosder et al. [22] applied NIPC to the three-dimensional wing flow considering the variations of free-stream Mach number and angle of attack. Zhao et al. [23] used an adjusted NIPC method to quantify the uncertainties of Mach number and lift coefficient in the robust design of the natural-laminar-flow (NLF) airfoil. Gopinathrao et al. [24] studied the propagation of inlet total pressure uncertainty for a transonic compressor rotor using NIPC. By using NIPC, Ghisu et al. [25] studied the performance uncertainties of a compression system based on an one-dimensional design procedure, by which the robust design optimization was then carried out. Ma et al. [26] studied the performance impact of manufacturing variability for a compressor blade using NIPC and then performed the robust design optimization. Generally, most of the present UQ studies concerned only the effects of either geometric or flow variations. Considering the simultaneous effects of geometric and flow variations, NIPC-based statistical evaluation of performance change requires more samples used for model construction and thus increased computational cost. Moreover, the uncertainty effects of these two kinds of uncertain parameters require to be further studied and demonstrated. To overcome the curse of dimensionality for UQ involving a large number of uncertain parameters, the adaptive NIPC method has been studied and applied to UQ in aerodynamics [15,27]. It should be noticed that the aforementioned open literatures [15,27] concerns the recursive approximation based on an iterative adaptation of polynomial basis to the solution. In the present study, an adaptive sparse grid technique is employed in the adaptive NIPC to reduce the required computational cost with improved response accuracy of the model. Generally, there are two main sampling approaches for NIPC, random and deterministic samplings. By the former sampling method, the coefficients of the polynomials vary as the number of random samples changes. In the study, the deterministic sampling method with Galerkin projection is used. The adaptive NIPC is investigated through function experiments to illustrate the improvements on response performance. Then the performance impact of geometric and flow variations for a turbine 4

blade is separately investigated (separated cases). Moreover, the simultaneous effects of flow and geometric variations to the change of aerodynamic performance are statistically evaluated. The results are presented and compared with those obtained from the separated cases in detail. 2. Uncertainty Quantification 2.1. Flow Uncertainty In the open literatures [9,12,22], two different description methods for flow variations have been widely used, the Gaussian and uniform distributions. The former has been used in the natural and social sciences to represent the real-valued random variables whose distributions are not clear [28]. Thus it is employed to describe the flow variations for turbomachinery blades in the present study. The perturbation of outlet back pressure of a three-dimensional turbine blade is assumed to agree the following normal distribution ⎧ 1 ξ2 ⎪ ⎪ ⎨ √ exp(− ), ξ ∈ [−E, E] 2 2π f (ξ) = (1) ⎪ ⎪ ⎩ 0, otherwise where E = 2.0, ξ is defined as the pressure perturbation ξ=

p − pb σp

(2)

p and pb are the perturbed and reference back pressure, respectively, σp is the standard deviation of pressure perturbation. Although only the back pressure is regarded as the uncertain parameter, ξ can be a vector to represent a series of uncertain flow variables. In this study, it is assumed that the outlet pressure and geometrical perturbation are independent of each other. 2.2. Geometric Uncertainty Principal component analysis (PCA) has been widely used to extract the geometric variation characteristics [17-19]. It is thus employed in the study to determine the basis modes of geometric variations for a turbine blade. A scatter matrix S can be firstly 5

constructed as following. S = XT X, X = (xij )n×p = (˜ xij − xj )n×p n  1 xj = x ˜ij , x ˜ij = x ˆij − x0ij n i=1

(3)

where n and p represent the numbers of manufactured blades and scanning points for each, respectively, x ˆij and x0ij are the coordinates of the scanning points on the manufactured and nominal blades, respectively, xj indicates the mean geometric variation at the j-th scanning point. Decomposing the scatter matrix by the singular value decomposition (SVD), we have S = Q(ΛT Λ)QT

(4)

where Q and Λ are the eigenvector and eigenvalue matrices of X. The basis modes, Ψ can be determined by Ψ = ΛQT

(5)

All the characteristics of geometric variations are included in Ψ = {ψ1 , · · · , ψp }. Then the geometric variations of the manufactured blades from the same machine can be described by the basis modes in a weighted summation form. x=x+

p 

s i ψi

(6)

i=1

where si is the weight of ψi . It is well known that the geometric variations are randomly produced, all the weights appear in Eq. (6) are random variables to describe the geometric uncertainty. Without a large number of manufactured blades, a widely used modeling method [17-20] is employed in the present study to model the geometric variations, which can be written as x = x + xp

(7)

where xp denotes the zero-mean Gaussian random perturbation, the corresponding co-

6

variance function is given as Cov(xi , xj ) = σi σj exp(−

xi − xj 2 ) l2

(8)

where σ and l are the geometric standard deviation and user defined spatial correlation length, respectively; xi − xj  is the Euclidean distance between the scanning points, xi and xj . Decomposing the covariance function, rather than S, the basis modes can be determined. To illustrate the procedures of the geometric uncertainty quantification by using the modeling method, the spatial correlation length l = 0.4 and an uniform geometric standard deviation σ = 0.2 millimeter (mm) are used for a three-dimensional turbine blade, which has been studied in previous work [17,18]. The blade configuration is shown in Figure 1. The flow on the outer spans is transonic, while it is supersonic on the inner spans. The aerodynamic performance parameters are obtained by solving the Euler equations on a grid with the resolutions 97, 33, 33 in the streamwise, pitchwise and spanwise, respectively. Figure 2 shows the histograms of the energy ratio for the first six basis mode, Pk and the accumulative energy ratio, APk , the definitions of which are

Casing

Outlet Blade

Inlet

Hub

Figure 1: Configuration of the turbine blade

7

P

Figure 2: Histograms of Pk and APk

given as

k λi , APk = i=1 p λ i i=1 i=1 λi

λk Pk =  p

(9)

where λi denotes the i-th eigenvalue. For any singular value system, the state vector in the system can be described by a series of primary basis modes with sufficient accuracy. From Figure 2, it is known that the accumulative energy ratio of the first five basis modes is already over 99.5%, demonstrating that the geometric variations can be approximately described by the first five basis modes, as shown by Eq. (6). Figure 3 shows the geometric variations for the first four basis modes. The blade is cut open at the trailing edge and then stretched in the axial direction. S.S., LE and P.S. in the figures represent the suction side, leading edge and pressure side, respectively. In the first basis mode, the blade is thinned on the whole span and the maximum change of blade thickness is about 0.17 mm. In the second basis mode, the blade is thickened on the inner spans and thinned on the outer spans. In the third and fourth basis modes, the blade is thickened on the middle spans and after 50% axial chord, respectively. Figure 4 presents the relative variations of isentropic Mach number on the S2 streamsurface for the 8

.1 4

4

-0

-0. 1

8

.1 1

-0.16

7 -0.0

-0.16

-0

.1 1

8

-0 .0

0 -0.

-0

.1 3

-0 .1 3

-0.07

0.055

0.055

0.1 -0

S.S.

(a)

P.S.

LE

0.1

S.S.

0.075

-0 1 .1

-0.1

P.S.

LE

(b)

-0.05

0.065

0.07

-0.04

-0.04

0.065 -0.015

S.S.

-0 . 00

5

0.07

(c)

P.S.

LE

S.S.

P.S.

LE

(d)

Figure 3: Geometric variations for: (a) mode 1; (b) mode 2; (c) mode 3; (d) mode 4

2

0.40

0.00

0.60

-0.20

0 -0.

0.28

0.20

0.40

-0.01

0.2 0.15

0. 06

0.13

0

-0.20

0.20

-0.20

0

(b)

0.0

0. 20

-0.40

(a)

0.00

1.00

-1.00

20 -0 .

(c)

0.06

(d)

Figure 4: Relative variations of isentropic Mach number on the S2 streamsurface: (a) mode 1; (b) mode 2; (c) mode 3; (d) mode 4

corresponding blades. Generally, the changes of isentropic Mach number are significant. For the first blade, the maximum change is over 1% in the supersonic flow region, which inevitably influence the aerodynamic performance of the turbine blade.

9

3. Adaptive Non-Intrusive Polynomial Chaos 3.1. Sparse Grid 3.1.1. Principle of Sparse Grid As presented in the reference [29], the sparse grid interpolation starts from the available univariate interpolation formula. U i (f ) =

mi 

f (xij ) · αji

(10)

j=1

where xij and αji are the one-dimensional interpolation points and the basis functions such as the piecewise linear splines and Lagrange polynomials, respectively. In the multivariate case with d > 1, define the tensor product as (U

i1

mi1 

⊗ ··· ⊗ U ) = id

mid

···

j1 =1



jd =1

f (xij11 , · · · , xijdd ) · (αji11 ⊗ · · · ⊗ αjidd )

(11)

where ⊗ indicates bilinear composition. Smoljak [30] proposed the Smoljak formula which are the linear combinations of the tensor products in the cases: only the products with a relatively small number of knots are used; the linear combination is chosen in such a way that an interpolation property for d = 1 is preserved for d > 1 [29]. U0 = 0 Δi = U i − U i−1

(12)

|i| = i1 + · · · + id ,

i ∈ Nd+

In such cases, the Smoljak formula can be written as A(q, d)(f ) =



(Δi1 ⊗ · · · ⊗ Δid )

(13)

|i|≤q

For the k-level sparse grid [31,32], Eq. (13) can also be written as A(q, d)(f ) =



 (−1)q−|i| ·

k+1≤|i|≤q

10

d−1 · (U i1 ⊗ · · · ⊗ U id )(f ) q − |i|

(14)

where k = |i| − d. One only needs to evaluate the function value on the sparse grid

H=

(Xi1 × · · · × Xid )

(15)

k+1≤|i|≤q

where Xi = {xi1 , · · · , ximi } denotes the points for U i . The sparse grid interpolation accuracy is O(N −2 · (log N )d−1 ) with respect to the maximum-norm where N is the number of points used. And the discretizations on sparse grids involve O(N · (log N )d−1 ) degrees of freedom only rather than the exponential dependence O(N d ) of conventional approaches, which can overcome the curse of dimensionality to some extent [33]. In the present study, the Newton-Cotes grid and piecewise linear basis functions are employed in the adaptive sparse grid for the easy implementations. 3.1.2. Adaptive Sparse Grid Since the adaptive sparse grid has already been introduced in the open literatures [34,35], the principles are conceptually introduced herein. First we define the hierarchical surplus [36,37] on the k-th grid level wjk,i = f (xij ) − A(k + d − 1, d)(xij )

(16)

where j denotes the number of the necessary basis nodes on the k-th grid level. As k tends to infinity, the hierarchical surplus tends to zero for the continuous functions. Because at the points xij , A(k + d − 1, d) is corrected by f (xij ) to approach the exact value. If the function has steep gradients, the hierarchical surplus should be large. So the the hierarchical surplus can be thought of to be a natural candidate for error estimation and tolerance control in the adaptive sparse grid. Following is the procedure of adaptive sparse grid. Step 1: Define the maximum grid level, kmax ; produce a number of initial nodes on the grid level k = 1. Step 2: Produce a number of actnum nodes on the grid level k = 2, which are stored in actgrid. Step 3: If k < kmax and actnum > 0, move all the nodes in actgrid to oldgrid and 11

then calculate the hierarchical surplus (w) for the elements in oldgrid; or else end the interpolation process. Step 4: If w is larger than the error tolerance, , produce 2d additional nodes and store in actgrid, d represents the dimension. Step 5: k = k + 1, goto Step 3. The additional nodes produced in Step 4 are on the grid level k + 1. They are closest to the corresponding basis nodes. In the meanwhile, the repetitive nodes are eliminated. Moreover, the nodes with the hierarchical surpluses lower than the tolerance are maintained to keep the adaptive sparse grid balanced [34]. To illustrate the improvements of response performance for the adaptive sparse grid, a function experiment is presented. The function is −0.0001 ((x1 − 0.3)2 + 0.010)((x2 − 0.3)2 + 0.010) −0.0001 + ((x1 − 0.7)2 + 0.015)((x2 − 0.7)2 + 0.015)

f=

(17)

where x1 , x2 ∈ [0, 1]. The function reaches the extremes at (0.3,0.3) and (0.7,0.7). Table 1 gives the numbers of grid points, N , interpolation extreme value at (0.3,0.3) and the maximum deviation, max for both the static and adaptive sparse grid methods. The maximum deviation between the interpolation and the exact function is defined as max = max|fi − Aq,2 (fi )|,

i = 1, · · · , N

(18)

where f and Aq,2 indicate the exact function and interpolation result, respectively. When the node numbers are in the same order of magnitude, the response accuracy and the interpolated extreme value are significantly improved by using the adaptive sparse grid. Figure 5 gives the grid evolution for the adaptive sparse grid. From the figures it can be

Table 1: Comparisons of static and adaptive sparse grid methods

static adaptive

N 1537 1493

Aq,2 (0.3, 0.3) -0.9718 -0.9993 12

max 2.9E-002 3.7E-003

(a)

(b)

(c)

(d)

Figure 5: Grid evolution of the adaptive sparse grid with tolerance  = 10−3

found that most of the nodes are distributed around the extremes, demonstrating that the proposed adaptive sparse grid technique is efficient. 3.2. Adaptive Non-Intrusive Polynomial Chaos NIPC has already been introduced and widely used for UQ studies, hereby the principles are also briefly introduced in the following. A general second-order random process X(θ) regarded as a function of θ, i.e., the random event, can be represented in the form [38,39] X(θ) = a0 H0 +

∞ 

ai1 H1 (ξi1 (θ)) +

i1 =1

+

i1  i2 ∞  

i1 ∞  

ai1 i2 H2 (ξi1 , ξi2 )

i1 =1 i2 =1

ai1 i2 i3 H3 (ξi1 (θ), ξi2 (θ), ξi3 (θ)) + · · ·

i1 =1 i2 =1 i3 =1

13

(19)

where Hn (ξi1 , · · · , ξin ) denotes the polynomial chaos basis of order n in terms of the multi-dimensional independent random variables, ξ = (ξi1 , · · · ξin ). Different random distributions correspond to different polynomial chaos[40]. Hermite polynomials can be used when Gaussian distribution is specified. On the other hand, X(θ) can also be described as X(θ) =

∞ 

uj Ψj (ξ)

(20)

j=0

where Ψj (ξ) is the random basis function corresponding to the j th mode. There is a one-to-one correspondence between Hn (ξi1 , · · · , ξin ) and Ψj (ξ), and also between the coefficients ai1 ,··· ,ir and uj , which is the deterministic component and the amplitude of the j th fluctuation. The polynomial chaos forms a complete orthogonal basis in the L2 space of the Gaussian random variables.  Ψi Ψj  = Ψ2i δij

(21)

·, · means the inner product in the L2 space. In practice, X(θ) can be described in a summation form of p basis functions by the truncation theory. X(θ) ≈

P 

uj Ψj (ξ), P =

j=0

(n + p)! −1 n!p!

(22)

Exponential convergence of the above description versus p can be obtained for the random process X with the same density as that of ξ [41]. Orthogonality of the PC basis enables the evaluation of the truncated PC representation, as shown by Eq. (22) by projecting the random process X onto the PC basis. X(θ) =

P 

uj Ψj (ξ), uj =

j=0

XΨj  Ψ2j 

(23)

There are two methods to get the NIPC model coefficients uj , regression and Galerkin projection. By the regression method, the model coefficients vary when the number of random collocation points changes [22]. In such case, the Galerkin projection combined 14

with the adaptive sparse grid is used. Based on the adaptive sparse grid, the corresponding adaptive NIPC can be implemented. Firstly, the function in the space can be approximately represented by using a series of adaptive sparse grids. Then the model coefficients of NIPC can be determined by using the Galerkin projection, as shown by Eq. (23). Figure 6 presents the procedure of adaptive NIPC.

Figure 6: Procedure for the adaptive NIPC

Another function experiment is used to study the response performance of adaptive NIPC. The function is f = 0.1 ∗ ln(x1 + 5) ∗ cos(x2 − x1 ), x1 , x2 ∈ [−2, 2]

(24)

where x1 , x2 are independent Gaussian random variables. The average error (avr ) and

15

the standard deviation (std ) are calculated by the following equation 1  |f (xi1 , xi2 ) − fˆ(xi1 , xi2 )| N i=1

 N  2  1  = f (xi1 , xi2 ) − fˆ(xi1 , xi2 ) N − 1 i=1 N

avr =

std

(25)

where N is the number of randomly produced variables. f and fˆ are exact solutions and the model responded values, respectively. Table 2 shows the statistical results by using the eight-order static and adaptive NIPC methods. The error tolerance (max ) of the adaptive NIPC is 1.0E − 04, 1.0E − 03 and 1.0E − 02 for Case A, Case B and Case C, respectively. Case D is the static NIPC. The results demonstrate that as the tolerance reduces, the model response accuracy improves; moreover, the statistical results approach closer to those of MCS. Furthermore, it is obvious that higher accuracy needs more function evaluations. Compared with Case B, the statistical results of Case D are more deviated from those of MCS and the model response accuracy is lower, indicating that the response performance of the present adaptive NIPC is improved. Figure 7 compares the probability density function (PDF) distributions of Case A and MCS. The adaptive NIPC based PDFs are almost the duplicates of MCS ones, further demonstrating the effectiveness and accuracy of the adaptive NIPC on function response. Table 2: Function experiment by the adaptive NIPC

μ σ avr std N

MCS 0.0719 0.0893 – – 10000

Case A 0.0719 0.0892 3.1E-05 1.2E-05 1989

Case B 0.0717 0.0890 2.7E-04 1.2E-04 474

Case C 0.0705 0.0872 2.3E-03 1.2E-03 121

Case D 0.0717 0.0888 4.5E-04 2.6E-04 501

The three-dimensional turbine blade in Figure 1 is studied for the purpose of further validating the response performance of the adaptive NIPC. Four order adaptive NIPC is used to propagate the geometric uncertainties. The geometric variations are constructed 16

0.09

0.09

PDF

0.12

PDF

0.12

0.06

0.03

0.00

0.06

0.03

-0.16

-0.08

0.00

0.08

0.00

0.16

-0.16

f

-0.08

0.00

0.08

0.16

f

(a) NIPC of Case A

(b) MCS

Figure 7: Comparisons of PDFs

by the first five basis modes extracted by PCA from the covariance function, the spatial correlation length l = 0.4 and the uniform geometric standard deviation σ = 0.2mm. The back pressure at the turbine outlet is pb /P0 = 0.35, where P0 is the inlet total pressure. A number of five thousand random samples are used for MCS. The results of adaptive NIPC are presented and compared in Table 3, where the subscript η indicates the adiabatic efficiency, μ, σ and skew indicate the statistic mean value, standard deviation and skewness, respectively, for the adiabatic efficiency change. The definition of skewness is given as skew =

N 1  fi − μ 3 ) ( N i=1 σ

(26)

Case A and Case B are adaptive NIPC studies with the tolerance max = 1.0E − 04 and max = 1.0E − 03, respectively, while Case C is the static NIPC. The results of Case A and Case B are more closer to those of MCS than Case C, demonstrating that the

Table 3: Adaptive NIPC validation for the turbine blade

μη /% ση skewη avr std N

MCS -0.0841 3.10E-03 -0.860 – – 5000

Case A -0.0874 3.09E-03 -0.865 3.68E-05 3.07E-05 879 17

Case B -0.0987 3.08E-03 -0.847 1.43E-04 7.14E-05 126

Case C -0.1138 3.00E-03 -0.920 1.95E-04 1.65E-04 182

adaptive NIPC is practicable to quantify the performance uncertainty for the turbine blade. Furthermore, even with less samples, the response error for Case B is significantly lower than Case C. Through the UQ studies of both function experiment and geometric variations for the turbine blade, the effectiveness of the adaptive NIPC is validated when compared with the direct MCS and the static NIPC. Moreover, the response accuracy of the adaptive NIPC is convinced through comparing with the direct MCS results. 4. Statistical Evaluations In the following, the effects of geometric variations, outlet back pressure perturbations and those of these two on the adiabatic efficiency change are statistically evaluated for the three-dimensional turbine blade by using the adaptive NIPC. 4.1. Performance Impact of Geometric Variations The modelling method as shown by Eq. (8) is used, in which the spatial correlation length l = 0.4 and the uniform geometric standard deviation σ = 0.2mm are given. Each zero mean random variable agrees the following standard normal distribution ⎧ 1 s2 ⎪ ⎪ ⎨ √ exp(− i ), si ∈ [−E, E] 2 2π f (si ) = ⎪ ⎪ ⎩ 0, otherwise

(27)

where E = 2.0. The inlet incidence angle α = 0 and the outlet back pressure is pb /P0 = 0.35. By solving the Euler equations, the nominal adiabatic efficiency is η0 = 0.9583. The first five basis modes of geometric variations are used to construct the geometric deviations of manufactured blades. Four order Hermite PCs are used to quantify the

18

adiabatic efficiency change and p = 5 in the present study. η(ξ) = a0 +

p 

ai H1 (ξi ) +

i=1

+

p 

ci H2 (ξi ) +

+

bij H1 (ξi )H1 (ξj )

i=1 j=i+1 p p  

dij H2 (ξi )H1 (ξj )

i=1 j=i

i=1 p p  

p p  

p 

eijk H1 (ξi )H1 (ξj )H1 (ξk ) +

i=1 j=i+1 k=j+1

+

p p  

p 

p 

fi H3 (ξi )

i=1 p 

(28)

hijkm H1 (ξi )H1 (ξj )H1 (ξk )H1 (ξm )

i=1 j=i+1 k=j+1 m=k+1

+

+

p  p  p 

gijk H2 (ξi )H1 (ξj )H1 (ξk )

i=1 j=i k=i p p  

p p  

i=1 j=i+1

i=1 j=i

qij H2 (ξi )H2 (ξj ) +

Qij H3 (ξi )H1 (ξj ) +

p 

pi H4 (ξi )

i=1

The statistical results have already been given in Table 3. The MCS results indicate that the absolute adiabatic efficiency drops by 0.0841% compared with the nominal value. Compared with the MCS results, the absolute deviations of μη and ση for Case A are about 0.0033% and 1.0E − 05, respectively, while those for Case B are about 0.0156% and 9.0E −05, respectively. Although the results of Case A are more accurate, it requires more time-consuming CFD computations. In fact, the results of Case B are comparatively accurate enough and most importantly, the CFD computations are dramatically reduced. Figure 8 presents the PDF and CDF distributions of adiabatic efficiency change for both MCS and NIPC. From the figures, it is known that the adiabatic efficiency change is nonlinearly dependent on the geometric variations, which can also be convinced by skewness given in Table 3. The PDF and CDF obtained from the adaptive NIPC are quite close to those of MCS, demonstrating the sufficient accuracy of the adaptive NIPC on statistically evaluating the performance change due to the geometric variations. To further evaluate the performance change and validate the response performance of the adaptive NIPC, performance impact of different geometric standard deviations for the turbine blade are investigated. The geometric standard deviation varies from 19

0.16

0.16

0.12

0.12

PDF

0.20

PDF

0.20

0.08

0.08

0.04

0.04

0.00

-1.03

-0.63

-0.23

0.17

0.00

0.57

-1.03

-0.63

dh(%)

-0.23

0.17

0.57

dh(%)

(a) Case A

(b) Case B

0.20

MCS

1.0

Case B Case A 0.16

0.8

0.6

PDF

CDF

0.12

0.4

0.08

0.2 0.04

0.0 0.00

-1.03

-0.63

-0.23

0.17

-1.03

0.57

-0.63

-0.23

0.17

0.57

dh(%)

dh(%)

(c) MCS

(d) CDF

Figure 8: PDF and CDF distributions of adiabatic efficiency change due to the geometric variations

0.15mm to 0.70mm with an uniform increment of 0.05mm. The geometric variations are also constructed by the first five basis modes extracted by PCA from the corresponding covariance function with the spatial correlation length l = 0.4. It should be noticed that the maximum geometric standard deviation is much lower than the threshold value σmax = 2.1mm, beyond which the blade geometry is dramatically deformed and the CFD computations fail to converge. Similarly, the four order Hermite PCs are used for UQ of adiabatic efficiency change. The tolerance of the adaptive NIPC is 1.0E − 03. A number of five thousand samples are used for MCS. Figure 9 compares the statistic mean values, standard deviations and skewness of adiabatic efficiency change between the direct MCS and adaptive NIPC. As the geometric standard deviation increases, the negative change of adiabatic efficiency increases fast, 20

MCS NIPC

-1.4

1.2

-1.0

1.0

m(%)

s(%)

-1.2

-0.8

0.8

-0.6

0.6

-0.4

0.4

-0.2

0.2

0.0

MCS NIPC

1.4

1

2

3

4

5

6

Case

7

8

9

10

11

0.0

12

1

(a) Mean value 0.0

3

4

5

6

Case

7

8

9

10

11

12

11

12

(b) Standard deviation 2000

MCS NIPC

1500

CFD evaluation

-0.5

Skew

2

1000

-1.0

500 -1.5

1

2

3

4

5

6

Case

7

8

9

10

11

0

12

(c) Skewness

1

2

3

4

5

6

Case

7

8

9

10

(d) CFD Evaluation

Figure 9: Comparisons of statistics between MCS and adaptive NIPC

demonstrating more significant effects of geometric variations. Moreover, the standard deviation of adiabatic efficiency change performs linear dependence on the geometric standard deviation. In fact, Luo et al. [17,18] has already found the linear dependence in the study of performance impact by using sensitivity-based method within a reduced tolerance space. The skewness increases as the geometric standard deviation increases, demonstrating more intensive nonlinear dependence of performance change on the increasing geometric variations. The statistics based on the four-order adaptive NIPC match well with the MCS results when the geometric standard deviation is lower than 0.45mm. When the geometric standard deviation is larger than 0.45mm, the four-order adaptive NIPC-based skewness evidently deviates from the results of MCS. However, in the study it is found that when the seven-order adaptive NIPC is used, the deviations can be significantly reduced. In the open literatures [18], statistics obtained by 21

sensitivity-based method obviously deviate from the MCS ones when the geometric standard deviation reaches 0.3mm. In such situations, compared with the sensitivity-based method, the statistic results of adaptive NIPC are more accurate for the large-scale UQ problems. Besides, the computational cost of adaptive NIPC is also investigated. The number of CFD evaluations increases rapidly with the increasing geometric standard deviation. So the advantage of the adaptive NIPC over the direct MCS will disappears when the geometric standard deviation is too large. For example, although the sevenorder adaptive NIPC can evaluate the adiabatic efficiency change with sufficient accuracy when σ = 1.0mm, the number of CFD evaluations reaches up to 4165. 4.2. Performance Impact of Back Pressure Perturbations For the study of performance impact due to the back pressure perturbation, the inlet incidence angle is α = 0 and the reference back pressure is pb /P0 = 0.35. The zero mean pressure perturbation agrees a standard normal distribution, as shown by Eq. (1). The standard deviation is σp /P0 = 0.05. The maximum deviation of back pressure is about 10% of the inlet total pressure. Four-order Hermite PCs are used to quantify the adiabatic efficiency change. η(ξ) = a0 + a1 H1 (ξ) + a2 H2 (ξ) + a3 H3 (ξ) + a4 H4 (ξ)

(29)

Table 4 presents the statistics of both MCS and NIPC. The tolerance, max for Case A and Case B are 1.0E − 04 and 1.0E − 03, respectively. The MCS results show that the absolute adiabatic efficiency drops by about 0.040% from the nominal. Compared with

Table 4: Performance impact of back pressure perturbation

μη /% ση skewη avr std N

MCS -0.0401 1.79E-03 -1.051 – – 5000

Case A -0.0413 1.78E-03 -1.051 1.19E-05 7.50E-06 15 22

Case B -0.0454 1.78E-03 -1.077 5.03E-05 1.65E-05 7

geometric variations, the effects of back pressure perturbation to performance change are weaker for the turbine blade. Besides, compared with MCS, the absolute deviations of μη and ση are 0.0053% and 1.0E − 05, respectively, even for Case B with only seven samples, demonstrating sufficient response accuracy of the adaptive NIPC. Moreover, as the sample number increases, the response accuracy improves. Compared with Table 3, the relative deviations of statistics and the response accuracy of the adaptive NIPC for the present UQ study with only one uncertain parameter are significantly improved. That is because for any model method, the responded outputs are more accurate with less inputs. In fact, the four-order adaptive NIPC model with tolerance 1.0E − 03 still performs well even if the outlet pressure varies in a wider range [0.25, 0.65], for which the number of CFD evaluation increases to 11 and the model errors are avr = 8.69E − 05 and std = 1.33E − 04. Figure 10 presents the PDF and CDF distributions of adiabatic efficiency change for both MCS and NIPC. The left-skewed PDF illustrates the nonlinear dependence of performance change on the back pressure perturbation. The PDF and CDF obtained from NIPC Case A are quite close to those of MCS, while only slight deviations are visible in the region of large adiabatic efficiency drops for NIPC Case B, which can be found from the CDF distributions. 4.3. Geometric and Back Pressure Variations The effects of geometric variations and back pressure perturbation are separately investigated previously. It is more important to concern on the simultaneous effects of geometric and flow variations to the performance change for turbomachinery blades. The uniform geometric standard deviation is σg = 0.2mm and the standard deviation of pressure perturbation is σp /P0 = 0.05. All the samples used for the UQ studies of geometric variations and back pressure perturbation are used for the present UQ study. The four-order Hermite PCs, as shown by Eq. (28) with p = 6, are used to quantify the adiabatic efficiency change. There are totally six uncertain parameters, five from the geometric variations and one from back pressure perturbation. 23

0.16

0.12

0.12

PDF

PDF

0.16

0.08

0.04

0.00

0.08

0.04

-0.63

-0.23

0.00

0.17

-0.63

-0.23

dh(%)

0.17

dh(%)

(a) Case A

(b) Case B

0.16

MCS

1.0

Case B Case A 0.8

0.12

CDF

PDF

0.6 0.08

0.4

0.04

0.2

0.0 0.00

-0.63

-0.23

-0.63

0.17

-0.23

0.17

dh(%)

dh(%)

(c) MCS

(d) CDF

Figure 10: PDF and CDF distributions of adiabatic efficiency change due to back pressure perturbation

Table 5 gives the statistics of both MCS and NIPC. The tolerance, max of Case A and Case B are 1.0E − 04 and 1.0E − 03, respectively. Compared with MCS, the absolute deviation of μη is 0.0051% and ση keeps unchanged for Case A, while the absolute deviations of μη and ση are about 0.0186% and 1.0E − 05, respectively, for Case B. The results demonstrate that the response accuracy of the adaptive NIPC is sufficient

Table 5: Performance impact of geometric variations coupled with back pressure perturbation

μη /% ση skewη avr std N

MCS -0.1237 3.59E-03 -0.523 – – 5000

Case A -0.1287 3.59E-03 -0.522 5.25E-05 3.71E-05 1878 24

Case B -0.1423 3.58E-03 -0.513 1.82E-04 8.17E-05 186

for the present compounded UQ study. Figure 11 shows the PDF and CDF distributions of adiabatic efficiency change. The NIPC-based PDFs and CDFs agree well with the MCS ones.

0.08

0.08

PDF

0.12

PDF

0.12

0.04

0.00

0.04

-1.43

-1.03

-0.63

-0.23

0.17

0.00

0.57

-1.43

-1.03

-0.63

dh(%)

-0.23

0.17

0.57

dh(%)

(a) Case A

(b) Case B

0.12

MCS

1.0

Case B Case A 0.8

0.08

PDF

CDF

0.6

0.4 0.04 0.2

0.0 0.00

-1.43

-1.03

-0.63

-0.23

0.17

-1.43

0.57

-1.03

-0.63

-0.23

0.17

0.57

dh(%)

dh(%)

(c) MCS

(d) CDF

Figure 11: PDF and CDF distributions of adiabatic efficiency change due to geometric variations and back pressure perturbation

For the compounded UQ study, it is quite interesting to investigate whether the effects of these two kinds of uncertain parameters are linearly or nonlinearly coupled. Table 6 presents the summated and compounded statistic mean values, standard deviations and skewness of adiabatic efficiency change obtained from the direct MCS, where G and P indicate geometric variations and back pressure perturbation, respectively, while GP indicates the geometric variations coupled with back pressure perturbation and G+P indicates the summation of G and P . It is known that the summated statistic mean value of adiabatic efficiency change is slightly deviated from the compounded one. It seems 25

that the effects of back pressure perturbation and geometric variations are linearly imposed. However, the compounded standard deviation and skewness of adiabatic efficiency change are less than the summated ones. Table 6: Comparisons of the compounded and uncoupled effects

G+P GP

μη /% -0.1242 -0.1237

ση 4.89E-03 3.59E-03

skew -1.911 -0.523

In order to investigate whether the previous results are occasionally obtained at one coupled case or not, eight additional compounded cases as shown in Table 7 are studied. The second row gives the reference back pressure and the third row gives the geometric standard deviation. Both the pressure perturbation standard deviation, σp /P0 = 0.05, and the geometric standard deviation are much lower than the threshold ones. The tolerance and polynomial order for adaptive NIPC are 1E − 03 and four, respectively. The statistics of adiabatic efficiency change for the separate and compounded cases and also the linear summation are displayed in Figure 12. The compounded results obtained from the adaptive NIPC are almost the same as those of MCS. Besides, the summated statistic mean values of adiabatic efficiency change are close to the compounded ones for all the nine points, while the linear summations of standard deviation and skew are larger than the compounded ones. The results reveal that the effects of back pressure perturbation and geometric variations on the adiabatic efficiency change are nonlinearly coupled. As the geometric standard deviation increases, the difference of CFD evaluations between the linear and the compounded cases becomes larger, demonstrating that the nonlinear interaction between the two kinds of uncertainty parameters is more significant as the geometric standard deviation increases.

Table 7: Coupled variations at multiple points

Case pb /p0 σ/mm

1 0.35 0.2

2 0.375 0.2

3 0.4 0.2

4 0.35 0.3 26

5 0.375 0.3

6 0.4 0.3

7 0.35 0.4

8 0.375 0.4

9 0.4 0.4

MCS

MCS

(a) Mean

(b) Standard deviation

GP MCS

(c) Skew

(d) CFD Evaluation

Figure 12: Comparisons of the compounded and uncoupled effects at multiple points

5. Conclusion To overcome the curse of dimensionality, the adaptive NIPC method based on an adaptive sparse grid technique is firstly introduced in the study. Compared with the direct MCS and the static NIPC, the response performance including the response accuracy and model construction efficiency are significantly improved for the adaptive NIPC, which are validated and convinced through the function experiments and the performance impact evaluation for a three-dimensional turbine blade. Then uncertainty quantification of aerodynamic performance change due to both the geometric and flow variations using the proposed method are presented. The Gaussian distribution is employed to quantify both the flow and geometric variations, the latter of which is constructed from a series of basis modes of geometric variations extracted by PCA. In order to ensure the 27

sufficiently converged flow solutions, the maximum geometric standard deviation in the present study is much lower than the threshold value (2.1mm) and the minimum back pressure is guarantied to be greater than the threshold value (pb /P0 = 0.25). The effects of two kinds of uncertain parameters, geometric variations and back pressure perturbation, to adiabatic efficiency change of the turbine blade are separately investigated by using the adaptive NIPC. The statistic results show that compared with the nominal blade, the adiabatic efficiency decreases when taking the uncertainty effects into account; moreover, the geometric variations induce a larger efficiency drop for the turbine blade in the study. The statistic skewness and asymmetric PDFs reveal that the adiabatic efficiency change exhibits obvious nonlinear dependence on the geometric and flow variations. Furthermore, statistic evaluations are performed by using the adaptive NIPC for a number of different geometric standard deviations. In opposite to the sensitivity-based evaluation method[17,18], the proposed method can evaluate the adiabatic efficiency change with sufficient accuracy even for significantly larger geometric standard deviation. When the geometric standard deviation is larger than σg = 0.45mm, the polynomial order of adaptive NIPC should be increased to propagate the uncertainty accurately. Meanwhile, the advantage of the adaptive NIPC over the direct MCS will disappear when the geometric standard deviation is too large. Compared with the performance of adaptive NIPC for the geometric UQ studies, the adaptive NIPC model still performs well even if the outlet pressure varies in a wider range. The simultaneous effects of geometric and flow variations to the adiabatic efficiency change are also investigated at multi-point variations. Compared with the linear summation results, although the statistic mean value is almost unchanged, the statistic standard deviation and skew exhibit significant deviations for the compounded case, demonstrating that the effects of these two kinds of uncertain parameters are nonlinearly dependent.

28

Acknowledgements The authors would like to thank The National Natural Science Foundation of China (Grant No. 51676003) and the Fundamental Research Funds for the Central Universities of China (Grant No. 2019QNA4058) for supporting the research work. References [1] A. Stenning, Inlet distortion effects in axial compressors, J. Fluid Eng. 102(1) (1980) 7-13. [2] P. Bry, P. Laval, G. Billet, Distorted flow field in compressor inlet channels, J. Eng. Gas Turbines Power 107(3) (1985) 782-791. [3] T.L. Butler, O.P. Sharma, H.D. Joslny, R.P. Dring, Redistribution of an inlet temperature distortion in an axial flow turbine stage, J. Propul. Power 5(1) (1989) 64-71. [4] D.J. Dorney, K.L. Gundy-Burlet, D.L. Sondak, A survey of hot streak experiments and simulations, Int. J. Turbo Jet Eng. 16(1) (1999) 1-15. [5] K. Bammert, H. Sandstede, Influence of manufacturing tolerances and surface roughness of blades on the performance of turbines, J. Eng. Gas Turbines Power 98(1) (1976) 29-36. [6] R.J. Roelke, J.E. Haas, The effects of rotor blade thickness and surface finish on the performance of a small axial flow turbine, J. Eng. Gas Turbines Power 105(2) (1983) 377-382. [7] K.L. Suder, R.V. Chima, A.J. Strazisar, W.B. Roberts, The effects of adding roughness and thickness to a transonic axial compressor rotor, J. Turbomach. 117(4) (1995) 491-505. [8] V.E. Garzon, D.L. Darmofal, Impact of geometric variability on axial compressor performance, J. Turbomach. 125(4) (2003) 692-703. [9] G.J.A. Loeven, J.A.S. Witteveen, H. Bijl, Probabilistic collocation: an efficient non-intrusive approach for arbitrayily distributed parametric uncertainties, in: 45th AIAA Aerospace Sciences Meeting and Exhibit, Conference Paper AIAA 2007-0317. [10] S. Bhattrai, J.H.S. de Baar, A.J. Neely, Efficient uncertainty quantification for a hypersonic trailingedge flap, using gradient-enhanced kriging, Aerosp. Sci. Technol. 80 (2018) 261-268. [11] Y. Zheng, Z. Qiu, Uncertainty propagation in aerodynamic forces and heating analysis for hypersonic vehicles with uncertain-but-bounded geometric parameters, Aerosp. Sci. Technol. 77 (2018) 11-24. [12] G.J.A. Loeven, H. Bijl, The application of the probabilistic collocation method to a transonic axial flow compressor, in: 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Conference Paper AIAA 2010-2923. [13] R. Pecnik, J.A.S. Witteveen, G. Iaccarino, Uncertainty quantification for laminar-turbulent transition prediction in RANS turbomachinery applications, in: 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Conference Paper AIAA 2011-0660. [14] G. Po¨ ette, D. Lucor, Non intrusive iterative stochastic spectral representation with application to compressible gas dynamics, J. COMPUT. PHYS. 231(9) (2012) 3587-3609. [15] D. Lucor, J. Witteveen, P. Constantine, D. Schiavazzi, G. Iaccarino, Comparison of adaptive uncertainty quantification approaches for shock wave-dominated flows, Summer Program Center for Turbulence Research, Stanford, CA, June, (2012) 219-228. [16] M.M. Putko, P.A. Newman, A.C. Taylor III, L.L. Green, Uncertainty propagation and robust design in CFD using sensitivity derivatives, J. Fluid Eng. 124(1) (2002) 60-69. [17] J. Luo, F. Liu, Performance impact of manufacturing tolerances for a turbine blade using second order sensitivities, in: ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, Conference Paper GT2018-75999. [18] J. Luo, F. Liu, Statistical evaluation of performance impact of manufacturing variability by an adjoint method, Aerosp. Sci. Technol. 77 (2018) 471-484. [19] C. Schillings, S. Schmidt, V. Schulz, Efficient shape optimization for certain and uncertain aerodynamic design, Comput. Fluids 46(1) (2011) 78-87. [20] D.I. Papadimitriou, C. Papadimitriou, Aerodynamic shape optimization for minimum robust drag and lift reliability constraint, Aerosp. Sci. Technol. 55 (2016) 24-33. [21] H.N. Najm. Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics, Annu. Rev. Fluid Mech. 41 (2009) 35-52.

29

[22] S. Hosder, R. Walters, M. Balch, Efficient sampling for non-intrusive polynomial chaos applications with multiple uncertain input variables, in: 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Conference Paper AIAA 2007-1939. [23] H. Zhao, Z. Gao, Y. Gao, C. Wang, Effective robust design of high lift NLF airfoil under multiparameter uncertainty, Aerosp. Sci. Technol. 68 (2017) 530-542. [24] N.P. Gopinathrao, D. Bagshaw, C. Mabilat, S. Alizadeh, Non-deterministic CFD simulation of a transonic compressor rotor, in: ASME Turbo Expo 2009: Power for Land, Sea, and Air, Conference Paper GT2009-60122. [25] T. Ghisu, G.T. Parks, J.P. Jarrett, P.J. Clarkson, Robust design optimization of gas turbine compression systems, J. Propul. Power 27(2) (2011) 282-295. [26] C. Ma, L. Gao, Y. Cai, R. Li, Robust optimization design of compressor blade considering machining error, in: ASME Turbo Expo 2017: Turbomachinery Technical Conference and Exposition, Conference Paper GT2017-63157. [27] T. Ghisu, G.T. Parks, J.P. Jarrett, P.J. Clarkson, Adaptive polynomial chaos for gas turbine compression systems performance analysis, AIAA J. 48(6) (2010) 1156-1170. [28] G. Casella, R.L. Berger, Statistical inference, second ed., Pacific Grove, Duxbury, 2002. [29] V. Barthelmann, E. Novak, K. Ritter, High dimensional polynomial interpolation on sparse grids, Adv. Comput. Math. 12(4) (2000) 273-288. [30] S.A. Smoljak, Quadrature and interpolation formulae on tensor products of certain function classes, Dokl.akad.nauk Sssr 4(5) (1963) 1042-1045. [31] G.W. Wasilkowski, H. Woniakowski, Weighted tensor product algorithms for linear multivariate problems, J. Complexity. 15(3) (1999) 402-447. [32] F. Xiong, S. Greene, W. Chen, Y. Xiong, S. Yang, A new sparse grid based method for uncertainty propagation, Struct. Multidiscip. O. 41(3) (2010) 335-349. [33] H.J. Bungartz, M. Griebel, Sparse grids, Acta Numerica. 13 (2004) 147-269. [34] X. Ma, N. Zabaras, An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations, J. Comput. Phys. 228(8) (2009) 3084-3113. [35] M.Griebel, Adaptive sparse grid multilevel methods for elliptic PDEs based on finite differences, Computing. 61(2) (1998) 151-179. [36] A. Klimke, B. Wohlmuth, Algorithm 847:Spinterp: piecewise multilinear hierarchical sparse grid interpolation in MATLAB, Acm T. Math. Software. 31(4) (2005) 561-579. [37] H.J. Bungartz, Finite Elements of Higher Order on Sparse Grids, Habilitation Tu Munchen Institut Fur Informatik. (1998). [38] N. Wiener, The Homogeneous Chaos, Am. J. Math. 60(4) (1938) 897-936. [39] R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag. 224 (1991). [40] D. Xiu, G.E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24(2) (2002) 61944. [41] J.P. Boyd, The Rate of Convergence of Hermite Function Series, Math. Comput. 35(152) (1980) 1309-1316.

30