University of Southern California Dissertations and Theses

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

(USC Thesis Other)

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DesignOptimizationunderUncertaintyfor RotorBladesofHorizontalAxisWindTurbines by Vahid Keshavarzzadeh Viterbi School of Engineering University of Southern California Los Angeles, California 90089-2531, USA vkeshava@usc.edu 1 ABSTRACT This dissertation is motivated by challenges involved in the design optimization under uncertainty for the structures undergoing fluid structure interaction phenomena. In particular, this study focuses on the aerodynamic simulation and optimization under uncertainty for rotor blades of Horizontal Axis Wind Turbines (HAWT). To that end, a reliability-based design framework integrating the aeroelastic computation based on the blade element method with a stochastic optimization, is developed. The wind model is a stochastic process which accounts for the variation in time and space. The aeroelastic calculation determines the aerodynamic forces on the blade. The uncertainty propagated in the aerodynamic forces stems from the uncertain characteristics of the wind. Another source of uncertainty is introduced in the material property of the blade which results in random structural parameters such as stiffness and capacity. These sources of uncer- tainty lead to a random objective function and constraints. The objective function is the generated power subjected to the constraints on stresses in the blade. The reliability is defined based on the probability of the safety in the system. The design variables are considered to be shape parameters, namely twist angles along the blade. In order to show the effect of uncertainty in the design, the result of the deterministic optimization is obtained and compared with the optimization results under uncertainty. Polynomial chaos expansion is used in various parts of this work as a mathematical tool for the uncertainty propagation. In particular, polynomial chaos expansion can be used to char- acterize the solution of stochastic differential equations such as the governing equation 2 of the motion of the blade. Significant computational challenges are associated with the polynomial chaos expansion approach for problems with high dimensional parameter space. This challenge is addressed in the present dissertation by developing a methodol- ogy for convergence acceleration of polynomial chaos-based stochastic Galerkin solu- tions. Specifically, Nonlinear sequence transformations are adapted to these expansions, viewed as a one-parameter family of functions with the parameter being the polynomial degree of the expansion. Stochastic Galerkin approach yields polynomial chaos repre- sentations that have the requisite analytical properties to ensure suitable convergence of these nonlinear sequence transformations. In particular, the properties of Shanks and Levin transformations are explored in the context of a stochastic initial value problem and a stochastic elliptic problem. Keywords: Aerodynamic optimization; Blade element method; Reliability based design optimization; Polynomial chaos expansion; Convergence acceleration; Stochastic dif- ferential equation; Nonlinear sequence transformation 3 Acknowledgements First and foremost, I would like to express my sincere gratitude to my advisor Prof. Roger Ghanem for his support, patience, insightful suggestions and critical comments throughout the course of my PhD. I would like to thank him for providing me the oppor- tunity to work in a motivational scientific environment where made me enthusiastic about academia. His high academic standards will always be a source of inspiration for me. I would like to thank my dissertation committee members, Prof. Sami Masri and Prof. Paul Newton for their valuable comments and suggestions for improvement of this work. I would also like to thank Prof. Jiin-Jen Lee and Prof. Carter Wellford for serving on my qualifying committee and providing insightful suggestions regarding my research. I gratefully acknowledge the financial support received from US Department of Energy, National Science foundation during the course of this work, in particular, I would like to express my appreciation to the Graduate School at USC for offering me the Provost Fellowship. I would like to extend my sincere thanks to my former and cur- rent colleagues, friends who supported directly or indirectly during this work. I would like to thank Miguel Ricardo Hernandez Garcia, Dr. Ramakrishna Tipireddy and Dr. Hadi Meidani for what I have learned from them and their continuous help and sup- port throughout the course of my study. I would also like to thank my other past and current fellow group members, Dr. Mohammadreza Jahanshahi, Dr. Arash Noshadra- van, Dr. Reza Jafarkhani, Dr. Ali Bolourchi, Dr. Hamed Haddadzadegan, Dr. Maarten Arnst, Dr. Evangelia Kalligiannaki, Dr. Bedrick Sousdik, Dr. Maud Comboul, Armen Derkevorkian, Iman Yadegaran, Charanraj Thimmisetty, Ruda Zhang, Nima Jabbari and Matthew Sexton for all the discussions and sharing thoughts on interesting problems. I am indebted to several people who extended their helps in different aspects through- out my study at USC. It is too numerous to enlist them. I am very grateful for their 4 invaluable friendship and support and send my warmest wishes to all of them. I am also grateful to all my teachers who taught me and inspired me throughout my student life. Lastly, I would like to reserve my special thanks for my family. Without them this work would not have been completed. My very deep gratitude is due to my father Ali who taught me how to be patient and hopeful, my mother Akram for her endless kindness and unconditional love, and my sister Raheleh for her support and motivations. Thank you. 5 Dedication To my parents 6 Contents ListofTables 9 ListofFigures 11 1 Introduction 14 2 UncertaintyAnalysis&ConvergenceAcceleration 19 2.1 Uncertainty Propagation . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1.1 Stochastic Initial Value Problem . . . . . . . . . . . . . . . . . 22 2.1.2 Stochastic Elliptic Problem . . . . . . . . . . . . . . . . . . . . 26 2.2 Sequence Transformation . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.1 Shanks Transformation . . . . . . . . . . . . . . . . . . . . . . 31 2.2.2 Levin Transformation . . . . . . . . . . . . . . . . . . . . . . . 35 2.2.3 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . 39 2.3 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3 WindTurbineModel 60 3.1 Wind model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2 Aeroelastic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.3 Structural model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.4 Uncertainty propagation . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.5 Design optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.5.1 Deterministic optimization . . . . . . . . . . . . . . . . . . . . 89 3.5.2 Optimization under uncertainty . . . . . . . . . . . . . . . . . 92 4 AnumericalstudyofdesignoptimizationofHAWTrotorblade 97 4.1 Wind Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.2 Aeroelastic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Structural Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.4 Convergence Acceleration . . . . . . . . . . . . . . . . . . . . . . . . 108 4.5 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.6 Design optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7 4.6.1 Deterministic Optimization . . . . . . . . . . . . . . . . . . . . 118 4.6.2 Optimization under Uncertainty . . . . . . . . . . . . . . . . . 122 5 SummaryandFuturework 128 ReferenceList 131 8 ListofTables 2.1 Stochastic Initial Value Problem: Relative error in the mean ofu(t;) . 42 2.2 Stochastic Initial Value Problem: Computational effort in terms of esti- mated number of operations . . . . . . . . . . . . . . . . . . . . . . . 47 2.3 Stochastic Initial Value Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with ( 2 6 1)( 2 8 1) applyinge 2 1 transformation . . . . . . . . . . . . . . . . . . 48 2.4 Stochastic Initial Value Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with ( 2 6 1)( 2 8 1) applyinge 1 transformation . . . . . . . . . . . . . . . . . . 48 2.5 Stochastic Initial Value Problem: Range of polynomial degrees in which computational efficiency is gained . . . . . . . . . . . . . . . . . . . . 49 2.6 Stochastic Elliptic Problem: Relative error in the mean ofu(x;) . . . . 52 2.7 Stochastic Elliptic Problem: Computational effort in terms of number of operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.8 Stochastic Elliptic Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with 1 2 (3 2 10 1) applyingv 1 transformation . . . . . . . . . . . . . . . . . . . . . . . 56 2.9 Stochastic Elliptic Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with 1 2 (3 2 10 1) applyinge 1 transformation . . . . . . . . . . . . . . . . . . . . . . . 56 4.1 Mean and standard deviation of the random coefficients obtained from the projection of stochastic modes on deteminisitc eigenmodes . . . . . 105 4.2 Mean and standard deviation of the random eigenvalues corresponding to the first three modes of the blade . . . . . . . . . . . . . . . . . . . . 107 4.3 Convergence acceleration results of the equation of motion of the blade: Relative error in the mean of the first mode solutionX (1) (t; m ; w ) . . 110 4.4 Computational effort in terms of estimated number of operation . . . . 110 4.5 Values of the root mean square of the loading parameters . . . . . . . . 113 4.6 Values of the sensitivity indices corresponding to different sources of uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 9 4.7 Deterministic optimization results . . . . . . . . . . . . . . . . . . . . 120 4.8 Optimization under uncertainty results . . . . . . . . . . . . . . . . . . 125 10 ListofFigures 2.1 Stochastic Initial Value Problem: Mean, standard deviation and coefficient of variation ofu(t;) fors = 10 . . . . . . . . . . . . . . . . . . . . . . . . 41 2.2 Stochastic Initial Value Problem: Relative error in the mean ofu(t;) . . . . 43 2.3 Stochastic Initial Value Problem: (top) Relative error in polynomial chaos coef- ficient associated with 10 , (bottom) Relative error in polynomial chaos coeffi- cient associated with 4 7 . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.4 Stochastic Initial Value Problem: (top) Relative error in polynomial chaos coef- ficient associated with ( 2 5 1) 7 , (bottom) Relative error in polynomial chaos coefficient associated with ( 2 6 1)( 2 8 1). . . . . . . . . . . . . . . . . 45 2.5 Stochastic Initial Value Problem: Relative error in the mean ofu(t;) versus computational effort defined as operation count . . . . . . . . . . . . . . . 48 2.6 Stochastic Elliptic Problem: Relative error in the mean ofu(x;) . . . . . . . 52 2.7 Stochastic Elliptic Problem: (top) Relative error in polynomial chaos coeffi- cient associated with 1 , (bottom) Relative error in polynomial chaos coeffi- cient associated with 1 2 (3 2 10 1). . . . . . . . . . . . . . . . . . . . . . 53 2.8 Stochastic Elliptic Problem: (top) Relative error in polynomial chaos coeffi- cient associated with 1 2 (3 2 5 1) 2 , (bottom) Relative error in polynomial chaos coefficient associated with 1 2 (5 3 4 3 4 ) 9 . . . . . . . . . . . . . . . 54 2.9 Stochastic Elliptic Problem: Relative error in the mean ofu(x;) versus com- putational effort defined as operation count . . . . . . . . . . . . . . . . . 55 3.1 Illustration of annular elements in BEM method . . . . . . . . . . . . . . . 67 3.2 Input and output of BEM method . . . . . . . . . . . . . . . . . . . . . . 68 3.3 Definition of lift and drag . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.4 Velocity and pressure up- and downstream of the rotor . . . . . . . . . . . . 70 3.5 Control volume around a wind turbine . . . . . . . . . . . . . . . . . . . . 70 3.6 Velocities at rotor plane . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.7 Normal and tangential forces . . . . . . . . . . . . . . . . . . . . . . . . 74 3.8 Reliability based design optimization of the HAWT blade . . . . . . . . . . 96 4.1 Eigenvalues of the covariance matrixR . . . . . . . . . . . . . . . . . . . 99 11 4.2 Spatial-temporal eigenmodes of the Karhunen-Loeve expansion of the wind process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Histograms of the random variables associated with the eigenmodes in Karhunen- Loeve expansion of the wind process . . . . . . . . . . . . . . . . . . . . 101 4.4 Distribution of the aerodynamic induction factors at 0.8 of the blade radius and t = 2:5sec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.5 Distribution of the aerodynamic forces on the section at 0.8 of the blade radius ndt = 2:5sec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.6 Distribution of the generated power . . . . . . . . . . . . . . . . . . . . . 103 4.7 Square of the natural periods of the blade corresponding to the mean value of the elastic modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.8 Histograms of the random coefficients obtained from the projection of stochas- tic modes on deteminisitc eigenmodes . . . . . . . . . . . . . . . . . . . . 106 4.9 Histograms of the random eigenvalues corresponding to the first three modes of the blade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.10 Histograms of the Generalized force in three modes . . . . . . . . . . . . . 108 4.11 Histograms of the solution of the governing equation in three modes . . . . . 109 4.12 Convergence acceleration results of the equation of motion of the blade: Rela- tive error in the mean of the first mode solutionX (1) (t; m ; w ) . . . . . . . 111 4.13 Relative error in the mean of the solution versus computational effort defined as operation count . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.14 Convergence acceleration results of the equation of motion of the blade: (a) Relative error in polynomial chaos coefficient associated with w 1 , (b) Rela- tive error in polynomial chaos coefficient associated with w 2 in the first mode solutionX (1) (t; m ; w ) . . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.15 Convergence acceleration results of the equation of motion of the blade: (a) Relative error in polynomial chaos coefficient associated with w 3 , (b) Rela- tive error in polynomial chaos coefficient associated with w 4 in the first mode solutionX (1) (t; m ; w ) . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.16 Convergence acceleration results of the equation of motion of the blade: (a) Relative error in polynomial chaos coefficient associated with w 5 , (b) Rela- tive error in polynomial chaos coefficient associated with m 1 in the first mode solutionX (1) (t; m ; w ) . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.17 Convergence acceleration results of the equation of motion of the blade: (a) Relative error in polynomial chaos coefficient associated with m 2 , (b) Rela- tive error in polynomial chaos coefficient associated with m 3 in the first mode solutionX (1) (t; m ; w ) . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.18 Convergence acceleration results of the equation of motion of the blade: (a) Relative error in polynomial chaos coefficient associated with 2 w 1 , (b) Relative error in polynomial chaos coefficient associated with w 2 m 1 in the first mode solutionX (1) (t; m ; w ) . . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.19 Distribution of the solution of the governing equation at the critical time . . . 118 12 4.20 Three dimensional blade model corresponding to the initial design variables . 119 4.21 Three dimensional blade model corresponding to the final deterministic design 122 4.22 Twist values of initial and final design corresponding to deterministic optimization123 4.23 Distribution of the generated power of the initial and final design . . . . . . . 124 4.24 Distribution of the stress and limit state function of the initial and final design . 125 4.25 Three dimensional blade model corresponding to the final design optimization under uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.26 Twist values of initial and final design corresponding to optimization under uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 13 Chapter1 Introduction The mechanical response of a wind turbine is based on the complex interaction of dif- ferent components. Wind turbines also operate in the hostile environment that should be designed for fluctuating loads due to the nature of the wind. Understanding the inter- action of the wind with the structure is the key parameter for a reliable and optimized design. The increased size of the rotor requires more comprehensive computational procedures to account for the dynamic response of the system as well as the effects of uncertainty. Efforts have been devoted to simulation of wind turbines considering the aeroelastic effects. Aeroelastic analysis has the aerodynamic part to describe the load- ing from the wind and a structural part to describe the dynamic response. There are methods with different levels of complexity for determining aerodynamic loads. The most common method for describing the aeroelastic behavior of the wind turbine is the Blade Element Momentum method (BEM) [36], [1], [35], [63]. However, there are computational works based on Navier-Stokes equations for simulation of the wind flow in wind turbine response. These models are basically used to determine the response of the wind turbine structure [7], [6], [37]. There are different ways to describe the structural model of the blade. The main purpose in structural model is to determine the structural load dynamically. One way is using the principle of virtual work and beam theory which provides an effective computational way to describe the dynamics of the blade [47], [16], [58], [52]. The other way is to employ full finite element approach which allows a more detailed state representation of the wind turbine blade [56]. In addition to the efforts for dynamic simulation of the wind turbine, the aerodynamic 14 design of the blade has also been studied. The aeroelastic design of rotors for horizontal axis wind turbines is a multi-disciplinary task that should integrate the dynamic analysis of the system for determining the loads and the optimization procedure for the maximum energy production. Considerable efforts have been devoted to aerodynamic and struc- tural optimization [43], [54], [53], [65], [23]. The deterministic optimization has been applied to reduce the cost associated with design and maximize the efficiency in these studies while the existence of the uncertainties in every engineering systems requires an optimization framework that incorporates uncertainty. A proper design framework should consider the uncertainties, as they may significantly change the performance and reliability of the final design obtained with deterministic design procedure. The applica- tion of design under uncertainty approaches ensures that the system will be performing with a desired reliability within certain limits. For optimization under uncertainty, two main approaches, reliability based design optimization and robust design optimization are taken. In the reliability based design optimization the basic objective is to design a system that satisfy certain performance during its lifetime [22], [25], [38], [48]. The objective performance is not necessarily deterministic due to the inherent random nature of different parameters in the system such as loading condition or structural parameters or condition of operation. Therefore, reliability based design optimization is concerned with quantitative estimation of the effects of uncertainty in terms of failure probabilities or expected values. The formulation of RBDO requires the definition of design variables and random parameters, the failure criterion of the system and the objective function. Another approach, robust design optimization determines the solution that has the least sensitivity to variation in the system parameters [10], [39], [45], [49]. The solution of the optimization problem under uncertainty can be sensitive to the variation from dif- ferent sources and the performance of the system could vary significantly. Therefore, the objective in robust optimization is to obtain a solution with the least variation in 15 parameters. In this study, the main objective is to incorporate the effect of uncertainty in the design optimization of the wind turbine blade. Since the safety and reliability of the system is the significant objective in the design of the wind turbine, the reliability based design optimization is considered in this study. A computational framework is developed for the design optimization under uncertainty. The mathematical framework utilizes polynomial chaos expansion for uncertainty modeling [30], [31]. The polynomial chaos expansion approach can be applied to the problems where the uncertain parameters are represented by random fields or stochastic processes. This approach provides a general method to characterize a solution for the problems with uncertain parameters. The characterized solution as an expansion can then be used to obtain different quantities that form the formulation of the optimization problem. Poly- nomial chaos expansion is an attractive approach due to its fast convergence property and its ability to provide a functional representation of the response of the system sub- jected to random input parameters. Coefficients in these expansions have been approx- imated based on various weighted residual criteria that have included Galerkin projec- tion [29, 62, 68], collocation [24, 33, 42], and least squares [9]. Variations on these methods have been categorized as intrusive or non-intrusive, depending on whether the specific implementation can use current deterministic analysis tools as black-boxes [26, 28, 33] or whether they require the solution of new adapted equations [29, 32]. One challenge with the polynomial chaos formulation pertains to the characterization of the solution as a polynomial in high-dimension and the ensuing curse of dimension- ality. The number of terms in these representations indeed grows factorially both as the order of polynomial approximation and as its dimension. Efficiency and error analy- sis issues related to these problems continue to be addressed by a growing number of researchers [4, 8, 19, 41, 50, 61, 62, 66]. The focus in many of these methods has been on adapting the approximating stochastic bases to the problem at hand by relying on 16 various error estimates. The present dissertation adopts a slightly different perspective on this challenge that augments current capabilities and pursues a methodology that permits gains in accuracy without modifications to the approximation basis. Specifi- cally, the focus is on nonlinear sequence transformations to significantly improve the accuracy of pre-computed approximations. This can be viewed as a nonlinear adapta- tion of Richardson Extrapolation [11, 12, 51, 67]. Thus, starting with two polynomial chaos approximations obtained at polynomial orders p 1 and p 2 > p 1 , respectively, a new expansion is synthesized that has comparable accuracy to an orderp > p 2 expan- sion. The suitability of various sequence transformations for convergence acceleration is closely tied to the convergence rate of polynomial chaos approximations, as func- tion of polynomial order, attainable with the stochastic Galerkin approach. In order to take advantage of this behavior, polynomial chaos approximations of a given stochastic equation are viewed as a one-parameter family of sequences, namely the coefficients in the expansion as function of the polynomial degree of the expansion. The approach is demonstrated using, specifically, the Shankse-transformation [57] and the Levinv- transformation [44] that are both adapted to the convergence rates associated with the stochastic Galerkin approach. The present dissertation is organized as follows. Chapter (2) introduces the uncer- tainty analysis approach in this study. The stochastic Galerkin formulation is intro- duced in the context of stochastic initial value and stochastic boundary value problems. The fundamental ideas underlying sequence transformations, focusing in particular on Shanks and Levin transformations are reviewed. The benefits of applying sequence transformations to stochastic problems are then demonstrated using two numerical exam- ples. The details of a wind turbine model consists of three main components, wind model, aeroelastic model and structural model are explained in Chapter (3). Follow- ing the details of wind turbine model, the design optimization framework is introduced 17 for both cases of deterministic optimization and optimization under uncertainty in this chapter. Numerical results of the optimization problem in this study for both deter- ministic optimization and optimization under uncertainty are presented in Chapter (4). Numerical results corresponding to different components of the wind turbine model are also presented in this chapter. Chapter (5) concludes this dissertation and discusses the future work. 18 Chapter2 UncertaintyAnalysis&Convergence Acceleration In this section, we explain the details of the uncertainty analysis method based on the polynomial chaos expansion used in the wind turbine model. The computational chal- lenges associated with this method is addressed by pursuing a new approach on conver- gence acceleration. The details of the new approach are presented in this chapter. The applicability of the approach is demonstrated on two numerical examples of a stochastic initial value problem and a stochastic boundary value problem. 2.1 UncertaintyPropagation Let an experiment that defines the physical context of the problem be described by the probability triple ( ;F;P). Real-valued random variables relevant to this experiment are measurable functions 7! R. Any second-order random variableU(!) : 7! R can be expressed as a convergent series in terms of Hermite polynomials in Gaussian variables [14], U(!) = u 0 H 0 + 1 X i 1 =1 u i 1 H 1 ( i 1 (!)) + 1 X i 1 =1 i 1 X i 2 =1 u i 1 i 2 H 2 ( i 1 (!); i 2 (!)) + 1 X i 1 =1 i 1 X i 2 =1 i 2 X i 3 =1 u i 1 i 2 i 3 H 3 ( i 1 (!); i 2 (); i 3 (!)) +:::; (2.1) 19 whereH n ( i 1 ;:::; in ) is a Hermite polynomial in terms of independent standard Gaus- sian random variables ( i 1 ;:::; in ). The equation above can be re-written in a more notationally convenient form as follows, U(!) = 1 X j=0 ^ u j j (!): (2.2) where j (!) = j ((!)) and(!) is the vector of independent standard Gaussian ran- dom variables ((!) = i 1 ;:::; in ) that is mentioned above. It should be noted that, there is a one to one correspondence between j (:) andH n (:) and also the coefficients associated with these polynomials. In practice, the above expansion is truncated and the random variable is expressed with finite number of polynomials. As mentioned earlier, any second-order random variable (random variable with finite variance) can be expressed with expansion given in equation (2.1). Therefore, this expansion can be applied for most physical processes. Moreover, an important property of these poly- nomials is their orthogonality with respect to the Gaussian density function that can be written as, E[ i j ] =E[ 2 i ] ij ; (2.3) where ij is a Kronecker delta. It should be noted that, if the underlying random source is non-Gaussian the polynomial chaos expansion is written with polynomials that are orthogonal with respect to the probability density of the underlying random variable. StochasticGalerkinProjection In this section the stochastic Galerkin projection is briefly described by highlighting underlying mathematical constructs that motivate the search for the computational effi- ciency provided by nonlinear sequence transformation methods. Two examples of stochas- tic initial value problem and stochastic boundary value problem are considered. The 20 existing results that substantiate the applicability of specific sequence transformation methods that are adapted to particular asymptotic behaviors of the error are also described. The deterministic domain over which the experiments (described by the probabil- ity triple ( ;F;P)) unfold consists of either space or time or space-time, and will be denoted byDR d whered denotes the dimension of the embedding space (1, 2, 3 or 4). A coefficient in the differential equation will be taken to be random, and denoted by k(x;!) with x2 D and !2 . Randomness in this parameter will be assumed to derive from other random variables, considered to be more basic and indeed to be the fundamental (also called “germ”) random variables [62] associated with this exper- iment. In other words, it is assumed thatF is generated by this germ, which is denoted by and assumed it has values inR s . It is further assumed that is square-integrable, and thus2 L 2 (( ;F;P);R s ). Thek(x;!) will henceforth be replaced withk(x;), ands is referred to as the stochastic dimension of the problem. In the construction of polynomial chaos expansion in this work, it is supposed that the random variablesf i g are statistically independent. Dependence of the parameterk onx2D signifies thatk is a stochastic process indexed byD. The solution,u(x;!), of the stochastic equation, is then a random process whose randomness also derives from and will thus be re-written asu(x;). A polynomial chaos decomposition ofu(x;) is then sought in the form of a multidimensional polynomial with respect to the random vector. The polynomial basis in this decomposition is denoted byf i ()g or simply byf i g. The polynomial chaos expansion is construed as a representation of a stochastic variable that is parameterized by the order of the polynomial approximation. The com- puted coefficients in this expansion will be considered as an approximation, converg- ing to their respective exact values as the polynomial order increases. In other words, 21 different values of a polynomial chaos coefficient corresponding to different polyno- mial orders in the expansion, form a sequence converging to an asymptotic value. The computational cost to estimate the polynomial chaos coefficients grows quickly with the stochastic dimension thus motivating the acceleration of their convergence that is addressed in this section. The next two sections describe the stochastic Galerkin pro- jection for stochastic initial value problem and stochastic boundary value problem that gives rise to a polynomial approximation to the solution of equations with random coef- ficients. The Galerkin scheme induces convergence rates in these approximations, as function of the polynomial order, that makes them suitable candidates for convergence acceleration via sequence transformations. 2.1.1 StochasticInitialValueProblem Consider a stochastic ordinary differential equation of the form, d dt u(t;) =k(t;)u(t;); u(0;) =u 0 ; 8t2D; a:s: in ; (2.4) whereD is a bounded domain and the coefficientk(t;) has the following decomposi- tion, k(t;) = L X i=0 k i (t) i (); (2.5) with known coefficientsk i (t). These coefficients are typically estimated using adapted statistical estimation techniques which may involve an inverse problem [17, 18, 27, 55, 60]. The polynomial chaos expansion of the solution can be expressed in a similar form as, u(t;) = N X i=0 u i (t) i () +" N ; (2.6) 22 where" N denotes a truncation error. The polynomial basisf i g is taken to be orthogo- nal with respect to the probability density of, which can be expressed as h i j i =E( i j ) = Z R s i (z) j (z)(z)dz =h 2 i i ij (2.7) whereh:i andE(:) denote the operator of mathematical expectation, is the probabil- ity density function of relative to Lebesgue measure, and ij is the Kronecker delta function. Substituting the expansion of the solution and the random parameter in equa- tion (2.4) results in, N X i=0 d dt u i (t) i = L X i=0 N X j=0 k i (t)u j (t) i j +e N ; (2.8) where an error is incurred due to the finite-dimensional representation of the solutionu. The stochastic Galerkin solution is obtained by requiring the above error to be orthogo- nal to the approximating space resulting in, N X i=0 d dt u i (t)h i k i = L X i=0 N X j=0 k i (t)u j (t)h i j k i; k = 0; 1;:::;N; (2.9) which is a deterministic equation that can be solved foru i (t); i = 0; ;N. This last equation can be re-written in the compact form, d dt U =A(t)U; 8t2D; (2.10) whereA(t)2R (N+1)(N+1) is a square matrix with elements given by, A kj (t) = 1 h 2 k i L X i=0 k i (t)h i j k i; j;k = 0; 1;:::;N: (2.11) 23 In many problems of interest, k(t;) is a linear form in, in which case the quantity h i j k i in equation (2.11) is simplified toh i j k i. This in turn significantly simpli- fies the sparsity structure of matrix A. Given the importance of this sparsity structure on the computational burden, we will make an explicit distinction in the remainder of the paper between linear and nonlinear parameter representations associated withh i j k i andh i j k i, respectively. As previously mentioned, equation (2.10) is a deterministic equation and is solved by a standard fourth-order Runge-Kutta method. We note that the exact analytical solution of this stochastic ordinary differential equation has the general form, u ex (t;) =u 0 e R t 0 k(;)d : (2.12) The exact coefficients in the polynomial chaos expansion can thus be expressed by pro- jecting this exact solution on the respective polynomials resulting in, u i (t) = 1 h 2 i i Z R s u ex (t;z) i (z)(z)dz i = 0; ;N: (2.13) The number of termsN + 1 in the polynomial decomposition in equation (2.6) depends both on the sizes of the stochastic germ, and on the polynomial degreep used in the approximation, with the relationship specified by the following equation, N + 1 = (p +s)! p!s! : (2.14) It is thus clear that for a fixed value of s, accuracy gains relative to p have signifi- cant computational consequences. The present work develops methods to accelerate the convergence of the polynomial chaos coefficients to their asymptotic values associated with p!1. The solution of equation (2.9) is thus viewed as a function of parame- ter p, which is considered the one-parameter family of such approximations. To each 24 coefficient in the polynomial chaos decomposition of the solution is thus associated a sequence that spans this family. In the one-dimensional case (s = 1), the number of coefficients grows in a one-to-one correspondence with the polynomial degree of the approximation, allowing a simple visualization shown in the following equation, where two such sequences have been individually marked, (2.15) For multiple stochastic dimensions, the mapping from p to N is one-to-many and the diagram is therefore more complicated. In other words, the increase in the number of coefficients as the polynomial order is increased from p to p + 1 is no longer one and indeed depends on the value ofp. In this case algorithmic book-keeping is required to keep track of the coefficients as they get created and updated as function of p. To each polynomial chaos coefficient, a sequence that describes the changes in that coef- ficient as p increases, is associated. The task then, in this work, will be to transform these sequences into new sequences that are closer to their respective asymptotic values. This process can be seen in the diagram shown below where only the first row of coef- ficients are obtained by solving the stochastic Galerkin equations, while the subsequent rows are obtained from the first row through the repetitive application of the sequence transformationsT , 25 u (1) 0 ;u (2) 0 ;:::;u (n) 0 + T u 0 (1) 0 ;u 0 (2) 0 ;:::;u 0 (n 0 ) 0 + T u 00 (1) 0 ;u 00 (2) 0 ;:::;u 00 (n 00 ) 0 . . . (2.16) 2.1.2 StochasticEllipticProblem Consider a stochastic elliptic boundary value problem consisting of finding a function u(x;) : D ! R, such that the equation below holds almost surely in [5], r(k(x;)ru(x;)) =f(x); x2 D u(x;) = 0; x2@D; (2.17) where D is an open bounded domain and the stochastic parameter satisfies, 0<k min k(x;)k max a:e: inx; a:s: in ; (2.18) for some deterministic scalars k min ;k max . Although the foregoing analysis can be extended to the setting wherek max depends on [46], it is assumed in this study that k max is independent ofx and. In addition, the right-hand side in equation (2.17) satis- fies, kfk L 2 (D) < +1: (2.19) Here again s-dimensional stochastic germ, with values in R s is considered. The numerical approximation of stochastic functions can concisely be formulated by using a tensorization of functions on D with functions on . Let the function space H 1 0 (D) 26 be the subspace of H 1 (D) consisting of functions which vanish at the boundary @D of D, and let L 2 ( ) be the space of random variables with finite variance defined on . Then the space H = H 1 0 (D) N L 2 ( ) is the completion of formal sumsu(x;) = P i;j v i (x)w j (), wherefv i g H 1 0 (D) andfw j g L 2 ( ). Considering the space H equipped with the inner product, (u;v) = Z Z D rurv(z)dxdz =E Z D rurvdx ; (2.20) let the bilinear formB : H H! R be defined as, B(u;v) =E Z D krurvdx 8u;v2 H; (2.21) and the linear formL : H! R be given by, L(v) =E Z D fvdx 8v2 H: (2.22) The variational formulation for equation (2.17) consists then of findingu2 H such that, B(u;v) =L(v) 8v2 H: (2.23) It should be noted that condition (2.18) ensures the continuity and coercivity of the bilinear form; hence by Lax-Milgram lemma [40], the solutionu2 H to the above vari- ational problem exists and is a unique solution. This solution can then be characterized as follows, u(x;) = M X i=1 N X j=0 u ij i (x) j (); (2.24) wheref i g is a piecewise linear continuous function, such as standard finite element shape functions constituting a basis in a finite dimensional subspace of H 1 0 (D) andf j g 27 is an orthogonal polynomial in L 2 ( ). With some abuse of notation, equation (2.24) will be referred to as a polynomial chaos decomposition ofu(x;). The Galerkin approxi- mation of the solution is obtained by employing test functions given by, v(x;) = k (x) l (); k = 1;:::;M; l = 0;:::;N: (2.25) Consequently, equation (2.23) yields a linear system of equations in the form, M X i=1 N X j=0 E j l Z D kr i r k dx u ij =E l Z D f k dx k = 1;:::;M; l = 0;:::;N; (2.26) the solution of which are the coefficients u ij . Expanding k(x;) according to the approximation, k(x;) = L X i=0 k i (x) i (); (2.27) equation (2.26) becomes, M X i=1 N X j=0 L X m=0 E [ m j l ] Z D k m (x)r i r k dx = E [ l ] Z D f k dx ; k = 1;:::;M; l = 0;:::;N: (2.28) Care should be taken to ensure that the approximation given in equation (2.27) satisfies condition (2.18). The system of equations above can be written in the following compact form, AU =b: (2.29) It is noted that for an s-dimensional stochastic dependence and a polynomial approx- imation of the solution of degree p, the value of N is given by equation (2.14). The coefficientsu ij in equation (2.24) converge to their asymptotic value asp!1. The 28 sequence of polynomial chaos coefficients is obtained by increasing the polynomial order in the expansion; therefore, as mentioned earlier, each element in the sequence corresponds to a polynomial order for which a solution has been computed. 2.2 SequenceTransformation LetfA n g =fA 1 ;A 2 ; ;g denote a sequence of real numbers, and letA denote the set of all such sequences. A sequence transformation,T is a mapping fromA ontoA that maps a given sequencefA n g into a new sequencefA 0 n g. The transformationT is said to accelerate the convergence of the sequencefA n g toA if, lim n!1 A 0 n B A n B = 0; A n 2fA n g; A 0 n 2fA 0 n g (2.30) wherefA 0 n g is the image offA n g underT andB is the asymptotic limit of the sequence fA n g. Sequence transformations are categorized as linear or nonlinear. Linear transfor- mations have the following general form, A 0 n = 1 X i=0 a ni A i ; (2.31) where the constant weightsa ni are independent of the sequence elementsA i . Among such linear transformations, the well known Richardson Extrapolation is utilized in Romberg’s method for improving the accuracy of the trapezoid rule [20,51]. Richardson Extrapolation was also studied in the analysis of the numerical solutions of ordinary dif- ferential equations [13,34]. It improves the accuracy by forming a linear combination of 29 two approximations corresponding to two different discretization schemes. Thus, con- sider an approximation of the exact quantityB in terms of a small parameterh (perhaps discretization increment) of the form, A n =B +C h r n p +O(h p 0 ) r> 1; (2.32) where the leading term in the approximation error is known to have the indicated form, and wherefr n ; n = 0;g induces a sequence of refined approximations,fA 0 n g, and the corresponding sequence of errors. In this equation,O(h p 0 ) represents higher order corrections with p 0 > p. The leading error term, C h r n p , can then be eliminated by linearly combining two such approximationsA n andA n+1 , resulting in a new sequence of approximants,fA 0 n g of the form, A 0 n = r p A n+1 A n r p 1 ; (2.33) whereA 0 n =B+O(h p 0 ) is the improved approximation. It should be noted that Richard- son Extrapolation accelerates convergence if r and p are known and the leading error term in the approximation of interest has the form of equation (2.32). In view of this, Richardson Extrapolation has been successfully applied for the acceleration of conver- gence with respect to spatial discretizations of partial differential equations. The key idea underlying Richardson Extrapolation can be generalized to other forms of error terms using nonlinear transformations that adapt to the specific error terms in a given problem. In the sequel, two well known classes of such transformations are reviewed which are particularly adapted to the form of the error induced by the stochas- tic Galerkin projection using polynomial chaos expansions. These are the Shanks and the Levin transformations. 30 2.2.1 ShanksTransformation One of the first attempts at convergence acceleration through nonlinear sequence trans- formation was developed by Shanks [57], where he construed a sequence as consisting of a “base” plus a transient that evolves with the sequence. If the transient decays to zero, then the sequence will converge to a value given by its base. If the transient diverges, then the sequence diverges. Clearly, the transient in equation (2.32) consists of a single term with exponential decay along the sequence. While a linear transfor- mation proved adequate for accelerating the convergence associated with this transient, adapted nonlinear transformations are required for more general forms of the transient. Transients, in this context, should be understood as referring to the error in approximat- ing the base. Improving the accuracy with which the transients are estimated clearly improves our estimate of the limit of the sequence, and sequence transformations can be viewed as a process through which the sequence of errors is iteratively estimated, based on some knowledge of error form. Shanks transformation extends the range of convergence acceleration methods to transients consisting of a finite sum of exponen- tially decaying terms. This form is significantly more robust than the exponential rate associated with Richardson extrapolation. Consider thus a sequence whose elements adhere to the following dynamics, A n =B + k X i=1 a i q n i ; q i = 2f0; 1g: (2.34) The objective here is to find the base, B, of the transient. A challenge is that, while the functional form of the transient can be known from theoretical considerations, the coefficients in the transient are generally not known, and the sequence is not observed sufficiently within its asymptotic regime to ensure that the transients are small enough. Depending on the values ofk,a i andq i , the sequence can be convergent or divergent, 31 monotonic or oscillatory, as function of n. If the conditionfjq i j < 18ig holds, then the sequence converges to B as n!1. Ifjq i j > 1 for some i2f1; ;kg, then fA n g diverges from B. In this latter case, B is called the anti-limit of the sequence fA n g. While efforts have gone into interpreting convergence of sub-sequence to anti- limits [57], we will not be concerned with that aspect of sequence transformations. Knowledge of two terms in the sequence is not sufficient, in the present case, to elim- inate the transient terms in a manner analogous to the procedure followed in developing Richardson Extrapolation. To address this problem, the sequencefA n g is regarded as a locally transient sequence characterized by 2k + 1 successive observations of its values centered aroundn, in the following form, A r =B kn + k X i=1 a in q r in ; nkrn +k; nk; q in = 2f0; 1g: (2.35) The local baseB kn and 2k quantitiesa in ;q in are thus constrained by the 2k+1 nonlinear equations above. These values are local because they correspond to the elementsA r with nkrn +k; nk. This procedure transforms the local sequenceA r into the local sequenceB kr having a faster decaying transient. The system of equations (2.35) is greatly simplified when expressed in terms of the difference of two consecutive values in the sequence (A n =A n+1 A n ) resulting in, A nr =B kn + k1 X s=0 s A nr+s ; 0rk; (2.36) with A nr+s = k X i=1 a in (q s+1 in q s in )q nr in ; k1 X s=0 s (q s+1 in q s in ) = 1; i = 1;:::;k: (2.37) 32 The new system of equations is linear with k + 1 unknowns which include s ; (s = 0;:::;k 1) andB kn . Consequently,B kn can be determined from the system of equa- tions (2.36) by Cramer’s rule, B kn = A nk ::: A n1 A n A nk ::: A n1 A n A nk+1 ::: A n A n+1 . . . . . . . . . A n1 ::: A n+k2 A n+k1 1 ::: 1 1 A nk ::: A n1 A n A nk+1 ::: A n A n+1 . . . . . . . . . A n1 ::: A n+k2 A n+k1 : (2.38) This procedure defines thee k transformation of Shanks as follows, e k (A) n =B kn nk; (2.39) requiring (2k + 1) observations of the sequence to bootstrap it. The iterated transforma- tione m k can be formed bym successive applications ofe k as follows, B kn =e k (A) n nk; e k (e m1 k (A)) n =e m k (A) n nmk: (2.40) 33 e 1 Transformation Of particular interest is the transformation e 1 as it is predicated on the least amount of information required to generate an accelerated sequence. For k = 1 in equa- tion (2.38), Shanks transformation reduces to the well known Aitken’s 2 process [2] which requires three consecutive elements in the sequence and is given as follows e 1 (A n ) = A n1 A n A n1 A n 1 1 A n1 A n = A n+1 A n1 A 2 n A n+1 +A n1 2A n : (2.41) Equation (2.35) evaluated at three consecutive values, corresponding tor =n1;n;n+ 1, yields A n1 =B 1n +a 1n q n1 1n ; A n =B 1n +a 1n q n 1n ; A n+1 =B 1n +a 1n q n+1 1n : (2.42) The relative error can be expressed as A n B 1n B 1n , (here, for relative error, we have r=n). An important point is the form of this error, a 1n B 1n q n 1n , which decreases exponentially withn for convergent sequences (jq 1n j< 1). This type of convergence is analogous to that which is expected in stochastic Galerkin approximations of polynomial chaos coefficients. It is noted that the error form associated with Shanks transformation, indicated in equation (2.34), generalizes the exponential error form in Richardson extrapolation shown in (2.32). The transient in equation (2.34) is quite general, and the development of the Shanks transform reviewed above does not require knowledge of the coefficients in the linear combination of the exponential; the mere knowledge of the suitability of 34 this form is sufficient. Further, a closer inspection of equation (2.42) shows that within any given locally transient sequence, the e 1 transformation presumes an exponential transient with constant decay rate. This is analogous to the error form associated with Richardson extrapolation as indicated in equation (2.32). The analytical significance of Shanks transformation follows from the transformation of the nonlinear equation (2.34) into the linear equation (2.36) by a change of variable fromq to A n . It should also be noted that there is no analytical expression in the form of equa- tion (2.42) that expresses the error in estimating the polynomial chaos coefficients as function of the polynomial order. Such a result only exists for the zeroth order coeffi- cient, and it is postulated that higher order coefficients behave in an analogous manner. 2.2.2 LevinTransformation As stated in the previous section, Shanks transformation is based on the model sequence (2.34) which is characterized as a sum of exponentials. Levin generalized this transformation for sequences with other rates of convergence [44]. By analogy with equation (2.35), k + 1 equations for the sequencefA n g can be written in the following form, A r =B kn +R k (r;n) nrn +k; (2.43) whereB kn is the approximation to the limit of the sequence andR k (r;n) is given by, R k (r;n) =R 1 (r) k1 X i=0 i (n)f i (r) k 1; (2.44) 35 where i (n) are constants andf i (r) are chosen based on convenience. Substituting the expression forR k (r;n) into equation (2.43) and rewritingR 1 (r) asR r results in, A r =B kn +R r k1 X i=0 i (n)f i (r) nrn +k: (2.45) The limitB kn can be obtained from the abovek + 1 equations using Cramer’s rule, B kn = A n A n+1 ::: A n+k R n f 0 (n) R n+1 f 0 (n + 1) ::: R n+k f 0 (n +k) R n f 1 (n) R n+1 f 1 (n + 1) ::: R n+k f 1 (n +k) . . . . . . . . . R n f k1 (n) R n+1 f k1 (n + 1) ::: R n+k f k1 (n +k) 1 1 ::: 1 R n f 0 (n) R n+1 f 0 (n + 1) ::: R n+k f 0 (n +k) R n f 1 (n) R n+1 f 1 (n + 1) ::: R n+k f 1 (n +k) . . . . . . . . . R n f k1 (n) R n+1 f k1 (n + 1) ::: R n+k f k1 (n +k) (2.46) For the particular case of f i (r) = r i , the determinants in B kn can be expressed in terms of Vandermonde determinants assumingR n 6= 0;8n. Dividing each column by R n ;R n+1 ;:::;R n+k and expanding the determinants by the first row results in the rep- resentation of Levin transformation as a ratio of two finite sums as follows, B kn = k X j=0 (1) j k j n +j n +k k1 A n+j R n+j k X j=0 (1) j k j n +j n +k k1 1 R n+j : (2.47) 36 The Levin transformation depends on the sequencefR n g, which is referred to as the remainders. If these remainders are directly obtained from the elements of the sequence fA n g, the Levin transformation is a nonlinear sequence transformation since each remain- derR n is explicitly dependent on at least one element of the sequencefA n g. Levin also suggested different remainder estimates that lead to different forms of Levin transfor- mation [44]. In particular, Levin’sv transformation, associated with a particular model for the remainder sequence, is utilized in the present paper. Rewriting elements of the sequence in terms of differences, in the form, A n = n X i=0 d i n = 0; 1;::: (2.48) whered n+1 = A n =A n+1 A n andd 0 =A 0 , the remainders for Levin’sv transfor- mation are given by, R n = d n d n+1 d n+1 d n : (2.49) With this model for the remainders, the sequencefB kn g is referred to as thev-transform of sequencefA n g. It should be noted that the remainder estimate in the v-transform relies implicitly on thee 1 transform of Shanks as an estimate of the limit of the sequence. To see that, letk = 1 in equation (2.43) and lete 1 (A r ) be an approximation to the limit B 1n , the expression forR r is then derived as follows, R r =R 1 (r) =A r e 1 (A r ) =A r A r1 d r+1 A r d r d r+1 d r = d r d r+1 d r+1 d r : (2.50) The Levinv 1 transformation utilized in this study is presented next. 37 v 1 Transformation Levin’sv 1 transformation is associated withk = 1 andf 0 (n) = 1 in equation (2.46). Specifically, it can be written as, v 1 (A) n = A n A n+1 R n R n+1 1 1 R n R n+1 : (2.51) Comparison with equation (2.41) shows that the difference in these two transformations is in using R n instead of A n . It is observed that while three consecutive values are required for the Shanks e 1 transformation, four consecutive values are needed for the Levinv 1 transformation. These four values are associated with equations of the form, A n1 =B 1n +a 1n q n1 1n ; A n =B 1n +a 1n q n 1n ; A n+1 =B 1n +a 1n q n+1 1n ; A n+2 =B 1n +a 1n q n+2 1n : (2.52) From these equations,d n andR n are evaluated as, d n =A n A n1 =a 1n q n1 (q 1); (2.53) R n = a 2 1n q 2n1 (q 1) 2 a 1n q n1 (q 1) 2 =a 1n q n ; (2.54) resulting in a value ofv 1 (A) n equal to, v 1 (A) n = A n a 1n q n+1 A n+1 a 1n q n a 1n q n+1 a 1n q n = B 1n q n+1 B 1n q n q n+1 q n =B 1n : (2.55) 38 The above equation shows that applying the Levinv 1 transformation to a sequence of numbers exhibiting convergence in the form of equation (2.52) recovers the exact limit for that sequence. 2.2.3 NumericalExamples The convergence acceleration of polynomial chaos approximation utilizing nonlinear sequence transformations is explored numerically in this section. In all cases, random- ness is introduced in the differential equation through a coefficient modeled as a random process. Results are presented for both cases of an initial value problem and a stochastic elliptic problem. For the sake of illustration, the procedure of transformation is shown as what follows. For instance for the mean value, in the case ofe 1 transformation which iterates over three elements, the transformation procedure is as follows, (u (1) 0 ;u (2) 0 ;u (3) 0 ) e 1 7!u 0 (1) 0 (u (2) 0 ;u (3) 0 ;u (4) 0 ) e 1 7!u 0 (2) 0 (u (3) 0 ;u (4) 0 ;u (5) 0 ) e 1 7!u 0 (3) 0 ! (u 0 (1) 0 ;u 0 (2) 0 ;u 0 (3) 0 ) e 1 7!u 00 (1) 0 (u (4) 0 ;u (5) 0 ;u (6) 0 ) e 1 7!u 0 (4) 0 ! (u 0 (2) 0 ;u 0 (3) 0 ;u 0 (4) 0 ) e 1 7!u 00 (2) 0 . . . (2.56) The values forming the second sequence (u 0 (1) 0 ;u 0 (2) 0 ;u 0 (3) 0 ;:::) can be transformed to form the third sequence (u 00 (1) 0 ;u 00 (2) 0 ;u 00 (3) 0 ;:::) as shown above. This procedure can 39 be continued until there are three available values in the sequence. In the case of v 1 transformation, four numbers are needed as shown below, (u (1) 0 ;u (2) 0 ;u (3) 0 ;u (4) 0 ) v 1 7!u 0 (1) 0 (u (2) 0 ;u (3) 0 ;u (4) 0 ;u (5) 0 ) v 1 7!u 0 (2) 0 (u (3) 0 ;u (4) 0 ;u (5) 0 ;u (6) 0 ) v 1 7!u 0 (3) 0 (u (4) 0 ;u (5) 0 ;u (6) 0 ;u (7) 0 ) v 1 7!u 0 (4) 0 ! (u 0 (1) 0 ;u 0 (2) 0 ;u 0 (3) 0 ;u 0 (4) 0 ) v 1 7!u 00 (1) 0 (u (5) 0 ;u (6) 0 ;u (7) 0 ;u (8) 0 ) v 1 7!u 0 (5) 0 ! (u 0 (2) 0 ;u 0 (3) 0 ;u 0 (4) 0 ;u 0 (5) 0 ) v 1 7!u 00 (2) 0 . . . (2.57) This procedure is applied independently to each coefficient in the polynomial chaos approximation of the solution. As explained in the previous section, the effective appli- cation of sequence transformations for the convergence acceleration of polynomial chaos solutions requires knowledge of the form of the error. While convergence of the mean solution of equation (2.17) with respect to the polynomial order is well established as having the form in equation (2.34) [4, 5], similar convergence for higher order polyno- mial chaos coefficients is observed in the numerical examples investigated in this paper. This form of the error motivates the application of nonlinear maps such as Shanks and Levin transformations that have underlying model sequences with similar error form. These maps are further justified by the generality of the functional form described by equation (2.34). StochasticInitialValueProblem In this section, a stochastic ordinary differential equation is expressed in the form, d dt u(t;) =k(t;)u(t;); u(0;) = 1;8t2D = [0; 1]; a:s: in ; (2.58) 40 wherek(t;) is given by, k(t;) = 10 X i=1 (0:3)( 2i 10 ) cos( 2it 10 ) i : (2.59) The vectorf i g 10 i=1 is a set of 10 independent standard Gaussian variables, resulting in a stochastic dimension of s = 10. The analytical solution to the above differential equation is, u(t;) =e 0 B B @ 10 X i=1 (0:3) sin( 2it 10 ) i 1 C C A : (2.60) The mean, standard deviation and coefficient of variation (CV) of the solutionu(t;) are shown in Fig. 2.1 as function of time. In order to obtain a sequence of numbers for each 0 0.2 0.4 0.6 0.8 1 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 Time(s) Mean 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 Time(s) Standard deviation 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time(s) Coefficient of variation Figure 2.1: Stochastic Initial Value Problem: Mean, standard deviation and coefficient of vari- ation ofu(t;) fors = 10 polynomial chaos coefficient of u(t;), successive polynomial orders are considered in its expansion. In this example each polynomial chaos coefficient is evaluated as a 41 function of time. Convergence is investigated relative to the following relative error measure which is defined for each polynomial chaos coefficientu (t), " = v u u t ( 100 X i=1 (u (t i )u ;exact (t i )) 2 )=( 100 X i=1 u 2 ;exact (t i )) 8; (2.61) wheret i =idt;dt = 0:01s. The relative errors in the mean of the solution correspond- ing to different polynomial orders and transformed sequences are listed in Table 2.1 and plotted in Fig. 2.2. For Shanks and Levin transformations, each number is obtained by combination of three and four consecutive values, respectively. For instance, the first number in e 1 column is obtained by combination of coefficients corresponding to p = 1; 2; 3, the second number by combination of coefficients corresponding to p = 2; 3; 4 and the procedure continues for higher orders. Similarly, the first number in v 1 column is obtained by combination of coefficients corresponding top = 1; 2; 3; 4 and the second number by combination of coefficients corresponding top = 2; 3; 4; 5 and so on for other coefficients. It should be noted that the values in the columns labelede 2 1 , v 2 1 ande 3 1 are obtained by iteratinge 1 twice,v 1 twice ande 1 three times, respectively. Table 2.1: Stochastic Initial Value Problem: Relative error in the mean ofu(t;) p PCE e 1 v 1 e 2 1 v 2 1 e 3 1 1 4.96e-02 2 1.24e-02 3 2.39e-03 5.45e-04 4 3.68e-04 7.15e-05 3.00e-05 5 4.68e-05 7.78e-06 2.91e-06 1.42e-06 6 4.06e-06 7.29e-07 2.38e-07 1.03e-07 7 4.76e-07 6.02e-08 1.72e-08 7.14e-09 2.18e-09 1.81e-09 8 3.96e-08 4.45e-09 1.12e-09 4.61e-10 5.02e-11 3.07e-11 9 2.95e-09 2.99e-10 6.74e-11 2.77e-11 1.70e-12 1.87e-12 The transformation is applied independently to the sequence of polynomial chaos coeffi- cients associated with each term in the expansion. The convergence acceleration results 42 1 2 3 4 5 6 7 8 9 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 p Relative error PCE e 1 v 1 e 1 2 v 1 2 e 1 3 Figure 2.2: Stochastic Initial Value Problem: Relative error in the mean ofu(t;) of some polynomial chaos coefficients associated with higher order polynomials are plotted in Fig. 2.3 and Fig. 2.4. Note that the ending point in all plots is p = 9, but the starting point varies depending on the order of the polynomial corresponding to the coefficient of interest (see equation (2.15)). In the case of the mean value, the start- ing point is p = 1 instead of p = 0 since the lowest order expansion in this study is p = 1. In order to illustrate the effect of applying the transformation in reducing the computational effort, the number of operations in each case is listed in Table 2.2. The Runge-Kutta algorithm for solving the stochastic initial value problem described by equation (2.10) consists of a sequence of matrix-vector multiplication (MATVEC). The number of operations for each such MATVEC is of ordern 2 , wheren = N + 1. It should be noted that, due to the sparsity of the matrix in this problem, the number of operations for each MATVEC is of order significantly less thann 2 . For a sparse matrix, the number of operations will be function ofn, sayf(n), that depends on the number of non-zero elements in the matrix, which is problem dependent. For the present problem, 43 1 2 3 4 5 6 7 8 9 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 p Relative error PCE e 1 v 1 e 1 2 v 1 2 e 1 3 2 3 4 5 6 7 8 9 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 p Relative error PCE e 1 v 1 e 1 2 v 1 2 e 1 3 Figure 2.3: Stochastic Initial Value Problem: (top) Relative error in polynomial chaos coef- ficient associated with 10 , (bottom) Relative error in polynomial chaos coefficient associated with 4 7 . the number of non-zero elements in each row of the matrix is constant and does not grow with the polynomial order. This is due to the linear representation of the random 44 3 4 5 6 7 8 9 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 p Relative error PCE e 1 v 1 e 1 2 4 5 6 7 8 9 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 p Relative error PCE e 1 v 1 e 1 2 Figure 2.4: Stochastic Initial Value Problem: (top) Relative error in polynomial chaos coeffi- cient associated with ( 2 5 1) 7 , (bottom) Relative error in polynomial chaos coefficient asso- ciated with ( 2 6 1)( 2 8 1). parameters in this problem. In the stochastic Galerkin approach, the sparsity structure of the matrix is identical to the sparsity structure of the (N + 1N + 1) matrix with 45 (k;j) entry given by X i h i j k i. Therefore, it is useful to scrutinize the expression X i h i j k i for the present problem. Each polynomial j is the product of orthogo- nal one-dimensional polynomials in each of the ten variables l . Letting q denote the order of the one-dimensional i -polynomial within j , the product of this polynomial with i can be shown to be the sum of two Hermite polynomials, of order q 1 and q + 1, in accordance with the recursion formula of Hermite polynomials. It can then be concluded thath i j k i is zero unless the order of i in k is q 1 or q + 1 and the order of the remaining l ; (l6= i); is the same in k and j . It can thus be seen that, for givenj, the summation X i h i j k i contains at most a single non-vanishing term. Consequently, given a stochastic dimensions, there can be at most 2s non-zero entries in each row, irrespective of the polynomial order. Therefore, in this particular example it can be assumed that the computational cost of matrix-vector multiplication is of order f(n) =n. The number of operations can then be given byrn, wherer is the number of iterations required for convergence. The number of iterations can generally vary with the size of the problem (here the size of the problem increases asp increases) and the coef- ficient of variation. However, the number of iterations in the Runge-Kutta scheme only depends on the number of discrete points in the time domain. In the present problem, the number of points in the time domain is fixed for all polynomial ordersp, (r = 100). It should also be noted that, at each step (iteration) of Runge-Kutta, the matrix-vector multiplication is carried out four times (according to the algorithm) but this cost is not included in the estimate of the computational effort since it is incurred analogously, by all methods. The number of operations for each transformation is then the summation of the number of operations of each solution required for the transformation. For instance, the first value in the column labelede 1 is the summation of the number of operations for the expansions corresponding top = 1; 2; 3 resulting in 3.630e+04 operations. 46 Table 2.2: Stochastic Initial Value Problem: Computational effort in terms of estimated number of operations p No. of coeff. PCE e 1 v 1 e 2 1 v 2 1 e 3 1 1 11 1.100e+03 2 66 6.600e+03 3 286 2.860e+04 3.630e+04 4 1001 1.001e+05 1.353e+05 1.364e+05 5 3003 3.003e+05 4.290e+05 4.356e+05 4.367e+05 6 8008 8.008e+05 1.201e+06 1.230e+06 1.236e+06 7 19448 1.945e+06 3.046e+06 3.146e+06 3.175e+06 3.182e+06 3.182e+06 8 43758 4.376e+06 7.121e+06 7.422e+06 7.522e+06 7.557e+06 7.557e+06 9 92378 9.238e+06 1.556e+07 1.636e+07 1.666e+07 1.679e+07 1.679e+07 From Table 2.2, it can be seen that the computational effort in all transformations at orderp is remarkably less than the computational effort in PCE at orderp + 1. From Table 2.1 for the mean of u(t;), it can also be seen that applying the transformation at orderp in all cases excepte 1 transformation, provides polynomial chaos coefficients with higher accuracy than PCE at orderp+1. In other words, higher accuracy and more computational efficiency are gained in all transformations excepte 1 transformation. It should be noted that in the case ofe 1 transformation, while accuracy is improved over an order p approximation, it remains slightly below order p + 1 approximation. The comparison for the accuracy and computational cost is illustrated again in Fig. 2.5 where the relative error in the mean of the solution is plotted against the computational effort. In order to show the comparison of the computational cost and accuracy for higher order polynomial chaos coefficients, results corresponding to the coefficient associated with ( 2 6 1)( 2 8 1) are listed in Table 2.3 and Table 2.4 for e 2 1 and e 1 transformation respectively. It is worth noting that, in order to illustrate higher accuracy and increased computational efficiency, results in Table 2.3 are presented and compared between PCE atp = 9 ande 2 1 atp = 8. 47 10 3 10 4 10 5 10 6 10 7 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 Computational effort (Op. count) Relative error PCE e 1 v 1 e 1 2 v 1 2 e 1 3 Figure 2.5: Stochastic Initial Value Problem: Relative error in the mean ofu(t;) versus com- putational effort defined as operation count Table 2.3: Stochastic Initial Value Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with ( 2 6 1)( 2 8 1) applying e 2 1 transformation Computational effort Relative error PCE (p = 9) e 2 1 (p = 8) PCE (p = 9) e 2 1 (p = 8) 9.238e+06 7.522e+06 1.15e-04 1.44e-05 Table 2.4: Stochastic Initial Value Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with ( 2 6 1)( 2 8 1) applying e 1 transformation Computational effort Relative error PCE e 1 PCE e 1 8.008e+05 1.201e+06 3.37e-02 1.15e-02 1.945e+06 3.046e+06 6.28e-03 1.45e-03 4.376e+06 7.121e+06 9.34e-04 1.73e-04 9.238e+06 1.556e+07 1.15e-04 1.81e-05 It should be noted that results in this paper are presented for the stochastic dimen- sion,s = 10. However, as noted earlier, the computational effort increases significantly and nonlinearly with the polynomial degree and an assessment of gains in computational 48 efficiency must be described as function of the dimension and order of the PCE. While for problems in lower dimensions the increase in the computational effort may not be significant, for each stochastic dimensions, a range of polynomial degrees can be spec- ified where computational efficiency is improved by applying the proposed transforma- tions. Specifically, the total computational cost of implementing these transformations in these ranges is less than the computational cost of PCE with one order higher. These ranges of polynomial degrees are listed in Table 2.5 for each stochastic dimension and transformation. It should be noted that, ranges in Table 2.5 correspond to the case where the operation count of each MATVEC is assumed to be of ordern, resulting in a conser- vative estimate. As expected, it is observed that a wider range is covered as the stochastic dimension increases. This further justifies the application of these transformations on problems with significant growth in computational cost. Table 2.5: Stochastic Initial Value Problem: Range of polynomial degrees in which computational efficiency is gained s e 1 v 1 e 2 1 v 2 1 e 3 1 3 p = 3 4 3p 4 p = 4 5 3p 5 4p 5 p = 5 6 3p 6 4p 6 5p 6 7 3p 7 4p 7 5p 7 p = 7 p = 7 8 3p 9 4p 8 5p 8 7p 8 7p 8 9 3p 10 4p 9 5p 9 7p 9 7p 9 10 3p 11 4p 10 5p 10 7p 10 7p 10 StochasticEllipticProblem The stochastic elliptic boundary value problem in this section is specified as, d dx (a(x;) d dx u(x;)) =x; x2D = [0; 1]; a:s: in ; u(x;) = 0; x2@D; (2.62) 49 wherea(x;) is given by, a(x;) = 1 + 0:073 10 X i=1 i i (x) i ; (2.63) and, i = exp(k 2 =8) = 10 4 ; i = 2k 1 ori = 2k (2.64) i (x) = 8 > > < > > : sin((k 0:5)x) i = 2k 1; cos((k 0:5)x) i = 2k: (2.65) The vectorf i g 10 i=1 is a set of independent and identically distributed random variables with uniform distribution on the interval [ p 3; p 3]. Consequently, the random vari- ables have zero mean and unit variance. This construction of a(x;) results in a sta- tionary random process which satisfies condition (2.18) [5]. The mean, standard devi- ation, and coefficient of variation of a(x;) are 1; 0:163; 0:163, respectively. Similar to the previous section, successive polynomial orders are considered in the expansion to obtain a sequence of numbers. The exact analytical solution is not available in this case. To investigate convergence of the acceleration schemes, results corresponding to coefficients with the expansion order up top = 6 are compared to coefficients with a higher expansion order that provide enough accuracy. The highest expansion order con- sidered isp = 8. The relative error measure is then defined for each polynomial chaos coefficient as follows, " = v u u t ( 19 X i=1 (u (x i )u ;8 (x i )) 2 )=( 19 X i=1 u 2 ;8 (x i )) 8; (2.66) where is the index of the deterministic polynomial chaos coefficient andx i =idx;dx = 0:05. The relative errors in the mean of the solution of equation (2.62) are listed in 50 Table 2.6 and plotted in Fig. 2.6. Convergence acceleration results of polynomial chaos coefficients associated with higher order polynomials are plotted in Fig. 2.7 and Fig. 2.8. Similarly to the previous case, the starting point is the polynomial order corresponding to the coefficient of interest (except the case of the mean which is started fromp = 1) and remarkable improvement is observed. The computational effort in terms of the number of operations is listed in Table 2.7. According to equation (2.24), the number of coefficients in the expansion is the product of the number of terms in the polyno- mial expansionN + 1, and the number of degrees of freedom in the spatial domainM, (M = 19). The number of operations in this case is also dependent on the total num- ber of coefficients (N + 1)M and the number of iterations. Similarly to the previous problem, the random parameter is a linear form, as defined in section 3.1. The recur- sion relation for Legendre polynomials is analogous to that for Hermite polynomials and consequently, the same conclusion can be drawn that, the number of non-zero entries in matrixA in equation (2.29) remains constant asp increases. It should be noted that, in this problem the spatial domain of the problem is discretized and the number of non-zero elements is indeed dependent on the spatial discretization. However, this discretization is independent of the polynomial order. Therefore, in this particular example, it can also be assumed that the computational cost of matrix-vector multiplication is of order n, where n=(N+1)M. Consequently, the number of operations required for convergence is of orderrn. In this case, the iterative linear solver used is based on the conjugate gradi- ent algorithm. Unlike the previous example, the number of iterations in this algorithm is not fixed and is dependent on the size of the problem. It can be seen from Table 2.7 that, the number of iterations to converge to a solution increases with the size of the problem. Similarly to the previous results for the stochastic initial value problem, the computa- tional cost is here again decreased while the accuracy is increased for all transformations 51 Table 2.6: Stochastic Elliptic Problem: Relative error in the mean ofu(x;) p PCE e 1 v 1 e 2 1 1 9.07e-05 2 1.57e-06 3 3.47e-08 7.54e-09 4 9.07e-10 1.47e-10 1.91e-11 5 2.67e-11 3.86e-12 7.19e-13 5.90e-13 6 8.69e-13 1.60e-13 2.94e-14 3.09e-14 1 2 3 4 5 6 10 −14 10 −12 10 −10 10 −8 10 −6 10 −4 p Relative error PCE e 1 v 1 e 1 2 Figure 2.6: Stochastic Elliptic Problem: Relative error in the mean ofu(x;) Table 2.7: Stochastic Elliptic Problem: Computational effort in terms of number of operations p No. of coefficients r PCE e 1 v 1 e 2 1 1 209 39 8.151e+03 2 1254 59 7.399e+04 3 5434 115 6.249e+05 7.070e+05 4 19019 205 3.899e+06 4.598e+06 4.606e+06 5 57057 290 1.655e+07 2.107e+07 2.114e+07 2.115e+07 6 152152 501 7.623e+07 9.667e+07 9.730e+07 9.737e+07 excepte 1 . It is again noted that thee 1 transformation is still effective in accelerating con- vergence if three consecutive coefficients are combined and the result is compared with 52 1 2 3 4 5 6 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 p Relative error PCE e 1 v 1 e 1 2 2 3 4 5 6 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 p Relative error PCE e 1 v 1 e 1 2 Figure 2.7: Stochastic Elliptic Problem: (top) Relative error in polynomial chaos coefficient associated with 1 , (bottom) Relative error in polynomial chaos coefficient associated with 1 2 (3 2 10 1). 53 3 4 5 6 10 −8 10 −6 10 −4 10 −2 10 0 p Relative error PCE e 1 v 1 4 5 6 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 p Relative error PCE e 1 (a) (b) Figure 2.8: Stochastic Elliptic Problem: (top) Relative error in polynomial chaos coefficient associated with 1 2 (3 2 5 1) 2 , (bottom) Relative error in polynomial chaos coefficient associated with 1 2 (5 3 4 3 4 ) 9 . 54 the coefficient corresponding to the highest order among those three orders. The relative error in the mean of the solution is plotted against the computational effort in Fig. 2.9. The computational effort and accuracy results associated with applyingv 1 ande 1 to the term 1 2 (3 2 10 1) are shown in Tables 2.8 and 2.9, respectively. Similarly to the pre- vious example, results in Table 2.8 are presented and compared for different orders, PCE at p = 6 and v 1 transformation at p = 5. It should be noted that, the growth in computational effort in the stochastic elliptic problem is more significant than stochastic initial value problem since the number of iterations to converge to a solution is not fixed in the former case. The same analysis can be done for the stochastic elliptic problem to specify the range of polynomial orders in which the computational efficiency is gained. This analysis requires the knowledge of the number of iterations overp for each stochas- tic dimension. However, as a conservative assumption the number of iterations can be considered to be fixed over p for each stochastic dimension. Consequently, the same results will be obtained as those listed in Table 2.5. 10 4 10 5 10 6 10 7 10 8 10 −14 10 −12 10 −10 10 −8 10 −6 10 −4 Computational effort (Op. count) Relative error PCE e 1 v 1 e 1 2 Figure 2.9: Stochastic Elliptic Problem: Relative error in the mean ofu(x;) versus computa- tional effort defined as operation count 55 Table 2.8: Stochastic Elliptic Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with 1 2 (3 2 10 1) applyingv 1 transformation Computational effort Relative error PCE (p = 6) v 1 (p = 5) PCE (p = 6) v 1 (p = 5) 7.623e+07 2.114e+07 9.56e-09 2.23e-09 Table 2.9: Stochastic Elliptic Problem: Comparison of the computational effort and accuracy for the polynomial chaos coefficient associated with 1 2 (3 2 10 1) applyinge 1 transformation Computational effort Relative error PCE e 1 PCE e 1 3.899e+06 4.598e+06 9.74e-06 1.67e-06 1.655e+07 2.107e+07 2.88e-07 4.13e-08 7.623e+07 9.667e+07 9.56e-09 1.05e-09 2.3 Sensitivityanalysis The sensitivity analysis of a model, is an attempt to quantify the importance of each input parameter. One of the approaches for sensitivity analysis is a variance based approach that has been employed in this paper. In this method, the variance of the output is decomposed into a sum of the variance of each input random variable. For this pur- pose, Sobol introduced sensitivity indices that are usually evaluated using Monte Carlo simulations [59]. The estimation of these sensitivity indices are greatly facilitated using Polynomial Chaos Expansion. For a model represented byy = f(x) = f(x 1 ;:::;x n ), wherey is the output andx =f0x i 1;i = 1;:::;ng is the set of input parameters, the Sobol decomposition of the output function into summands of increasing dimension is written as follows [64], f(x) =f 0 + n X i=1 f i (x i ) + n X 1i<jn f ij (x i ;x j ) +::: +f 1;2;:::;n (x 1 ;:::;x n ); (2.67) 56 wheref 0 is a constant and the summands are orthogonal as follows, R f i 1 ;:::;is (x i 1 ;:::;x is )f j 1 ;:::;jt (x j 1 ;:::;x jt )dx i 1 :::dx is :dx j 1 :::dx jt = 0; fi 1 ;:::;i s g6=fj 1 ;:::;j t g: (2.68) Assuming the functionf is square integrable, each summand in equation (2.67) is square integrable consequently. Therefore, the following relation can be written due to the orthogonality given in equation (2.80). Z f 2 dxf 2 0 = n X s=1 X 1i 1 <:::<isn Z f 2 i 1 :::is dx i 1 :::dx is : (2.69) From the equation above, the total varianceD and partial variancesD i 1 :::is are defined as follows, D = Z f 2 dxf 2 0 ; D i 1 :::is = Z f 2 i 1 :::is dx i 1 :::dx is : (2.70) Accordingly, the global sensitivity indices are defined as follows, S i 1 :::is = D i 1 :::is D : (2.71) It should be noted thatS i 1 :::is are non negative and their sum is one as follows, n X s=1 X 1i 1 <:::<isn S i 1 :::is = 1: (2.72) As explained in the previous section, the polynomial chaos expansion of a function can be written as follows, U() = 1 X j=0 ^ u j j (): (2.73) 57 In practice, the above expansion is truncated and is represented with a finite number of term as follows, U() = M X j=0 ^ u j j (): (2.74) The number of terms in this expansionM + 1 depends on the number of random vari- ables, namely the stochastic dimension of the problemn and the degree of the expansion p as follows, M + 1 = (n +p)! n!p! : (2.75) Similarly to the equation (2.67), the polynomial chaos decomposition of a function can be written as follows, U() = ^ u 0 + n X i=1 X 2 i ^ u ( i ) + X 1i 1 <i 2 n X 2 i 1 ;i 2 ^ u ( i 1 ; i 2 ) +::: + X 1i 1 <:::<isn X 2 i 1 ;:::;is ^ u ( i 1 ;:::; is ) +::: + X 2 i 1 ;:::;is ^ u ( 1 ;:::; n ): (2.76) In the equation above, =f i ;i = 1;:::;ng is a set of integer numbers corresponding to degree of each one dimensional polynomial in the multivariate polynomial as follows, ( 1 ;:::; n ) = n Y i=1 K i ( i ); (2.77) whereK i is a one dimensional polynomial. It should be noted that, the degree of the expansionp in equation (2.75) can be given as follows, p = n X j=1 j : (2.78) 58 Accordingly, the set i 1 ;:::;is is defined as a set of’s such that, only the indices (i 1 ;:::;i s ) are nonzero. This definition is written as follows, i 1 ;:::;is = 8 > > > > > < > > > > > : j > 0 8j = 1;:::;n; j2 (i 1 ;:::;i s ) : j = 0 8j = 1;:::;n; j = 2 (i 1 ;:::;i s ) (2.79) According to equation(2.76), summands in the Sobol decomposition in eauation(2.67) can be written as follows, f i 1 ;:::;is ( i 1 ;:::; is ) = X 2 i 1 ;:::;is ^ u ( i 1 ;:::; is ): (2.80) Consequently the polynomial chaos based Sobol indices are denoted by, S i 1 :::is = X 2 i 1 ;:::;is ^ u 2 E[ 2 ]=D PCE ; (2.81) whereD PCE is the total variance and according to equation(2.74) is given by, D PCE = M X j=1 ^ u 2 j E[ 2 ()]: (2.82) Note that the summation above starts from (j = 1). The equation (2.81) shows the sensitivity of the output with respect to the joint uncertainty sources i 1 ;:::; is . Simi- larly, first order sensitivity indicesS i 1 can be calculated by only taking the single index i 1 = 1;:::;n. 59 Chapter3 WindTurbineModel In this section, we explain the details of the wind turbine model used in the optimization process. The wind turbine model consists of three main components, the wind model, the aeroelastic model and the structural model that are explained separately in what follows. 3.1 Windmodel The aerodynamic loading is the force of the wind flow passes through the structure. The wind field that is seen by the blade varies in space and time. Therefore, the wind simulation has become an important part in the wind turbine structural analysis. Due to the atmospheric turbulence for a mid to large size wind turbine spatial variation should be considered. Let V (x;t) 2 R 3 denote the multivariate wind speed, as a function of spatial position and timet2 [0;1]. ConsideringV (x;t) as a stochastic process, for each position in spacex2R 3 the power spectral density can be defined. The time series of the wind should then satisfy the prescribed power spectral density. For the sake of simplification in presenting the procedure, the random process is only considered with respect to time. Letv(t) denote a random process at a spatial pointx 0 . The covariance function of this random process is given by, C(t 1 ;t 2 ) =E[v(t 1 )v(t 2 )]: (3.1) 60 For a staionary random process the covariance function can be represented as a function of the time difference, =t 1 t 2 . The above equation can then be written as follows, C() =C(t 1 ;t 2 ): (3.2) The realizations of vectorv can be characterized by a linear combination of the eigen- vectors of the covariance function as follows, v(t) = 1 X i=1 x i i (t); (3.3) where i ’s are the eigenvectors of the covariance function and can be given as follows, C i = i i ; (3.4) andx i ’s satisfy the following relation, E[x i x j ] = i ij ; (3.5) where i ’s are the eigenvalues of the covariance function. Eachx i can then be written in terms of zero mean, unit variance uncorrelated Gaussian random variables i as follows, x i = p i i (3.6) Accordingly the random processv(t) can be expressed as, v(t) = 1 X i=1 p i i (t) i : (3.7) 61 For a periodic covariance function, the equation (3.4) can be written in the following form, Z T=2 T=2 C(t 1 ;t 2 ) i (t 2 )dt 2 = i i (t 1 ): (3.8) Therefore, the covariance function can be written as an expansion in terms of its eigen- vectors and eigenvalues as follows, C(t 1 ;t 2 ) = 1 X i=1 i (t 1 ) i (t 2 ) i (3.9) The covariance function can also be represented by the fourier transformation of the power spectral density as follows, C(t 1 ;t 2 ) = 1 X j=1 e i! j (t 1 t 2 ) s j ; (3.10) wheres j is given by, s j = 1 T Z T=2 T=2 e i! j C()d (3.11) Substituting equation (3.10) in equation (3.8), it can be concluded that the eigenvalues and eigenvectors have the following form, i =Ts i ; i (t) = 1 p T e i! k t : (3.12) Consequently, the random processv(t) in equation (3.3) can be written as, v(t) = Z 1 1 e i!t dz(!); (3.13) 62 wheredz(!) satisfies the following relation, E[dz(!)dz (! 0 )] = 8 > > < > > : S(!)d! ! =! 0 0 !6=! 0 (3.14) Similarly to the equation (3.9), the covariance function can be obtained from equations (3.13), (3.14). According to Wiener-Khintchine equation, the covariance function can be written as, C(t 1 ;t 2 ) = Z 1 1 e i!(t 1 t 2 ) S(!)d! (3.15) Comparing two equations (3.10) and (3.15), the following relation can be written, S(!)d! = 1 X k=1 s k (!! k ) (3.16) In practice, the random process can be expressed in terms of finite number of terms in equation (3.7). The random process can now be characterized in terms of the power spectral density as follows [15], v(t)' n d X j=n d e i! j t q S(! j )! j : (3.17) While this representation is similar to the Karhunen-Loeve expansion, there is a subtle difference between the number of terms that should be considered in the expansion. In KL expansion the process can be characterized with a relatively small number of terms since the few first eigenvalues are the largest eigenvalues and have the most contribution to the characterization of the process. But, the eigenvalues in (3.17) are not in descend- ing order and each term can have significant contribution to the characterization of the process. Therefore, in this case a relatively large number of terms should be considered 63 to generate the random process. It should be noted that the random process in (3.17) is only obtained for a point in space. In reality, the magnitude of the wind varies in space as well. Therefore, in the case of multivariate random processes, the spectral density matrix that includes cross spectral densities, is considered. The cross spectral densi- ties between two points are related to the spectral densities of those points through the coherence function as follows, S ij (!) = q S i (!)S j (!)Coh ij (!); (3.18) whereCoh i j(!) is the coherence function andS i (!) is the power spectral density at the spatial pointi. The spectral density matrix denoted by S(!) can also be decomposed as follows, S(!) = 1 X k=1 k k (!) k (!); (3.19) where k (!) and k (!) satisfy the eigenvalue problem, S(!) k (!) = k (!) k (!): (3.20) In this case, the multivariate random processv(t) can be written in terms of the decom- position of the spectral density matrix as follows, v(t)' n d X j=n d ns X k=1 e i! j t q k (! j )! k (! j ) jk ; (3.21) where jk are zero mean, unit variance uncorrelated Gaussian random variables. It should be noted that the random process above is a multivariate process since k (! j ) 64 is a vector with sizeM as the number of spatial points. The above equation can then be written in the following form, V (x i ;t)' n d X j=n d ns X k=1 e i! j t q k (! j )! k (x i ;! j ) jk i = 1;:::;M: (3.22) As mentioned earlier, in order to generate a stationary random process with the pre- scribed power spectral density, a large number of independent random variables should be considered. In the sequel, the generated random process is used in the optimization problem and it is desirable to reduce the number of random dimension of the problem. For this purpose, a Karhunen-Loeve expansion (KLE) can be applied to the generated random process and only first few significant modes can be considered in the charac- terization of the random process. Now, for the random process denoted byV (x;t) the physical domain can be considered to be a pair of space and time and can be defined as, z = 2 4 x t 3 5 : (3.23) Now, the windV (z) can be denoted as a random processV (z;) :D !R, where is a set of elementary events on which a probability space is constructed and2 . The covariance matrix can be formed from the generated random process as follows, R(z 1 ;z 2 ) =E[V (z 1 ;)V (z 2 ;)] =<V (z 1 ;);V (z 2 ;)>; (3.24) where<:;:> denoted the inner product of two vectors. The eigenvectors of the covari- ance matrix are a set of orthonormal basis in the Karhunen-Loeve decomposition and is the solution of the following eigen problem, R i = i i : (3.25) 65 Each eigenvector can also be considered as an eigenmode in the physical domain of the process. The projection of the random process on the eigenvectors is a random variable that is associated with the corresponding eigenmode. For example, the corresponding random variable to the eigenmodei can be given as follows, i () =<V (z;); i >= p i ; (3.26) where i is the eigenvalue associated with the modei and p i is a normalizing factor. Therefore, the random variables i are unit variance random variables. They are zero mean since the process is zero mean and also they are mutually orthonormal as follows, < i ; j >= ij (3.27) Finally the Karhunen-Loeve decomposition of a zero mean process can be written as follows, V (z;) = 1 X i=1 p i i (z) i (): (3.28) In practice, only first few largest eigenvalues are considered in the expansion, therefore the expansion is truncated and can be denoted as follows, V (z;)' n V X i=1 p i i (z) i (): (3.29) This reduced order model for the wind can then be used in the aeroelastic analysis of the blade. The details of the aeroelastic model are explained in the next section. 66 3.2 Aeroelasticmodel A wind turbine extracts the mechanical energy from the wind. Apart from the wind sim- ulation, aeroelastic analysis is needed to determine the wind loads on the blades. There are various levels of complexity in the models that describe the interaction of the fluid with the structure. These models range from the popular Blade Element Momentum method (BEM) to the models incorporating the Navier-Stokes equations. In this study, the BEM as the most common tool for evaluating the aerodynamic forces on the wind turbine rotor is used. In BEM, the rotor blade is divided to several annular elements as shown in Fig. 3.1. There is an assumption that all sections along the rotor are indepen- Figure 3.1: Illustration of annular elements in BEM method 67 dent and is analyzed separately. In addition, the method assumes that the forces from the blades on the flow is constant in each annular elements. This assumption is valid for a rotor with indefinite number of blades. In order to use BEM for a rotor with finite num- ber of blades a correction factor known as Prandtl’s tip loss is introduced in the sequel. At a given radial section, there is a difference between the wind speed from far upstream to deep in the wake that is considered as the induced velocity. The main objective in the BEM method is to evaluate the induced velocities knowing the input parmaeters such as incoming velocity, rotational velocity of the rotor blade, chord length and twist of the blade at each section, and airfoil data (lift and drag coefficients) at each section. The induced velocities in both directions are expressed with two parameters, namely aerodynamic induction factors denoted bya;a 0 . The input and output of BEM method is shown in Fig. 3.2. The computation of the aerodynamic forces and extracted power is then based on the aerodynamic induction factors. The details of the BEM model is described in the sequel. The blade of the wind turbine is a long and slender structure. Figure 3.2: Input and output of BEM method The flow at each section can be considered two dimensional and therefore, 2D airfoil data can be applied for the aeroelastic analysis. The reacting force from the flow is decomposed into a component parallel and perpendicular to the incoming velocity. The parallel component is drag forceD, and the perpendicular component is the lift forceL, as shown in Fig. 3.3. 68 Figure 3.3: Definition of lift and drag If the airfoil is designed for the airplane, the main objective is the maximization of theL=D ratio. The lift is the force to overcome gravity. The lift and drag coefficients are, C l = L 1=2V 0 c ; C d = D 1=2V 0 c : (3.30) wherec is the chord length and is the air density. For a slow moving aircraft such as a wind turbine, the lift and drag coefficients are functions of and Reynolds number. The angle is the angle of attack defined as the angle between the chordline and V 0 . The upstream flow with velocity V 0 is decreased to u on the rotor plane and to u 1 in downstream. With the ideal rotor assumption the relationship betweenV 0 ,u andu 1 can be derived as explained in what follows. Due to the change in the velocity, there should be a change in the pressure resulting in the thrust force on the rotor plane. The shaft power P is the absorbed energy by the turbine that can be obtained by the Bernoulli equation. The flow is incompressible and no external force acts on the fluid in the up 69 Figure 3.4: Velocity and pressure up- and downstream of the rotor or downstream. Therefore, the Bernoulli equation can be written for the flow from upstream to the rotor and from the rotor to downstream as written in Eq. (3.31). Figure 3.5: Control volume around a wind turbine p 0 + 1 2 V 2 0 =p + 1 2 u 2 ; p p + 1 2 u 2 =p 0 + 1 2 u 2 1 : (3.31) 70 From above equations, p = 1 2 (V 2 0 u 2 1 ): (3.32) Considering the axial momentum and an ideal rotor for the control volume in Fig. 3.5, a relation can be written as follows, u 2 1 A 1 +V 2 0 (A cv A 1 ) + _ m side V 0 V 2 0 A cv =T: (3.33) The term _ m side is found from the conservation of mass, A 1 u 1 +(A cv A 1 )V 0 + _ m side =A cv V 0 ; _ m side =A 1 (V 0 u 1 ): (3.34) The relationship betweenu andu 1 is also found from conservation of mass, uA =u 1 A 1 : (3.35) Combining the equations of conversation of momentum and mass, the equation below is found. T =uA(V 0 u 1 ): (3.36) The thrust and the pressure drop are related and denoted by, T = PA; (3.37) 71 whereA = R 2 is the area of the rotor. Combining the pressure drop equation (3.32), the thrust equation (3.36) and the relationship between the pressure drop and thrust, equation (3.37) it can be found that, u = 1 2 (V 0 +u 1 ): (3.38) The velocity in the rotor plane is the mean of the wind speedV 0 and the final velocity in the wakeu 1 . The velocity at the turbine is related to the velocity in the far upstream through the axial induction factora as follows, u = (1a)V 0 : (3.39) Therefore, from Eq. (3.38), the velocity in the wake is found, u 1 = (1 2a)V 0 ; (3.40) which can be introduced in Eq. (3.36) to obtain the thrustT . T = 2V 2 0 a(1a)A: (3.41) As mentioned earlier, similarly to the axial induction factor, there is an induction factor for the rotational velocity a 0 . The air at the wind turbine will rotate in the opposite direction of the blade. Therefore, the relative rotational velocity felt by the blade is (1 +a 0 )r! rot , where! rot is the angular velocity of the rotor. Moreover, considering the annular element and angular momentum equation, the torquedM is found by, dM =rC d _ m = 2r 2 uC dr; (3.42) 72 whereC is the rotational velocity in the wake. The absolute rotational velocity in the upstream is zero and is considered a 0 r! rot at the rotor. The rotational velocity in the wake is found 2a 0 r! rot according to the relationship for the axial velocities. The thrust and moment for the ring on the control volume with the area 2rdr can also be written as, dT = 2rdru(V 0 u 1 ) = 4rV 2 0 a(1a)dr; dM = 4r 3 V 0 ! rot (1a)a 0 dr: (3.43) As stated earlier, the relative velocity felt by a section of the blade is the combination of the axial velocityV 0 (1a) and the tangential velocityr! rot (1+a 0 ) as shown in Fig. 3.6. where is the local twist of the blade, is the angle between the plane of rotation and Figure 3.6: Velocities at rotor plane the relative velocityV rel . The local angle of attack is given by =, moreover it can be found that, tan = (1a)V 0 (1 +a 0 )r! rot : (3.44) 73 If the lift and drag coefficients are known the liftL and dragD force per length can be found from the equations below(Eq. (3.30)): L = 1 2 V 2 rel cC l ; D = 1 2 V 2 rel cC d : (3.45) The lift and drag forces are projected into the directions normal and tangential to the plane as shown in Fig. 3.7. The relationship for normal and tangential forces is as Figure 3.7: Normal and tangential forces follows, p N =L cos +D sin; p T =L sinD cos: (3.46) The above equations are normalized with respect to 1 2 V 2 rel c, C n = p N 1 2 V 2 rel c ; C t = p T 1 2 V 2 rel c : (3.47) 74 Further, a parameter in each section is defined as follows, (r) = c(r)B 2r ; (3.48) whereB denoted the number of blades andc(r) is the local chord length. The thrust and moment on each section can be written as, dT =Bp N dr; dM =rBp T dr: (3.49) From Fig. 3.6, it can also be written that, V rel sin =V 0 (1a); (3.50) V rel cos =r! rot (1 +a 0 ): (3.51) Combining the above equations for theV rel with the equations fordT anddM, it can be written that, dT = 1 2 B V 2 0 (1a) 2 sin 2 cC n dr; dM = 1 2 B V 0 (1a)r! rot (1 +a 0 ) sin cos cC t rdr: (3.52) 75 If the above equations are equalized with the equations (3.43) and the definition of(r) is used, then the expressions for the axial and tangential induction factors are obtained as follows, a = 1 4 sin 2 C n + 1 ; a 0 = 1 4 sin cos C t 1 : (3.53) All necessary equations for the BEM model are given. Each strip on the control volume is computed separately. The main objective is to find the induction factors and finding the velocity that is felt by the wind turbine. The thrust, power and loads in the blade can be computed accordingly. In order to obtain the induction factors, a procedure should be followed iteratively. The induction factors first initialized,a = a 0 = 0. The flow angle is then computed from Eq. (3.44). The local angle of attack is computed by having the twist at the corresponding section. The values of lift and drag coefficients at the obtained angle of attack are read from the table for the specific airfoil. The normalized normal and tangential forces are obtained from Eq. (3.47). Accordingly, the values for the induction factors can be obtained from Eq. (3.53). If the obtained values for a and a 0 have changed more than the tolerance, this procedure is continued until the convergence within the tolerance is obtained. There are empirical relations that are used to correct some assumptions in the aeroelastic analysis. Prandtl’s tip loss factor corrects the assumption of an infinite number of blades. In the real situation the number of blades are finite and therefore the vortex system in the wake is different from the system 76 with infinite number of blades. A correction factorF is introduced by Prandtl into the equations fordT anddM in equation (3.43), dT = 4rV 2 0 a(1a)Fdr; dM = 4r 3 V 0 ! rot (1a)a 0 Fdr; (3.54) In the equations above,F is given by, F = 2 cos 1 (e f ); (3.55) where f = B 2 Rr r sin : (3.56) If the prandtl’s tip loss is considered then the expressions for the induction factors are changed as follows, a = 1 4F sin 2 C n + 1 ; a 0 = 1 4F sin cos C t 1 : (3.57) Another correction is the Glauert correction when the axial induction factor is larger than a specific value. The relationship between the thrust coefficient anda can be found as following, C T = T 1 2 V 2 0 A = 4a(1a): (3.58) 77 Different empirical relations are found between C T and a to fit the measurement, for example, C T = 8 > > < > > : 4a(1a)F aa c 4(a 2 c + (1 2a c )a)F a>a c ; (3.59) where a c is 0:2 in this model [63] and F is the Prandtl’s tip loss factor. Considering this expression for the thrust coefficient and accordinglydT , the axial induction factor is corrected as follows, ifa<a c , a = 1 4F sin 2 C n + 1 ; (3.60) and ifa>a c , a = 1 2 h 2 +K(1 2a c ) p (K(1 2a c ) + 2) 2 + 4(Ka 2 c 1) i ; (3.61) where, K = 4F sin 2 C n : (3.62) As stated earlier, the aerodynamic induction factors are obtained from the aeroelastic model to compute the relative velocity that is seen by the blade. The normal load p N and tangential loadp T are given in equation (3.46). In order to compute the lift and drag forces, the relative velocity is needed and is given by, V rel = p [V 0 (1a)] 2 + [r! rot (1 +a 0 )] 2 : (3.63) It should be noted that the normal load and tangential load are random due to the ran- domness in the wind. In other words, for each realization of the wind, there is a realiza- tion of the aerodynamic force. Therefore similarly to the wind, the aerodynamic force 78 can be considered as a random process that varies in space and time. In the next section, the details of the reduced order model for the blade structure and construction of the generalized forces used in that model, are explained. 3.3 Structuralmodel In order to determine the temporal variation of the stresses in the blade, a structural model is required. The dynamic response is obtained from the model which is subjected to time dependent loads. The time dependent loads are then obtained through using an aerodynamic model that was described in the last section. In this section the details of the structural model is discussed. In order to get the stress state of the blade, a detailed finite element model of the blade is used. However, the dynamic response of the blade is obtained through a reduced order model based on modal analysis. The reduced order model is constructed from the detailed finite element model. The uncertainty is also incorporated in the reduced order model and consequently, the dynamic response of the blade is a random process . In addition to the randomness in the wind, one source of the randomness can be the uncertainty in the material properties of the blade that can be observed through the elastic modulus of the blade. In this work, the elastic modulus of the blade is considered to be random throughout the blade. The randomness in the elastic modulus will then be propagated in the generalized mass, stiffness and damping of the structure. The most common approach for such a problem incorporating uncertainty is to use Monte Carlo simulation. In Monte Carlo simulation, one can only obtain the realizations of the quantity of interest. However, in the current study each segment of the analysis is needed for the optimization problem. Therefore, it is desirable to have a functional form for the quantities of interest that can be used later in the optimization process. 79 In order to consider the randomness in the elastic modulus, the blade is divided to three different sections and the elastic modulus is considered to vary in each section through an independent uniform random variable. LetR denote the length of the blade, then the elastic modulus throughout the blade can be written as, E(r; m ) = 8 > > > > > < > > > > > : E 1 m 1 0<r< R 3 ; E 2 m 2 R 3 <r< 2R 3 ; E 3 m 3 2R 3 <r<R; (3.64) wherer is the radius in the spanwise direction,E 1 ;E 2 ;E 3 are mean values of the elastic modulus in three sections and m = ( m 1 ; m 2 ; m 3 ) is a vector of three independent and identically distributed uniform random variables varying around 1. The stiffness matrix of the blade is then a random matrix according to the randomness in the elastic modulus. In modal Analysis, the dynamic response of the blade is characterized in terms of eigenmodes. The eigenmodes are obtained from the eigenvalue problem as follows, L (i) = (! (i) ) 2 (i) ; (3.65) where L = M 1 K and M;K are mass and stiffness matrices that obtained from the finite element analysis. In the equation above, (i) ;! (i) are the eigenvector and the frequency associated with the modei respectively. As mentioned earlier, the matrixK is a random matrix due to the randomness in the elastic modulus of the blade. Therefore, the frequencies and eigenvectors are also random variables and random vectors. The 80 polynomial chaos expansion of the frequency and eigenvector of the mode i can be written as follows, ! (i) ( m )' n! X j=0 ! (i) j j ( m ); (i) ( m )' n X j=0 (i) j j ( m ); (3.66) where n ! and n are the finite numbers of terms considered in the expansion. The coefficients ! (i) j ; (i) j are deterministic coefficients and are obtained by projecting the random frequency and random vectors on the polynomials j ( m ). It should be noted that the polynomials are orthogonal with respect to the probability density function of the random variable in the problem. In this case, Legendre polynomials are used since the randomness in the eigenvalue problem is introduced through uniform random vari- ables. The projection of random frequency and random vectors on the polynomials is an integration that can be done based on the quadrature rule as follows, ! (i) j ' nq X q=1 ! (i) ( (q) m ) j ( (q) m )w (q) ; (i) j ' nq X q=1 (i) ( (q) m ) j ( (q) m )w (q) ; (3.67) where n q is the number of quadrature points and j ( (q) m );! (i) ( (q) m ); (i) ( (q) m ) are the values of the polynomials, frequency and eigenvector in modei evaluated at the quadra- ture point (q) m respectively. In other words, the equation (3.65) is solved at quadrature points and the values of the frequencies and eigenvectors are obtained on those quadra- ture points. There is a weight associated with each quadrature point that is denoted by w (q) . The coefficients (i) j are vectors that can also be characterized in terms of the 81 eigenvectors associated with the mean of the stiffness matrix. Let (i) denote the eigen- vector associated with the mean of the stiffness matrix K at modei. The (i) is obtained from the following eigenproblem, L (i) = ( ! (i) ) 2 (i) ; (3.68) where L = M 1 K. It should be noted that, eigenvectors (i) are not mutually orthog- onal. In order to characterize the polynomial chaos coefficients of the stochastic eigen- values efficiently, orthogonal (orthonormal) bases are preferred. For this purpose, orthonormalization process such as Gram-Schmit process can be carried out on the non- orthogonal eigenvectors (i) . Let the new orthonormal vectors be denoted by (i) GS . Consequently, the polynomial chaos coefficient (i) j can be written as follows, (i) j = n mode X k=1 < (i) j ; (k) GS > (k) GS ; (3.69) where < :;: > denotes the inner product and n mode is the finite number of modes. The generalized mass and stiffness matrices can also be characterized as a polynomial chaos expansion. LetM (i) andK (i) denote the generalized mass and stiffness of mode i that are diagonal elements in the generalized mass and stiffness matrices. They can be written as follows, M (i) = ( (i) ) T M (i) ; K (i) = ( (i) ) T K (i) : (3.70) 82 It should be noted that, these values are random values due to the randomness in the stiffness matrix and eigenvectors. It should also be noted that, due to the orthogonality relation between modes, the off diagonal values are zero as follows, ( (i) ) T M (j) = 0; ( (i) ) T K (j) = 0: (3.71) Similarly to the expansion for frequencies and eigenvectors, the polynomial chaos expan- sion for generalized mass and stiffness for modei can be given as follows, M (i) ( m )' n M X j=0 M (i) j j ( m ); K (i) ( m )' n K X j=0 K (i) j j ( m ); (3.72) whereM (i) j ;K (i) j are obtained as follows, M (i) j ' nq X q=1 M (i) ( (q) m ) j ( (q) m )w (q) ; K (i) j ' nq X q=1 K (i) ( (q) m ) j ( (q) m )w (q) ; (3.73) 83 andn M ;n K are finite numbers of terms in the expansion. Consequently the generalized mass and stiffness matrices can be formed as follows, M g ( m ) = 2 6 6 6 6 6 6 6 4 M (1) ( m ) 0 ::: 0 0 M (2) ( m ) 0 . . . . . . 0 . . . 0 0 ::: 0 M (n mode ) ( m ) 3 7 7 7 7 7 7 7 5 ; K g ( m ) = 2 6 6 6 6 6 6 6 4 K (1) ( m ) 0 ::: 0 0 K (2) ( m ) 0 . . . . . . 0 . . . 0 0 ::: 0 K (n mode ) ( m ) 3 7 7 7 7 7 7 7 5 : (3.74) In practice, first few modes are considered in the analysis. The damping matrix can also be assumed based on the mass and stiffness matrices (C g =(M g +K g ), = 0:05). Sim- ilarly to the mass and stiffness, the generalized force can be obtained by pre-multiplying the eigenvectors to the aerodynamic forces. It should be noted that both eigenvectors and aerodynamic forces are random values but the sources of randomness are differ- ent. Let w denote the randomness in the wind. Then w is a vector of random vari- ables identical to first few random variables in the Karhunen-Loeve expansion, equa- tion (3.29)( w = 1 ; 2 ;:::). The generalized force can then be given as follows, F (i) (t; m ; w ) = ( (i) ( m )) T F (t; w ); (3.75) 84 whereF (t; w ) is a vector of arodynamic forces given as, F (t; w ) = 2 6 6 6 6 6 6 6 4 f(x 1 ;t; w ) f(x 2 ;t; w ) . . . f(x n dof ;t; w ) 3 7 7 7 7 7 7 7 5 ; (3.76) andn dof is the total number of degrees of freedom in the structure. The generalized force can also be characterized through a polynomial chaos expansion. The eigenvectors are already expanded with polynomials j ( m ) according to equation (3.66). Therefore, the expansion of the generalized force can be given as follows, F (i) (t; m ; w ) = n X j=0 C (i) j (t; w ) j ( m ); (3.77) whereC (i) j (t; w ) are given as, C (i) j (t; w ) = ( (i) j ) T F (t; w ); (3.78) and can also be characterized as an expansion as follows, C (i) j (t; w ) = m X k=0 D (i) jk (t) k ( w ) (3.79) where 0 k s are orthogonal polynomials with respect to the probability density function of w . The deterministic coefficients D (i) jk (t) can be obtained based on the quadrature rule as follows, D (i) jk ' nq X q=1 C (i) j (t; (q) w ) k ( (q) w )w (q) : (3.80) 85 Consequently, the expansion of the generalized force can be given as a double sum as follows, F (i) (t; m ; w ) = n X j=0 m X k=0 D (i) jk (t) j ( m ) k ( w ): (3.81) The displacement of the blade is obtained through the solution of the governing equation given as follows, M g ( m ) X +C g ( m ) _ X +K g ( m )X =F g (t; m ; w ); (3.82) whereF g (t; m ; w ) is given by, F g (t; m ; w ) = 2 6 6 6 6 6 6 6 4 F (1) (t; m ; w ) F (2) (t; m ; w ) . . . F (n mode ) (t; m ; w ) 3 7 7 7 7 7 7 7 5 ; (3.83) The solution of the equation (3.82) is in the form ofX(;t), X g (;t) = 2 6 6 6 6 6 6 6 4 X (1) (;t) X (2) (;t) . . . X (n mode ) (;t) 3 7 7 7 7 7 7 7 5 ; (3.84) whereX (i) (;t) corresponds to the displacement of the modei and = ( m ; w ). The displacement of the blade can then be given as follows, x(;t) = n mode X i=1 X (i) (;t) (i) ( m ): (3.85) 86 From this analysis, the displacement of the blade can be obtained. The displacement can then be used as an input in the detailed finite element model to calculate the stresses in the blade. The details of the polynomial chaos expansion for the solution of the stochastic differential equation (3.72) are explained in the next section. 3.4 Uncertaintypropagation As shown in the previous section, the governing equation based on the modal analysis can be written in the following form, M g ( m ) X +C g ( m ) _ X +K g ( m )X =F g (t; m ; w ): (3.86) The numerical approximation of the stochastic functionX(;t), where = ( m ; w ), can be obtained through Galerkin projection [30]. We can now only focus on the mode i since the modes are decoupled. Therefore, the governing equation for the modei can be written as follows, M (i) ( m ) X (i) +C (i) ( m ) _ X (i) +K (i) ( m )X (i) =F (i) (t; m ; w ): (3.87) The polynomial chaos expansion of the solution can be written as follows, X (i) (;t) = n X j=0 X (i) j (t) j () +" N ; (3.88) 87 where" N is the truncation error. It should be noted that, the solution is a function of time and consequently the deterministic coefficientsX (i) j are functions of time. Substituting the expansion of solution in equation (4.10) results in, n X j=0 M (i) ( m ) d 2 X (i) j (t) dt 2 j () +C (i) ( m ) dX (i) j (t) dt j () +K (i) ( m )X (i) j (t) j () = F (i) (t; m ; w ) +e N ; (3.89) where the errore N is incurred due to the finite-dimensional representation of the solu- tionX. The stochastic Galerkin solution is obtained by requiring the above error to be orthogonal to the approximating space resulting in, n X j=0 d 2 X (i) j (t) dt 2 hM (i) ( m ) j () k ()i + dX (i) j (t) dt hC (i) ( m ) j () k ()i+ X (i) j (t)hK( m ) j () k ()i =hF (t; m ; w ) k ()i; j = 0; 1;:::;n; (3.90) whereh:i denotes the operator of mathematical expectation,h:i = E(:). The above equation is a deterministic equation that is solved for theX (i) j (t) by a standard fourth order Runge-Kutta method. Once theX (i) j (t) are obtained the stochastic representation of the solution is obtained in the form of equation (3.88). It should be noted that each point in time is a random variable, therefore the stress in the blade is a random variable at each time and is obtained from the dynamic displacement. 3.5 Designoptimization The design of aeroelastic structures requires the high reliability on the structural parts and performance of the system. The failure in parts results in an unacceptable cost. To assess the reliability of the system, the stochastic variations in the performance of 88 the system and structural parameters should be considered. The uncertainty in the tra- ditional design has been accounted for through safety factors. This design approach results in a conservative or on the other end unsafe design while a more realistic design is the one with considering uncertainty in the model. The high fidelity aerodynamic models can only reduce modeling uncertainty but in order to characterize the uncer- tainty in the model accurately, stochastic analysis techniques should be integrated into the design process. For design process with uncertainty, the probability distribution of the response of the system should be precisely obtained in order to evaluate the proba- bility of failure in the system. For this purpose the reliability based design optimization is employed [3], [21]. This approach is able to account for uncertainty from different parts in the system. The uncertainty can be considered in the structural model parame- ters and overall system performance that will be explained later in this section. In the sequel, first the deterministic optimization and consequently the details of optimization under uncertainty for the rotor blade of the HAWT are explained. 3.5.1 Deterministicoptimization A deterministic optimization problem in a standard form is as follows, min x z =f(x) subject to G i (x)> 0 i = 1;:::;N g x j min x j x jmax j = 1;:::;N x (3.91) 89 wherez is the cost or objective function andx is the design variable.G i is the constraint function to be satisfied while the design variable is within lower and upper bounds. The limit state function is defined as, G =RS; (3.92) whereR is the capacity specified for the structure andS is the magnitude of the response in a specific point. The objective function in this study is the generated power which is the result of BEM analysis. The instantaneous generated power is given as, P ins (t;) = 4! 2 rot Z R 0 V 0 (r;t)a 0 (r;t;)(1a(r;t;))r 3 dr; (3.93) The equation above is the generated power at each time step and is dependent on the twist of the blade at each section in the span wise direction. These twists are considered as the design variables along the blade. In the optimization process, as design variable changes, different values are obtained for aerodynamic induction factors that will affect the value of power and the structural response of the blade accordingly. It should also be noted that, the value that is used in the optimization process is the average of the power in time as follows, P () = 1 T Z T 0 P ins (t;)dt; (3.94) whereT can be assumed according to the duration of the dynamic analysis. As men- tioned earlier, In order to obtain the stress in the blade , first the reduced order model of the blade based on the modal analysis given in equation (3.95) is constructed and 90 solved. The response of this reduced order model is then used in the detailed finite element model to obtain the stress in the blade. M g X +C g _ X +K g X =F (t): (3.95) The solution of this equation is the dynamic displacement of the structure that can be given as follows, x(t) = n mode X i=1 X (i) (t) (i) : (3.96) It can be seen that the solution of the equation is a deterministic solution and is differ- ent from the equation (3.85). However, the response is still characterized as the linear combination of different modes. The stress in the blade is time dependent. Therefore the reliability assessment should be carried out for the most critical point in the time domain that has the least safety. The stress in the blade can be obtained from the eigen modes since the structure remains in the elastic state. LetS (i) denote the stress in the blade associated with the eigenvector at mode i. The dynamic stress in the blade can then be written as follows, s(t) = n mode X i=1 X (i) (t)S (i) : (3.97) The response s(t) is obtained for all elements in the structure. Therefore, the most critical point is the point that has the largest magnitude of stress in time among all elements. This value is then considered asS in the limit state function in equation (3.92). It should be noted that this value changes from time to time and element to element in different optimization iterations. Therefore, at each iteration step the maximum value of the stress should be searched both in time and in space among all elements. In the next section, details of the optimization under uncertainty are explained. 91 3.5.2 Optimizationunderuncertainty For a system with uncertainty, a reliability based design optimization is maximizing the probability of the benefit while satisfying probabilistic and deterministic constraints by varying the design variables within lower and upper bounds. An optimization problem under uncertainty can be written as follows, max x z =Pr(f(x)> 0) subject to Pr(G i (x)> 0)> Pr i i = 1;:::;N g x j min x j x jmax j = 1;:::;N x (3.98) In the equations above, the objective functionf(x) and limit state functionsG i are random. As mentioned earlier, there might be deterministic constraints in addition to probabilistic constraints in the problem. The objective is to maximize the probability of the objective function exceeding a threshold. The constraint is to maintain the reliabil- ity level. In other words, the probability of safety remains above a threshhold, where the safety is denoted by G i (x) > 0 and Pr i is the threshhold for the ith constraint. Similarly to the deterministic optimization, the design variablesx j varies in a range of [x j min ;x jmax ]. As mentioned in the previous section, the randomness can be considered in the structural parameters and the wind. The objective function, the generated power is then random, through the randomness in the wind. The randomness in the wind prop- agates also in the aerodynamic induction factors. Consequently the random form of the instantaneous power can be written as follows, P ins (t;; w ) = 4! 2 rot Z R 0 V 0 (r;t; w )a 0 (r;t;; w )(1a(r;t;; w ))r 3 dr; (3.99) 92 The instantneous power can then be averaged in time and given by, P (; w ) = 1 T Z T 0 P ins (t;; w )dt: (3.100) As mentioned earlier, this averaged generated power is dependent on the values of the twist along the blade. In this case, it is a random variable and the distribution changes as design variables change. The objective function can then be defined as the probability level, z, that should be maximized. The probability level is written as follows, z =Pr(P ()> P ); (3.101) where P is the specified value for the generated power. In order to evaluate the reliability the response of the blade is needed. As mentioned in the previous section, the governing equation of the motion of the blade can be written as follows, M g ( m ) X +C g ( m ) _ X +K g ( m )X =F g (t; m ; w ): (3.102) The response of the blade can then be characterized as an expansion with the following form, x(;t) = n mode X i=1 X (i) (;t) (i) ( m ): (3.103) In this case, in addition to dependence on time, the stress state of the blade is random due to the randomness in the response of the blade. Therefore, in order to search for the most critical point in time, the mean values of the stress are considered. After the critical point is found, then the full representation of the random stress for that specific point is considered in the optimization process. The mean of stress can be obtained from the mean of the responseX (i) 0 (t) and the mean of the eigenmodes (i) 0 . It should be noted that the mean of eigenmodei, (i) 0 is different from the eigenmodei associated 93 with the mean of the stiffness matrix, (i) , introduced in equation (3.68). Similarly to the deterministic case, letS (i) 0 denote the stress in the blade associated with the mean of the eigenmodei. The mean of stress in the blade can then be given as follows, s 0 (t) = n mode X i=1 X (i) 0 (t)S (i) 0 : (3.104) The mean value of the stress is obtained for all elements and is a function of time. The most critical point can then be searched in all elements and in time. It should again be noted that, this value is not used in the optimization process. The sole purpose to obtain the mean value is to search for the critical point. It should be noted thatS (i) 0 is a vector and its size is the number of elements. LetS cr (i) j denote the value corresponding to the critical point throughout the elements in the vector S (i) j . The vector S (i) j is the stress in the blade associated with the polynomial chaos coefficients (i) j of the random eigenvectors (i) ( m ) in equation (3.66). The polynomial chaos expansion ofS cr (i) ( m ) can then be given as, S cr (i) ( m )' n X j=0 S cr (i) j j ( m ); (3.105) In addition, let the critical point in time be denoted byt cr . The random response of the governing equation from equation (3.88) can then be denoted byX cr (i) () and given by (neglecting the truncation error), X cr (i) () =X (i) (;t cr )' n X j=0 X (i) j (t cr ) j (): (3.106) Finally, the random representation of the stress at the critical point is given by, S cr () = n mode X i=1 X cr (i) ()S cr (i) ( m ): (3.107) 94 Again, it should be noted that the randomness in the response of the governing equa- tion is from both material and wind sources, = ( m ; w ). But, the randomness in the eigenvectors and the stresses resulted form the deformations associated with these eigenvectors are from material properties m . The complete representation is therefore dependent on both sources of randomness = ( m ; w ). The limit state function for the reliability based optimization can be written as follows, G() =R( m )S(); (3.108) whereG is the limit state function andR is the capacity that is considered to be random and dependent on the randomness in the material properties. For the structure to be safe the probability of safety should be higher than a target value. Therefore, the constraint in the optimization problem under uncertainty can be considered in the following form, Pr(G> 0)> Pr; (3.109) where Pr is the target probability of safety. One way to formulate the reliability based optimization is employing First-order (FORM) or Second-order (SORM) reliability meth- ods. Both approaches search for the most probable point (MPP) on the limit state func- tion. FORM uses a first order approximation of the failure function. In other words, it is based on a first order Taylor series expansion of the performance (failure) function linearized at the mean value of the random variables. Therefore, FORM can be accurate if the limit state function at the MPP is not highly nonlinear. For the cases that fail- ure surface in the uncertainty space exhibits nonlinearity, SORM approach can be more accurate since it uses first and second order derivatives to build a failure surface at MPP. Accordingly, the computational cost for the SORM is higher. In this work the reliability 95 assessment is done based on the polynomial chaos expansion of the limit state func- tion. This representation allows us to precisely compute the probability of failure from the distribution of the limit state function. It also enables us to get higher accuracy by employing higher order representation in the solution. The flowchart shown in Fig. 3.8 depicts the reliability based design optimization framework in this study. The details of the optimization of the blade design is explained and illustrated through a numerical example in the next section. Figure 3.8: Reliability based design optimization of the HAWT blade 96 Chapter4 Anumericalstudyofdesign optimizationofHAWTrotorblade In this section, different steps of the dynamic analysis and design optimization of the turbine blade is presented. The proposed aeroelastic analysis and optimization under uncertainty is illustrated by a numerical example on a 3D model of a blade. The values of lift and drag coefficients are in different radius along the span are form the AWT(advance wind turbine) blade data of National Renewable Energy Laboratory (NREL). As the first step, the generation of the wind is illustrated in the next section. 4.1 WindModel In order to generate the wind, V on Karman power spectral density function is considered as follows, S(!) = 4(5:7u 2 )L x =U(x) 1:339(1 + 39:48(!L x =U(x)) 2 ) 5 6 : (4.1) In the above equation,U(x) is the mean speed at pointx.andL x is the integral length scale and considered to be 600m. The mean velocity at different points can be obtained relative to the mean velocity at heighth as follows, U(x) U(h) = ln(x=z 0 ) ln(h=z 0 ) ; (4.2) 97 wherez 0 is called the roughness length and is considered 0:03m. For the wind model, the blade is divided to ten sections with lengthdr' 4m. The mean speed at the hub height,h = 100m is considered to be 22m=s. In equation (4.1),u is the shear velocity and given by, u = 0:4U(x) ln(x=z 0 ) ; (4.3) With the given specifications, the value of the shear velocity can then be obtained u = 1:2739. The value of the power spectral density at each spatial point is obtained from equation (4.1). As mentioned earlier, the values of cross spectral densities can be obtained through the assumed coherence function as follows, Coh ij (!) = exp(12(!4r ij =U ij )); (4.4) where4r ij is the spacing between two points and U ij is the average of mean speed between two points. The full spectral density matrix can then be obtained having the values of the spectral density at each point and cross spectral densities. The wind can be generated from equation (3.22). Ten modes are considered in the expansion of the spec- tral density matrix, (n s = 10 in equation (3.22)). The number of discrete points in the frequency domain is 201 which corresponds ton d = 100. This construction, results in 2010 independent standard normal random variables denoted by jk in equation (3.22) that can be seen as the number of random dimensions in the process. As described ear- lier, the number of random dimensions in the generated process can be reduced through a Karhunen-Loeve expansion. In this expansion, only first few significant modes are considered. The number of modes, n in equation (3.29) is chosen such that, the ratio P n i=1 i = P N i=1 i is close to 1, whereN >>n. The first 100 eigenvalues, i are plotted in the Fig. 4.1. 98 0 10 20 30 40 50 60 70 80 90 100 10 −1 10 0 10 1 10 2 10 3 10 4 Mode indices Eigenvalues Figure 4.1: Eigenvalues of the covariance matrixR The above ratio forn = 5 andN = 100 is 0:81. Consequently, the number of modes in the Karhunen-Loeve expansion is considered to ben = 5 in equation (3.29). It should be noted, in the figure, only eigenvalues for the first 100 modes out of 2000 modes are plotted. The first five eigenmodes associated with physical domain (denoted by i (z) in equation (3.29)) of the process are plotted in Fig. 4.2. The histograms of the random variables associated with these eigenmodes (denoted by i (z) in equation (3.29)) are plotted in Fig. 4.3. As expected, these random variables are obtained as standard normal random variables since the original random process is a gaussian random process. 99 0 20 40 0 5 10 0.015 0.02 0.025 0.03 x(m) t(s) γ 1 0 20 40 0 5 10 −0.04 −0.02 0 0.02 0.04 x(m) t(s) γ 2 0 20 40 0 5 10 −0.04 −0.02 0 0.02 0.04 x(m) t(s) γ 3 0 20 40 0 5 10 −0.04 −0.02 0 0.02 0.04 x(m) t(s) γ 4 0 20 40 0 5 10 −0.06 −0.04 −0.02 0 0.02 0.04 x(m) t(s) γ 5 Figure 4.2: Spatial-temporal eigenmodes of the Karhunen-Loeve expansion of the wind process 4.2 AeroelasticModel The Aeroelastic analysis is done to determine the aerodynamic forces on the structure. As described earlier, the aerodynamic forces are obtained through aerodynamic induc- tion factorsa anda 0 . Inputs for this analysis at this stage are the drag and lift coefficients, the twist and chord distribution as well as the upstream velocity at each time. The rota- tional velocity of the blade is considered(rad=s). The polynomial chaos expansion of 100 −5 0 5 0 100 200 300 400 η 1 Frequency −5 0 5 0 100 200 300 400 η 2 Frequency −5 0 5 0 100 200 300 400 η 3 Frequency −5 0 5 0 100 200 300 400 η 4 Frequency −5 0 5 0 100 200 300 400 η 5 Frequency Figure 4.3: Histograms of the random variables associated with the eigenmodes in Karhunen- Loeve expansion of the wind process the aerodynamic induction factors can be obtained similarly to the any random quantity such as generalized mass and stiffness. These expansions can be written as follows, a(r;t; w )' na X j=0 a j (r;t) j ( w ); a 0 (r;t; w )' n a 0 X j=0 a 0 j (r;t) j ( w ); (4.5) 101 wherea j (r;t);a 0 j (r;t) are obtained as follows, a j (r;t)' nq X q=1 a(r;t; (q) w ) j ( (q) w )w (q) ; a 0 j (r;t)' nq X q=1 a 0 (r;t; (q) w ) j ( (q) w )w (q) ; (4.6) The only quantity needed to obtain the polynomial chaos coefficients according to equa- tion (4.6) is the value of the induction factors on the quadrature points denoted by a(r;t; (q) w );a(r;t; (q) w ). The same approach can be followed to obtain the random repre- sentation of the aerodynamic forces and generated power. The distributions of the aero- dynamic factors,a anda 0 at location 0:8R and timet = 2:5sec are plotted in Fig. 4.4. Similarly, the distributions of the normal and tangential forces are plotted in Fig. 4.5. 0.04 0.06 0.08 0.1 0.12 0 50 100 150 200 250 300 350 400 a Frequency 2.4 2.6 2.8 3 3.2 3.4 x 10 −3 0 100 200 300 400 500 600 700 800 900 1000 a’ Frequency Figure 4.4: Distribution of the aerodynamic induction factors at 0.8 of the blade radius and t = 2:5sec The distribution of the generated power corresponding to the initial design variables (twist angles) is plotted in Fig. 4.6. 102 7000 7500 8000 8500 9000 0 50 100 150 200 250 300 350 400 P N /m Frequency 500 1000 1500 2000 0 50 100 150 200 250 300 350 400 450 P T /m Frequency Figure 4.5: Distribution of the aerodynamic forces on the section at 0.8 of the blade radius nd t = 2:5sec 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 x 10 6 0 50 100 150 200 250 300 350 400 450 Power(W) Frequency Figure 4.6: Distribution of the generated power 103 4.3 StructuralModel The dynamic response of the model is obtained through a reduced order model as described earlier. In order to show the significance of first modes, the fast decay of the square of the natural periods are plotted in Fig. 4.7. It should be noted that, the first significant modes are corresponding to the modes with smallest natural frequen- cies, equivalently largest natural periods. The natural periods are expressed in terms of the natural frequencies! (i) (the eigenvalues in equation (3.65)) as follows, T (i) = 2 ! (i) : (4.7) It should also be noted that, results in Fig. 4.7 correspond to the mean value of the 0 5 10 15 20 25 30 10 −3 10 −2 10 −1 10 0 10 1 Modes T 2 Figure 4.7: Square of the natural periods of the blade corresponding to the mean value of the elastic modulus elastic modulus and show the first 30 modes of the blade. However, in practice smaller number of modes are considered in the response analysis. As a quantitative measure, 104 the ratio P n i=1 (T (i) ) 2 = P 30 i=1 (T (i) ) 2 wheren<< 30 is considered to choose the number of significant modes. In this analysis, first three eigenmodes are considered since the above ratio is 0:97 forn = 3. The random form of the generalized mass and stiffness in the first mode according to the equation (3.72) is as follows, M (1) ( m )' 2:7538 0:0181 m 1 + 0:0132 m 2 + 0:0048 m 3 ; K (1) ( m )' 18:616 + 0:5014 m 1 + 0:2976 m 2 + 0:1323 m 3 : (4.8) These expansions are obtained through performing stochastic collocation as mentioned earlier. The values of the generalized mass and stiffness are obtained on quadrature points and consequently the polynomial chaos coefficients are obtained according to equation (3.73). According to equations (3.66) and (3.69), the stochastic eigenmodes can also be written in terms of the deterministic eigenmodes as, (i) ( m ) = n mode X j=1 ij ( m ) (j) GS (4.9) The histogram of each random variable ij in equation (4.9) is plotted in Fig. 4.8. The mean and standard deviation of these random variables are listed in Table 4.1. Table 4.1: Mean and standard deviation of the random coefficients obtained from the projection of stochastic modes on deteminisitc eigenmodes Mean Standard deviation 11 0.9999 2.4659e-05 12 -1.3887e-05 0.0061 13 -2.7928e-05 0.0018 21 -0.2494 0.0044 22 0.9683 0.0011 23 -1.9942e-05 0.0057 31 -0.3885 0.0035 32 0.0525 0.0060 33 0.9198 0.0018 105 0.9998 0.9999 1 0 0.5 1 1.5 2 x 10 4 α 11 Frequency −0.02 0 0.02 0 500 1000 1500 2000 α 12 Frequency −10 −5 0 5 x 10 −3 0 1000 2000 3000 α 13 Frequency −0.26 −0.25 −0.24 −0.23 0 1000 2000 3000 α 21 Frequency 0.965 0.97 0.975 0.98 0 500 1000 1500 2000 α 22 Frequency −0.02 0 0.02 0 500 1000 1500 2000 α 23 Frequency −0.4 −0.38 −0.36 0 1000 2000 3000 α 31 Frequency 0.02 0.04 0.06 0.08 0 500 1000 1500 2000 α 32 Frequency 0.91 0.92 0.93 0 1000 2000 3000 α 33 Frequency Figure 4.8: Histograms of the random coefficients obtained from the projection of stochastic modes on deteminisitc eigenmodes In addition to eigenmodes, polynomial chaos expansion of random eigenvalues can be obtained as discussed earlier. In figure4.8, the distribution of each random eigenvalue (square of the natural frequency of the structure, i = (! (i) ) 2 ) is plotted. Table 4.2 shows the mean and standard deviation of these random eigenvalues. The eigenvalues of the deterministic blade model associated with the mean of the elastic modulus are obtained as 1 = 6:7628; 2 = 27:216; 3 = 106:85. From Table4.2, it can be seen that these values are very close to the mean of the random eigen- values of the blade. The distribution of the generalized force (F (i) g in equation (3.75)) of 106 6 6.5 7 7.5 0 200 400 600 800 1000 1200 1400 1600 λ 1 Frequency 24 26 28 30 0 200 400 600 800 1000 1200 1400 1600 1800 λ 2 Frequency 100 105 110 115 0 500 1000 1500 2000 2500 λ 3 Frequency Figure 4.9: Histograms of the random eigenvalues corresponding to the first three modes of the blade Table 4.2: Mean and standard deviation of the random eigenvalues corresponding to the first three modes of the blade Mean Standard deviation 1 6.7612 0.1396 2 27.209 0.5163 3 106.83 1.8543 different modes att = 7:5sec is plotted in Fig. 4.10. The histogram of the solution of the stochastic governing equation,X (i) (;t) in equation (3.85), att = 7:5sec is plotted in Fig. 4.11. 107 −2000 −1500 −1000 0 50 100 150 200 250 300 350 400 450 F g (1) Frequency 1000 1500 2000 2500 0 50 100 150 200 250 300 350 400 450 500 F g (2) Frequency −3000 −2500 −2000 0 50 100 150 200 250 300 350 400 F g (3) Frequency Figure 4.10: Histograms of the Generalized force in three modes 4.4 ConvergenceAcceleration In this section, results of the convergence acceleration of the equation of motion for the blade are presented. As described earlier, the reduced order model describing the motion of the blade can be written as follows, M (i) ( m ) X (i) +C (i) ( m ) _ X (i) +K (i) ( m )X (i) =F (i) (t; m ; w ): (4.10) The solution of the above stochastic ordinary differential equation is characterized as a polynomial chaos expansion. The stochastic Galerkin projection is used to obtain this solution. As described earlier, in the stochastic Galerkin solutions, the accuracy of the solution is improved as the order of the expansion increases. However, more accurate 108 −80 −75 −70 −65 0 50 100 150 200 250 300 350 X (1) Frequency 20 25 30 0 50 100 150 200 250 300 X (2) Frequency −7 −6 −5 −4 0 50 100 150 200 250 300 350 X (3) Frequency Figure 4.11: Histograms of the solution of the governing equation in three modes solutions can be obtained from the lower order expansions through numerical transfor- mations which are described in this dissertation. Similarly to the numerical examples of the stochastic initial value problem and stochastic boundary value problem, convergence is investigated relative to the error measure which is defined for each polynomial chaos coefficient of the solutionX (i) (t) as follows, " = v u u t ( 200 X j=1 (X (i) (t j )X (i) ;6 (t j )) 2 )=( 200 X j=1 (X (i) ;7 (t j )) 2 ) 8; (4.11) wheret j =j dt;dt = 0:05s. Results in this section are presented for the initial design variables. The accuracy is measured relative to the highest order expansionp = 6 since the exact solution is not available. The relative errors in the mean of the first mode 109 solution corresponding to different polynomial orders and transformed solutions using e 1 transformation are listed in Table 4.3 and plotted in Fig. 4.12. The application of the e 1 transformation results in a higher accuracy as shown in this section and in the numerical results in Chapter 2. In order to illustrate the computational cost associated with this transformation, the computational effort in terms of the number of operations based on the similar analysis in Chapter 2 is listed in Table 4.4. It should be noted that the stochastic dimension in this problem iss = 8. The relative error in the mean of the solution is plotted against the computational effort in Fig. 4.13 for further illustration. Table 4.3: Convergence acceleration results of the equation of motion of the blade: Relative error in the mean of the first mode solutionX (1) (t; m ; w ) p PCE e 1 1 1.38e-02 2 2.91e-05 3 4.90e-08 4.12e-09 4 5.71e-11 4.39e-12 Table 4.4: Computational effort in terms of estimated number of operation p No. of coefficients r PCE e 1 1 9 200 1.800e+03 2 45 200 9.000e+03 3 165 200 3.300e+04 4.380e+04 4 495 200 9.900e+04 1.410e+05 5 1287 200 2.574e+05 3.894e+05 The convergence acceleration results of polynomial chaos coefficients associated with the linear terms w 1 ;:::; m 3 are plotted in Figures 4.14,4.15,4.16,4.17. Similarly, the convergence acceleration results of polynomial chaos coefficients asso- ciated with higher order polynomials 2 w 1 and w 2 m 1 are plotted in Fig. 4.18. 110 1 2 3 4 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error PCE e 1 Figure 4.12: Convergence acceleration results of the equation of motion of the blade: Relative error in the mean of the first mode solutionX (1) (t; m ; w ) 4.5 SensitivityAnalysis The results of this section is the sensitivity of the governing equation of the blade (equa- tion (3.82)) with respect to the uncertainty sources. Results correspond to the first mode and the initial design variables. The governing equation of the first mode is written as follows, M (1) ( m ) X (1) +C (1) ( m ) _ X (1) +K (1) ( m )X (1) =F (1) (t; m ; w ): (4.12) The solution of the above differential equation at the critical time is considered for the sensitivity analysis purposes. Prior to carrying out the sensitivity analysis, the knowl- edge of the magnitude of the input parameters is also useful. The parameters for the 111 10 3 10 4 10 5 10 6 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 Computational effort (Op. count) Relative Error PCE e 1 Figure 4.13: Relative error in the mean of the solution versus computational effort defined as operation count mass and stiffness are given in equation.(4.8). The RHS of the equation can be written as follows, F (1) (t; m ; w ) = F (1) 0 (t) +F (1) 1 (t) w 1 +::: +F (1) 5 (t) w 5 + F (1) 6 (t) m 1 +F (1) 7 (t) m 2 +F (1) 8 (t) m 3 +::: (4.13) Since the input parameters in the RHS are in the form of time histories, the root mean square (RMS) of each parameter is obtained in order to be able to compare the intensity of input parameters. The RMS of the second parameter can be written as follows, RMS 1 = v u u u t n X i=1 (F (1) 1 (t i )) 2 n ; (4.14) 112 1 2 3 4 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error 1 2 3 4 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error PCE e 1 PCE e 1 (a) (b) Figure 4.14: Convergence acceleration results of the equation of motion of the blade: (a) Rel- ative error in polynomial chaos coefficient associated with w 1 , (b) Relative error in polynomial chaos coefficient associated with w 2 in the first mode solutionX (1) (t; m ; w ) where n is the number of points in the time history (n = 200). The RMS values are given in the table below, The distribution of the solution at the critical point is plotted in Table 4.5: Values of the root mean square of the loading parameters RMS 1 RMS 2 RMS 3 RMS 4 RMS 5 RMS 6 RMS 7 RMS 8 38.824 2.737 9.390 1.030 9.997 6.389 6.519 0.132 Fig. 4.19. The first order Sobol indices denoted bySI 1 ;SI 2 ;:::;SI 8 , are the sensitivity of the output with respect to uncertainty sources w 1 ;:::; w 5 ; m 1 ; m 2 ; m 3 respectively. It should be noted that, the variation in the response of the system is not directly depen- dent on the magnitude of the parameters on the RHS of the equation. The variation in the material property results in uncertain parameters in the LHS of the governing equation 113 1 2 3 4 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error 1 2 3 4 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 P Relative Error PCE e 1 PCE e 1 (a) (b) Figure 4.15: Convergence acceleration results of the equation of motion of the blade: (a) Rel- ative error in polynomial chaos coefficient associated with w 3 , (b) Relative error in polynomial chaos coefficient associated with w 4 in the first mode solutionX (1) (t; m ; w ) that also contributes to the total system response variation. Therefore, the importance of each uncertainty source can not be assessed by only ranking the magnitude of the RHS coefficients and a formal sensitivity analysis is required. The first order sobol indices are listed in Table 4.6 for two cases of critical time corresponding to the initial design variables (t cr = 0:7sec) and an arbitrary timet = 7:5sec. It can be seen from Table 4.6 that, the variation in the first uncertainty source w 1 in the wind has the most significant contribution in the variation of the output response. In addition, the first uncertainty source in the material property m 1 , has a considerable importance. 114 1 2 3 4 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error 1 2 3 4 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error PCE e 1 PCE e 1 (a) (b) Figure 4.16: Convergence acceleration results of the equation of motion of the blade: (a) Rel- ative error in polynomial chaos coefficient associated with w 5 , (b) Relative error in polynomial chaos coefficient associated with m 1 in the first mode solutionX (1) (t; m ; w ) Table 4.6: Values of the sensitivity indices corresponding to different sources of uncer- tainty SI 1 SI 2 SI 3 SI 4 SI 5 SI 6 SI 7 SI 8 t =t cr 0.4787 0.0051 0.0562 8.2377e-04 0.0599 0.2218 0.1572 0.0200 t = 7:5sec 0.6295 0.0013 0.0302 3.9019e-05 0.0207 0.2980 0.0083 1.0845e-04 4.6 Designoptimization Results for design optimization are presented for two scenarios, deterministic optimiza- tion and optimization under uncertainty. In each case, the design variables are consid- ered to be the twist angle along the blade in different sections. For this purpose, the 115 1 2 3 4 4 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error 1 2 3 4 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error PCE e 1 PCE e 1 (a) (b) Figure 4.17: Convergence acceleration results of the equation of motion of the blade: (a) Rel- ative error in polynomial chaos coefficient associated with m 2 , (b) Relative error in polynomial chaos coefficient associated with m 3 in the first mode solutionX (1) (t; m ; w ) blade is divided in ten equally spaced sections and a twist angle is assigned to each section. 116 2 3 4 10 −8 10 −6 10 −4 10 −2 10 0 P Relative Error 2 3 4 10 −10 10 −8 10 −6 10 −4 10 −2 P Relative Error PCE e 1 PCE e 1 (a) (b) Figure 4.18: Convergence acceleration results of the equation of motion of the blade: (a) Rel- ative error in polynomial chaos coefficient associated with 2 w 1 , (b) Relative error in polynomial chaos coefficient associated with w 2 m 1 in the first mode solutionX (1) (t; m ; w ) 117 −185 −180 −175 −170 −165 −160 −155 −150 −145 0 50 100 150 200 250 300 350 X (1) (t cr ) Frequency Figure 4.19: Distribution of the solution of the governing equation at the critical time 4.6.1 DeterministicOptimization The initial design variables i (twists of the blades in different radiuses) and the allow- able ranges are given as follows, Normalized radius = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0:5 0:15 0:25 0:35 0:45 0:55 0:65 0:75 0:85 0:95 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; i = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0:6380 0:3261 0:2360 0:1331 0:1179 0:0850 0:0564 0:0286 0:0103 0:0013 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; [ min ; max ] = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0:63; 1:00 0:32; 0:63 0:23; 0:32 0:13; 0:23 0:11; 0:13 0:08; 0:11 0:05; 0:08 0:02; 0:05 0:01; 0:02 0:00; 0:01 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : 118 The three dimensional blade model is shown in Fig. 4.20. The objective function in the 0 10 20 30 40 50 −2 0 2 −1 0 1 2 3 4 X Y Z Figure 4.20: Three dimensional blade model corresponding to the initial design variables optimization is the averaged generated power given by equation (3.94). The duration of the analysis,T is considered to be 10 seconds. In this case, the wind values are from one 119 realization of the generated random wind in equation (3.29). The optimization problem can be stated as follows, max P () subject to RS()> 0 min max 2R 10 (4.15) The maximum allowable stress (capacity), R is set to be 20 Mpa. Table 4.7 is the summary of the results for the deterministic optimization. After the optimization, the Table 4.7: Deterministic optimization results Critical Stress(Mpa) Power(KW) Initial design 19.64 1049.84 Final design 19.25 1102.58 120 magnitude of the stress value stays within the limit and the power is increased. The values of the twist angle corresponding to the final design are as follows, f = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0:7808 0:4321 0:2797 0:1703 0:1100 0:0859 0:0500 0:0311 0:0100 0:0000 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : (4.16) The blade model corresponding to the final deterministic design is shown in Fig. 4.21. The initial and final twist values corresponding to deterministic optimization, are plotted in Fig. 4.22. 121 0 10 20 30 40 50 −2 0 2 −1 0 1 2 3 4 X Y Z Figure 4.21: Three dimensional blade model corresponding to the final deterministic design 4.6.2 OptimizationunderUncertainty The initial design variables and the allowable ranges are identical to the deterministic optimization problem. The optimization under uncertainty can be written as follows, max Pr(P ()> P ) subject to Pr(G> 0)> Pr min max 2R 10 (4.17) 122 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Normalized radius Twist Initial Final (Det. optimization) Figure 4.22: Twist values of initial and final design corresponding to deterministic optimization whereP () and G are given by equations (3.100) and (3.108) respectively. The maxi- mum allowable stress is considered to be a random variable given as follows, R( m ) = 7( m 1 + m 2 + m 3 ) Mpa; (4.18) where m i ’s are uniform random variables with the range [0:95; 1:05] equivalent to the random variables considered in the material properties of the blade. The distribution of power for initial and final design is plotted in Fig. 4.23. The distribution of the stress at the critical point and the limit state function corresponding to the initial and final design is shown in Fig. 4.24. The value of P is the mean of the generated power of the initial design, P = 1040:15KW . The value of Pr is set to be 0:99. The initial values of design 123 variables are identical to the previous case. The final values of design variables are as follows, 5 6 7 8 9 10 11 12 13 14 x 10 5 0 200 400 600 Power initial (W) Frequency 5 6 7 8 9 10 11 12 13 14 x 10 5 0 200 400 600 Power final (W) Frequency Figure 4.23: Distribution of the generated power of the initial and final design f = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0:7894 0:4341 0:2706 0:1577 0:1183 0:0909 0:0552 0:0339 0:0159 0:0039 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : (4.19) 124 1.6 1.8 2 2.2 x 10 7 0 100 200 300 400 S i Frequency −2 0 2 4 6 x 10 6 0 100 200 300 400 G i Frequency 1.6 1.8 2 2.2 x 10 7 0 100 200 300 400 S f Frequency −2 0 2 4 6 x 10 6 0 100 200 300 400 G f Frequency Figure 4.24: Distribution of the stress and limit state function of the initial and final design The blade model corresponding to the final design optimization under uncertainty is shown in Fig. 4.25. The initial and final twist values corresponding to optimization under uncertainty, are plotted in Fig. 4.26. The summary of the results of the optimiza- tion under uncertainty is listed in Table 4.8. Table 4.8: Optimization under uncertainty results Pr(G> 0) Pr(P > P ) Initial design 0.9979 0.5182 Final design 0.9973 0.7564 As mentioned earlier, in order to have a reliable design the constraint is defined such that the probability of safety exceeds a target valuePr(G > 0) > 0:99. It can be seen that, this condition is satisfied in the final design while the generated power is increased. The mean value of the generated power corresponding to the initial and final design is 125 0 10 20 30 40 50 −2 0 2 −1 0 1 2 3 4 X Y Z Figure 4.25: Three dimensional blade model corresponding to the final design optimization under uncertainty 1040:15KW and 1115:85KW respectively, that further illustrates the maximization of the generated power. 126 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Normalized radius Twist Initial Final (Opt. under uncertainty) Figure 4.26: Twist values of initial and final design corresponding to optimization under uncer- tainty 127 Chapter5 SummaryandFuturework A computational framework for design optimization of the horizontal-axis wind tur- bine rotor blade is developed. The framework incorporates the effects of uncertainty in the optimization process. Specifically, the generated power of the blade is maxi- mized which is subjected to the reliability constraints based on the allowable stress in the blade. The aerodynamic shape of the blade, namely twist angle along the blade is optimized. An accurate reliability assessment is developed employing polynomial chaos expansion. The reliability assessment is based on the dynamic response of the blade sub- ject to the aerodynamic forces. The response of the blade and consequently, the stress state of the blade structure is obtained by constructing the reduced order model of the structure based on modal analysis. A probabilistic model is developed that propagates the randomness inherent in the wind to the aerodynamic induction factors and gener- ated power as well as the randomness in the material property to the parameters of the reduced order structural model. These probabilistic models are achieved without mod- ifying the underlying forward deterministic model through a non-intrusive polynomial chaos approach. The output response of the probabilistic reduced order model of the blade is then obtained by an intrusive approach namely, stochastic Galerkin. The opti- mization under uncertainty as well as deterministic optimization are illustrated through a numerical example of a wind turbine model. The proposed design optimization under uncertainty framework can be extended to a more advanced optimization framework by employing a more detailed fluid structure interaction model. For this purpose, a detailed wind flow analysis based on Navier-Stokes equations should be coupled with a detailed 128 structural model. Consequently, these developments require the application of high per- formance computing that needs to be addressed. In this dissertation, the challenge of high computational cost associated with intrusive polynomial chaos approaches in high dimensional parameter space is addressed by pursuing an innovative methodology for convergence acceleration of approximations to polynomial chaos coefficients using non- linear sequence transformations. Two classes of linear stochastic differential equations, a stochastic IVP and a stochastic BVP, are used to exemplify the method. In each case, the random coefficient in the differential equation is specified as a random process and consequently, for each polynomial chaos coefficient, a convergent sequence is obtained with respect to different polynomial orders in the expansion. The computational cost grows very fast in problems with a large stochastic dimension as the expansion order increases. The methodology pursued in this paper, is shown to be an effective way of increasing the maximum accuracy attainable using available computational resources. Nonlinear transformations which only utilize elements of the available sequence are well suited for improving accuracy while no further information is required. While Shanks and Levin transformations have provable performance only under conditions of exponential convergence rate with respect to polynomial order, they are observed to pro- vide excellent results for conditions where the exact convergence behavior is unknown. In most cases, it is shown that applying the transformation provides a more accurate result together with less computational cost, however in all cases, improvement in the accuracy is demonstrated. The procedure is not limited to the type of equations that have been studied in this dissertation and can be implemented on various problems that incorporate polynomial chaos expansion. It is expected, however that extensions into nonlinear regimes of behaviors would necessitate the development of new sequence transformations that are adapted to the particular rates of convergence. It is also noted that the one-parameter representation based on the polynomial order does not apply to 129 non-intrusive polynomial chaos solutions, since the accuracy of the approximation in that case depends solely on the quadrature rule used and not on the order of polynomial chaos representation. In this case, an appropriate parametrization may be achieved by considering transients induced by a refinement in the order of numerical quadrature. 130 ReferenceList [1] Andres Ahlstrom. Aeroelastic simulation of wind turbine dynamics. Doctoral thesis in structural mechanics, 2005. [2] Alexander Aitken. On Bernoulli’s numerical solution of algebraic equations. Pro- ceedings of the Royal Society of Edinburgh, 46:289–305, 1926. [3] M. Allen and K Maute. Reliability-based shape optimization of structures undergo- ing fluid-structure interaction phenomena. Computer Methods in Applied Mechan- ics and Engineering, 194(30):3472–3495, 2004. [4] I. Babuska, R. Tempone, and GE Zouraris. Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM Journal on Numerical Analysis, 42(2):800–825, 2004. [5] Ivo Babuska, Ral Tempone, and Georgios Zouraris. Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation. Computer Methods in Applied Mechanics and Engineer- ing, 194(12-16):1251 – 1294, 2005. [6] Yr Bazilevs, MC Hsu, Jc Kiendl, Rc Wuechner, and KU Bletzinger. 3d simulation of wind turbine rotors at full scale. part ii: Fluid-structure interaction. International Journal of Numerical Methods in Fluids, (65):236 – 253, 2011. [7] Yrr Bazilevs, MC Hsu, Ic Akkerman, Sc Wright, Kc Takizawa, Bc Henicke, Tc Spielman, and TE Tezduyar. 3d simulation of wind turbine rotors at full scale. part i: Geometry modeling and aerodynamics. International Journal of Numerical Methods in Fluids, (65):207 – 235, 2011. [8] F. Benth and J. Gjerde. Convergence rates for finite element approximations of stochastic partial differential equations. Stochastics and stochastics reports, 63:313–326, 1998. [9] M. Berveiller, B. Sudret, and M. Lemaire. Stochastic finite element: a non intrusive approach by regression. European Journal of Computational Mechanics, 15(1-3), 2006. 131 [10] HG Beyer and BC Sendhoff. Robust optimization - a comprehensive survey. Com- put Methods Appl Mech Eng, 196(33-34):31903218, 2005. [11] C. Brezinski. Extrapolation algorithms and pad approximations: a historical sur- vey. Applied Numerical Mathematics, 20(3):299 – 318, 1996. [12] C. Brezinski. Convergence acceleration during the 20th century. J. Comput. Appl. Math., 122:1–21, October 2000. [13] R. Bulirsch and J. Stoer. Numerical treatment of ordinary differential equations by extrapolation methods. Numer. Math., 8(1):1 – 13, 1966. [14] RH Cameron and WT Martin. The orthogonal development of nonlinear function- als in series of fourier hermite functionals. Annalsof Mathematics, (48):385 – 392, 1947. [15] L. Carassale, G. Solari, and F. Tubino. Proper orthogonal decomposition in wind engineering, part 2: Theoretical aspects and some applications. Wind and Struc- tures, 10(2):177–208, 2007. [16] PK Chaviaropoulos. Development of a state-of-the-art aeroelastic simulator for horizontal axis wind turbines. part 1. structural aspects. Wind Engineering, 25(6):405 – 422, 1996. [17] S. Das, R. Ghanem, and J. Spall. Asymptotic sampling distribution for polyno- mial chaos representation of data : A maximum-entropy and fisher information approach. SIAM Journal on Scientific Computing, in press, 2008. [18] C. Descelliers, R. Ghanem, and C. Soize. Maximum likelihood estimation of stochastic chaos representation from experimental data. International Journal for Numerical Methods in Engineering, 66(6):978–1001, 2006. [19] A. Doostan, R. Ghanem, and J. Red-Horse. Stochastic model reduction for chaos representations. Computer Methods in Applied Mechanics and Engineer- ing, 196(37-40):3951–3966, 2007. [20] J. Dutka. Richardson extrapolation and romberg integration. Historia Mathemat- ica, 11(1):3 – 21, 1984. [21] M. Eldred. Design under uncertainty employing stochastic expansion methods. International Journal for Uncertainty Quantification, 1(2):119–146, 2011. [22] Ic Enevoldsen and JD Srensen. Reliability-based optimization in structural engi- neering. Structural Safety, 15(3):169 – 196, 1994. 132 [23] P. Fuglsang and HA Madsen. Optimization method for wind turbine rotors. Jour- nal of Wind Engineering and Industrial Aerodynamics, 80(2):191–206, 1998. [24] B. Ganapathysubramanian and N. Zabaras. Sparse grid collocation schemes for stochastic natural convection problems. Journal of Computational Physics, 225(1):652–685, 2007. [25] Markus Gasserand and Gerhart Schueller. Reliability-based optimization of struc- tural systems. Mathematical Methods of Operations Research, 46(3):287–307, 1997. [26] R. Ghanem and J. Abras. A general purpose library for stochastic finite element computations. In Second MIT Conference on Computational Mechanics, Cam- bridge, MA, June 17-20 2003. [27] R. Ghanem and R. Doostan. On the construction and analysis of stochastic mod- els: Characterization and propagation of the errors associated with limited data. Journal of Computational Physics, 217(1):63–81, 2006. [28] R. Ghanem and D. Ghiocel. A new implementation of the spectral stochastic finite element method for stochastic constitutive relations. In ASCE 12th Engineering Mechanics Conference, La Jolla, CA, May 17-20 1998. [29] R. Ghanem and D. Spanos. Stochastic Finite Elements: A Spectral Approach. Dover, New York, 2002. [30] RG Ghanem and PD Spanos. Stochastic finite elements: A spectral approach. Dover publications, 2002. [31] Roger Ghanem. Ingredients for a general purpose stochastic finite elements imple- mentation. Computer Methods in Applied Mechanics and Engineering, 168:19 – 34, 1999. [32] Roger Ghanem and John Red-Horse. Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach. Phys. D, 133:137–144, September 1999. [33] D. Ghiocel and R. Ghanem. Stochastic finite element analysis of seismic soil- structure interaction. ASCE, Journal of Engineering Mechanics, 128(1):66–77, 2002. [34] William Gragg. On extrapolation algorithms for ordinary initial value problems. Journal of the Society for Industrial and Applied Mathematics, 2(3):384–403, 1965. 133 [35] AC Hansen and CP Butterfield. Aerodynamics of horizontal-axis wind turbines. Annual Review of Fluid Mechanics, 25(115), 1993. [36] Martin Hansen. Aerodynamics of wind turbines. Earthscan publications, 2008. [37] MC Hsu, Ic Akkerman, and Yc Bazilevs. High-performance computing of wind turbine aerodynamics using isogeometric analysis. Computers and Fluids, 2012. [38] Hector A. Jensen. Design and sensitivity analysis of dynamical systems subjected to stochastic loading. Computers and Structures, 83(14):1062 – 1075, 2005. [39] Monu Kalsi, Kurt Hacker, and Keper Lewis. A comprehensive robust design approach for decision trade-offs in complex systems design. Journal of Mechani- cal Design, 123(1):1–10, 2001. [40] P. Lax and A. Milgram. Parabolic equations. Ann. Math. Studies, 33:167–190, 1954. [41] O. Le Maitre, H. Najm, R. Ghanem, and O. Knio. Multi-resolution analysis of wiener-type uncertainty propagation schemes. J. Comput. Phys., 197:502–531, July 2004. [42] Olivier Le Maˆ ıtre and Omar Knio. Spectral methods for uncertainty quantification: with applications to computational fluid dynamics. Springer, 2010. [43] D. J. Lekou and T. P. Philippidis. Pre- and post-thin: a tool for the probabilistic design and analysis of composite rotor blade strength. Wind Energy, 12(7):676– 691, 2009. [44] D. Levin. Development of nonlinear transformations for improving convergence of sequences. Int. J. Comput. Math., 3:371–388, 1973. [45] John M. Mulvey, Robert J. Vanderbei, and Stavros A. Zenios. Robust optimization of large-scale systems. Operation Research, 43(2):264–281, 1995. [46] A. Nouy and C. Soize. Random fields representations for stochastic ellip- tic boundary value problems and statistical inverse problems. Archived(2013), http://arxiv.org/abs/1304.2902. [47] SS Oye. Flex 4 simulation of wind turbine dynamics. Proceedings of 28th IEA meeting of experts concerning Satate of the Art of Aeroelastic Codes for wind turbine calculations, 1996. [48] MN Papadrakakis, ND Lagaros, and VB Plevris. Design optimization of steel structures considering uncertainties. Engineering Structures, 27(9):1408 – 1418, 2005. 134 [49] Gyung-Jin Park, Kwang-Hyeon Hwang, Tae-Hee Lee, and Kwon Hee Lee. Robust design: An overview. AIAA Journal, 44(1):181–191, 2006. [50] M. Pellissetti and R. Ghanem. Iterative solution of systems of linear equations arising in the context of the stochastic fem. Journal of Advances in Engineering Software, 31:607–616, 2000. [51] L. Richardson. The approximate arithmetical solution by finite differences of phys- ical problems involving differential equations, with an application to the stresses in a masonry dam. Philosophical Transactions of the Royal Society of London., 210:307–357, 1911. [52] V A Riziotis and SG V outsinas. Gast: A general aerodynamic and structural pre- diction tool for wind turbines. Proceedings of the ECWEC, Dublin, Ireland, 1997. [53] Knut O Ronold and Carl J Christensen. Optimization of a design code for wind- turbine rotor blades in fatigue. Engineering Structures, 23(8):993 – 1004, 2001. [54] Knut O. Ronold, Jakob Wedel-Heinen, and Carl J. Christensen. Reliability-based fatigue design of wind-turbine rotor blades. Engineering Structures, 21(12):1101 – 1114, 1999. [55] G. Saad and R. Ghanem. Characterization of reservoir simulation models using a polynomial chaos-based ensemble Kalman filter. Water Resources Research, 45(Article Number: W04417), 2009. [56] JG Schepers. Verification of european wind turbine design codes, vewtdc: final report. Technical report ECNC01055, Netherlands energy research foundation ECN, Petten, 2002. [57] D. Shanks. Non-linear transformations of divergent and slowly convergent sequences. J. Math. Phys., 34:1–42, 1955. [58] HH Snel and CC Liendenburg. Aeroelastic rotor system code for horizontal axis wind turbine. Proceedings of the ECWEC, Madrid, Spain, 1990. [59] Ii Sobol. Sensitivity analysis for non-linear mathematical models. Mathematical Modeling and Computational Experiment, (1):407–414, 1993. [60] C. Soize. A computational inverse method for identification of non-gaussian ran- dom fields using the bayesian approach in very high dimension. Computer Meth- ods in Applied Mechanics and Engineering, 200(4546):3083 – 3099, 2011. [61] C. Soize and C. Descelliers. Computational aspects for constructing realizations of polynomial chaos in high dimension. SIAM Journal on Scientific Computing, 32(5):2820–2831, 2010. 135 [62] C. Soize and R. Ghanem. Physical systems with random uncertainties: Chaos rep- resentations with arbitrary probability measure. SIAM Journal of Scientific Com- puting, 26(2):395–410, 2004. [63] DA Spera. Wind turbine technology. 1994. [64] Bruno Sudret. Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering and System Safety, (93):964 – 979, 2008. [65] Henrik Stensgaard Toft and John Dalsgaard Srensen. Reliability-based design of wind turbine blades. Structural Safety, 33(6):333 – 342, 2011. [66] X. Wan and G. Karniadakis. Multi-element generalized polynomial chaos for arbi- trary probability measures. SIAM Journal on Scientific Computing, 28(3):901–928, 2006. [67] Ernst Weniger. Nonlinear sequence transformations for the acceleration of conver- gence and the summation of divergent series. Computer Physics Reports, 10:189 – 371, 1989. [68] D. Xiu and G. Karniadakis. The wiener-askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24:619–644, 2002. 136

