UFDC Home  myUFDC Home  Help 



Full Text  
INVESTIGATION OF NAVIERSTOKES CODE VERIFICATION AND DESIGN OPTIMIZATION By RAJKUMAR VAIDYANATHAN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 Copyright 2004 by Raj kumar Vaidyanathan To my parents and sisters. ACKNOWLEDGMENTS I express my sincere and heartfelt gratitude to all the people in my life who have influenced and helped in shaping it to date. I thank Dr. Wei Shyy for believing in my capabilities and his support, both academically and financially. He has been patient in guiding and encouraging me in my research when I needed it the most. I also thank Drs. Raphael Haftka and Marc Garbey, my committee members, for guiding me in my research. I would like to thank my other committee members, Drs. Renwei Mei, Norman FitzCoy and Samim Anghaie, for reviewing my work and offering their suggestions. I wish to acknowledge Mr. Kevin Tucker and NASA Marshall Space Flight Center for collaborating with me and supporting this work. I also wish to acknowledge Drs. Nilay Papila, Nestor Queipo and Siddharth Thakur for their timely help and guidance. I am thankful to my mother, Meenakshi Vaidyanathan, my father, Vaidyanathan Balasubramanian, and my three sisters, Rajeswari, Vidya and Vijaya, for being supportive throughout my academic career. I wish to thank all my "special" friends back home in India and here in the USA, who have been there for me whenever I needed them. My family and the "special" friends have been responsible for keeping my motivation up with a personal touch and also the reason for maintaining a balance between my professional and personal life. Last but not least, I appreciate all my colleagues, some of whom I have collaborated with, for providing a friendly work atmosphere and lighter moments that have helped me relax at tense times. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST O F TA BLE S ............................................ .. .............. .... .............. vii LIST O F FIG U R E S ..................... ................................... ......... .. .. ........... viii A B S T R A C T ............... ................................................................................. ......... ..... x i CHAPTER 1 IN TR OD U CTION AN D SCOPE .................................................................. ..........1 In tro d u ctio n ............................................................................... 1 S c o p e ...................................................................................... . 5 2 NAVIERSTOKES CODE VERIFICATION...............................................................9 L iteratu re R ev iew ............................................................................ .......... .. .. ...9 S c o p e ........................................................................................ 14 A approach ................ ..... ................ ...... ........... ............... 15 Richardson Extrapolation (RE) .............................................. ...............15 Least Square Extrapolation (LSE) ............ .............................................16 Error versus Residual ..................................................................... ......19 Flow Solver ........................................................ 21 A lgorithm for L SE .............................................. ............................. 23 LSE for pressure equation only.........................................................23 LSE for coupled pressurevelocity................................ ... ..................24 Results and Discussion.................................................25 TwoDimensional Turning Point Problem ............................................. 26 Lam inar LidDriven Cavity Flow ........................................ ............... 29 LSE for pressure only .............................. ........................... 33 LSE for coupled pressurevelocity computations ..............................41 3 GLOBAL SENSITIVITY AND OPTIMIZATION TECHNIQUES ...........................44 Response Surface Methodology (RSM)................................................................44 Design of Experiments (DOE) ...........................................................................49 Sensitivity A analysis .. .. .... ........................................................ ............. 5 1 O ptim ization A lgorithm ........................................ ............................................55 M ultiObjective Genetic Algorithm ...................................................................... 58 H ierarchical Clustering M ethod..................................... ............................ ........ 60 V isualization U sing B ox Plot........................................................ ............... 60 4 DESIGN OPTIMIZATION SINGLE ELEMENT ROCKET INJECTOR .................62 L iteratu re R ev iew .......................................................................... ..................... 62 S c o p e ........................................................................................6 6 Inj sector model .................................... .......................... ..... ......... 67 N um erical Procedure .................. ........................... .... .. ................. 70 D design Space ........................................................................7 1 R results and D discussion .......... ..... ...................................................... .... .... ....... 72 C FD A analysis .................................... ............... ......... 72 Grid Sensitivity Investigation...................................................................... 76 Response Surface Generation...................................................................... 78 O ptim ization P process .............. .... ...................................... .......... .. .... 8 1 Singleobjective optimization ............................................................82 M ultiobjective optim ization...................................... ............... 86 Prelim inary O observations ........................................ .......................... 93 Sensitivity and TradeOff Analysis ..................................... .....................95 Global analysis ...............................................95 Pareto front analyses ..................................................... 99 5 SUMMARY, CONCLUSION AND FUTURE WORK ......... ................................107 N avierStokes Code V erification ........................................ ....................... 107 D esign Optim ization Study ............................................................................. 108 Future w ork ............................................. ............. ............... .111 NavierStokes Code Verification........................................................112 D esign O ptim ization ......................................................... ............. ..112 L IST O F R E FE R E N C E S ................................................................. .... ...................... 114 BIOGRAPHICAL SKETCH ............................................................. ...............121 LIST OF TABLES Table page 4.1: Range of design variables (a is an acute angle in degrees and x is the thickness of O P T T in in ch es) .................................................... ................ 7 0 4.2: Fitting designs (unacceptable designs are marked in bold). .....................................73 4 .3 : T testing designs.. .........................................................................74 4.4: Independent and dependent variable (objectives) for Cases X and Y from CFD computations (nonnormalized values shown). .................................. .................74 4.5: Performance of full and reduced quadratic RSAs for the four objective functions (nonnorm alized values used) ............................................................... .... ........... 78 4.6: Performance of RSAs for the TTmax. Reduced cubic RSA has 21 coefficients (non norm alized values used). ................................................ ................................ 80 4.7: Minimizing individual objectives. Value in parenthesis (1) indicates which objective function is minimized (normalized values shown)............... ................................... 83 4.8: Study of effect of life/survivability on performance. Value in parenthesis indicates the weighting on the objective functions (normalized values shown)........................87 4.9: Study of influence of life/survivability and performance objectives on each other. Value in parenthesis indicates the weighting on the objective functions (normalized values show n). .........................................................................91 4.10: List of essential (V) and nonessential (x) variables for each objective and the mean errors between the modified RSA and original RSA. An essential variable accounts for at least 5% of the objective's variability. ........................................................96 4.11: Coefficient associated with the different terms in the linear RSA and comparison of R2 values of linear RSA with that of nonlinear RSA. ..........................................98 4.12: Objective function and design variables of nine (9) representative designs from the Pareto optim al solution set. ............................................ ............................ 102 4.13: Partial correlation coefficients (Ror) of design variables vs. objectives for different clusters along the PO F .................................................... ..... .. .. ...... ...... 106 LIST OF FIGURES Figure pge 2.1: Collocated grid and notation for 2D grid .......................... ..... .......................... 22 2.2: 2D turning point problem (dim tensions are r x ) ................................. ..................... 26 2.3: Solution (u) contours for the 2D turning point problem, grid = 81 x 81 ..................27 2.4: Application: turning point problem with e= 0.1, xaxis is the number of grid points in each direction for the fine grid on which extrapolation is done...........................28 2.5: Two dimensional liddriven cavity flow (dimensions 1 x 1). ...................................30 2.6: uvelocity contour for grid 121 x 121 and Re = 5.33. ..............................................31 2.7: uvelocity contour for grid 171 x 171 and Re = 1000. ....................... ...............31 2.8: Pressure along the lid for grid 121 x 121 and Re = 5.33 ........................................32 2.9: Pressure along the lid for grid 171 x 171 and Re = 1000........................................32 2.10: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 w ith variable lidvelocity. ........................................ .......................... 34 2.11: Least square L2 norm error and error in extrapolated pressure on grid 101 x 101 for different a for Re = 5.33 and variable lidvelocity. ............................................35 2.12: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 and variable lidvelocity ......... .............................................. ............... 36 2.13: Least square L2 norm error and error in extrapolated pressure on grid 101 x 101 for different a for Re = 5.33 and constant lidvelocity (full domain)..........................37 2.14: The shaded domain is used for LSE. A region between y = 0 0.95 for all x is u sed for ex trap olation ....................................................................... ..................37 2.15: Least square L2 norm error and error in extrapolated pressure on grid 101 x 101 for different a for Re = 5.33 and constant lidvelocity (reduced domain). ..................38 2.16: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with constant lidvelocity. ............................. ...... ........................... 39 2.17: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 with constant lidvelocity. ................................. ............... 40 2.18: Extrapolated pressure and velocity compared with CFD solution on grid 121 x 121 (M,), Re = 5.33 and variable lidvelocity. .........................................42 2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 x 171 (M,), Re = 1000 and variable lidvelocity.....................................43 3.1: Schematic of the procedure for global design optimization (Papila, 2001).............46 3.2: Full factorial 3level design (dark circles refer to points in the foreground and light circles refer to points in the background). ..........................................50 3.3: Two phases of the optimization process, where Phase 1 deals with RSA and Phase 2 deals with optimization...... ... ....... .......................56 3.4: Desirability function (d) for various weight factors, s. (Note: B < A) (Myers and M ontgomery, 1995). ............................................... ..... ...58 3.5: Schem atic of a box plot ............................................................... 61 4.1: Schematic of the G02/GH2 impinging and coaxial injector elements..................68 4.2: Schematic of Hybrid Boeing Element (U. S. Patent 6253539) ...............................68 4.3: Design variables and objectives of the single element rocket injector...................69 4.4: Simulation domain and boundary conditions ................................... .....71 4.5: Temperature field for Cases X and Y. .......................................... 74 4.6: Nearinjector temperature field for Cases X and Y..............................................75 4.7: Large recirculation zone in the combustion chamber.......................................75 4.8: Comparing the unrefined (336x81) {thick lines} and refined (430x81) {thin lines} grids. .................................................77 4.9: Comparison between the best fit possible and as predicted by quadratic response surface. Optimum refers to the case when RSA and CFD values are identical. RS CFD represents the value as for the current case (normalized values shown). ..........79 4.10: Minimization of different objectives. (Case 1: TFmax, Case 2: TW4, Case 3: TTmax, Case 4: X,,. Normalized values shown). .......................................84 urrent case (normalized values shown). ..........79 4.10: Minimization of different objectives. (Case 1: TFmax, Case 2: TW4, Case 3: TTmax, Case 4: X c,. Normalized values shown). ........................................ ............... 84 4.11: Variation of objectives with respect to AOA for fixed a, AHA and OPTT. (normalized values shown and D is the desirability function. All objectives are equally weighed in the composite function). ................................... .................86 4.12: Composite minimization of objectives with different weightings.(Case 4: (0,0,0,1), Case 5: (0,0,1,1), Case 6: (1,0,1,1), Case 7: (1,1,1,1). The values in parenthesis indicate weights for (TFmax, TW4, TTmax, Xcc). Normalized values shown) ...........88 4.13: Variation of objectives with respect to AHA or AOA with other design variables fixed. (normalized values shown and D is the desirability function. All objectives are equally weighed in the composite function)........................................... 89 4.14: Composite minimization of objectives with different weightings. (Case 8: (1,1,1,0.5), Case 9: (5,5,5,0.1), Case 10: (0.5,0.5,0.5,1), Case 11: (0.1,0.1,0.1,5). The values in parenthesis indicate weights for (TFmax, TW4, TTmax, Xcc. normalized values show n). ........................................................................92 4.15: Variation of TFmax, TW4, TTmax and Xcc, with respect to AOA for a=0, AHA=1 and OPTT=0 (normalized values shown and D is the desirability function. Case 9: Temperatures weighted over performance in the ratio of 50:1) .............. ..............93 4.16: Main factor (Si,) influence on objective variability. ....................... ...............96 4.17: Twoobjective Pareto optimal front ............................ ..... ........... .... 100 4.18: Pareto optimal solution set (+) and the multiobjective optima listed in Tables 4.7 4 .9 (o ). ..............................................................................10 1 4.19: Pareto optimal solution set and nine (9) representative solutions from the same set. The circles identify the representative of a specific cluster. ................................ 103 4.20: Box plots for the design variables in clusters 1, 3, 6 and 9. ..................................103 4.21: Box plots for the objectives in clusters 1, 3, 6 and 9. ..................... ...............105 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy INVESTIGATION OF NAVIERSTOKES CODE VERIFICATION AND DESIGN OPTIMIZATION By Raj kumar Vaidyanathan August, 2004 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering With rapid progress made in employing computational techniques for various complex NavierStokes fluid flow problems, design optimization problems traditionally based on empirical formulations and experiments are now being addressed with the aid of computational fluid dynamics (CFD). To be able to carry out an effective CFDbased optimization study, it is essential that the uncertainty and appropriate confidence limits of the CFD solutions be quantified over the chosen design space. The present dissertation investigates the issues related to code verification, surrogate modelbased optimization and sensitivity evaluation. For NavierStokes (NS) CFD code verification a least square extrapolation (LSE) method is assessed. This method projects numerically computed NS solutions from multiple, coarser base grids onto a finer grid and improves solution accuracy by minimizing the residual of the discretized NS equations over the projected grid. In this dissertation, the finite volume (FV) formulation is focused on. The interplay between the concepts and the outcome of LSE, and the effects of solution gradients and singularities, nonlinear physics, and coupling of flow variables on the effectiveness of LSE are investigated. A CFDbased design optimization of a single element liquid rocket injector is conducted with surrogate models developed using response surface methodology (RSM) based on CFD solutions. The computational model consists of the NS equations, finite rate chemistry, and the keturbulence closure. With the aid of these surrogate models, sensitivity and tradeoff analyses are carried out for the injector design whose geometry (hydrogen flow angle, hydrogen and oxygen flow areas and oxygen post tip thickness) is optimized to attain desirable goals in performance (combustion length) and life/survivability (the maximum temperatures on the oxidizer post tip and injector face and a combustion chamber wall temperature). A preliminary multiobjective optimization study is carried out using a geometric mean approach. Following this, sensitivity analyses with the aid of variancebased nonparametric approach and partial correlation coefficients are conducted using data available from surrogate models of the objectives and the multiobjective optima to identify the contribution of the design variables to the objective variability and to analyze the variability of the design variables and the objectives. In summary the present dissertation offers insight into an improved coarse to fine grid extrapolation technique for NavierStokes computations and also suggests tools for a designer to conduct design optimization study and related sensitivity analyses for a given design problem. CHAPTER 1 INTRODUCTION AND SCOPE Introduction With the advancements in computational technologies, different complex design problems can now be analyzed with reduced economical and computational cost as compared to an experiment. Additionally, computational models can be used as a substitute to study environments that are too hazardous to conduct physical experiments. These aspects have motivated the development of different tools in the areas of computational fluid dynamics (CFD) and multidisciplinary optimization (MDO) which in turn promotes CFDbased design optimization studies. Computational fluid dynamicsbased design studies can be broadly divided into two parts. The first part is the accurate modeling of the fluid flow problem and the second is the sensitivity estimation of the CFD solutions to the design variables and exploration of optimum designs. Solving fluid flow problems based on NavierStokes (NS) equations and other related models introduces challenges that include physical modeling uncertainty, geometric complexities, nonuniform and nonorthogonal meshing, disparate length scales and characteristics of the flow variables (such as velocity and pressure), and acceptable run time for engineering analysis and design. To obtain an accurate numerical solution after addressing all these issues is an elaborate topic of research. As pointed out in the review by Oberkampf and Trucano (2002), there are 5 sources of error in a CFD analysis, namely, insufficient spatial discretization convergence, insufficient temporal discretization convergence, insufficient convergence of an iterative procedure, computer roundoff and computer programming error. Much work has been done to estimate the first three errors. The remaining two are related to the field of computer programming and usually involve the verification of the code or the software. In order to ascertain the overall accuracy and identify uncertainties associated with physical modeling there is a strong need to develop techniques for accuracy assessment of a given numerical approach and mesh resolution. A priori and aposteriori error estimates are widely used to estimate errors and verify CFD solutions. A priori error estimation uses the information available before a solution is obtained based on the partial differential equation (PDE) that has to be solved and the initial and boundary conditions. In CFD analysis these estimates aid in the estimation of the stability and existence of a solution. On the other hand, an aposteriori error estimate uses the information of the computed solution. This allows the performance of a given numerical approach to be judged by estimating the errors arising out of the discretization of the PDEs. The issue of mesh resolution is of importance because a coarse grid is computationally economical but maybe inadequate in capturing the detail flow features whereas a fine grid, although it provides a wellresolved solution, increases the computational cost. Hence it is of interest to use a coarse grid and obtain an accurate solution. A well known method in this category is the Richardson extrapolation (RE). Details of this method can be found in Roache (1997) and Shyy et al. (2002). The Richardson extrapolation is based on eliminating the leading order error term based on the local truncation error analysis, implemented via solutions obtained on two or more base grid levels. An attractive feature of this method is that it can be applied to any partial differential equation (PDE) without direct information about the solution procedure or the equation itself. On the other hand, the disadvantage of this method is that it assumes the order of the solution a priori, whereas generally the order of convergence of a CFD solution is space and solution profile dependent, and may not be consistent with that indicated by the local error analysis. As pointed out by Shyy et al. (2002), the flow field with high gradient in flow properties can deteriorate the performance of this method. Typically, RE is usually conducted on two base grids with one having twice the number of nodes as the other in each direction. For practical problems with complex geometries, it is often difficult to satisfy such a requirement. A very coarse grid may not efficiently capture the flow features for extrapolation and a very fine grid might increase the computational cost, thereby spoiling the whole purpose of extrapolation. Recently, Garbey and Shyy (2003) have proposed a least square extrapolation (LSE) technique. This method addresses some of the disadvantages noticed in RE. Unlike RE, LSE directly uses the discretized PDE model under consideration to estimate the degree of convergence. In this method residual is estimated on the fine grid onto which the base grid solutions are extrapolated based on the discretized equation and the specific computational scheme. Additionally the space dependency of the solution's convergence order can be accounted for by assuming spatially dependent weight functions. Garbey and Shyy (2003) have shown that LSE is effective even when the base grids are not based on grid doubling or halving. Details of the LSE method will be discussed in chapter 2. In the area of MDO tools like design of experiment (DOE) (Myers and Montgomery, 1995) and response surface methodology (RSM) (Myers and Montgomery, 1995) have been shown to work well for a wide range of design problems. Hence it is of interest to use these tools to design devices that have complex geometries and flow features. Sensitivity analysis and optimization studies related to these devices require numerous function evaluations and the computational cost of CFD simulations restricts its use in such scenario. For such situation, RSM has been found to be very useful. A set of designs can be identified using DOE techniques and the solution for the designs obtained using CFD analyses. These solutions are then represented with response surface approximations (RSAs) (polynomials or any other analytical function), which can then be used as the function evaluator. Sensitivity analysis helps the designer by providing the knowledge of the influence of design variables on the essential flow features. Global sensitivity analysis (Sobol, 1993) is a method which has the potential for providing such information. The global sensitivity indices can be used to rank the design variables in order of their importance and also to identify the interaction between them. Based on these observations, nonessential variables can be fixed and the dimensionality of the design problem reduced. Most of the design studies require that more than one aspect of the solution (objective) be improved upon at the same time. There are different methods available for solving such multiobjective problems. Multiobjective optimization problems usually have numerous optimal solutions, known as Pareto optimal solutions (Miettinen, 1999). These solutions are such that no improvement is possible in any objective without sacrificing at least one of the remaining objectives. Hence these solutions are non dominated. On the other hand if for a given solution there are other solutions where improvement in any objective is possible without sacrificing the remaining objectives, that solution is said to be dominated. The function space of all the nondominated solutions in the Paretooptimal set is termed the Paretooptimal front (POF). This set provides the designer with a clear idea of the tradeoffs involved in a design study. Evolutionary algorithms (EAs) are popular global optimizers that have been used to find multiple Pareto optimal solutions in a single simulation run. Different multiobjective evolutionary algorithms (MOEA) have been proposed in literature (Deb, 2001). Insight into the tradeoffs between different objectives can be obtained from the POF. In this dissertation different aspects of NS code verification and CFDbased design process are addressed and methods explored with the help of individual case studies. The LSE used for NS code verification is implemented on a liddriven cavity flow with different boundary conditions and Reynolds numbers. The goal is to assess the effectiveness of the LSE technique for finite volume (FV) NS computations. The CFDMDO coupled design study is addressed with the aid of a single element liquid rocket injector. The study is conducted using response surface methodology along with global sensitivity analyses and genetic algorithms. The overall goal is to estimate the global sensitivity and tradeoffs involved in the design of a rocket injector. Scope The scope of the dissertation can be broadly divided into two, namely, * Investigation ofNS CFD code verification. * CFDbased design optimization. The NS CFD code verification deals with the demonstration of LSE for numerical accuracy improvement and computational cost reduction with the aid of few widely used benchmark problems. The implementation of the method for complex flows involved in injector design is an issue of future study. The goal of the LSE related study is to * Assess the effectiveness of the LSE technique for FV computations. * Evaluate the implication of solution gradients and singularities on the performance of LSE. * Address the issue of coupling between the pressure and velocity in the NS equations. The CFDbased design optimization study is addressed by modeling a single element liquid rocket injector. The goal is to integrate CFD with the optimization tools and obtain better design methodologies than those used in the past. While the ultimate goal is to analyze multielement injectors, much of the detailed work in injector design can be done, or at least initiated, at the single element level. Design variables governing the life/survivability and performance of the injectors are identified, namely, H2 flow angle, H2 and 02 flow areas with fixed mass flow rates of fuel and oxidizer and 02 post tip thickness. The design objectives that are life/survivability indicators are the maximum temperature on the oxidizer post tip, the maximum temperature on the injector face and the combustion chamber wall temperature taken three inches from the injector face. The performance indicator is the length of the combustion zone. To facilitate the development of the present methodology, a baseline element design is needed. This baseline concept is generated by an empirical design methodology based on a specific set of propellant flow rates, mixture ratio and chamber pressure. The selected design variables are then varied based on this baseline design and the design space populated with the aid of a DOE technique. The prescribed CFD cases are executed and post processed to extract the required dependent variable data. This data are then used to generate a RSA for each objective in terms of the independent design variables. Sensitivity and tradeoff analyses for the design are then carried out using the data obtained from surrogate models of the design objectives, and the POF (multiobjective optima) generated by Goel et al. (2004) with the aid of a multiobjective genetic algorithm (MOGA) and a local search method (e constraint strategy). The regions of the POF that represent different tradeoffs among the objectives are obtained through a hierarchical clustering algorithm which will be discussed in chapter 3. This study is broadly divided into the following two parts. 1. A sensitivity study is carried over the whole design space. The contribution of the design variables to the objectives variability is calculated using a variancebased nonparametric approach and correlations between objectives are investigated. 2. Sensitivity analyses are conducted on clusters along the POF. Box plots are used to highlight the variability of the design variables and the objectives within a cluster. Additionally, the linear relationships between the design variables and the objectives are explored with the aid of partial correlation coefficients. The remaining chapters of this dissertation address each of the issues mentioned above in detail. In chapter 2 the different aspects ofNS CFD code verification are given. In this chapter a literature survey is presented initially followed by the details of Richardson extrapolation (RE) and least square extrapolation (LSE) techniques. These techniques are implemented on a twodimensional turning point problem and a laminar liddriven cavity flow. The performances of the RE and LSE are evaluated and issues related to their implementation are identified. In chapter 3 different tools used for the design optimization study are presented in detail. Response surface methodology, DOE techniques, MOGA, optimization algorithm and sensitivity analyses tools are presented. Additionally a clustering method and a visualization tool are presented which along with other tools mentioned above are used in chapter 4 to study the design problem in detail. In chapter 4 the past and present studies related to the design of a single element rocket injector are presented. The design variables and objectives of the current injector design are identified and optimization studies carried out for different design goals using surrogate models developed based on CFD solutions. Sensitivity analyses are done with the aid of surrogate models and sensitivity indices. The design trends are observed with the aid of results obtained from these studies and visualization techniques. In chapter 5 the dissertation is summarized and conclusions drawn. The unresolved issues are identified and scopes for future studies identified. CHAPTER 2 NAVIERSTOKES CODE VERIFICATION In this chapter, the literature survey related to the work done on computational fluid dynamics (CFD) code verification is presented. Following this, details of Richardson extrapolation (RE) and least square extrapolation (LSE) techniques are presented along with case studies and results. Literature Review In the area of CFD code verification, there has been considerable research done to address the issues of grid sensitivity and quantification of uncertainties in the computations (Roache, 1997; Pelletier et al., 2001; Gu et al., 2001; Turgeon et al., 2001). The governing partial differential equations (PDEs) of fluid mechanics are either linear or nonlinear depending on the problem solved. The most frequently used computational methods for solving these PDEs are finite element (FE), finite difference (FD) and finite volume (FV) methods. To obtain an accurate solution using any of these discretization schemes, adequate grid resolution is essential since for a consistent scheme the computed solution approaches the exact continuous solution as the grid spacing tends to zero (Shyy, 1994). The error arising out of inadequate grid distribution falls under the category of domain discretization error. At the same time the discretization operators used for the derivatives in the PDEs should be according to the order of variation of the solution over the domain. If an improper variation is assumed, then it gives rise to the equation discretization error. Hence the goal is to achieve a fine grid distribution over the domain and a valid discretization of the derivatives in the PDEs without violating issues related to the stability of the computational model. Additionally, to obtain a solution for these PDEs valid boundary and initial conditions need to be specified which will define the flow field. The error estimation can be done before or after the CFD analysis is carried out. The former is called apriori error estimate and the later is referred to as aposteriori error estimate. In this dissertation the aim is to implement the aposteriori error estimates to a NavierStokes (NS) CFD code. Hence a literature survey of the available aposteriori error estimates is presented. As pointed out by Oden (2001), there are two broad classes of error estimation methods, namely, residual methods and recovery methods. In the residual method, the residual estimates the degree with which the approximate solution fails to satisfy the equations of the original problem (Babuska and Rheinbolt, 1978; Demkowicz et al., 1984). The residual estimates are used to define error norms which are largely used for controlling adaptive meshing procedures. The other type, referred to as the recovery method, attempts to enhance the computed solution during the post processing step, like Richardson extrapolation (Roache, 1997; Shyy et al., 2002) and least square extrapolation (Garbey and Shyy, 2003). There is a vast amount of literature available with respect to aposteriori error estimators for finite element approximation. A recent book by Ainsworth and Oden (2000) provides a survey of the available literature on such error estimators. Babuska and Rheinbolt (1978), Demkowicz et al. (1984), Zienkiewicz and Zhu (1987) are some of the earliest works in this area. These works concentrate on error estimation in finite element problems. Zienkiewicz and Zhu (1992a, 1992b) in a twopart paper introduced a recovery technique which was an improvement on their previous effort. They achieved higher order accuracy by using this approach and also found it to be cost effective. Their approach uses a single and continuous polynomial expansion for the nodal values of the solutions (displacements, stresses) over a patch of elements surrounding the node of interest. A least square method or an L2 projection is used to obtain this polynomial. The error in the computed solution is then measured by comparing it with the recovered solution using different forms of error norms. Ainsworth and Oden (1992, 1993a, 1993b) have presented an element residual method for finite element computations. This method has later been implemented by Jasak and Gosman (2001) for finite volume approach. Oden et al. (1994) extended the application of element residual methods to NS equations using a finite element approach. They designed error norms to measure the error in velocity and pressure for incompressible flow problems. Use of a posteriori estimator in finite volume approach has been of interest only over the past few years. Ilinca et al. (2000) have compared three different error estimators for finite volume solutions. They have compared a technique using Richardson extrapolation, a technique based on the difference of the computed solution, and a higher order reconstruction obtained using a least square method and a technique which solves an error equation. Richardson extrapolation (Roache, 1994) uses the solutions on different grids with different levels of refinement, one finer than the other. Using the Taylor series representation of the discrete solution and combining the available multiple solutions, a higher order approximation of the desired variable is obtained. The drawback is to generate an adequately resolved solution on multiple grids for complex problems and the a priori assumption of the solution accuracy order. Solution reconstruction (Barth, 1993; Olliviergooch, 1996) based on an approximation of the derivatives in the Taylor series expansion using a weighted least square method improves solution accuracy. This is then used to compare with the initial solution to obtain an error measure. The third method proposed by Zhang et al. (1997) estimates error using an error equation. In this method the error equation is derived directly from the governing equation and has a source term, which needs to be computed. This source term is evaluated using the flux differences at the cell interface. These methods were tested on subsonic, transonic and supersonic flows. Richardson extrapolation and the error equation method are shown to perform reasonably well as compared to the solution reconstruction method. Richardson extrapolation is based on eliminating the leading order error term in the assumed Taylor series expansion of the solution. This is done by combining solutions obtained on two different grids with (uniform) discrete spacings. The attractive feature about this method is that it can be implemented on a solution for any PDE without bothering about the solution procedure or the equation itself. On the other hand, the disadvantage of this method is based on the fact that it assumes a fixed order of accuracy (when using two grids) of the solution all over the computational domain. In practical fluid problems, due to numerical schemes and fluid physics involved, the order of accuracy is not uniform over the computational domain. Hence, the local truncation error based on Taylor series expansion may not represent the global accuracy of the solution. To address this, three grids of different resolution are required to estimate the order of accuracy of the solution. Additionally, problems involving turbulent flow features will result in sharp gradients in flow parameters in different regions of the domain. It has been noticed by Shyy et al. (2002) that flow field with high gradient in flow properties can deteriorate the performance of this method. Although the most common use of RE is with two grids, one twice as fine as the other, this is not always essential. But for coarse grids, the assumption of monotonic truncation error convergence may not be valid thereby requiring fine base grids. For practical problems with complex geometries, it is tough to obtain two grids that satisfy such requirements. A very coarse grid may not efficiently capture the flow features for extrapolation, and a very fine grid might increase the computational cost, thereby spoiling the whole purpose of extrapolating. Another important disadvantage of RE is that the extrapolation does not preserve the conservativeness of the flow properties. The book by Roache (1998) has more information on different approaches to address some of these issues. But in its simplest form, RE is based on an apriori assumed order of convergence of the continuous solution. In the work of Garbey and Shyy (2003), basic properties of RE and its computational implications are presented in detail. They have presented complementary views that are asymptotic expansion for continuous function in a normed vector space and numerical approximation for discrete functions defined in a mesh. Comparing these views they have pointed out that it is difficult to achieve asymptotic order of convergence unless the numerical perturbations (arising out of imperfect convergence of the fluid solver) in the computations are small. Unless a grid is very fine, it is hard to reduce the perturbation error. Based on their study on the convergence order approximation and RE in relation to CFD problems, they have found that convergence order is space and grid dependent. Hence, they have concluded that there is no way to justify the efficiency of RE which uses a constant weight formulation. Additionally, RE is easy to implement only when the grid spacing is uniform and the subsequent grids used in the extrapolation are obtained by grid doubling or halving. Scope The aim of this work is to investigate the performance of LSE method in NS computations (Vaidyanathan et al., 2004a). In particular, FV formulation is considered, which is different from the FD approach adopted by Garbey and Shyy (2003). The lid driven cavity flow with different boundary conditions and Reynolds numbers (Re) is adopted as the main test problem. The goal of this study is to address the following issues: * To assess the effectiveness of the LSE technique for FV computations. * To evaluate the implication of solution gradients and singularities on the performance of LSE. * To address the issue of coupling between the pressure and velocity in the NS equations. It is noted that the geometric complexity plays a major role in practical flow computations. This aspect can have substantial impact on the performance of any extrapolation techniques. This issue will be addressed in future studies. In the following sections, the salient features of both RE and LSE techniques will be first presented. Following this, a brief discussion of the flow solver and the algorithms for the implementation of the LSE in relation to the NS computations will be discussed. Finally, their implementation on test cases and a discussion based on the results obtained will be provided. Approach In the following, RE is reviewed first to help motivate the main topic of interest, namely the least square extrapolation. Richardson Extrapolation (RE) In essence, RE is based on an apriori assumed order of accuracy of the continuous solution. Considering a flowfield output, say a velocity component or pressure, of a CFD simulation on a grid size, h, and assuming secondorder convergence apriori, the Taylor series expansion can be written as U(x) = u(x; h)+a2h2 + ah3 +... (2.1) where U(x) is the exact solution and u(x;h) is the numerical solution based on h. Similarly the solution on a grid size, h/2 can be written as U(x)=u x; +a2,+a+... (2.2) ( 22 4 8 Solving for U(x), using Eqs. (2.1) and (2.2) and eliminating O(h2) term and neglecting the higherorder terms leads to U(x)= I4u x; hu(x;h)] +O(h3) (2.3) This is a secondorder RE which results in a gain of order one. Similarly if the solution is assumed to have firstorder convergence apriori and O(h) term is eliminated, it results in a firstorder RE with a gain of order one. U(x) = 2ux; hu(x; h)++0(h2) (2.4) When in a practical fluids problem, the degree of convergence is not uniform over the computational domain, the order of extrapolation used can either improve the order of accuracy or not affect it at all. If Eq. (2.4) was based on assumptions made in Eqs. (2.1) and (2.2), there would be no gain in the degree of convergence. Equations (2.3) and (2.4) can be generalized for two grids with spacing hi and h2, resulting in the following respective extrapolation scheme. U(x) (= h, (2.5) U(x) = (x ) (x ) (2.6) (4h2) where i7(\v.1 ) represents the interpolated values from two coarse grid solutions, which are not necessarily based on grid doubling or halving. In implementing the extrapolation techniques, it is required that the coarse grid solutions be interpolated to the locations of the fine grid on which the extrapolated solution is obtained. The order of interpolation has to be such that it does not deteriorate the order of the computed solution. Least Square Extrapolation (LSE) In this dissertation the least square extrapolation method is used which, most importantly, estimates the order of convergence as the solution to a least square problem instead of assuming it a priori. Additionally, it accounts for the dependence of the order of accuracy on spatial coordinates by using spatially dependent extrapolating weights. The details of this method are presented below. In this approach, two coarse grid solutions, ui and u2, not necessarily based on grid doubling or halving, are combined linearly using a weighting function, a, which can be spatially dependent. The extrapolated solution is given as U=aiu,+ 1a ui (2.7) where 1i and 52 represent the coarse grid solutions interpolated onto a fine grid, Mo. Following Garbey and Shyy (2003), a modified Fourier expansion has been used for the weighting function, a. The weighting function, ais given by M a(x)= o +a, cos(xxr)+ C a sin((j1) xr) (2.8) J=1 with x,j = 0,... ,Mreals. The weighting function, ais dependent on the spatial co ordinate, x for a onedimensional problem. For two dimensional problem it is a function of spatial coordinates, x and y such that a(x, y) = a(x) a(y) (2.9) This modified Fourier expansion is ideal for rectangular domains. Different function may have to be used for more complex domains. The coefficients, c, are estimated by solving a least squares problem. For a given linear PDE, say L[u]= f (2.10) with u (Ea, ) and f (Eb, b), its numerical approximation can be written as Lh [uh]= f (2.11) with uh e (E,, a) and fh e (Eh, ). Based on this, the least square problem can be formulated as Pa: Find aE= A(0,1) c L such that Lh [a1, +(la)2,]fh (2.12) is minimum in L2(Mo), where Mo represents the fine grid to which the coarse grid solutions are interpolated. For a onedimensional problem using two Fourier modes, Eq. (2.12) that has to be minimized can be rewritten as Lh [a +(1ao)2 ] +Lh[a cos (xr) +(la cos (xr))2]fh (2.13) This approach can be generalized for nonlinear PDE problems via a Newtonlike loop. For a given nonlinear problem, say, N[u= f (2.14) linearization results in J(u)[v]= g (2.15) The least square problem then becomes h k 1l +2( k+1f+ ( k+1 2 h (2.16) which is minimum in L2(Mo). The iteration is started from the initial condition of a = 0, until a k+1 a is less than some tolerance value. This method requires that the initial guess should not be too far from the final solution for convergence; i.e., the solution 92 should be close to the grid solution on Mo. In situations where J(u) is not available, a nonlinear set in ais obtained, which can be solved using nonlinear solvers. This is an issue of further study. In this initial implementation, a linearization is adapted to simplify the problem in hand. A point to be noted is that LSE uses interpolated functions on a fine grid, and the differential order, say n, of the PDE that is being solved impacts the choice of the interpolating scheme. The scheme used should give a result that should be smooth enough to be ntimes differentiable. In the present study, cubic spline is used for interpolation of the solutions. The LSE technique is based on the minimization of the residual. However, the real goal of LSE is to maximize the solution accuracy, or equivalently, minimize the solution error on the extrapolated grid. However, minimization of the residual is not necessarily the same as the minimization of the solution error. This aspect is explained in the following subsection. Error versus Residual To see the above mentioned point, let us consider only the linear problem, assuming the following form: AU=b (2.17) where A is a symmetric coefficient matrix of size n x n and U is the n x 1 exact solution vector of the linear system. Assuming u to be a nonexact solution of size n x 1, we have the following relationship: Au = b + r (2.18) where r represents the n x 1 residual vector of the system of equations. Defining E as the difference between u and U and subtracting Eq. (2.18) from Eq. (2.17) gives: As = r (2.19) Denoting the eigenvectors of matrix A by a, i = 1, ..., n such that Ia1l = 1.0 and the eigenvalues as ,, i = 1, ..., n the following relationships can be obtained: n r = aa 1=1 where ao, i = 1, ..., n are linear weights associated with each eigenvector and E = A'r which can be rewritten as (2.20) (2.21) (a AJa (2.22) The residual, F, is then defined as 7=1 F Z< (2.23) The least square approach aims at minimizing the L2 norm of the residual, F; i.e., min (F) = min = min a2 (2.24) Now the adequacy of the nonexact solution, u, can be measure by the L2 norm error measured between the u and U, which can be written as: error = E = (2.25) 7=1 Minimizing the error, gives min(error)=min ( /2)1 (2.26) From Eqs. (2.24) and (2.26) it is clear that minimization ofF is equivalent to minimization of L2 error norm (error) only when A's are equal. Therefore it is clear that properties of the matrix A define the properties of i's and the relationship between F and L2 error norm (error). The implication of the residual minimization in accuracy control will be examined along with the test problems. In the context of the current study, the solution error measure is defined as follows 2 En N =[ 1iZ ( ,,),, n ] (2.27) =1 =1I where n indicates the source of the nonexact solution, m indicates the source of the solution with respect to which the error is measured; i.e., the reference grid, (in )n represents the extrapolated or interpolated solution (nonexact solution) onto the fine grid, Mo, at the (ij)th node, (i)m represents the reference grid solution, at the (ij)th node and N is the number of nodes of the fine grid (Mo). The reference grid can be either the fine grid (Mo) solution or a fixed grid solution (Mn) fine enough to be considered as the benchmark solution. Flow Solver The governing equations are solved using the inhouse code STREAM which is based on the pressurebased algorithm and finite volume approach (Thakur et al., 2002). The equations are solved on a multiblock structured curvilinear grid. The convective terms in the momentum equations are discretized using the secondorder upwind scheme (Shyy, 1994) and the diffusion terms are discretized using a central difference scheme. The continuity and momentum equations in Cartesian coordinates are presented below. Steady state computations are carried out. Continuity: (pu)+(p = 0 (2.28) ax ax Momentum: (p ) + ( )= + i P SP (2.29) ax ay ax ax ay ay where p is the density, u and v are the velocity components in Cartesian coordinates, u is the laminar viscosity and Sp represents the pressure gradient term. In Eq. (2.29), 0 represents u or v component of the velocity. Figure 2.1 shows the collocated grid system. All the flow characteristics, namely pressure, velocity, density and viscosity are stored at the cell center (P). Carrying out the volume integral on Eq. (2.29) and using the notations shown in Figure 2.1 for the nodal references, the following discretized equation is obtained for the momentum equations. app = aEE + a(,,(1, + aNON + aSS SP + S (2.30) where the a's represent the convective and diffusive fluxes at the respective grid locations, Sp is the source term arising out of the pressure gradient term and Sc is the source term arising out of higher order derivatives in the convective fluxes. More details of this equation can be found in Shyy (1994) and Shyy et al. (1997). The residual in the context of LSE is defined by subtracting the LHS from the RHS of Eq. (2.30). The flow solver is based on the SIMPLEC (Van Doormal and Raithby, 1984) algorithm where the pressure is updated by solving a pressure correction equation formulated out of combining momentum and continuity equations. For the current study, uniformly spaced Cartesian mesh is used since the case studies involve a square domain. SNE \ N \ V NW \ (\ ^j fl \ l"  Figure 2.1: Collocated grid and notation for 2D grid. In the implementation of LSE, a distinct equation is needed for each of the flow parameters that have to be extrapolated. Hence for pressure, the pressure equation is used which is obtained by substituting the velocity components, obtained from the momentum s obtained by substituting the velocity components, obtained from the momentum equations, into the discretized (in the FV sense) continuity equation. The obtained pressure equation is given as Appp = AEE + Ap, + ANpN + Aps + S (2.31) where the S is a function of the velocity components and the A's are obtained by modifying the momentum and continuity equations (note that A's are different than the a's in the momentum equation). The details of the pressure equation can be found in Patankar (1980). In summary, there are three equations to estimate the residuals required for LSE, namely Eq. (2.30) for u and vcomponent of the velocity and Eq. (2.31) for the pressure. The residuals in the context of LSE are defined by subtracting the LHS from the RHS of Eqs. (2.30) and (2.31). Algorithm for LSE As already mentioned there are issues that need to be addressed while implementing the LSE technique, namely, the nonlinearity in the momentum equations (a's have velocity components in them which make the equation nonlinear) and the coupling of pressurevelocity so as to preserve the conservativeness of the flow properties in each cell. The pressure equation is linear in nature and therefore is simpler to work on provided that the source term is adequately represented. To proceed systematically, the NavierStokes problem is tackled in two steps. In the first step, only the pressure is extrapolated and in the second step the pressurevelocity coupled extrapolation is carried out. The algorithms are presented in the following sections. LSE for pressure equation only This step is to test the efficiency of the approach in extrapolating a given solution accurately; therefore the velocities obtained from the CFD solution on the finer of the two coarse grids used for extrapolation are used for the source term in Eq. (2.31). The extrapolated pressure, PLSE, is defined as pLSE = ap + (1 a) P2 (2.32) where p, and p2 are the interpolated pressure from the coarse grids to the fine grid using spline. The residual for the LSE method is then obtained by using the pressure equation, (Eq. (2.31)). For the preliminary analysis of the NavierStokes problem, only one mode (Eq. (2.8)) is used to obtain the pressure field. The algorithm for the LSE of pressure can be stated as * DefinepjE =pa + (1a)2 * InputpLsE, f2, v2 into Eq. (2.31). * Find that minimizes the residual of Eq. (2.31). * Calculate PLE. LSE for coupled pressurevelocity The momentum equation (Eq. (2.30)) is nonlinear in nature as the coefficients contain the velocity components as well. Hence, to linearize the system, a Picardlike iterative scheme (Ferziger and Peric, 1999) is adapted. The equation used for the velocity extrapolations comes from Eq. (2.30). n+ _=a nn+l nn+l +a n n+1 SP+Sn' (2.33) P LSE E LSE W WLSE N VLSE S SISE where Sc"+l is defined based on the new velocity components. Now the velocity components and pressure are extrapolated sequentially. Again, only one mode (Eq. (2.8)) is used for the extrapolation of all the three parameters. The algorithm of the approach is given below * START * Define pLE = a, + (la) p2 * Input pulse, f2, v2 into Eq. (2.31) (ULSE and VESE are used from the 2nd loop). * Find that minimizes the residual of Eq. (2.31). * Calculate PLE. * Define ULSE =Y7; + (1) f12. * Input pulse, ULSE, V2 into Eq. (2.33) with 0 replaced by u (VLSE is used from the 2nd loop). * Find that minimizes the residual. * Calculate ULSE. * Define VLSE =J8 V + (1f3)i. * Input pulse, ULSE, VLSE into Eq. (2.33) with 0 replaced by v. * Find / that minimizes the residual of the equation. * Calculate VLSE. * Go back to START for the next iteration. Continue the loop until a, y and / converge within the specified tolerances. These algorithms are so designed that a direct access to the coding of the CFD software is necessary. The governing equations have to be modified to accommodate the estimation of the weights. Results and Discussion The least square extrapolation method is tested by implementing it on the following cases: * Linear, twodimensional turning point problem (Figure 2.2) and * Laminar, liddriven cavity flow (Figure 2.5). 1. variable lidvelocity 2. constant lidvelocity TwoDimensional Turning Point Problem The linear equation that is solved is given below. eAu+ a(x,y) = 0,xe (0,r) (2.34) ax where e= 0.1 and a(x,y)=x +0.3 jy (2.35) u=O u =yu yy y. u=0 Sx Figure 2.2: 2D turning point problem (dimensions are r x r). The boundary conditions of the problem are defined such that xy e [0,7r]; u =y(r y) at x = 0 and u = y(ny) at x = r; u = 0 aty = 0 and 7r (Figure 2.2). There is no velocity component along the y direction. This case has been tested by Garbey and Shyy (2003) where they have used a finite difference approach. In this work a finite volume approach is used. The problem is defined so that the transition layer of thickness e (where the velocity changes direction) centered on the curve a(x,y) = 0 is not parallel to x ory axis thereby making the problem twodimensional. A sample solution is shown in Figure 2.3. In the finite volume implementation, central difference is used for the diffusion term and firstorder upwind is used for the convection term. Two modest base grids MI and M2 of sizes 23 x 23 and 29 x 29, respectively are chosen such that there is at least one or two grid points in the transition layer for the solution on grid M2. Both firstorder and second order RE are evaluated and compared with the performance of LSE. For this case study LSE is implemented with 4 Fourier modes for ax in each spatial variable (instead of one mode as in the case ofNavierStokes problem). The extrapolation is done onto fine grids, Mo, ranging from size 41 x 41 to 101 x 101 in steps of 10. In Figure 2.4, subscripts, RE1 and RE2 refers to RE of first and secondorder, respectively. The errors are represented as E which denotes the error in solution obtained using a certain method (RE or LSE) or CFD simulation on a given grid (represented as a) as compared to the solution obtained through CFD simulation on grid denoted by b (e.g. Em" is the RMS error in u estimated using all the points in the computational domain, of M, as compared to the interpolated solution from a finer grid (M,) of spacing h/2). grid = 81 x 81. Figure 2.3: Em. vs N 1 E MO23 1.5 MOREo 2 E M29 E Mn EMO ) E Mo 2 2.5  3 M0 EMo LSE 3.5 30 40 50 60 70 80 90 100 110 Grid Points (N) Figure 2.4: Application: turning point problem with e= 0.1, xaxis is the number of grid points in each direction for the fine grid on which extrapolation is done. From Figure 2.4 it is seen that RE2 gives better results than REI for grid with poor resolution. This is because the transition layer is underresolved and the dominant error comes from the second order diffusion term. As the grid gets finer, performance of REI is better than RE2 as the first order convective error dominates. For some intermediate grid (51 x 51), there maybe a cancellation between the convection and diffusion errors, which might result in large improvements for RE as seen in Figure 2.4. Richardson extrapolation fails as the grid Mo gets finer than N = 61. In all the cases, LSE predicts the fine grid solution with an error less than that of the fine grid, Mo, which is an approximation of the exact solution for N as large as 101. With the modest grid sizes of 23 x 23 and 29 x 29, there are only one or two grid points in the transition layer which has a thickness of e= 0.1. Still, the LSE gives a solution, which is more accurate than the fine grid solution on which the extrapolation is done resulting in a gain of more than one order of accuracy. This problem is relatively easy to solve and is used to confirm the potential ofLSE for FV implementation. Use of 4 modes showed better performance of LSE than using one mode. Increasing the number of modes beyond four provided marginal improvement in the extrapolated solution. Least square extrapolation is found effective for FV computations. Additionally it shows that LSE works well with flows having transition layers. For this case the resolution of the coarse grids is adequate. There are still issues like nonlinear equations, coupling of flow parameters, sharp gradients and singularities in the flow that needs to be addressed. To do so, liddriven cavity flow is used and the details are given in the following section. Laminar LidDriven Cavity Flow Garbey and Shyy (2003) have addressed the liddriven cavity flow problem in a vorticitystreamfunction formulation with a finite difference implementation. In the current study the complete NS equations are solved in 2D for the liddriven cavity flow and the pressure along with u and vcomponents of the velocity field are extrapolated. A laminar incompressible flow computation is carried out in a square cavity of dimensions 1 x 1 for Reynolds numbers 5.33 and 1000. This problem addresses the motion of the fluid in a square container induced due to the lidvelocity uhd in the positive xdirection (Figure 2.5). The remaining walls of the container have the no slip condition. Two different scenarios are considered. In the first case, the singularity at the end points of the sliding lid of the unit square cavity is avoided by choosing a variable speed as follows: Uhd = 16x2 (1 x)2 (2.36) where uhd is the ucomponent of the velocity at the lid. The vcomponent at the lid is zero and x varies from 0 to 1. In the second case, the lid velocity is kept constant, uhd = 1.0. This results in singularities in uvelocity at the lid covers where there is an abrupt change of uhd from one to zero. U'bhd v=0 u=O u=0 v=O v=0 u=O v=0 Figure 2.5: Two dimensional liddriven cavity flow (dimensions 1 x 1). The Reynolds number for the flow is defined based on the mean lidvelocity. The uvelocity contours are shown in Figures 2.6 and 2.7 for Reynolds numbers 5.33 and 1000, respectively. Some difference is noticed in the uvelocity contours at the upper corner regions of the cavity for both the Reynolds numbers. The negative variable uhd results in a uvelocity profile which looks like the mirror image of the positive constant uhd for Reynolds number 1000. The constant distribution of velocity results in the presence of singularity at the corners of the lid for the pressure as noticed from Figures 2.8B and 2.9B. The magnitude of pressure near the corners for the case with variable uhd is about 10 times more than that for the case with constant Uhd (Figures 2.8 and 2.9). The presence of singularity will be an issue while implementing the LSE technique. In the later sections, the effect will be highlighted by means of the obtained results. 0.7 0.7 A  0.4 044  0.3 03  0.2 0.2 0.1 0\ 1  0 05 1 0 05 X x A B Figure 2.6: uvelocity contour for grid 121 x 121 and Re = 5.33. A) ud = 16x2(lx)2. B) uhzd 1.0. I 0.9 0.8 08 07 0.7 06 0.6 > 0.5 5  > 0.5 . 04 0.4 /0.4 03 0.3 0.2 0.2 01 0.1  0 2 0 05 1 0 05  x X A B Figure 2.7: uvelocity contour for grid 171 x 171 and Re = 1000. A) uhd = 16x2(x)2. B) U1d = 1.0. The positive and negative values of pressure identifies whether the pressure is acting on the lid or away from it, respectively. essurealne thelid ulid 10 Pressure along the lid, uhd= 16x(lx) 0 5 0 025 05 075 1 0 025 05 075 1 X X A B Figure 2.8: Pressure along the lid for grid 121 x 121 and Re = 5.33. A) ud = 16x2 (x)2. B) uld = 1.0. Pressure along the lid, ulid = 16x2(1x)2 0 0.25 0.5 X 0 0.25 0.5 x A B Figure 2.9: Pressure along the lid for grid 171 x 171 and Re = 1000. A) uhd = 16x2(1)2. B) Uhd = 1.0. The results for the following different scenarios are presented below: * LSE for Pressure only (Re = 5.33 and 1000) 1. Uhd is variable. 2. uhzd is constant. * LSE for PressureVelocity (Re = 5.33 and 1000) 0.09 0.08 0.07 0.06 0.05 0.02 0.03 IL 0.02 0.01 0 0.01 0.02 Pressure along the lid, ulid = 1.0 0.75 1 \I I I I I 1. ui1d is variable. LSE for pressure only The results are presented in two parts. In part (a) the case with variable lidvelocity is presented and in part (b) the case with constant lidvelocity is presented. In both the parts the results for Re = 5.33 are presented first followed by the results for Re = 1000. For Re = 5.33, the two selected base grids which captures the flow features adequately are 61 x 61 and 71 x 71. The extrapolated pressure fields are compared to the solution obtained using a fixed fine grid (Mn) of size 121 x 121. The extrapolated fine grids (Mo) vary between sizes 81 x 81 and 121 x 121. For Re = 1000, the two coarse grids are 111 x 111 and 121 x 121 and the fixed fine grid (Mn) for comparison is 171 x 171. The extrapolated grids (Mo) are between 131 x 131 and 171 x 171. Only one Fourier mode (Eq. (2.8)) is used for the extrapolation of pressure. The results are compared with results obtained from Richardson extrapolation. Variable lidvelocity (Re = 5.33). In Figure 2.10A the comparison of extrapolated solutions on grid M (81 x 81, ...) with the CFD solutions on the same grid (Mo) show that LSE approach performs better as expected. The error increases as the grid gets finer since resolution of the base grid solutions is not adequate for extrapolating onto very fine grids. From Figure 2.10B, the true error can be estimated where the extrapolated solutions on any grid M, is compared with the fixed fine grid (Mn) solution which represents the "correct" solution, obtained on a substantially finer grid, for the given Reynolds number. As expected it is seen that as the grid gets finer, LSE outperforms other methods and attains an asymptotic convergence limit. The improved performance of LSE for finer grid is expected as the minimization of the residual should satisfy the governing equations more accurately. Emn vs N Emn vs N Re = 5.33 Re =5.33 2.8 2.9 2.9 EMORE1 3.1 3  E EMn"R EEM6 E RE271 E 3.73 3.2 _ _ EMO E 4E 3.9 E LSE 3.8 34 E MnLSE 4.3 4.5 4.2 80 90 100 110 120 130 80 90 100 110 120 130 Grid Points (N) Grid Points (N) A B Figure 2.10: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with variable lidvelocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (Mo). B) Extrapolated pressure compared with CFD solution on fixed fine grid (Mn), 121 x 121. (E'": n refers to the source of the solution which is compared with the reference solution identified by m). To follow up on the previous discussion, the variation of the L2 norm of the residual, F, and the L2 norm error of the extrapolated solution on grid Mo as compared to the fixed grid 121 x 121, with respect to ais examined. Figure 2.11 shows that the trend of both the accuracy measures are similar but based on minimization ofF (Figure 2.11A) LSE picks a value of 1.48 for abut based on L2 norm error (Figure 2.11B) a value of  1.6 for a is a better choice. Although it does not identify similar ds that minimize both measures, they are close to each other. More assessment will be presented for the flow with a constant lid velocity. Variable lidvelocity (Re = 1000). The performance of LSE is compared with other methods for Re = 1000 and the obtained results are shown in Figure 2.12. The base grids used for the extrapolation are 111 x 111 and 121 x 121. It clearly shows in Figure 2.12A that for a given fine grid (Mo) LSE performs better. Similar to the low Reynolds number case, the LSE solution improves as the grid gets finer in terms of the true error measure (Figure 2.12B). 9.90 9.91 9.92 9.93 9.94 9.95 a vs F (Mo = 101) UoRE2 = 2.77, loglo(F) = 9.945 oLSE = 1.48, loglo(F) = 9.976 3.20 3.40 3.60 3.60 ' 9.96 9.97 9.98 3.50 2.50 1.50 0.50 0.50 1.50 3.80 a vs Emn (Mo = 101) UoRE2 = 2.77, loglo(EM"RE2)= 3.588 LLSE = 1.48, loglo(EM"LSE)= 3.865 4.00 .I 3.50 2.50 1.50 0.50 0.50 1.50 Figure 2.11: Least square L2 norm error and error in extrapolated pressure on grid 101 x 101 for different ao for Re = 5.33 and variable lidvelocity. A) Least square L2 norm error (F). B) L2 norm error (En") (comparing extrapolated pressure with CFD solution on grid 121 x 121 (Mn)). This exercise with variable lidvelocity and two different Reynolds numbers gives us incite into few important aspect of LSE. * The resolution of the coarse grid has to be adequate for efficient extrapolation. Hence for higher Reynolds number, resolution of the base grids is higher. * LSE performs better for a considerable range of Reynolds numbers (5.33 and 1000), thereby showing its efficiency for a wide range of laminar and incompressible flow regime. * LSE performs better as the grid gets finer since it satisfies the governing equation more accurately over the computational domain. * The minimization of the residual does not guarantee the minimization of the solution error. Emn vs N Emn vs N Re = 1000 2.8 Re = 1000 2.9 EMORE1 3 EMnRE1 3.1 E 3.5 EMol21 E l EMn1 3.4 A EMn 121 3.7 121 3.9 EMLSE 3.6 MnRE2 4.1 3.8 EMnLSE 4.3 4.5 4 130 140 150 160 170 180 130 140 150 160 170 180 Grid Points (N) Grid Points (N) A B Figure 2.12: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 and variable lidvelocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (Mo). B) Extrapolated pressure compared with CFD solution on grid 171 x 171. (E,7: n refers to the source of the solution which is compared with the reference solution identified by m) Constant lidvelocity (Re = 5.33). At this point it would be interesting to take a look at a case where the lidvelocity is uniform to evaluate the efficiency of LSE when singularities are present in the solution domain. Initially the LSE is implemented on the complete solution domain. In Figure 2.13, the sensitivity of L2 norm of the residual, F and the L2 norm of the solution error, E,7, with respect to a is shown. LSE identifies a = 0.48 as the optimum value, based on the minimization of F (Figure 2.13A). But as seen from Figure 2.13B this value of a is far from the value required to obtain an accurate solution which is about 4.0. To investigate the impact of the singularities, the LSE is carried out over a reduced solution domain. Figure 2.14 shows the reduced domain over which the LSE is implemented. The extrapolation is carried over a domain with x = 01 and y = 0~0.95. a vs F (Mo = 101) URE2 = 2.77, Ioglo(F) = 6.282 ULSE = 0.48, loglo(F) = 6.637 1.80 2.00 E 2.20 0 2.40 2.60 6.5 5.5 4.5 3.5 2.5 1.5 0.5 0.5 1.5 2.5 a vs Emn (Mo = 101) URE2 = 2.77, Ioglo(EMnRE2) = 2.365 ULSE = 0.48, loglo(EMnLsE) = 2.054 6.5 5.5 4.5 3.5 2.5 1.5 0.5 0.5 1.5 2.5 Figure 2.13: Least square L2 norm error and error in extrapolated pressure on grid 101 x 101 for different (a for Re = 5.33 and constant lidvelocity (full domain). A) Least square L2 norm error (F). B) L2 norm error (E,7) (comparing extrapolated pressure with CFD solution on grid 121 x 121 (M,)) . Figure 2.14: The shaded domain is used for LSE. A region between is used for extrapolation. 0 0.95 for all x As shown before in Figure 2.13 for the whole domain, the sensitivity of least square error norm, F and the L2 norm error (E ) in comparison to the CFD solution with respect to a for the reduced domain, of grid 101 x 101, is shown in Figure 2.15. It is seen 5.60 5.80 6.00 S6.20 0 _oC 6.40 6.60 6.80 . 1 1 1 1 I I I I that LSE identifies a = 0.48 as the optimum value for a, based on the minimization ofF (Figure 2.15A). But as seen from Figure 2.15B this value of a is not as far from the optimum value of a based on the minimization of E, The suggested value of a is approximately 1.3 as compared to the previous approximate value of4.0. This suggests that a subdomain which is adequately far from the region of singularity improves the performance of LSE. A point to be noted is that during interpolation of data from one grid to the other grid, the corer singularity tends to get smeared onto the final grid. Some method of treating or accounting for this might have to be looked into. a vs F (Mo = 101) a vs Emn (Mo = 101) URE2 = 2.77, Ioglo(F) = 6.559 URE2 = 2.77, Ioglo(EMnRE2) = 2.267 . LSE = 0.48, loglo(F) = 6.911 ULSE = 0.48, loglo(EMnLsE) = 2.374 5.80 6.00 1.90 6.20 E 2.10 6.40 6.60 2.30 6.80 7.00 2.50 6.5 5.5 4.5 3.5 2.5 1.5 0.5 0.5 1.5 2.5 6.5 5.5 4.5 3.5 2.5 1.5 0.5 0.5 1.5 2.5 A B Figure 2.15: Least square L2 norm error and error in extrapolated pressure on grid 101 x 101 for different ac for Re = 5.33 and constant lidvelocity (reduced domain). A) Least square L2 norm error (F). B) L2 norm error (Em ) (comparing extrapolated pressure with CFD solution on grid 121 x 121). Figure 2.16 shows the comparisons of different extrapolation techniques. The errors shown are based on the reduced domain. There are oscillations noticed in the convergence rates as the grid gets finer. This is due to the fact that the region for different grids need not always lie between y = 0 0.95, since the nearest grid point may be below or above. This results in different levels of improvement. The present domain reduction strategy suggests that when a domain farther away from the singularities is considered, the performance of LSE seems better. Of course, this is a very simplistic solution and may not necessarily be able to handle other flow problems with more complexities. In Figure 2.16A, the extrapolated solutions are compared to the fine grid (Mo) solutions, which show that LSE outperforms the others. In Figure 2.16B, the comparison of the extrapolated solution is made with a fixed fine grid solution of size 121 x 121. The solution on this grid is considered as the benchmark solution. Clearly LSE performs well as compared to RE. Emn vs N Emn vs N Re = 5.33 Re = 5.33 1.6 1.6  1.7 1.8 EMORE1 EMn 1.8 EMo61 2 9 E RE1 E EM61 U2.2 EMRE2 2 EM071 2.1 EMn61 S2.3 E) E p Mn 2.6 2.4 MnLSE 2.8 2.5 80 90 100 110 120 130 80 90 100 110 120 130 Grid Points (N) Grid Points (N) A B Figure 2.16: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with constant lidvelocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (M,). B) Extrapolated pressure compared with CFD solution on grid 121 x 121 (Mn). (EZ: n refers to the source of the solution which is compared with the reference solution identified by m) Constant lidvelocity (Re = 1000). LSE is implemented in the same reduced domain as shown in Figure 2.14. The coarse grids used are 111 x 111 and 121 x 121 and the fixed fine grid (Mn), which is the representative of the correct solution, is 171 x 171. Figure 2.17, shows the comparison of different extrapolation schemes. Similar trends are noticed for the reduced domain as before. Emn vs N Emn vs N Re = 1000 Re = 1000 2.8 2.8 3 m_ ) __ EMORE1 3 3.2 3.2 X A X X )K EMnRE1 3.23.2 E 3.4 EMORE2 EM 111 EM0121 I3.4 O  EMn111 S3.6 3 Mn 3. EM E 3. 6 LSE 41 3.8 4.2 4 130 140 150 160 170 180 130 140 150 160 170 180 Grid Points (N) Grid Points (N) A B Figure 2.17: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 with constant lidvelocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (Mo). B) Extrapolated pressure compared with CFD solution on grid 171 x 171. (E,": n refers to the source of the solution which is compared with the reference solution identified by m) This exercise helps gain an idea of shortcomings related to LSE technique. It shows that when singularities are present, the minimization of the L2 norm of the discretized governing equation residual, F, need not minimize the solution accuracy measure, E" It depends largely on the matrix of equations used for the least square implementation. When the singularity is removed and the LSE is carried out on the reduced domain, the minimization of F tends towards the minimization of E" and the performance of LSE is improved. LSE for coupled pressurevelocity computations The coupled algorithm is implement on a liddriven cavity flow with Re = 5.33 and 1000 and the variable lidvelocity as given in Eq. (2.36) is used. The geometry is same as before (Figure 2.5). The coarse grids (MI, M2), fine grids (Mo) and the fixed grids (Mn) are all the same as those used for pressure extrapolation. It is noticed that about 34 iterations of this algorithm results in converged values of a yand f with about 34 inner iterations for the velocity components. Only the error compared to the finest grid (Mn) is shown since the goal is to obtain a solution as accurate as the exact solution. For clarity, only the RE2 and LSE solutions are shown in Figures 2.18 and 2.19 since they performed better compared to the rest. Re = 5.33. The En measure for pressure, uvelocity and vvelocity are shown in Figures 2.18A, 2.18B and 2.18C, respectively. Single mode (constant, ro) is used for the extrapolation of all the flow parameters. The LSE performs well in extrapolating all the 3 flow parameters onto the finest grid. Re = 1000. In Figure 2.19, the LSE solution is compared with the solution of second order RE. There is considerable difference in the trends as compared to the low Reynolds number case. But as the grid gets finer, the performance of LSE is better than RE, suggesting that LSE gives a solution closer to the exact solution than RE. This exercise gives the following information: * The coupled algorithm is effective in treating all flow variables: pressure, u and v velocity, adequately, and works well for a wide range of Reynolds number. * It is surprising to notice that for the coupled algorithm, the improvement noticed in LSE as compared to RE2 seems to be more for the higher Reynolds number case than the lower. Additionally the trend of uvelocity is reversed for Re = 1000 and as the grid gets finer, the performance of LSE is poorer in terms of accuracy. This requires further study. Emn vs N Re = 5.33 pressure 80 90 100 110 Grid Points (N) 120 130 Emn vs N Re = 5.33 uvelocity 80 90 100 110 Grid Points (N) Emn norm error vs N Re = 5.33 vvelocity 4.0 4.1 4.2 X EMnRE2 ) o 4.3 EMnLSE 4.4 4.5 80 90 100 110 120 130 Grid Points (N) C Figure 2.18: Extrapolated pressure and velocity compared with CFD solution on grid 121 x 121 (M,), Re = 5.33 and variable lidvelocity. A) Pressure. B) uvelocity. C) vvelocity. 3.4 3.5 3.6 0 3.7 3.8 120 130 Emn vs N Re = 1000 pressure 3.4 3.5 3.6 o 6 3.7 3.8 3.9 2.7 2.8 2.9 E 0g3.0 o 3.1 3.2 3.3 Emn vs N Re = 1000 uvelocity 130 140 150 160 170 180 130 140 150 160 Grid Points (N) Grid Points (N) 170 180 Emn vs N Re = 1000 vvelocity 130 140 150 160 170 180 Grid Points (N) C Figure 2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 x 171 (M,), Re = 1000 and variable lidvelocity. A) Pressure. B) uvelocity. C) vvelocity. CHAPTER 3 GLOBAL SENSITIVITY AND OPTIMIZATION TECHNIQUES In this chapter detail of the surrogate modeling method namely response surface methodology (RSM), global sensitivity approach and multidisciplinary optimization (MDO) tools like design of experiments (DOE) and genetic algorithm are presented. These tools are applied in the design of a single element injector of a liquid rocket engine and the obtained results are presented in chapter 4. Response Surface Methodology (RSM) The approach of RSM is to perform a series of experiments, based on numerical analyses or physical experiments, for a prescribed set of design points (using DOE), and to construct a global approximation of the measured quantity (objective) over the design space. Figure 3.1 shows the three steps involved for a singleobjective optimization study which is dependent on two design variables. This is a general case and can be extended for multiobjective optimization study as explained later in this chapter. Step I (Figure 3.1) identifies the designs and the corresponding values of the objective for each design over the selected design space. In step II (Figure 3.1) a response surface approximation (RSA) of the objective is generated, which is a function of the two design variables. Finally step III (Figure 3.1) involves the search for the optimum design using the generated RSA for function evaluations. In this dissertation regression analysis is used to construct polynomialbased RSAs of assumed order and unknown coefficients. The number of coefficients to be evaluated depends on the order of the polynomial and the number of design parameters involved. For instance, a secondorder polynomial ofN design variables has (N+1)(N+2)/(2!) coefficients. A cubic model has (N+1)(N+2)(N+3)/(3!) coefficients. In this dissertation, RSAs are constructed using JMP (2002), a statistical analysis software that provides a variety of statistical features in an interactive format. In the practical application of RSM, it is necessary to develop an approximate model (RSA) for the true response or the objective. The second order (quadratic) RSA is used most frequently as it is the most economical one. Such an approximation for a response variable (objective) y with k regressors can be written as k k k1 k Y = A + A,x, + ,x2+ x, x + c (3.1) I=1 I=1 i=1 j=2 The above equation can be written in matrix notation as follows y = Xf + e (3.2) where y is the (nxl) vector of observations, X: (nxn,) matrix of the levels of the independent variables, fl: (npxl) vector of the regression coefficients, e: (nxl) vector of random error, n: the number of observations, and nr: the number of terms in the RSA. The purpose is to find the vector of least square estimators b that minimizes Error = _,2 = E e = (y Xf)T (y Xf) (3.3) I=l which yields the least squares estimator of f b=(XTX) XTy (3.4) Objective Function S Deslgn Vanable 2 step I: Define and populate the space step II: Generate RSA. step III: Search for the optimum Figure 3.1: Schematic of the procedure for global design optimization (Papila, 2001). The quality of different RSAs can be evaluated by comparing the adjusted rms error oa (Myers and Montgomogery, 1995) value that is defined as: S= e2 (3.5) nn where e, is the error at ith point of the data used for fitting. The accuracy of the RSA in representing the design space is measured by comparing their predictions to the actual values of the objective at test design points, different from those used to generate the fit. The prediction rmserror ofor the test set is given by: S=J (3.6) In this equation c, is the error at the ith test point and m is the number of test points. The test points can be selected using the crossvalidation technique (Papila, 2001), which is an established method for estimating the prediction accuracy. It is usually performed using either a number of random test/train partitions of the data, or a kfold cross validation (Mullin and Sukthankar, 2000). In kfold crossvalidation, the data points are divided in k equally sized subsets. One of the subset is kept for testing and the remaining k1 are used to fit the RSA. This is done for all the k subsets and the average error of all the k RSAs is computed. The fitting k1 sets which has the minimum prediction rms error, o, is selected and the remaining subset used as the testing set. In certain cases, especially when higher order polynomials are used, the data available for RSA may not be enough to spare some for testing the fitted model. Hence, an alternate method to estimate the performance of the RSA is to compute the PRESS rmserror. This method was proposed by Allen (1971, 1974) and it computes a sum of squares of the residuals. The residual is obtained by fitting a RSA over the design space after dropping one design point from the fitting set and then comparing the predicted value of the RSA for that point with the expected value. This is done for every point in the design space. The PRESS rmserror is given by PRESS,, = =l (3.7) 7) where y, is the expected value, ), is the value predicted by the RSA for the ith point, which is excluded while generating the RSA and n is the number of design points. If this value is close to the adjusted rmseror, oa, then the model performs well. A measure of the variability in an objective accounted by its RSA is given by the coefficient of multiple determination (Myers and Montgomogery, 1995), R2, given as R = 1 ES 33.8) SSyy n where SSE is the sum of squares of the residuals (= y (y ),2 ) where j is the predicted 7=1 value by the RSA. SSyy is the total sum of squares about the mean given by SS = SSE +SSR = (y)2 (3.9) i=1 where y is an overall average ofyz. SSR is the sum of squares due to regression n (= Y (, y)2 ) where y is the overall average ofy,. Since R2 increases as more terms are 7=1 added in the RSA, the adjusted coefficient of multiple determination R2, a normalized measure is usually preferred. It accounts for the degrees of freedom in the model and is given by S SS, (nn ) n (3.1 R=I1 = (IR (3.10) SS, /(n 1) nnP For a good fit, the value of R2 should be closer to one. The polynomialbased RSA is effective in representing the global characteristics of the design space. It can filter the noise associated with design data. The linearity of the polynomialbased RSA allows us to use statistical techniques known as design of experiments (DOE) to find efficient fitting sets. On the other hand, depending on the order of polynomial employed and the shape of the actual response surface, the RSA can introduce a substantial error in certain region of the design space. Design of Experiments (DOE) It is well established that the predictive capability of RSM is greatly influenced by the distribution of sampling points in design space (Unal et al., 1997, 1998). Design of experiment is required to select designs for RSA that minimizes the effect of noise. There are different types of DOE techniques in the literature as reported by Haftka et al. (1998). One of the most conservative DOE techniques available in literature is the full factorial design (JMP, 2002). Unal et al. (1997, 1998) studied response surface modeling using orthogonal arrays (OA) in computer experiments for reusable launch vehicles and illustrated that using this technique can minimize design, development, test and evaluation cost. Unal and Dean (1995) studied a robust design method based on the Taguchi method (Unal and Dean, 1991; Dean, 1996) to determine the optimum configuration of design parameters for performance, quality and cost. They demonstrated that using such a robust design method for selection of design points is a systematic and efficient approach for determining the optimum configuration. The full factorial design and the OA are explained below. In full factorial design the range of each design variable is divided into one or more intervals, which mark the number of levels. These intervals are usually, evenly spaced. All the possible combinations of the levels of all the design variables give the design points in the design space. For example, for three variables with three levels each there are totally 27 design points (Figure 3.2). Figure 3.2: Full factorial 3level design (dark circles refer to points in the foreground and light circles refer to points in the background). An OA is a fractional factorial matrix that assures a balanced comparison of levels of any factor or interaction of factors. Consider A, a matrix with elements of aj where denotes the row ( = 1,2... nr) and i denotes the column (i = 1,2.. .n) that a belongs to, supposing that each a e Q = {0,1...q1 }. If the columns representing each design variable are mutually orthogonal, matrix A is called an orthogonal array. It is of strength t nc if in each nrrowbytcolumn submatrix, the qt possible distinct rows occur 2 (index of the array) times. Such an array is denoted by OA(nr,nc,q,t) by Owen (1992). Since the points are not necessarily at the vertices, the OA can be more robust than the full factorial design in interior design space and is less likely to fail the analysis tool. Based on the DOE theory, OA can significantly reduce the number of experimental configurations (Papila, 2001). In this dissertation OA is used to obtain the design for the optimization study. Sensitivity Analysis Unlike local sensitivity where the partial derivatives are used to locally estimate the sensitivity of an objective to a specific design variable, global sensitivity allows the study of overall model behavior. To understand the concept, assume a surrogate model (in this dissertation a RSA) of an objective, f(X) as a function of the design variables, X, where X is the vector of the design variables x1, x2, ..., x,, scaled between zero and one. This surrogate model (in this dissertation, polynomials of nth order) represented in Eq. (3.11) has to be square integrable. f(X)= fo + If (x, )+ f,,(x,,x,)+...+ ,(x,,x,,...,x,) (3.11) An analysis of variance (ANOVA) (Archer et al., 1997) representation of Eq. (3.11) can be obtained by imposing the following condition f1 ,dk = 0 (3.12) 0 for k = i, ..., i,. This condition is essential for the uniqueness of Eq. (3.11). This results in a set of summands of different dimensions. Each of the components f/ is responsible for the interactive contribution of design variables x, ,..., x, to the variability of the objectivef(X) over the ndimensional unit hypercube design space. Integrating Eq. (3.11) over all the design variables gives If(x)ndxk fo (3.13) Integrating Eq. (3.11) over all the design variables except x, gives f (x)ldrk =0 o+fx) (3.14) from which(x) is obtained. from which f (x,) is obtained. Integrating Eq. (3.11) over all the design variables except x, and x, gives Jf(x)dxk f0 f, )+ f, (X )+ ) (3.15) k:z,j from which, (x,,x,) is obtained. This can similarly be extended to obtain the remaining summands off(X). The imposed condition ensures that the summands are orthogonal in nature and sincef(X) is square integrable the summands are too. Therefore the partial variances can be calculated as D, = f,2dx..., (3.16) and the total variances is D= ff2(x)dx f2 (3.17) Therefore it can be shown (Sobol, 1993) that D=3 ZDY (3.18) S=1 1,< Partial variance gives a measure of the spread of the function due to one of the independent variables. Total variance gives a measure of the spread of the dependent data due to all the independent variables. Sensitivity is a measure of the contribution of an independent variable to the total variance of the dependent data. Sobol (1993) had proposed a variancebased nonparametric approach to estimate the global sensitivity using Monte Carlo methods. To calculate the total sensitivity of any design variable, say xj, the design variable vector X is divided into two complementary subsets, x, and Z where Z is a vector containing x2, x3, ..., x,. The purpose of using these subsets is to isolate the influence of x, onf(X) variability from the influence of the remaining design variables included in Z. The total sensitivity index for x, is then defined as S toa = D toal ID (3.19) where Dotal = D + D (3.20) DxI is the partial variance associated with x, and Dx^,z is the sum of all partial variances associated with any combination of the remaining variables representing the interactions between x, and Z. Similarly the partial variance associated with Z can be defined as Dz. Therefore the total objective variability can be written as D= Dx, +D, +D xz (3.21) In cases where theJ(X) is an analytical function, the multidimensional integrals for computing the partial variances can be evaluated analytically. However, the Monte Carlo approach is applicable under general conditions (e.g., any model, design under uncertainty) and has been adopted in this dissertation. The designs for the Monte Carlo approach are selected without any preference of one design over the other from the uniformly distributed design space of unit sides. Hence the design variables are uncorrelated which is essential for the effective implementation of the method proposed by Sobol (1993). The variance estimates can be obtained using the following expressions: IN If (x,,Z) f (3.22) wherefo is the mean of the objective. N f(x,Z)f/(x,,Z') D, +f02 (3.23) f2(1, Z) D + f02 (3.24) N1= J N f f(x,,Z)f(x', Z) Dz + f02 (3.25) N (3.25) where the terms (xi, Z) and (x ', Z') represent random sample designs. Equations. (3.23), (3.24) and (3.25) give the estimates for Dx D and Dz, respectively and Dx,, can be calculated using Eq. (3.21). Once these estimates are known the total sensitivity index of the objective variability with respect to a given design variable can be obtained using Eq. (3.19). The influence of a design variable to an objective variability without accounting for any of its interactions with other variables is denoted as a main factor index and given as Sx = D, /D (3.26) Each pair of random samples requires three different objective function evaluations (e.g.,J(xl, Z),f(xi, Z') andf(x ', Z)). The mean and the total variance of an objective need to be estimated only once, and only two Monte Carlo integrals per design variable are necessary to compute the main factor and total sensitivity indices. Hence the increase in computational cost is linear with the increase in the number of design variables. These indices can be used to understand the influence of design variables on the variability of any given objective over the chosen design space. This method proposed by Sobol (1993) is effective when the design variables are uncorrelated. In a multiobjective optimization study multiple optima are obtained which collectively give the Pareto optimal front (POF). The design variables at different regions along the POF share some similar features. Therefore, the design variables are correlated and this correlation has to be accounted for before estimating the influence of a design variable on an objective. The total sensitivity and main factor indices, which require uncorrelated data, cannot be used for such a situation. Hence, for this purpose, a partial correlation coefficient (JMP, 2002) was calculated. Estimation of partial correlation coefficient (JMP, 2002) is a variancebased parametric approach that gives a measure of the linear relationship between the variances of a design variable, say xi and an objectiveJ(X), after the influence of other variables have been filtered out. It gives a measure of expected change inf(X) per unit change in xi which in other words gives the sensitivity of the function to the design variable. Linear approximations are obtained for xi andJ(X) as a function of the components of Z and the residuals (say, ri and r2, respectively) measured as the difference between the data used for the approximations and the predicted value of the approximations are estimated. A partial correlation coefficient is the correlation between these two residuals and is given by Rcor E (r )(r2 2)] (3.27) Optimization Algorithm The optimization process can be divided into two parts: 1. RSA phase, 2. Optimizer phase. In the first phase, RSA are generated with the available data set. In the second phase the optimizer uses the RSA during the gradientbased search for the optimum (the maximum or minimum of the objective). The initial design for the search is randomly selected from within the design space. The flowchart of the process is shown in Figure 3.3. The optimization problem at hand can be formulated as min(fX)}subject to lb X K ub, where lb is the lower boundary vector and ub is the upper boundary vector of the design variables vector X. If the goal is to maximize the objective function thenf(X) can be written as g(X), where g(X) is the objective function. Additional linear or nonlinear constraints can be incorporated if required. Solver (2002), an optimization tool available as part of Microsoft Excel package is an ideal tool for such studies. This tool uses the Generalized Reduced Gradient (GRG2) nonlinear optimization code developed by Lasdon et al. (1978). __________________ Initial Design Variables Second Phase First Phase Data set Perturb Design variables Generating RS RS generated Evaluate No Converged? Yes gradient /I Optimum < Soln. Figure 3.3: Two phases of the optimization process, where Phase 1 deals with RSA and Phase 2 deals with optimization. In most optimization studies it is desirable to simultaneously optimize more than one response/objective. One method is to build, from the individual objectives, a composite objective known as the desirability function. The method allows for the designers' priorities to be addressed during the optimization procedure. The first step is to define the desirability function, d for each objective. Desirability function is a normalizing of the objective such that its' maximum and minimum over the design space lie between zero and one. In the case where an objective should be maximized, say R1, the desirability takes the form: cd, = (3.28) BA where B is the target value and A is the lowest acceptable value such that d, = 1 for any R1 > B and d = 0 for R1 < A. The weighting factor t is set according to the designer's goal and the compromise he wants to achieve between different objectives. In the case where a response is to be minimized, say R2, the desirability takes on the form: d2 / E' (3.29) SCE where C is the target value and E is the highest acceptable value such that d2 = 1 for any R2 < C and d2 = 0 for R2 > E. Choices for A, B, C, and E are made according to the designer's priorities or simply as the maximum and minimum values of the objective over the design space. Values of s and t are set based on the priority of the objective. The sensitivity of the parameter, s, illustrated in Figure 3.4 can be instructive (Myers and Montgomery, 1995). Desirability functions with s << 1 imply that a product need not be close to the response target value, C, to be quite acceptable. However, s = 8, implies that the product is nearly unacceptable unless the response is close to C. A single composite desirability D is defined based on the geometric mean of the desirabilities of the individual objectives. D=(d d2 d3....dm,)' (3.30) This is then submitted to an optimization toolbox to be maximized. 58 A B 1.0 1.0 0.8/ 0.8 0.6 S /04 0.6 I 0.4 0.4 0.2 0.2 s8 0.0 0.0 A B Response Value Figure 3.4: Desirability function (d) for various weight factors, s. (Note: B < A) (Myers and Montgomery, 1995). MultiObjective Genetic Algorithm One of the other optimization approaches is to use Multiobjective genetic algorithm (MOGA). It works on the principle of survival of the fittest. Genetic operators like reproduction, crossover and mutation are employed to find the optimal solution. In a multiobjective optimization scenario, when the objectives are conflicting in nature, many representative optimal solutions can be obtained. All these solutions comprise a set called the Pareto optimal set (Miettinen, 1999). The solutions in this set are considered non dominated as they are not inferior to any other solution in the entire design space without having a bias towards some of the objectives. The function space of all the solutions in the Pareto optimal set is termed the Pareto optimal front (POF) (Miettinen, 1999). When there are two objectives, the POF is a curve. When there are three objectives, the POF is represented by a surface. If there are more than three objectives, it is represented by hypersurfaces. One of the MOGA that has been shown to be effective in finding the Pareto optimal solutions is the elitist nondominated sorting genetic algorithm (NSGAII) (Deb et al., 2000). The algorithm can be described as: 1. Randomly initialize population (designs in the design space) of size N. 2. Compute objectives for each design. 3. Rank the population using nondomination criteria. 4. Compute crowding distance (this distance finds the relative closeness of the solution with other solutions in the objective space.) 5. Employ genetic operators selection, crossover & mutation create new population. 6. Evaluate objectives for the new population. 7. Combine the two populations, rank them and compute the crowding distance. 8. Select N best individuals. 9. Go to step 3 and repeat till termination criteria is reached, which in the current study is chosen to be the number of generations. It is shown in a number of studies that using a combination of MOGA and local search (also known as hybrid GA), helps achieve faster convergence to the global Pareto optimal solution set (Deb and Goel, 2000, 2001; Goel, 2001; Goel and Deb, 2002; Ishibuchi and Yoshida, 2002). The a posteriori hybrid method (Deb and Goel, 2001) assumes the set of solutions generated by MOGA simulations as a starting point for the local search. Most of the local search methods are very efficient only for single objective optimization problems. Hence, a e constraint strategy (Deb and Goel, 2001) helps reduce the dimensionality of the problem. Sequential quadratic programming, available in Matlab (2002) can be used as the local search optimizer. This gives a set of solutions from which the dominated and duplicate solutions are removed to obtain the global Pareto optimal solution set and the POF. Hierarchical Clustering Method The POF can be divided into a number of clusters using a hierarchical clustering algorithm (Jain and Dubes, 1988) to assist the designer in selecting the optimal solution of choice. The clustering algorithm can be summarized as: 1. Start by assuming all the solutions as individual clusters. 2. Find the mean of each cluster. 3. Find the distance between clusters. 4. Merge closest clusters. 5. Go to step 2 till the number of clusters is equal to some prefixed value. 6. Find the member of each cluster closest to the mean of the cluster. This is the representative element. Visualization Using Box Plot Box plot (JMP, 2002) is a visualization tool, which can help understand the variability in the design variables and objectives within the clusters of the POF. It can also assist in identifying the design variable that should be fixed when analyzing a given cluster. A schematic of a box plot with the yaxis representing the range of the scaled design variable is shown in Figure 3.5. Figure 3.5 shows that 25%, 50% and 75% of the data lie below the lower hinge, median and upper hinge, respectively. The difference between the upper and lower hinges is known as the "Hspread". The inner fence is located 1 step beyond the hinges, which is equal to 1.5 times the Hspread. The upper adjacent value is identified as the largest value below the upper inner fence. The lower inner fence and lower adjacent value can be 61 similarly determined. These box plots will be used to visualize the variability along the POF. The spread between the upper and lower adjacent values gives the range of variation of a variable in a given set. Any variable lying beyond this spread is a potential outlier of the set. Small range gives a tight bound on that particular variable. Sometimes the box plot is known to collapse to a single point, which suggests that the variable should be fixed at that value. 1 Outside value S Inner fence Upper Adjacent Value Upper hinge S Median Lower hinge Lower Adjacent Value 0 Figure 3.5: Schematic of a box plot CHAPTER 4 DESIGN OPTIMIZATION SINGLE ELEMENT ROCKET INJECTOR In this chapter the sensitivity and optimization studies of a single element liquid rocket injector design are presented. Initially a literature survey of the past injector design studies is presented. Following this the injector model is described. Then the results related to the sensitivity analyses and optimization studies are presented. Literature Review A critical goal for space propulsion design is to make the device safer, more affordable and more reliable. The design of combustion devices, namely, injectors, chambers and nozzles, will help in meeting these goals. The characteristics of the injector design are a key factor for both performance and thrust chamber environments. The thrust chamber performance is estimated by the rate and the extent to which mixing and resultant combustion occurs. The location of the mixing and resultant combustion determines the injector and thrust chamber thermal environments. These environments include temperatures on the combustor wall, the injector face and, for coaxial injectors, the oxidizer post tip. The difficulty encountered in designing injectors that perform well and have manageable environments is that the factors that promote performance often lead to increased heating of the solid surfaces of the injector and combustor thereby reducing the life expectancy or survivability of these components. Current injector design tools have been in use for 30 years or more and are largely empiricalbased (Rupe, 1965; Dickerson et al., 1968; Pavli, 1979; Nurick, 1990; Pieper, 1991). Extensive sub and fullscale hotfire test programs often guided these injector design methodologies. The experimental databases, and thus the tools developed from them, are limited, in terms of design space, to specific element configurations that have been tested (Calhoon et al., 1973). In terms of scope, the design tools typically focus on performance, with the environment being a secondary consideration. The limited amount of environmental information available from these tools is usually onedimensional and not functionally related to details of the injector design. It is very doubtful that application of these traditional design tools will result in robust future propulsion devices. Over time space propulsion programs with compressed schedules, lower budgets and more stringent requirements have resulted in the development of broader and more efficient injector design methodologies. The initial work by Tucker et al. (1998) and Vaidyanathan et al. (2000) focused on the optimization of a shear coaxial injector element with gaseous oxygen (GO2) and gaseous hydrogen (GH2) propellants. In this study, the design data was generated using an empirical design methodology developed by Calhoon et al. (1973). Calhoon et al. (1973) conducted a large number of coldflow and hotfire tests over a range of propellant mixture ratios, propellant velocity ratios and chamber pressure for shear coaxial, swirl coaxial, impinging, and premixed elements. The data were correlated directly with injector/chamber design parameters, which were recognized from both theoretical and empirical standpoints as the controlling variables. For the shear coaxial element, performance, as measured by energy release efficiency, ERE, was obtained using correlations taking into account combustor length, Lcomb (length from injector to throat) and the propellant velocity ratio, Vf/Vo. The nominal chamber wall heat flux at a point just downstream of the injector, Qnom, was calculated using a modified Bartz equation and was correlated with propellant mixture ratio, O/F, and propellant velocity ratio, VfVo to yield the actual chamber wall heat flux, Q. The objective was to maximize injector performance while minimizing chamber wall heat flux (lower heat fluxes reduce cooling requirements and increase chamber life/survivability) and chamber length (shorter chambers lower engine weight). In Tucker et al. (2000) the designs of an impinging injector element and a swirl co axial injector element have been carried out. For the impinging injector element, the empirical design methodology of Calhoon et al. (1973) used the oxidizer pressure drop, APo, fuel pressure drop, APf, combustor length, Lcomb, and the impingement halfangle, a as design variables. Objectives included ERE (a measure of element performance), wall heat flux, Qw, injector heat flux, Q,y, relative combustor weight, Wrei, and relative injector cost, Crei. The gaseous propellants were injected at a temperature of 540R. The empirical design methodology used to characterize the ERE and Qw used a quantity called the normalized injection momentum ratio to correlate the mixing at the different design points for the triplet element. The swirl coaxial element has been used sparingly in USA, but has been widely used in Russia because of its reported ability to perform well over a large throttle range (Gill and Nurick, 1976). This has sparked the interest in exploring its possibilities. The empirical design methodology of Calhoon et al. (1973) used the oxidizer pressure drop, APo, fuel pressure drop, APf, combustor length, Lcomb, and the full cone swirl angle, 0, as design variables. The objectives modeled were ERE (a measure of element performance), wall heat flux, Qw, injector heat flux, Q,,, relative combustor weight, Wrei, and relative injector cost, Crei. The advent and advancement of CFD in the last 20 years have produced a capability that has shown potential as an improved design tool in many areas of rocket propulsion. The application of CFD to injector design has lagged behind other areas such as turbomachinery because the physical models are more complicated for multiphase, turbulent reacting flows. New models that efficiently account for some of the complex processes (Cheng and Farmer, 2002) and thus increase the solution fidelity have recently become available. However, the threedimensional geometry of multielement injectors and the complex physical processes inherent in the flows that issue from them create major obstacles in validating the solution and generating significant amount of parametric solutions. The harsh high pressure and temperature environments typical of injector flows create significant difficulty in obtaining experimental data of satisfactory quality to validate and guide further development of computational models. Finally, solving the equations, for multiphase reacting flows, with high resolution typically requires lengthy computational times. However, continuing increases in computer speed and progress in parallel processing have begun to mitigate this turnaround problem. It has long been known that small changes in injector geometry can have significant impact on performance, as well as on environmental variables such as combustion chamber wall and injector face temperatures and heat fluxes (Gill and Nurick, 1976). Several geometric design variables can be accounted over appropriate ranges by calculating the performance and environmental variables of the injector using CFD. Global approximation method like RSM can provide a continuous representation of the objectives over the design space thereby reducing the number of numerical computations required for sensitivity analyses. This global approximation can also guide the optimization study. Additionally, a global approximation method can also at times address issues related to CFD analysis, for e.g. grid resolution, thereby increasing the fidelity of the computations. The present work is first of its kind where a CFDbased design optimization methodology is proposed for the design of an injector (Vaidyanathan et al., 2003a; Vaidyanathan et al., 2004b). Scope The focus is on sensitivity and tradeoff analyses for the design of a gaseous injector for liquid rocket propulsion. The data for the analyses is obtained from surrogate models (RSAs) of the objectives, and the multiobjective optima (Pareto optimal front, POF) generated by Goel et al. (2004) with the aid of multiobjective genetic algorithm (MOGA) and a local search method. The regions of the POF that represent different tradeoffs among the objectives are obtained through a hierarchical clustering algorithm. The initial study concentrates on the generation of response surface approximations (RSAs) and preliminary optimization studies. Then it is followed by an elaborate sensitivity and tradeoff analyses. The later part of the study is broadly divided into two parts. Firstly a sensitivity study is carried over the whole design space. The contribution of the design variables to the objective variability is calculated using a variancebased nonparametric approach (Sobol, 1993) and correlations between objectives are investigated. Secondly, sensitivity analyses are conducted on clusters along the POF. Box plots are used to highlight the variability of the design variables and the objectives within a given cluster. Additionally, the linear relationships between the variance of the design variables and the objectives are explored with the aid of partial correlation coefficients. Injector model Liquid rocket propulsion injector elements can be categorized into two basic types based on propellant mixing. The first type is an impinging element (Figure 4.1A) where mixing occurs by direct impingement of the propellant streams at an acute angle. The impingement enhances mixing by headon interaction between the oxidizer and fuel (Gill and Nurick, 1976). The second type of injector consists of nonimpinging elements where the propellant streams flow in parallel, typically in coaxial fashion (Figure 4.1B). Here, mixing is accomplished through a shearmixing process (Calhoon et al., 1973). From a design standpoint, both element types have appealing characteristics. However both also have undesirable characteristics. For instance, if the impinging element has an FOF arrangement (Figure 4.1 A), the mixing occurs rapidly, which can yield high performance. However, since the combustion zone is close to the injector face, the potential for high levels of injector face heating must be considered. If the nonimpinging element is assumed to be a shear coaxial element, mixing across the shear layer is relatively slow, requiring longer chambers to allow complete combustion. However, since the combustion zone is spread over a longer axial distance, the injector face is generally exposed to less severe thermal environments. Important design parameters for the impinging element (assuming fixed mass flow rates and constant propellant inlet conditions) include relative orifice size (or, relative stream momentum ratio), impingement angle and orifice spacing. Important parameters for the shear coaxial element (assuming fixed mass flow rates and constant propellant inlet conditions) include the area ratio of the two concentric tubes (or velocity ratio) and the shear area between the two propellant streams (i.e., the oxidizer post tip thickness) (Calhoon et al., 1973). 68 002 I ^T GH2 s_ GH2 G GH,2 A B Figure 4.1: Schematic of the GO2/GH2 impinging and coaxial injector elements. A) FO F impinging element. B) Coaxial element. One can combine the abovementioned characteristics of injectar types to develop hybrid concepts. For example, it has been noted (Gill and Nurick, 1976) that performance improvement in the shear coaxial element can be realized by directing the fuel toward the oxidizer stream rather than parallel to it and thinning the oxidizer post wall. The first modification causes the shear coaxial element to take on some of the aspects of the FOF impinging element. These notions lead to the hybrid element shown in Figure 4.2 that has been developed by The Boeing Company (U. S. Patent 6253539). Figure 4.2: Schematic of Hybrid Boeing Element (U. S. Patent 6253539) The four independent design variables chosen for the element used in this study are shown in Figure 4.3A. They are the angle at which the H2 is directed toward the oxidizer, a, the change in H2 flow area from a baseline area, AHA, the change in 02 flow area from a baseline area, AOA, and the oxidizer post tip thickness, OPTT. The fuel and oxidizer 69 flow rates are held constant. The independent variable ranges considered in this exercise are shown in Table 4.1. V TT max Injector detail TFmax TW4 Figure 4.3: Design variables and objectives of the single element rocket injector. A) Design variables and their range. B) Objective functions. / Xcc   I Table 4.1: Range of design variables (a is an acute angle in degrees and x is the thickness of OPTT in inches) Variable Minimum Maximum Actual Scaled Actual Scaled a 00 0 a 1 AHA Baseline 0 Baseline+25% 1 AOA Baseline40% 0 Baseline 1 OPTT x in 0 2x in 1 The dependent variables chosen for design evaluation are shown in Figure 4.3B. They are the maximum injector face temperature, TFmax, the wall temperature at a distance three inches from the injector face, TW4, the maximum oxidizer post tip temperature, TTmax and centerline axial location where the combustion is 99% complete, (hereafter referred to as the combustion length) Xc (Figure 4.3B). The three temperatures (calculated as adiabatic wall temperatures in this study) were chosen as indicators of local environments. Lower temperatures would indicate a design that had longer life/survivability due to decreased thermal strain on the part. The combustion length, Xc, was chosen as a measure of performance. Shorter combustion lengths indicate better performance. Numerical Procedure A pressurebased, finite difference, NavierStokes solver, FDNS500CVS (Chen, 1989; Wang and Chen, 1990; Chen and Farmer, 1991) is used in this study. The Navier Stokes equations, the twoequation turbulence model, and kinetic equations are solved. Convection terms are discretized using either a second order upwind, third order upwind or a central difference scheme, with adaptively added second order and fourth order dissipation terms. For the viscous and source terms, second order central difference scheme is used. First order upwind scheme is used for scalar quantities, like turbulence kinetic energy and species mass fractions, to ensure positive values. Steady state is assumed and an implicit Euler time marching scheme is used for computational efficiency. The chemical species continuity equations represent the H2 02 chemistry. It is represented with the aid of a 7species and 9reaction set (Chen, 1989; Wang and Chen, 1990; Chen and Farmer, 1991). The simulation domain and the boundary conditions used in all the CFD cases are shown in Figure 4.4. Because of the very large aspect ratio, both the injector and chamber have been shortened (at the cross hatched areas) for clarity. Both fuel and oxidizer flow in through the west boundary where the mass flow rate is fixed for both streams. The nozzle exit, at the east boundary, is modeled by an outlet boundary condition. The south boundary is modeled with the symmetry condition. All walls (both sides of the oxygen post, the outside of the fuel annulus, the outside chamber wall, and the faceplate) are modeled with the no slip adiabatic wall boundary condition. Each CFD analysis was done on a 200 CPU PC cluster with an AMD Athlon MP 1800 (1.8 GHz) chip and 1 GB RAM. All the cases were run concurrently and took about 5 days. a diabaic waM injector face injectorface oxygen pO st p hte symm etry Figure 4.4: Simulation domain and boundary conditions Design Space In this study OA (with n, = 4, q = 3, t = 3 and 2 = 2) is used to generate 54 designs for RSA and testing of the approximation. To test the RSA, 14 of the 54 designs are selected using crossvalidation techniques (Papila, 2001). During the CFD simulations, two of the 40 designs (fitting set) were found to be unacceptable because they exhibited unsteady behaviors while the numerical algorithm used was a steady state model. Hence, the final information included 38 designs for fitting the RSA and 14 to test their predictive capabilities. All the design variables are scaled between zero and one based on their upper and lower bounds. All the objectives obtained from the CFD solutions of the 52 valid designs are scaled between zero and one based on the maximum and minimum values and polynomialbased RSAs generated. Once the RSAs are generated, the scaled objectives are normalized between zero and one based on the maximum and minimum of the generated RSAs, which will be referred to as the normalized objectives in the text. The scaled values are the nonnormalized values of the objectives. The fitting and testing designs are shown in Tables 4.2 and 4.3, respectively. It should be noted that when the AOA is 1 or 0, the 02 flow area is reduced by 0% or 40%, respectively, as compared to the baseline area. Results and Discussion CFD Analysis Comparison of two of the evaluated designs serves to illustrate the motivation for combining CFD analyses and an efficient global approximation model during the design process. The scaled independent design variables are shown for the two cases in Table 4.4. In terms of the design space evaluated, these two designs are seen to be quite different. The normalized dependent variables are also shown in Table 4.4 with the temperatures shown in contour plots (Figures 4.5 and 4.6). Table 4.2: Fitting designs (unacceptable designs are marked in bold). Case # a AHA AOA OPTT Case 1: 0 0 1 0 Case 2: 0 0 1 0.5 Case 3: 0 0 1 1 Case 4: 0.5 0.5 0 1 Case 5: 1 1 0.5 0 Case 6: 1 1 0.5 0.5 Case 7: 1 1 0.5 1 Case 8: 0 0.5 0.5 0 Case 9: 0.5 1 1 0 Case 10: 0.5 1 1 0.5 Case 11: 0.5 1 1 1 Case 12: 1 0 0 0 Case 13: 1 0 0 0.5 Case 14: 1 0 0 1 Case 15: 0 1 0 0 Case 16: 0 1 0 0.5 Case 17: 0 1 1 1 Case 18: 0.5 0 0.5 0 Case 19: 0.5 0 0.5 1 Case 20: 1 0.5 1 0 Case 21: 1 0.5 1 0.5 Case 22: 1 0.5 1 1 Case 23: 0 1 0.5 0 Case 24: 0 1 0.5 1 Case 25: 0.5 0 1 0 Case 26: 0.5 0 1 1 Case 27: 1 0.5 0 0 Case 28: 1 0.5 0 1 Case 29: 0 0 0 0 Case 30: 0 0 0 0.5 Case 31: 0 0 1 1 Case 32: 0.5 0.5 0.5 0.5 Case 33: 1 1 1 0 Case 34: 1 1 1 1 Case 35: 0 0.5 1 0 Case 36: 0 0.5 1 1 Case 37: 0.5 1 0 0 Case 38: 0.5 1 0 1 Case 39: 1 0 0.5 0 Case 40: 1 0 0.5 1 Table 4.3: Testing designs.. Case# a AHA AOA OPTT Case 41: 0.5 0.5 0 0 Case 42: 0.5 0.5 0 0.5 Case 43: 0 0.5 0.5 0.5 Case 44: 0 0.5 0.5 1 Case 45: 0.5 0 0.5 0.5 Case 46: 0 1 0.5 0.5 Case 47: 0.5 0 1 0.5 Case 48: 1 0.5 0 0.5 Case 49: 0.5 0.5 0.5 0 Case 50: 0.5 0.5 0.5 1 Case 51: 1 1 1 0.5 Case 52: 0 0.5 1 0.5 Case 53: 0.5 1 0 0.5 Case 54: 1 0 0.5 0.5 Table 4.4: Independent and dependent variable (objectives) for Cases X and Y from CFD computations (nonnormalized values shown). Case a AHA AOA OPTT TFmax TW4 TTmax Xcc X 1 0 0 0 0.998 0.927 0.128 0.004 Y 0 0.5 0.5 1 0.285 0.395 0.923 0.567 . TF = 0.285 TW,= 0.395  TF,,, = 0.998 .TW.= 0.927 Figure 4.5: Temperature field for Cases X and Y. The chamber wall and injector face temperatures for Case Y (as seen in Figure 4.5) are low or moderately low, while for Case X, they are high. Figure 4.7 shows a large recirculation zone located between the injector and the chamber wall. This recirculating flow strips hot gases from the flame and causes them to flow back along the chamber wall and injector face. This phenomenon regulates the chamber wall temperature, TW4 and the injector face temperature, TFmax. Figure 4.6 shows that the other life/survivability indicating variable, the maximum oxidizer post tip temperature, TTmax, has essentially the opposite trend as compared to the other two temperatures. Figure 4.6: Nearinjector temperature field for Cases X and Y. Injector face 02 flow Figure 4.7: Large recirculation zone in the combustion chamber. The performance indicator, combustion length, Xc, is seen (Table 4.4) to be at the minimum for Case X (shorter combustion lengths indicate better mixing elements) and at a moderate level for Case Y. Given these observations, it is clear that the dependent variables exhibit competing trends such that no design will produce the "best" values for all the dependent variables as desired in this study. These comparisons confirm the past observation that changes in the injector design details have major effects on injector performance and injectorgenerated environments. Injector designs which address space propulsion goals must be produced by a tool that accounts for performance and multiple, spatially resolved environmental variables. Efficient, validated CFD codes that model sufficient injector physics are necessary to meet these requirements. Generation of large amounts of complex information by these codes produces the need for a means to manage the data. The injector designer must be able to confidently and efficiently sort through this database to locate an acceptable compromise design. Global sensitivity and approximation tools can guide the designer effectively. Grid Sensitivity Investigation Initially the 54 cases identified by DOE were computed on an axisymmetric geometry with 336x81 nodes. Only 33 out of the 40 fitting cases gave valid results. Results of the remaining seven cases contained unsteady features, which do not represent solutions of the steady state model employed. The RSAs generated with these data for TTmax and Xc, had Ra2 values of 0.961 and 0.976, respectively, suggesting a less than desirable fit. On checking the grid distribution in the combustion zone for the 33 cases used, it was determined that the grid resolution was insufficient. After a series of tests involving addition of grid points in the axial direction in the combustion zone, a 430x81 grid was found to be appropriate and used for the second run of the optimization study. To highlight the grid refinement, a comparison of the grid distributions is shown in Figure 4.8. The final grid was the product of tripling the axial node density in the combustion zone. The thick lines show the initial grid density, while the thin lines show final grid density. Note that, for clarity, only every sixth jline is shown. New RSAs were generated from new solutions obtained on the fine grid. The new fits for TTmax and Xc, had Ra2 values of 0.989 and 0.995, respectively, representing a considerable improvement over the RSA performance based on the first, coarser grid. This time, only two out of the 40 designs failed to produce valid results. It is to be noted that the Ra2 values for TFmax and TW4 are 0.999 for both the initial and final grid. This experience indicates that in addition to facilitating design optimization, the RSM can also help address the adequacy of the CFD solution accuracy. It offers insight into potential problems, based on the statistical regressions, from which the computations can be refined, thus improving the fidelity of the individual and collective databases. While this approach does not guarantee universally satisfactory outcomes, it does suggest clear directions to assess the critical area of data quality. ,, , Figure 4.8: Comparing the unrefined (336x81) {thick lines} and refined (430x81) {thin lines) grids. Response Surface Generation In the following, data for each design variable from the 38 acceptable fitting cases were used to generate the RSAs. Both full and reduced quadratic polynomials are generated. Table 4.5 identifies oa and ofor the four objective functions and different polynomials generated. The full quadratic model is consistent in performance in terms of both a and a. The reduced quadratic models either have a poorer value of oor offer only marginal improvement over the full response surface model. Since there is no appreciable improvement by reducing the fits, there may be noise in the data or the quadratic models (both full and reduced) do not sufficiently represent the data. Table 4.5: Performance of full and reduced quadratic RSAs for the four objective functions (nonnormalized values used). Full quadratic Reduced quadratic TFmax Number of 38 38 observations Ga 0.00566 0.00546 C (14 points) 0.00460 0.00463 Mean 0.495 0.495 TW4 Number of 38 38 observations Ga 0.00803 0.00795 C (14 points) 0.00669 0.00799 Mean 0.514 0.514 TTmax Number of 38 38 observations Ga 0.0413 0.0401 0 (14 points) 0.0396 0.0382 Mean 0.560 0.560 Xcc Number of 38 38 observations Ga 0.0205 0.0197 0 (14 points) 0.0178 0.0186 Mean 0.497 0.497 Comparing the full quadratic model predictions for the fitting cases to the CFD results of the various objectives, the variations for TFmax, TW4 and Xcc were found to be 79 negligible (Figures 4.9A, 4.9B and 4.9D), suggesting no need for further improvement. However, a large number of points lie away from the best fit in the plot of TTmax (Figure 4.9C). TFmaxQuadratic RS 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 RS TTmaxQuadratic RS 0 0.2 0.4 0.6 RS TW4Quadratic RS RSCFD Optimum 0 0.2 0.4 0.6 0.8 1 RS XccQuadratic RS 1 0.8 0.6 0.4 0.2 0.8 1 0 0.2 0.4 0.6 RS Figure 4.9: Comparison between the best fit possible and as predicted by quadratic response surface. Optimum refers to the case when RSA and CFD values are identical. RSCFD represents the value as for the current case (normalized values shown). A) TFmax. B) TW4. C) TTmax. D) Xcc. Based on this observation, a cubic model is generated for TTmax. A full cubic fit in a 4design variable model requires a minimum of 35 designs points. To obtain a good fit, 0.8 0.6 0.4 0.2 1 0.8 0.6 0.4 0.2 0.8 1 0 the number of design points should be considerably larger than the required number. Since there are only 38 design points available from the fitting set, the testing set is also included in the fitting set. The PRESSrms is used to estimate the performance of the generated RSA along with the o. Using the 52 design points, a full cubic was generated. The values of oa and PRESSrms were 0.0348 and 0.0598, respectively. The difference between the two measures of error is noticeable. Hence, a reduced cubic model is generated with the available design points. This improves the fit with values for oa and PRESSrms being 0.0303 and 0.0388, respectively. Table 4.6 compares the performance of the full quadratic and the reduced cubic models for TTmax. The reduced cubic model is seen to perform better statistically than the quadratic model. Hence, this cubic model is used for TTmax in the optimization studies that follow. The RSAs for all four objectives are shown as Eqs. (4.1)(4.4). These RSAs can then be used for sensitivity and optimization study. The RSAs are based on the scaled values of the design variables and obj ectives. Table 4.6: Performance of RSAs for the TTmax. Reduced cubic RSA has 21 coefficients (nonnormalized values used). Full quadratic Reduced cubic TTax Number of observations 38 52 Ga 0.0413 0.0303 PRESSrms 0.0521 0.0388 Mean 0.560 0.591 TFmax = 0.692 + 0.477(a) 0.687(AHA) 0.080(AOA) 0.0650(OPTT) 0.167(a)2 0.0129(AHA)(a) + 0.0796(AHA)2 0.0634(AOA)(a) 0.0257(AOA)(AHA) + 0.0877(AOA)2 0.0521(OPTT)(a) + 0.00156(OPTT)(AHA) + 0.00198(OPTT)(AOA) + 0.0184(OPTT)2. (4.1) TW4 = 0.758 + 0.358(a) 0.807(AHA) + 0.0925(AOA) 0.0468(OPTT) 0.172(a)2 + 0.0106(AHA)(a) + 0.0697(AHA)2 0.146(AOA)(a) 0.0416(AOA)(AHA) + 0.102(AOA)2 0.0694(OPTT)(a) 0.00503(OPTT)(AHA) + 0.0151(OPTT)(AOA) + 0.0173(OPTT)2. (4.2) TTma = 0.370 0.205(a) + 0.0307(AHA) + 0.108(AOA) + 1.019(OPTT) 0.135(a)2 + 0.0141(AHA)(a) + 0.0998(AHA)2 + 0.208(AOA)(a) 0.0301(AOA)(AHA) 0.226(AOA)2 + 0.353(OPTT)(a) 0.0497(OPTT)(AOA) 0.423(OPTT)2 + 0.202(AHA)(a)2 0.281(AOA)(a)2 0.342(AHA)2(a) 0.245(AHA)2(AOA) + 0.281(AOA)2(AHA) 0.184(OPTT)2(a) 0.281(AHA)( a)(AOA) (4.3) Xcc = 0.153 0.322(a) + 0.396(AHA) + 0.424(AOA) + 0.0226(OPTT) + 0.175(a)2 + 0.0185(AHA)(a) 0.0701(AHA)2 0.251(AOA)(a) + 0.179(AOA)(AHA) + 0.0150(AOA)2 + 0.0134(OPTT)(a) + 0.0296(OPTT)(AHA) + 0.0752(OPTT)(AOA) + 0.0192(OPTT)2 (4.4) Optimization Process The generated RSA are used to conduct an optimization exercise and to study the relationship between the design variables and the objectives that are indicators of life/survivability and performance. Also, the ability to accommodate design features that promote extended life/survivability in a design while maintaining reasonable performance is explored. Three separate optimization studies, using the approach of desirability functions as explained in chapter 3, are presented below. First, single objective minimizations are shown. Secondly, multiobjective optimizations are performed with equal weights. Finally, the multiobjective optimizations are conducted with variable weights. The obtained optimum solutions are compared with CFD computations. Since the objectives are scaled to 0(1), the error measures have to be scaled accordingly to estimate the accuracy of the obtained solutions. In terms of the actual values the error for an objective is defined as error = YCF YR (4.5) YCFD where YCFD is the solution obtained from the CFD analyses and yRp is the prediction of the RSA. Using simple mathematics, not shown here, the error in the scaled variables can be written as YCFD YRS error = (4.6) CFD +K where the bar represents the scaled values, and K is defined as K = Ymn (4.7) Ymax Y,, Here yn,, and ym, are the actual minimum and maximum values, respectively, based on the available set of fitting and testing data for that objective from the CFD analyses. In the dissertation, no bar is used for the notations used to represent the scaled design variables and objectives to avoid any confusion. Nonnormalized values of the objectives refer to the scaled values and the normalized values refer to the rescaled values of the objectives based on the maximum and minimum of the RSAs. This is to make some of the plots in the result section more meaningful. Singleobjective optimization The purpose of this effort is twofold. First, the goal is to verify the performance of the optimization methodology in locating the minimum values for single objectives in the chosen design space. This verification is straightforward and necessary, but not sufficient to conclude that the technique is useful for injector design. Table 4.7 shows the results of optimizing (here, minimizing) each objective separately. The numbers in parentheses in Table 4.7 indicate the weights applied during the optimization process. Here, a weight of one means that objective was included while a weight of zero means that objective was excluded. Based on the normalized values of the objectives the minimum value for each is necessarily 0. The optimizer found the correct value for each of the four cases. Table 4.7: Minimizing individual objectives. Value in parenthesis (1) indicates which objective function is minimized (normalized values shown). Opt a AHA AOA OPTT TFmax TW4 TTmax Xcc Case 1 0 1 0.592 1 0.0 0.0725 0.914 0.769 (1) (0) (0) (0) CFD 0.00207 0.0656 0.936 0.758 Error 0.11 0.08 0.70 0.25 (%) 2 0 1 0 1 0.0309 0.0 1.0 0.440 (0) (1) (0) (0) CFD 0.091 0.0461 0.911 0.568 Error 0.90 0.57 2.84 2.93 (%) 3 1 0 1 0 0.944 0.976 0.0 0.153 (0) (0) (1) (0) CFD 0.943 0.969 0.103 0.158 Error 0.01 0.08 4.46 0.12 (%) 4 0.917 0 0 0 0.987 0.926 0.182 0.0 (0) (0) (0) (1) CFD 0.981 0.919 0.119 0.004 Error 0.08 0.08 2.69 0.12 (%) Secondly, this process is helpful in understanding the injector operational influences via trend identification and variable groupings. The results from Table 4.7 are shown graphically in Figure 4.10 to facilitate the discussion. It was shown earlier that the flow entrained in the large recirculation zone (Figure 4.7) regulated both the maximum temperature on the injector face, TFmax and the chamber wall temperature, TW4. The results from OptCases 1 and 2 support this conclusion. Minimization of TFma very nearly minimizes TW4 and vice versa. Not surprisingly, the designs for the two cases are also quite similar. Both designs are shear coaxial (a= 0) with the hydrogen flow area, AHA, and the oxidizer post tip thickness, OPTT, at their maximum values. Table 4.7 shows that AOA is the only inconsistent variable in the two designs. In terms of the other objectives, when either TFmax or TW4 is minimized, the maximum oxidizer post tip temperature, TTmax, is high. The resulting moderatetolong combustion lengths, Xcc, are consistent with the relatively slow mixing expected from a shear coaxial element. Optimization Study Optimization Study 0.8 0.8 D TFmax o 00.6 A TW4 0.6 AHA So TTmax ~ A OA 0.x x Xcc x OPTT '0.4 0.4 0.2 0.2 x 0 A1 $ 0 1 2 3 4 1 2 3 4 OptCase OptCase A B Figure 4.10: Minimization of different objectives. (Case 1: TFmax, Case 2: TW4, Case 3: TTmax, Case 4: Xcc. Normalized values shown). A) Objectives. B) Design variables. Reference to OptCases 3 and 4 in Table 4.7 and Figure 4.10 also suggests a correlation between TTmax and Xcc, although not as tight as the one for TFmax and TW4. When TTmax is minimized, Xcc is low and vice versa. Further investigation needs to be done to completely understand the physics that underlies this correlation. Here, both designs are impinginglike with the hydrogen flow angle, a, at or near its maximum value and AHA and OPTT at their minimum values in the chosen design space. Again, Table 4.7 shows that AOA is the only inconsistent variable between the two cases. Both minimizations result in very high values for TFmax and TW4. Not surprisingly, the impinging elements represented by the two designs yield significantly shorter combustion lengths than the shear coaxial designs. The trends for a, AHA and OPTT are consistent between the two pairs of dependent variables but AOA varies among the four cases. With the other design variables set at the optimum design levels, Figure 4.11A shows the variation of TFmax and TW4 with AOA. A similar plot for TTmax and Xc, is shown in Figure 4.1 1B. The maximum injector face temperature, TFmax, is least sensitive to AOA. It is also the only objective to exhibit a minimum value in the interior of the design space. This finding suggests that expanding the design space could result in more robust designs for the multiobjective optimization. The trends for TW4 and Xc, are similar to each other, with AOA zero for both cases. The trend for TTmax relative to AOA is opposite. Table 4.7 and Figure 4.10 show the resulting design is different for each of the four singleobjective optimizations. These results indicate that designing robust injectors demands consideration of all the independent variables together. Figure 4.10 clearly shows the competing design trends. Thus, meeting a set of design requirements will require a compromise design. 86 Response of TFmax and TW4 Response Functions 0.2 1 0.9  0.16  0.8  S0.12 TW4 0.7  0.6  S0.08  0.08 0.4  0.3  0.04  0.2 TTmax Xcc TFmax.1  0 0.2 0.4 0.6 0.8 1 0 0 0.2 0.4 0.6 0.8 1 AOA AOA A B Figure 4.11: Variation of objectives with respect to AOA for fixed a, AHA and OPTT. (normalized values shown and D is the desirability function. All objectives are equally weighed in the composite function). A) TFmax and TW4 for a=0, AHA=1 and OPTT=1. B) TTmax and X,, for a=1, AHA=0 and OPTT=0. Multiobjective optimization To concurrently evaluate component life/survivability and performance considerations, a multiobjective optimization study is carried out. Starting with performance, X,,, as most of the injector design tools do, the objectives influencing thermal environments are added one at a time to study the effect on the resulting optimum designs. These optimum designs are presented in Table 4.8 and Figure 4.12 where OptCase 4, with X,, minimized is repeated as the starting point. The design for this case is an impinging element (a = 0.917) with minimum flow areas and the thinnest oxidizer post tip. The consequences of this design are a minimum Xc, a low TTmax and very high values of TFmax and TW4. When TTmax is minimized along with X,, in OptCase 5, the optimum design does not change significantly from that in OptCase 4. The hydrogen flow angle increases form 0.917 to 1.0. The values of AHA, AOA and OPTT remain unchanged from Opt Case 4. The combustion length, X,,, is unchanged as compared to OptCase 4, whereas TTmax is marginally improved from 0.181 to 0.152. The other thermal objectives remain at very high levels. Recall from Table 4.7 that AOA was 1.0 for the singleobjective optimization of TTmax. Table 4.8: Study of effect of life/survivability on performance. Value in parenthesis indicates the weighting on the objective functions (normalized values shown). Opt a AHA AOA OPTT TFmax TW4 TTmax Xcc Case 4 0.917 0 0 0 0.987 0.926 0.182 0.0 (0) (0) (0) (1) CFD 0.981 0.919 0.119 0.004 Error 0.08 0.08 2.69 0.12 (%) 5 1 0 0 0 1.0 0.928 0.151 0.00097 (0) (0) (1) (1) CFD 0.998 0.927 0.128 0.004 Error 0.03 0.01 1.00 0.14 (%) 6 1 1 0 0 0.376 0.224 0.155 0.279 (1) (0) (1) (1) CFD 0.369 0.214 0.252 0.264 Error 0.10 0.12 3.93 0.38 (%) 7 1 1 0 0 0.376 0.224 0.155 0.279 (1) (1) (1) (1) CFD 0.369 0.214 0.252 0.264 Error 0.10 0.12 3.93 0.38 (%) More insight regarding the optimum design can be gained by visualizing Figure 4.1 IB. Although the individual TTmax and Xcc drive the design towards the opposite ends of the range of AOA, the composite desirability function, D, has marginal variation with respect to AOA. The optimizer picks the value of AOA that gives the maximum value of the desirability function, D. The optimum design at this point is affected mostly by the 88 other design variables. This shows the benefits of using a global optimization procedure, which helps identify the trend of the composite function over the design space even when the individual objectives have opposing trends. Optimization Study Optimization Study 1 in  a1 aI e  0.8 0.8 D TFmax A TW4 0.6 Ta 0.6 o AHA o TTmax N x Xcc AOA '0.4 0.4 OPT x 0.2 A 0.2 0o 0 x: 0  4 5 6 7 4 5 6 7 OptCase OptCase A B Figure 4.12: Composite minimization of objectives with different weightings.(Case 4: (0,0,0,1), Case 5: (0,0,1,1), Case 6: (1,0,1,1), Case 7: (1,1,1,1). The values in parenthesis indicate weights for (TFmax, TW4, TTmax, Xce). Normalized values shown). A) Objectives. B) Design variables. OptCase 6 minimizes TFmax simultaneously with TTmax and Xc,. In terms of the design variables, AHA has now shifted from 0 to 1, while the other three remain at the same values as in OptCase 5. This shift in AHA is consistent with the singleobjective minimization of TFmax shown as OptCase 1 in Table 4.7. This single change in the design produces dramatic changes in the objectives. As desired, the value of TFmax has decreased significantly from 1.0 to 0.376. The concurrent large drop in the value of TW4 is consistent with the earlier conclusion linking TFmax and TW4. The oxygen post tip temperature, TTmax, remains essentially unchanged. The increase in Xc, from 0 to 0.279 is a negative impact of the design change. Figure 4.13 sheds more light on the situation. ion. 