Abstract (if available)

Abstract This dissertation is motivated by challenges involved in the design optimization under uncertainty for the structures undergoing fluid structure interaction phenomena. In particular, this study focuses on the aerodynamic simulation and optimization under uncertainty for rotor blades of Horizontal Axis Wind Turbines (HAWT). To that end, a reliability‐based design framework integrating the aeroelastic computation based on the blade element method with a stochastic optimization, is developed. The wind model is a stochastic process which accounts for the variation in time and space. The aeroelastic calculation determines the aerodynamic forces on the blade. The uncertainty propagated in the aerodynamic forces stems from the uncertain characteristics of the wind. Another source of uncertainty is introduced in the material property of the blade which results in random structural parameters such as stiffness and capacity. These sources of uncertainty lead to a random objective function and constraints. The objective function is the generated power subjected to the constraints on stresses in the blade. The reliability is defined based on the probability of the safety in the system. The design variables are considered to be shape parameters, namely twist angles along the blade. In order to show the effect of uncertainty in the design, the result of the deterministic optimization is obtained and compared with the optimization results under uncertainty. Polynomial chaos expansion is used in various parts of this work as a mathematical tool for the uncertainty propagation. In particular, polynomial chaos expansion can be used to characterize the solution of stochastic differential equations such as the governing equation of the motion of the blade. Significant computational challenges are associated with the polynomial chaos expansion approach for problems with high dimensional parameter space. This challenge is addressed in the present dissertation by developing a methodology for convergence acceleration of polynomial chaos‐based stochastic Galerkin solutions. Specifically, Nonlinear sequence transformations are adapted to these expansions, viewed as a one‐parameter family of functions with the parameter being the polynomial degree of the expansion. Stochastic Galerkin approach yields polynomial chaos representations that have the requisite analytical properties to ensure suitable convergence of these nonlinear sequence transformations. In particular, the properties of Shanks and Levin transformations are explored in the context of a stochastic initial value problem and a stochastic elliptic problem.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Comprehensive uncertainty quantification in composites manufacturing processes

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

A polynomial chaos formalism for uncertainty budget assessment

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Optimal clipped linear strategies for controllable damping

PDF

Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Stochastic and multiscale models for urban and natural ecology

PDF

Studies into computational intelligence approaches for the identification of complex nonlinear systems

PDF

On practical network optimization: convergence, finite buffers, and load balancing

PDF

Iterative path integral stochastic optimal control: theory and applications to motor control

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

PDF

The next generation of power-system operations: modeling and optimization innovations to mitigate renewable uncertainty

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

Stochastic optimization in high dimension

PDF

A novel hybrid probabilistic framework for model validation

PDF

Robustness of gradient methods for data-driven decision making

Asset Metadata

Creator Keshavarzzadeh, Vahid (author)

Core Title Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering (Structural Engineering)

Publication Date 06/23/2016

Defense Date 04/24/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag aerodynamic optimization,blade element method,convergence acceleration,nonlinear sequence transformation,OAI-PMH Harvest,polynomial chaos expansion,reliability based design optimization,stochastic differential equation

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ghanem, Roger G. (committee chair), Masri, Sami F. (committee member), Newton, Paul K. (committee member)

Creator Email vkeshava@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-422612

Unique identifier UC11285965

Identifier etd-Keshavarzz-2577.pdf (filename),usctheses-c3-422612 (legacy record id)

Legacy Identifier etd-Keshavarzz-2577.pdf

Dmrecord 422612

Document Type Dissertation

Format application/pdf (imt)

Rights Keshavarzzadeh, Vahid

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA