<%BANNER%>

Investigation of Navier-Stokes Code Verification and Design Optimization


PAGE 1

INVESTIGATION OF NAVIER–STOKES 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

PAGE 2

Copyright 2004 by Rajkumar Vaidyanathan

PAGE 3

To my parents and sisters.

PAGE 4

iv 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 Fitz-Coy 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.

PAGE 5

v TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES............................................................................................................vii LIST OF FIGURES.........................................................................................................viii ABSTRACT....................................................................................................................... xi CHAPTER 1 INTRODUCTION AND SCOPE.....................................................................................1 Introduction ............................................................................................................1 Scope ............................................................................................................5 2 NAVIER-STOKES CODE VERIFICATION..................................................................9 Literature Review.......................................................................................................9 Scope ..........................................................................................................14 Approach ..........................................................................................................15 Richardson Extrapolation (RE).......................................................................15 Least Square Extrapolation (LSE)..................................................................16 Error versus Residual......................................................................................19 Flow Solver .......................................................................................21 Algorithm for LSE .......................................................................................23 LSE for pressure equation only.............................................................23 LSE for coupled pressure-velocity........................................................24 Results and Discussion.............................................................................................25 Two-Dimensional Turning Point Problem.....................................................26 Laminar Lid-Driven Cavity Flow...................................................................29 LSE for pressure only ....................................................................33 LSE for coupled pressure-velocity computations.................................41

PAGE 6

vi 3 GLOBAL SENSITIVITY AND OPTIMIZATION TECHNIQUES.............................44 Response Surface Methodology (RSM)...................................................................44 Design of Experiments (DOE).................................................................................49 Sensitivity Analysis..................................................................................................51 Optimization Algorithm...........................................................................................55 Multi-Objective Genetic Algorithm.........................................................................58 Hierarchical Clustering Method...............................................................................60 Visualization Using Box Plot...................................................................................60 4 DESIGN OPTIMIZATION – SINGLE ELEMENT ROCKET INJECTOR.................62 Literature Review.....................................................................................................62 Scope ..........................................................................................................66 Injector model ..........................................................................................................67 Numerical Procedure................................................................................................70 Design Space ..........................................................................................................71 Results and Discussion.............................................................................................72 CFD Analysis .......................................................................................72 Grid Sensitivity Investigation.........................................................................76 Response Surface Generation.........................................................................78 Optimization Process......................................................................................81 Single-objective optimization...............................................................82 Multi-objective optimization.................................................................86 Preliminary Observations...............................................................................93 Sensitivity and Trade-Off Analysis................................................................95 Global analysis ....................................................................95 Pareto front analyses ....................................................................99 5 SUMMARY, CONCLUSION AND FUTURE WORK...............................................107 Navier-Stokes Code Verification...........................................................................107 Design Optimization Study....................................................................................108 Future work ........................................................................................................111 Navier-Stokes Code Verification..................................................................112 Design Optimization .....................................................................................112 LIST OF REFERENCES.................................................................................................114 BIOGRAPHICAL SKETCH...........................................................................................121

PAGE 7

vii LIST OF TABLES Table page 4.1: Range of design variables ( is an acute angle in degrees and x is the thickness of OPTT in inches).......................................................................................................70 4.2: Fitting designs (unacceptable designs are marked in bold)........................................73 4.3: Testing designs........................................................................................................... 74 4.4: Independent and dependent variable (objectives) for Cases X and Y from CFD computations (non-normalized values shown)...........................................................74 4.5: Performance of full and reduced quadratic RSAs for the four objective functions (non-normalized values used).....................................................................................78 4.6: Performance of RSAs for the TTmax. Reduced cubic RSA has 21 coefficients (nonnormalized 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 shown).............................................................................................................91 4.10: List of essential () and non-essential () 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 optimal solution set.....................................................................................102 4.13: Partial correlation coefficients (Rcorr) of design variables vs. objectives for different clusters along the POF............................................................................................106

PAGE 8

viii LIST OF FIGURES Figure page 2.1: Collocated grid and notation for 2D grid....................................................................22 2.2: 2D turning point problem (dimensions are )......................................................26 2.3: Solution ( u ) contours for the 2D turning point problem, grid = 81 81....................27 2.4: Application: turning point problem with = 0.1, x -axis is the number of grid points in each direction for the fine grid on which extrapolation is done.............................28 2.5: Two dimensional lid-driven cavity flow (dimensions 1 1).....................................30 2.6: u -velocity contour for grid 121 121 and Re = 5.33.................................................31 2.7: u -velocity contour for grid 171 171 and Re = 1000................................................31 2.8: Pressure along the lid for grid 121 121 and Re = 5.33............................................32 2.9: Pressure along the lid for grid 171 171 and Re = 1000...........................................32 2.10: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with variable lid-velocity.................................................................................34 2.11: Least square L2 norm error and error in extrapolated pressure on grid 101 101 for different for Re = 5.33 and variable lid-velocity..................................................35 2.12: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 and variable lid-velocity..................................................................................36 2.13: Least square L2 norm error and error in extrapolated pressure on grid 101 101 for different for Re = 5.33 and constant lid-velocity (full domain)............................37 2.14: The shaded domain is used for LSE. A region between y = 0 0.95 for all x is used for extrapolation...............................................................................................37 2.15: Least square L2 norm error and error in extrapolated pressure on grid 101 101 for different for Re = 5.33 and constant lid-velocity (reduced domain)....................38

PAGE 9

ix 2.16: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with constant lid-velocity.................................................................................39 2.17: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 with constant lid-velocity................................................................................40 2.18: Extrapolated pressure and velocity compared with CFD solution on grid 121 121 ( Mn), Re = 5.33 and variable lid-velocity.................................................................42 2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 171 ( Mn), Re = 1000 and variable lid-velocity................................................................43 3.1: Schematic of the procedure for global design optimization (Papila, 2001)................46 3.2: Full factorial 3-level 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 Montgomery, 1995)....................................................................................................58 3.5: Schematic of a box plot..............................................................................................61 4.1: Schematic of the GO2/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: Near-injector temperature field for Cases X and Y....................................................75 4.7: Large recirculation zone in the combustion chamber.................................................75 4.8: Comparing the unrefined (336 81) {thick lines} and refined (430 81) {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. RSCFD 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: Xcc. Normalized values shown)...................................................................84

PAGE 10

x 4.11: Variation of objectives with respect to OA for fixed HA 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 HA or OA 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 shown)...........................................................................................................92 4.15: Variation of TFmax, TW4, TTmax and Xcc, with respect to OA for =0, HA=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: Two-objective Pareto optimal front........................................................................100 4.18: Pareto optimal solution set ( + ) and the multi-objective optima listed in Tables 4.74.9 (o).....................................................................................................................101 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

PAGE 11

xi 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 NAVIER-STOKES CODE VERIFICATION AND DESIGN OPTIMIZATION By Rajkumar Vaidyanathan August, 2004 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering With rapid progress made in employing computational techniques for various complex Navier-Stokes 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 CFD-based 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 model-based optimization and sensitivity evaluation. For Navier-Stokes (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

PAGE 12

xii 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 CFD-based 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 k turbulence closure. With the aid of these surrogate models, sensitivity and trade-off 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 multi-objective optimization study is carried out using a geometric mean approach. Following this, sensitivity analyses with the aid of variance-based non-parametric approach and partial correlation coefficients are conducted using data available from surrogate models of the objectives and the multi-objective 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 Navier-Stokes computations and also suggests tools for a designer to conduct design optimization study and related sensitivity analyses for a given design problem.

PAGE 13

1 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 multi-disciplinary optimization (MDO) which in turn promotes CFD-based design optimization studies. Computational fluid dynamics-based 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 Navier-Stokes (NS) equations and other related models introduces challenges that include physical modeling uncertainty, geometric complexities, non-uniform and non-orthogonal 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

PAGE 14

2 discretization convergence, insufficient convergence of an iterative procedure, computer round-off 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 a posteriori 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 a posteriori 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 well-resolved 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

PAGE 15

3 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,

PAGE 16

4 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 multi-objective problems. Multi-objective 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 nondominated. 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,

PAGE 17

5 that solution is said to be dominated. The function space of all the non-dominated solutions in the Pareto-optimal set is termed the Pareto-optimal front (POF). This set provides the designer with a clear idea of the trade-offs 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 multi-objective evolutionary algorithms (MOEA) have been proposed in literature (Deb, 2001). Insight into the trade-offs between different objectives can be obtained from the POF. In this dissertation different aspects of NS code verification and CFD-based 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 lid-driven 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 CFD-MDO 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 trade-offs involved in the design of a rocket injector. Scope The scope of the dissertation can be broadly divided into two, namely, Investigation of NS CFD code verification. CFD-based 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

PAGE 18

6 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 CFD-based 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 multi-element 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 O2 flow areas with fixed mass flow rates of fuel and oxidizer and O2 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.

PAGE 19

7 Sensitivity and trade-off analyses for the design are then carried out using the data obtained from surrogate models of the design objectives, and the POF (multi-objective optima) generated by Goel et al. (2004) with the aid of a multi-objective genetic algorithm (MOGA) and a local search method ( -constraint strategy). The regions of the POF that represent different trade-offs 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 variance-based non-parametric 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 of NS 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 two-dimensional turning point problem and a laminar lid-driven 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

PAGE 20

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

PAGE 21

9 CHAPTER 2 NAVIER-STOKES 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

PAGE 22

10 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 a priori error estimate and the later is referred to as a posteriori error estimate. In this dissertation the aim is to implement the a posteriori error estimates to a Navier-Stokes (NS) CFD code. Hence a literature survey of the available a posteriori 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 postprocessing 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 a posteriori 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 two-part paper introduced a recovery

PAGE 23

11 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 higherorder 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

PAGE 24

12 (Barth, 1993; Ollivier-gooch, 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

PAGE 25

13 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 a priori 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

PAGE 26

14 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 liddriven 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.

PAGE 27

15 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 a priori assumed order of accuracy of the continuous solution. Considering a flow-field output, say a velocity component or pressure, of a CFD simulation on a grid size, h and assuming second-order convergence a priori the Taylor series expansion can be written as 23 23()(;)... Uxuxhahah (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 23 23();... 248 hhh Uxuxaa !" #$ %& (2.2) Solving for U(x) using Eqs. (2.1) and (2.2) and eliminating O(h2) term and neglecting the higher-order terms leads to 31 ()4;; 32 h UxuxuxhOh'( !"#$ )* %& +, (2.3) This is a second-order RE which results in a gain of order one. Similarly if the solution is assumed to have first-order convergence a priori and O(h) term is eliminated, it results in a first-order RE with a gain of order one. 2()2;; 2 h UxuxuxhOh !" #$ %& (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

PAGE 28

16 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 h1 and h2, resulting in the following respective extrapolation scheme. 22 1221 22 12,, () huxhhuxh Ux hh !! (2.5) 1221 12,, () huxhhuxh Ux hh !! (2.6) where (x,hi) 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, u1 and u2, not necessarily based on grid doubling or halving, are combined linearly using a weighting function, which can be spatially dependent. The extrapolated solution is given as 121 Uuu !! (2.7)

PAGE 29

17 where !1 and !2 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, The weighting function, is given by 01 1cossin1M j j x xjx(2.8) with j, j = 0,…, M reals. The weighting function, is dependent on the spatial coordinate, x for a one-dimensional problem. For two dimensional problem it is a function of spatial co-ordinates, x and y such that x yxy (2.9) This modified Fourier expansion is ideal for rectangular domains. Different function may have to be used for more complex domains. The coefficients, j, are estimated by solving a least squares problem. For a given linear PDE, say L uf (2.10) with ,a auE and ,b bfE its numerical approximation can be written as hhh L uf (2.11) with ,h ha auE and ,h hb bfE Based on this, the least square problem can be formulated as P: Find (0,1) L such that 121hh L uuf '( +, !! (2.12) is minimum in L2(Mo) where Mo represents the fine grid to which the coarse grid solutions are interpolated. For a one-dimensional problem using two Fourier modes, Eq. (2.12) that has to be minimized can be rewritten as

PAGE 30

18 010211121cos1coshhh L uuLxuxuf'( '( +, +, !!!! (2.13) This approach can be generalized for nonlinear PDE problems via a Newton-like loop. For a given nonlinear problem, say, N uf (2.14) linearization results in J uvg (2.15) The least square problem then becomes 11 121211kkkk hh J uuuug'(+, !!!! (2.16) which is minimum in L2(Mo) The iteration is started from the initial condition of 0 = 0, until 1kk 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 !2 should be close to the grid solution on Mo. In situations where J ( u ) is not available, a nonlinear set in is 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 n -times 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

PAGE 31

19 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 n and U is the n 1 exact solution vector of the linear system. Assuming u to be a non-exact solution of size n 1, we have the following relationship: Au = b + r (2.18) where r represents the n 1 residual vector of the system of equations. Defining as the difference between u and U and subtracting Eq. (2.18) from Eq. (2.17) gives: A = -r (2.19) Denoting the eigenvectors of matrix A by ai, i = 1, …, n such that || ai|| = 1.0 and the eigenvalues as i, i = 1, …, n the following relationships can be obtained: 1n ii ia-r (2.20) where i, i = 1, …, n are linear weights associated with each eigenvector and 1 Ar (2.21) which can be rewritten as n iii iia-! (2.22) The residual, F, is then defined as

PAGE 32

20 2 1n i i F r(2.23) The least square approach aims at minimizing the L2 norm of the residual, F ; i.e., 22 11minminminnn ii iiFr!"!"#$#$ %&%&-(2.24) Now the adequacy of the non-exact solution, u, can be measure by the L2 norm error measured between the u and U, which can be written as: 2 1n i ierror -E" (2.25) Minimizing the error, gives 2 1minminn ii ierror!"#$ #$ %&(2.26) From Eqs. (2.24) and (2.26) it is clear that minimization of F is equivalent to minimization of L2 error norm ( error ) only when i’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 12 2 2 111NN m nijij nm jiEuu N'( )* +,--! (2.27) where n indicates the source of the non-exact solution, m indicates the source of the solution with respect to which the error is measured; i.e., the reference grid, ij nu represents the extrapolated or interpolated solution (non-exact solution) onto the fine

PAGE 33

21 grid, Mo, at the (i,j)th node, ij mu represents the reference grid solution, at the (i,j)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 in-house code STREAM which is based on the pressure-based 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 second-order 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: 0 uv x x (2.28) Momentum: PuvS xyxxyy !" !" #$ #$ %& %& (2.29) where is the density, u and v are the velocity components in Cartesian coordinates, is the laminar viscosity and SP represents the pressure gradient term. In Eq. (2.29), 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

PAGE 34

22 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. PPEEWWNNSSPCaaaaaSS (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. 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

PAGE 35

23 equations, into the discretized (in the FV sense) continuity equation. The obtained pressure equation is given as PPEEWWNNSSApApApApApS (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 v-component 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 non-linearity in the momentum equations (a’s have velocity components in them which make the equation non-linear) and the coupling of pressure-velocity 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 Navier-Stokes problem is tackled in two steps. In the first step, only the pressure is extrapolated and in the second step the pressure-velocity 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

PAGE 36

24 coarse grids used for extrapolation are used for the source term in Eq. (2.31). The extrapolated pressure, pLSE, is defined as 121LSE p pp!! (2.32) where 1 p !and 2 p 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 Navier-Stokes 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 Define pLSE =1 p + (1-)2 p !. Input pLSE, !2, 2v into Eq. (2.31). Find that minimizes the residual of Eq. (2.31). Calculate pLSE. LSE for coupled pressure-velocity 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 Picard-like iterative scheme (Ferziger and Peric, 1999) is adapted. The equation used for the velocity extrapolations comes from Eq. (2.30). 111111PPEEWWNNSSC LSELSELSELSELSEnnnnnnnnnnn PaaaaaSS (2.33) where SC n+1 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

PAGE 37

25 Define pLSE =1 p + (1-)2 p !. Input pLSE, !2, 2v !into Eq. (2.31) (uLSE and vLSE are used from the 2nd loop). Find that minimizes the residual of Eq. (2.31). Calculate pLSE. Define uLSE = !1 + (1-) !2. Input pLSE, uLSE, 2v !into Eq. (2.33) with replaced by u (vLSE is used from the 2nd loop). Find that minimizes the residual. Calculate uLSE. Define vLSE = 1v + (1-)2v Input pLSE, uLSE, vLSE into Eq. (2.33) with 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 ", # 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, two-dimensional turning point problem (Figure 2.2) and Laminar, lid-driven cavity flow (Figure 2.5). 1. variable lid-velocity 2. constant lid-velocity

PAGE 38

26 Two-Dimensional Turning Point Problem The linear equation that is solved is given below. 2,0,0, u uaxyx x (2.34) where = 0.1 and ,0.3 22axyxy!" !" #$ #$ %& %& (2.35) u = 0 u = 0 x y u = y( y) u = y( y) u = 0 u = 0 x y u = y( y) u = y( y) Figure 2.2: 2D turning point problem (dimensions are ). The boundary conditions of the problem are defined such that x,y [0, ]; u = y( -y) at x = 0 and u = -y( -y) at x = ; u = 0 at y = 0 and (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 (where the velocity changes direction) centered on the curve a(x,y) = 0 is not parallel to x or y axis thereby making the problem two-dimensional. A sample solution is shown in Figure 2.3. In the finite volume implementation, central difference is used for the diffusion term and

PAGE 39

27 first-order upwind is used for the convection term. Two modest base grids M1 and M2 of sizes 23 23 and 29 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 first-order and secondorder RE are evaluated and compared with the performance of LSE. For this case study LSE is implemented with 4 Fourier modes for in each spatial variable (instead of one mode as in the case of Navier-Stokes problem). The extrapolation is done onto fine grids, Mo, ranging from size 41 41 to 101 101 in steps of 10. In Figure 2.4, subscripts, RE1 and RE2 refers to RE of firstand second-order, respectively. The errors are represented as b a 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. 0 nM M E is the RMS error in u estimated using all the points in the computational domain, of Mo, as compared to the interpolated solution from a finer grid (Mn) of spacing h/2). 2 1 5 9 1 8 5 0 1 5 4 2 1 2 3 4 0 9 2 5 0 6 1 7 -0 6 1 7 0 3 0 8 0 3 0 8 0 0 0 0 0 0 0 0 0 3 0 8 0 3 0 8 0 6 1 7 0 6 1 7 0 9 2 5 0 9 2 5 1 2 3 4 1 2 3 4 1 5 4 2 1 8 5 0 1 5 9 X Y 0 1 2 3 0 0.5 1 1.5 2 2.5 3 Figure 2.3: Solution (u) contours for the 2D turning point problem, grid = 81 81.

PAGE 40

28 Em n vs N -3.5 -3 -2.5 -2 -1.5 -1 30405060708090100110 Grid Points (N)log10(Em n) EMo LSEEMo RE1EMo RE2EMo 29EMo 23EMn Mo Em n vs N -3.5 -3 -2.5 -2 -1.5 -1 30405060708090100110 Grid Points (N)log10(Em n) EMo LSEEMo RE1EMo RE2EMo 29EMo 23EMn Mo Figure 2.4: Application: turning point problem with = 0.1, x-axis 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 RE1 for grid with poor resolution. This is because the transition layer is under-resolved and the dominant error comes from the second order diffusion term. As the grid gets finer, performance of RE1 is better than RE2 as the first order convective error dominates. For some intermediate grid (51 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 23 and 29 29, there are only one or two grid points in the transition layer which has a thickness of = 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.

PAGE 41

29 This problem is relatively easy to solve and is used to confirm the potential of LSE 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, lid-driven cavity flow is used and the details are given in the following section. Laminar Lid-Driven Cavity Flow Garbey and Shyy (2003) have addressed the lid-driven cavity flow problem in a vorticity-streamfunction formulation with a finite difference implementation. In the current study the complete NS equations are solved in 2D for the lid-driven cavity flow and the pressure along with u and v-components of the velocity field are extrapolated. A laminar incompressible flow computation is carried out in a square cavity of dimensions 1 1 for Reynolds numbers 5.33 and 1000. This problem addresses the motion of the fluid in a square container induced due to the lid-velocity ulid in the positive x-direction (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: 2 2161liduxx (2.36)

PAGE 42

30 where ulid is the u-component of the velocity at the lid. The v-component at the lid is zero and x varies from 0 to 1. In the second case, the lid velocity is kept constant, ulid = 1.0. This results in singularities in u-velocity at the lid corners where there is an abrupt change of ulid from one to zero. u=0 v=0 u=ulid v=0 x y u=0 v=0 u=0 v=0 Figure 2.5: Two dimensional lid-driven cavity flow (dimensions 1 1). The Reynolds number for the flow is defined based on the mean lid-velocity. The u-velocity contours are shown in Figures 2.6 and 2.7 for Reynolds numbers 5.33 and 1000, respectively. Some difference is noticed in the u-velocity contours at the upper corner regions of the cavity for both the Reynolds numbers. The negative variable ulid results in a u-velocity profile which looks like the mirror image of the positive constant ulid 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 ulid is about 10 times more than that for the case with constant ulid (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.

PAGE 43

31 A B Figure 2.6: u-velocity contour for grid 121 121 and Re = 5.33. A) ulid = -16x2(1-x)2. B) ulid = 1.0. A B Figure 2.7: u-velocity contour for grid 171 171 and Re = 1000. A) ulid = -16x2(1-x)2. B) ulid = 1.0. The positive and negative values of pressure identifies whether the pressure is acting on the lid or away from it, respectively.

PAGE 44

32 A B Figure 2.8: Pressure along the lid for grid 121 121 and Re = 5.33. A) ulid = -16x2(1-x)2. B) ulid = 1.0. A B Figure 2.9: Pressure along the lid for grid 171 171 and Re = 1000. A) ulid = -16x2(1-x)2. B) ulid = 1.0. The results for the following different scenarios are presented below: LSE for Pressure only (Re = 5.33 and 1000) 1. ulid is variable. 2. ulid is constant. LSE for Pressure-Velocity (Re = 5.33 and 1000)

PAGE 45

33 1. ulid is variable. LSE for pressure only The results are presented in two parts. In part (a) the case with variable lid-velocity is presented and in part (b) the case with constant lid-velocity 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 61 and 71 71. The extrapolated pressure fields are compared to the solution obtained using a fixed fine grid (Mn) of size 121 121. The extrapolated fine grids (Mo) vary between sizes 81 81 and 121 121. For Re = 1000, the two coarse grids are 111 111 and 121 121 and the fixed fine grid (Mn) for comparison is 171 171. The extrapolated grids (Mo) are between 131 131 and 171 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 lid-velocity (Re = 5.33). In Figure 2.10A the comparison of extrapolated solutions on grid Mo (81 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 Mo 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.

PAGE 46

34 Em n vs N Re = 5.33-4.5 -4.3 -4.1 -3.9 -3.7 -3.5 -3.3 -3.1 -2.9 8090100110120130 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 71EMo 61 Em n vs N Re = 5.33-4.5 -4.3 -4.1 -3.9 -3.7 -3.5 -3.3 -3.1 -2.9 8090100110120130 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 71EMo 61 Em n vs N Re = 5.33 -4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 8090100110120130 Grid Points (N)log10(Em n) EMn RE1EMn 61EMn 71EMn RE2EMn LSE Em n vs N Re = 5.33 -4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 8090100110120130 Grid Points (N)log10(Em n) EMn RE1EMn 61EMn 71EMn RE2EMn LSE A B Figure 2.10: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with variable lid-velocity. 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 121. (m n 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 121, with respect to is examined. Figure 2.11 shows that the trend of both the accuracy measures are similar but based on minimization of F (Figure 2.11A) LSE picks a value of -1.48 for but based on L2 norm error (Figure 2.11B) a value of 1.6 for is a better choice. Although it does not identify similar 's that minimize both measures, they are close to each other. More assessment will be presented for the flow with a constant lid velocity. Variable lid-velocity (Re = 000). 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 111 and 121 121. It clearly shows in Figure

PAGE 47

35 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). vs F (Mo = 101) RE2 = -2.77, log10(F) = -9.945 LSE = -1.48, log10(F) = -9.976 -9.98 -9.97 -9.96 -9.95 -9.94 -9.93 -9.92 -9.91 -9.90 -3.50-2.50-1.50-0.500.501.50 log10(F) vs F (Mo = 101) RE2 = -2.77, log10(F) = -9.945 LSE = -1.48, log10(F) = -9.976 -9.98 -9.97 -9.96 -9.95 -9.94 -9.93 -9.92 -9.91 -9.90 -3.50-2.50-1.50-0.500.501.50 log10(F) vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -3.588 LSE = -1.48, log10(EMn LSE) = -3.865 -4.00 -3.80 -3.60 -3.40 -3.20 -3.50-2.50-1.50-0.500.501.50 log10(Em n) vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -3.588 LSE = -1.48, log10(EMn LSE) = -3.865 -4.00 -3.80 -3.60 -3.40 -3.20 -3.50-2.50-1.50-0.500.501.50 log10(Em n) A B Figure 2.11: Least square L2 norm error and error in extrapolated pressure on grid 101 101 for different for Re = 5.33 and variable lid-velocity. A) Least square L2 norm error (F). B) L2 norm error (m n E ) (comparing extrapolated pressure with CFD solution on grid 121 121 (Mn)). This exercise with variable lid-velocity 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.

PAGE 48

36 Em n vs N Re = 1000 -4.5 -4.3 -4.1 -3.9 -3.7 -3.5 -3.3 -3.1 -2.9 130140150160170180 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 121EMo 111 Em n vs N Re = 1000 -4.5 -4.3 -4.1 -3.9 -3.7 -3.5 -3.3 -3.1 -2.9 130140150160170180 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 121EMo 111 Em n vs N Re = 1000 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 130140150160170180 Grid Points (N)log10(Em n) EMn RE1EMn 111EMn 121EMn RE2EMn LSE Em n vs N Re = 1000 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 130140150160170180 Grid Points (N)log10(Em n) EMn RE1EMn 111EMn 121EMn RE2EMn LSE A B Figure 2.12: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 and variable lid-velocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (Mo). B) Extrapolated pressure compared with CFD solution on grid 171 171. (m n E : n refers to the source of the solution which is compared with the reference solution identified by m) Constant lid-velocity (Re = 5.33). At this point it would be interesting to take a look at a case where the lid-velocity 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, m n E with respect to is shown. LSE identifies = -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 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 = 0!1 and y = 0!~0.95.

PAGE 49

37 vs F (Mo = 101) RE2 = -2.77, log10(F) = -6.282 LSE = -0.48, log10(F) = -6.637-6.80 -6.60 -6.40 -6.20 -6.00 -5.80 -5.60 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(F) vs F (Mo = 101) RE2 = -2.77, log10(F) = -6.282 LSE = -0.48, log10(F) = -6.637-6.80 -6.60 -6.40 -6.20 -6.00 -5.80 -5.60 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(F) vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -2.365 LSE = -0.48, log10(EMn LSE) = -2.054 -2.60 -2.40 -2.20 -2.00 -1.80 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(Em n) vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -2.365 LSE = -0.48, log10(EMn LSE) = -2.054 -2.60 -2.40 -2.20 -2.00 -1.80 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(Em n) vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -2.365 LSE = -0.48, log10(EMn LSE) = -2.054 -2.60 -2.40 -2.20 -2.00 -1.80 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(Em n) A B Figure 2.13: Least square L2 norm error and error in extrapolated pressure on grid 101 101 for different for Re = 5.33 and constant lid-velocity (full domain). A) Least square L2 norm error (F). B) L2 norm error (m n E ) (comparing extrapolated pressure with CFD solution on grid 121 121 (Mn)) Figure 2.14: The shaded domain is used for LSE. A region between y = 0 0.95 for all x is used for extrapolation. As shown before in Figure 2.13 for the whole domain, the sensitivity of least square error norm, F and the L2 norm error (m n E ) in comparison to the CFD solution with respect to for the reduced domain, of grid 101 101, is shown in Figure 2.15. It is seen

PAGE 50

38 that LSE identifies = -0.48 as the optimum value for ", based on the minimization of F (Figure 2.15A). But as seen from Figure 2.15B this value of is not as far from the optimum value of based on the minimization of m n E The suggested value of is approximately -1.3 as compared to the previous approximate value of -4.0. This suggests that a sub-domain 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 corner singularity tends to get smeared onto the final grid. Some method of treating or accounting for this might have to be looked into. vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -2.267 LSE = -0.48, log10(EMn LSE) = -2.374 -2.50 -2.30 -2.10 -1.90 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(Em n) vs Em n (Mo = 101) RE2 = -2.77, log10(EMn RE2) = -2.267 LSE = -0.48, log10(EMn LSE) = -2.374 -2.50 -2.30 -2.10 -1.90 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(Em n) vs F (Mo = 101) RE2 = -2.77, log10(F) = -6.559 LSE = -0.48, log10(F) = -6.911 -7.00 -6.80 -6.60 -6.40 -6.20 -6.00 -5.80 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(F) vs F (Mo = 101) RE2 = -2.77, log10(F) = -6.559 LSE = -0.48, log10(F) = -6.911 -7.00 -6.80 -6.60 -6.40 -6.20 -6.00 -5.80 -6.5-5.5-4.5-3.5-2.5-1.5-0.50.51.52.5 log10(F) A B Figure 2.15: Least square L2 norm error and error in extrapolated pressure on grid 101 101 for different for Re = 5.33 and constant lid-velocity (reduced domain). A) Least square L2 norm error (F). B) L2 norm error (m n E ) (comparing extrapolated pressure with CFD solution on grid 121 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

PAGE 51

39 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 121. The solution on this grid is considered as the benchmark solution. Clearly LSE performs well as compared to RE. Em n vs N Re = 5.33 -2.5 -2.4 -2.3 -2.2 -2.1 -2 -1.9 -1.8 -1.7 -1.6 8090100110120130 Grid Points (N)log10(Em n) EMn RE1EMn 61EMn 71EMn RE2EMn LSE Em n vs N Re = 5.33 -2.5 -2.4 -2.3 -2.2 -2.1 -2 -1.9 -1.8 -1.7 -1.6 8090100110120130 Grid Points (N)log10(Em n) EMn RE1EMn 61EMn 71EMn RE2EMn LSE Em n vs N Re = 5.33-2.8 -2.6 -2.4 -2.2 -2 -1.8 -1.6 8090100110120130 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 71EMo 61 Em n vs N Re = 5.33-2.8 -2.6 -2.4 -2.2 -2 -1.8 -1.6 8090100110120130 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 71EMo 61 A B Figure 2.16: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 5.33 with constant lid-velocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (Mo). B) Extrapolated pressure compared with CFD solution on grid 121 121 (Mn). (m n E : n refers to the source of the solution which is compared with the reference solution identified by m) Constant lid-velocity (Re = 000). LSE is implemented in the same reduced domain as shown in Figure 2.14. The coarse grids used are 111 111 and 121 121 and

PAGE 52

40 the fixed fine grid (Mn), which is the representative of the correct solution, is 171 171. Figure 2.17, shows the comparison of different extrapolation schemes. Similar trends are noticed for the reduced domain as before. Em n vs N Re = 1000-4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 130140150160170180 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 121EMo 111 Em n vs N Re = 1000-4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 130140150160170180 Grid Points (N)log10(Em n) EMo RE1EMo LSEEMo RE2EMo 121EMo 111 Em n vs N Re = 1000 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 130140150160170180 Grid Points (N)log10(Em n) EMn RE1EMn 111EMn 121EMn RE2EMn LSE Em n vs N Re = 1000 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 130140150160170180 Grid Points (N)log10(Em n) EMn RE1EMn 111EMn 121EMn RE2EMn LSE A B Figure 2.17: Comparison of interpolated and extrapolated pressure with CFD solutions for Re = 1000 with constant lid-velocity. A) Extrapolated pressure compared with CFD solution on the extrapolated grid (Mo). B) Extrapolated pressure compared with CFD solution on grid 171 171. (m n 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, m n 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 m n E and the performance of LSE is improved.

PAGE 53

41 LSE for coupled pressure-velocity computations The coupled algorithm is implement on a lid-driven cavity flow with Re = 5.33 and 1000 and the variable lid-velocity as given in Eq. (2.36) is used. The geometry is same as before (Figure 2.5). The coarse grids (M1, M2) fine grids (Mo) and the fixed grids (Mn) are all the same as those used for pressure extrapolation. It is noticed that about 3-4 iterations of this algorithm results in converged values of and with about 3-4 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 m n E measure for pressure, u-velocity and v-velocity are shown in Figures 2.18A, 2.18B and 2.18C, respectively. Single mode (constant, 0) 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 = 000. 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, uand vvelocity, 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 u-velocity 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.

PAGE 54

42 Em n vs N Re = 5.33 pressure -3.8 -3.7 -3.6 -3.5 -3.4 8090100110120130 Grid Points (N)log10(Em n) EMn RE2EMn LSE Em n vs N Re = 5.33 pressure -3.8 -3.7 -3.6 -3.5 -3.4 8090100110120130 Grid Points (N)log10(Em n) EMn RE2EMn LSE Em n vs N Re = 5.33 u-velocity -4.2 -4.1 -4.0 -3.9 8090100110120130 Grid Points (N)log10(Em n) EMn RE2EMn LSE Em n vs N Re = 5.33 u-velocity -4.2 -4.1 -4.0 -3.9 8090100110120130 Grid Points (N)log10(Em n) EMn RE2EMn LSE A B Em n norm error vs N Re = 5.33 v-velocity-4.5 -4.4 -4.3 -4.2 -4.1 -4.0 8090100110120130 Grid Points (N)log10(Emn) EMn RE2EMn LSE Em n norm error vs N Re = 5.33 v-velocity-4.5 -4.4 -4.3 -4.2 -4.1 -4.0 8090100110120130 Grid Points (N)log10(Emn) EMn RE2EMn LSE C Figure 2.18: Extrapolated pressure and velocity compared with CFD solution on grid 121 121 (Mn), Re = 5.33 and variable lid-velocity. A) Pressure. B) u-velocity. C) v-velocity.

PAGE 55

43 Em n vs N Re = 1000 pressure -3.9 -3.8 -3.7 -3.6 -3.5 -3.4 130140150160170180 Grid Points (N)log10(Em n) EMn RE2EMn LSE Em n vs N Re = 1000 pressure -3.9 -3.8 -3.7 -3.6 -3.5 -3.4 130140150160170180 Grid Points (N)log10(Em n) EMn RE2EMn LSE Em n vs N Re = 1000 u-velocity-3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 130140150160170180 Grid Points (N)log10(Emn) EMn RE2EMn LSE Em n vs N Re = 1000 u-velocity-3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 130140150160170180 Grid Points (N)log10(Emn) EMn RE2EMn LSE A B Em n vs N Re = 1000 v-velocity-3.7 -3.6 -3.5 -3.4 -3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 130140150160170180 Grid Points (N)log10(Emn) EMn RE2EMn LSE Em n vs N Re = 1000 v-velocity-3.7 -3.6 -3.5 -3.4 -3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 130140150160170180 Grid Points (N)log10(Emn) EMn RE2EMn LSE Em n vs N Re = 1000 v-velocity-3.7 -3.6 -3.5 -3.4 -3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 130140150160170180 Grid Points (N)log10(Emn) EMn RE2EMn LSE C Figure 2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 171 (Mn), Re = 1000 and variable lid-velocity. A) Pressure. B) u-velocity. C) v-velocity.

PAGE 56

44 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 multi-disciplinary 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 single-objective optimization study which is dependent on two design variables. This is a general case and can be extended for multi-objective 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 polynomial-based RSAs of assumed order and unknown coefficients. The number of coefficients to be evaluated depends on the order of the polynomial and the

PAGE 57

45 number of design parameters involved. For instance, a second-order polynomial of N 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 -1 12 1 2 1 0 k i k j j i ij k i i ii k i i ix x x x y (3.1) The above equation can be written in matrix notation as follows X y (3.2) where y is the (n 1) vector of observations, X: (nnp) matrix of the levels of the independent variables, : (np 1) vector of the regression coefficients, !: (n 1) vector of random error, n: the number of observations, and np: the number of terms in the RSA. The purpose is to find the vector of least square estimators b that minimizes 2 1()()in TT iErroryXyX (3.3) which yields the least squares estimator of 1()TTbXXXy (3.4)

PAGE 58

46 Design Variable 2 Design Variable 1 Objective Function 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 rmserror a (Myers and Montgomogery, 1995) value that is defined as: 2 i a pe nn = "# (3.5) where ei 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,

PAGE 59

47 different from those used to generate the fit. The prediction rms-error for the test set is given by: 2 im (3.6) In this equation i is the error at the ith test point and m is the number of test points. The test points can be selected using the cross-validation 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 k-fold crossvalidation (Mullin and Sukthankar, 2000). In k-fold cross-validation, the data points are divided in k equally sized subsets. One of the subset is kept for testing and the remaining k-1 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 k-1 sets which has the minimum prediction rmserror, 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 rms-error. 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 rms-error is given by 2 1ˆn ii i rmsyy PRESS n (3.7)

PAGE 60

48 where yi is the expected value, ˆiy 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 rms-eror, a, 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 21E y ySS R SS (3.8) where SSE is the sum of squares of the residuals (=2 1ˆ ()n ii iyy-) where ˆ yis the predicted value by the RSA. SSyy is the total sum of squares about the mean given by n 2 yyERi i1SSSSSS(yy)(3.9) where yis an overall average of yi. SSR is the sum of squares due to regression (= 2 1ˆ ()n i iyy-) where y is the overall average of yi. Since R2 increases as more terms are added in the RSA, the adjusted coefficient of multiple determination 2 a R a normalized measure is usually preferred. It accounts for the degrees of freedom in the model and is given by 22/() 1 111 /(1)Ep a yypSSnn n R R SSnnn!" #$ #$%& (3.10) For a good fit, the value of 2 a R should be closer to one. The polynomial-based 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 polynomial-based RSA allows us to use statistical techniques known as design of

PAGE 61

49 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).

PAGE 62

50 Figure 3.2: Full factorial 3-level 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 j ia where j denotes the row (j = 1,2… nr) and i denotes the column (i = 1,2…nc) that j iabelongs to, supposing that each j iaQ = {0,1…q-1}. 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 nr-row-by-t-column sub-matrix, the qt possible distinct rows occur (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. X Y z X Y z X Y z X Y z

PAGE 63

51 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, …, xn, 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. j i n n j i ij i i ix x x f x x f x f f X f ,..., ... ,2 1 ... 12 0 (3.11) An analysis of variance (ANOVA) (Archer et al., 1997) representation of Eq. (3.11) can be obtained by imposing the following condition .1 0 ...01k i idx fs (3.12) for k = i1, …, is. 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 si if...1is responsible for the interactive contribution of design variables si ix x ,...,1to the variability of the objective f (X) over the n–dimensional unit hypercube design space. Integrating Eq. (3.11) over all the design variables gives 0 k f xdxf (3.13) Integrating Eq. (3.11) over all the design variables except xi gives i k i i kx f f dx x f0 (3.14) from which fi( xi) is obtained.

PAGE 64

52 Integrating Eq. (3.11) over all the design variables except xi and xj gives j i k j i j i j j i i kx x f x f x f f dx x f, 0, (3.15) from which fi,j( xi,xj) is obtained. This can similarly be extended to obtain the remaining summands of f (X). The imposed condition ensures that the summands are orthogonal in nature and since f (X) is square integrable the summands are too. Therefore the partial variances can be calculated as .s s si i i i i idx dx f D ...1 1 12 ... ... (3.16) and the total variances is 2 0 2f dx x f D (3.17) Therefore it can be shown (Sobol, 1993) that -n s n i i i is sD D1... ...1 1 (3.18) 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 variance-based non-parametric approach to estimate the global sensitivity using Monte Carlo methods. To calculate the total sensitivity of any design variable, say x1, the design variable vector X is divided into two complementary subsets, x1 and Z where Z is a vector containing x2, x3, …, xn. The purpose of using these subsets is to isolate the influence of x1 on f (X) variability from the influence of the remaining design variables included in Z. The total sensitivity index for x1 is then defined as

PAGE 65

53 D D Stotal x total x1 1 (3.19) where Z x x total xD D D,1 1 1 (3.20) ixD is the partial variance associated with x1 and Z xD,1 is the sum of all partial variances associated with any combination of the remaining variables representing the interactions between x1 and Z. Similarly the partial variance associated with Z can be defined as DZ. Therefore the total objective variability can be written as Z x Z xi iD D D D, (3.21) In cases where the f (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: -!N j jf Z x f N1 0 1, 1 (3.22) where f0 is the mean of the objective. "N j x jf D Z x f Z x f N1 2 0 1 11, 1 (3.23) !N j jf D Z x f N1 2 0 1 2, 1 (3.24)

PAGE 66

54 "N j Z jf D Z x f Z x f N1 2 0 1 1, 1 (3.25) where the terms ( x1, Z) and ( x1’ Z ’) represent random sample designs. Equations. (3.23), (3.24) and (3.25) give the estimates for 1xD, D and DZ, respectively and Z xiD, 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 D D Si ix x (3.26) Each pair of random samples requires three different objective function evaluations (e.g., f ( x1, Z), f ( x1, Z ’) and f ( x1’ 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 multi-objective 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

PAGE 67

55 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 variance-based parametric approach that gives a measure of the linear relationship between the variances of a design variable, say x1 and an objective f (X), after the influence of other variables have been filtered out. It gives a measure of expected change in f (X) per unit change in x1 which in other words gives the sensitivity of the function to the design variable. Linear approximations are obtained for x1 and f (X) as a function of the components of Z and the residuals (say, r1 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 121122 corr rr E rrrr R '( +, (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 gradient-based 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 { f (X)}subject to lb X ub, where lb is the lower boundary vector and ub is the upper boundary vector of the

PAGE 68

56 design variables vector X. If the goal is to maximize the objective function then f(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). RS generated Evaluate gradient Yes No Initial Design Variables Perturb Design variables Optimum Soln. Converged? Data set Generating RS First Phase Second Phase 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

PAGE 69

57 lie between zero and one. In the case where an objective should be maximized, say R1, the desirability takes the form: 1 1t R A d B A !"#$%& (3.28) where B is the target value and A is the lowest acceptable value such that d1 = 1 for any R1 > B and d1 = 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: 2 2 s R E d CE !"#$%& (3.29) 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. m 1 m 3 2 1d .... d d d D # # (3.30) This is then submitted to an optimization toolbox to be maximized.

PAGE 70

58 A A B BResponse Value 0.00.0 0.20.2 0.40.4 0.60.6 0.80.8 1.01.0Desirability, dDesirability, d s=0.2 s=0.4 s=1 s=2 s=8 A A B BResponse Value 0.00.0 0.20.2 0.40.4 0.60.6 0.80.8 1.01.0Desirability, dDesirability, d s=0.2 s=0.4 s=1 s=2 s=8 Figure 3.4: Desirability function (d) for various weight factors, s (Note: B < A) (Myers and Montgomery, 1995). Multi-Objective Genetic Algorithm One of the other optimization approaches is to use Multi-objective 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 multi-objective 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 nondominated 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

PAGE 71

59 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 hyper-surfaces. One of the MOGA that has been shown to be effective in finding the Paretooptimal solutions is the elitist non-dominated sorting genetic algorithm (NSGA-II) (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 non-domination 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 -constraint strategy (Deb and Goel, 2001) helps reduce the dimensionality of the problem. Sequential quadratic programming, available in

PAGE 72

60 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 y-axis 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 “H-spread”. The inner fence is located 1 step beyond the hinges, which is equal to 1.5 times the H-spread. 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

PAGE 73

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. 0 1 Inner f ence Upper Adjacent Value Upper hinge Median Lower hinge O utside value Lower Adjacent Value Figure 3.5: Schematic of a box plot

PAGE 74

62 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 empirical-based (Rupe, 1965; Dickerson et al., 1968; Pavli, 1979; Nurick, 1990; Pieper,

PAGE 75

63 1991). Extensive suband full-scale hot-fire 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 one-dimensional 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 cold-flow and hot-fire 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

PAGE 76

64 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, Vf/Vo 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 coaxial 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, Po, fuel pressure drop, Pf, combustor length, Lcomb, and the impingement half-angle, as design variables. Objectives included ERE (a measure of element performance), wall heat flux, Qw, injector heat flux, Qinj, relative combustor weight, Wrel, and relative injector cost, Crel. 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, Po, fuel pressure drop, Pf, combustor length, Lcomb, and the full cone swirl angle, $, as design variables. The objectives modeled were ERE (a measure of element performance),

PAGE 77

65 wall heat flux, Qw, injector heat flux, Qinj, relative combustor weight, Wrel, and relative injector cost, Crel. 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 three-dimensional geometry of multi-element 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

PAGE 78

66 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 CFD-based 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 trade-off 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 multi-objective optima (Pareto optimal front, POF) generated by Goel et al. (2004) with the aid of multi-objective genetic algorithm (MOGA) and a local search method. The regions of the POF that represent different trade-offs 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 trade-off 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 variance-based non-parametric 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.

PAGE 79

67 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 head-on interaction between the oxidizer and fuel (Gill and Nurick, 1976). The second type of injector consists of non-impinging elements where the propellant streams flow in parallel, typically in coaxial fashion (Figure 4.1B). Here, mixing is accomplished through a shear-mixing 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 F-O-F arrangement (Figure 4.1A), 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 non-impinging 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).

PAGE 80

68 dodf GO2 GH2GH2 s dodf GO2 GH2GH2 s GO2GH2GH2 do dfidfo GO2GH2GH2 do dfidfo GO2GH2GH2 do dfidfo A B Figure 4.1: Schematic of the GO2/GH2 impinging and coaxial injector elements. A) F-OF impinging element. B) Coaxial element. One can combine the above-mentioned characteristics of injector 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 F-O-F 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, the change in H2 flow area from a baseline area, HA, the change in O2 flow area from a baseline area, OA, and the oxidizer post tip thickness, OPTT. The fuel and oxidizer

PAGE 81

69 flow rates are held constant. The independent variable ranges considered in this exercise are shown in Table 4.1. A B Figure 4.3: Design variables and objectives of the single element rocket injector. A) Design variables and their range. B) Objective functions. O2 Flow Area ( OA=0--40%) H2 Flow Area ( HA=0-25%) O2 Post Tip Thickness (OPTT=0.01-0.02 in.) H2 Flow Angle ( =0o-20o) H2 O2 H2 O2 Flow Area ( OA=0--40%) H2 Flow Area ( HA=0-25%) O2 Post Tip Thickness (OPTT=0.01-0.02 in.) H2 Flow Angle ( =0o-20o) H2 O2 H2 H2 Flow Angle (000) O2 Post Tip Thickness (OPTT=x-2x in) Injector detail XccTTmaxTF max TW4

PAGE 82

70 Table 4.1: Range of design variables ( is an acute angle in degrees and x is the thickness of OPTT in inches) Minimum Maximum Variable Actual Scaled Actual Scaled 0o 0 o 1 HA Baseline 0 Baseline+25% 1 OA Baseline-40% 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) Xcc (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, Xcc, was chosen as a measure of performance. Shorter combustion lengths indicate better performance. Numerical Procedure A pressure-based, finite difference, Navier-Stokes solver, FDNS500-CVS (Chen, 1989; Wang and Chen, 1990; Chen and Farmer, 1991) is used in this study. The NavierStokes equations, the two-equation 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

PAGE 83

71 assumed and an implicit Euler time marching scheme is used for computational efficiency. The chemical species continuity equations represent the H2 O2 chemistry. It is represented with the aid of a 7-species and 9-reaction 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. Figure 4.4: Simulation domain and boundary conditions Design Space In this study OA (with nc = 4, q = 3, t = 3 and = 2) is used to generate 54 designs for RSA and testing of the approximation. To test the RSA, 14 of the 54 designs are

PAGE 84

72 selected using cross-validation 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 polynomial-based 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 non-normalized 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 OA is 1 or 0, the O2 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).

PAGE 85

73 Table 4.2: Fitting designs (unacceptable designs are marked in bold). Case # HA OA 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 7: 0 " 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 3 : 0 0 " 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

PAGE 86

74 Table 4.3: Testing designs.. Case # HA OA 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 (non-normalized values shown). Case HA OA 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 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

PAGE 87

75 indicating variable, the maximum oxidizer post tip temperature, TTmax, has essentially the opposite trend as compared to the other two temperatures. CaseYCaseXTTmax = 0.923 TTmax = 0.128 Figure 4.6: Near-injector temperature field for Cases X and Y. H2 flow O2 p ost ti p wall O2flow Injector face Large recirculation zone Figure 4.7: Large recirculation zone in the combustion chamber.

PAGE 88

76 The performance indicator, combustion length, Xcc 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 injector-generated 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 axi-symmetric geometry with 336 81 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 Xcc had Ra 2 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

PAGE 89

77 involving addition of grid points in the axial direction in the combustion zone, a 430 81 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 j-line is shown. New RSAs were generated from new solutions obtained on the fine grid. The new fits for TTmax and Xcc had Ra 2 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 Ra 2 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 (336 81) {thick lines} and refined (430 81) {thin lines} grids.

PAGE 90

78 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 a and for the four objective functions and different polynomials generated. The full quadratic model is consistent in performance in terms of both a and The reduced quadratic models either have a poorer value of or 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 (non-normalized values used). Full quadratic Reduced quadratic Number of observations 38 38 a 0.00566 0.00546 (14 points) 0.00460 0.00463 TFmax Mean 0.495 0.495 Number of observations 38 38 a 0.00803 0.00795 (14 points) 0.00669 0.00799 TW4 Mean 0.514 0.514 Number of observations 38 38 a 0.0413 0.0401 (14 points) 0.0396 0.0382 TTmax Mean 0.560 0.560 Number of observations 38 38 a 0.0205 0.0197 (14 points) 0.0178 0.0186 Xcc 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

PAGE 91

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). TFmax-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimum TFmax-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimum TW4-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimum TW4-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimum A B Xcc-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimu m Xcc-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimu m Xcc-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimu m TTmax-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimu m TTmax-Quadratic RS 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 RSCFD RS-CFD Optimu m C D 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. RS-CFD 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 4-design variable model requires a minimum of 35 designs points. To obtain a good fit,

PAGE 92

80 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 a. Using the 52 design points, a full cubic was generated. The values of a 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 a 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 objectives. Table 4.6: Performance of RSAs for the TTmax. Reduced cubic RSA has 21 coefficients (non-normalized values used). Full quadratic Reduced cubic Number of observations 38 52 a 0.0413 0.0303 PRESSrms 0.0521 0.0388 TTmax Mean 0.560 0.591 TFmax = 0.692 + 0.477() – 0.687( HA) – 0.080( OA) – 0.0650(OPTT) – 0.167()2 – 0.0129( HA)() + 0.0796( HA)2 – 0.0634( OA)() – 0.0257( OA)( HA) + 0.0877( OA)2 – 0.0521(OPTT)() + 0.00156(OPTT)( HA) + 0.00198(OPTT)( OA) + 0.0184(OPTT)2. (4.1)

PAGE 93

81 TW4 = 0.758 + 0.358() – 0.807( HA) + 0.0925( OA) – 0.0468(OPTT) – 0.172()2 + 0.0106( HA)() + 0.0697( HA)2 – 0.146( OA)() – 0.0416( OA)( HA) + 0.102( OA)2 – 0.0694(OPTT)() – 0.00503(OPTT)( HA) + 0.0151(OPTT)( OA) + 0.0173(OPTT)2. (4.2) TTmax = 0.370 – 0.205() + 0.0307( HA) + 0.108( OA) + 1.019(OPTT) – 0.135()2 + 0.0141( HA)() + 0.0998( HA)2 + 0.208( OA)() – 0.0301( OA)( HA) – 0.226( OA)2 + 0.353(OPTT)() – 0.0497(OPTT)( OA) – 0.423(OPTT)2 + 0.202( HA)()2 – 0.281( OA)()2 – 0.342( HA)2() – 0.245( HA)2( OA) + 0.281( OA)2( HA) – 0.184(OPTT)2() – 0.281( HA)( )( OA) (4.3) Xcc = 0.153 0.322() + 0.396( HA) + 0.424( OA) + 0.0226(OPTT) + 0.175()2 + 0.0185( HA)() – 0.0701( HA)2 – 0.251( OA)() + 0.179( OA)( HA) + 0.0150( OA)2 + 0.0134(OPTT)() + 0.0296(OPTT)( HA) + 0.0752(OPTT)( OA) + 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, singleobjective minimizations are shown. Secondly, multi-objective optimizations are performed with equal weights. Finally, the multi-objective optimizations are conducted with variable weights. The obtained optimum solutions are compared with CFD

PAGE 94

82 computations. Since the objectives are scaled to O(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 CFDRS CFDyy error y (4.5) where yCFD is the solution obtained from the CFD analyses and yRS is the prediction of the RSA. Using simple mathematics, not shown here, the error in the scaled variables can be written as CFDRS CFDyy error yK (4.6) where the bar represents the scaled values, and K is defined as min maxminy K yy (4.7) Here ymin and ymax 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. Non-normalized values of the objectives refer to the scaled values and the normalized values refer to the re-scaled 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. Single-objective 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

PAGE 95

83 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). OptCase HA OA OPTT TFmax TW4 TTmax Xcc 1 0 1 0.592 1 0.0 (1) 0.0725 (0) 0.914 (0) 0.769 (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.0 (1) 1.0 (0) 0.440 (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) 0.976 (0) 0.0 (1) 0.153 (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) 0.926 (0) 0.182 (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

PAGE 96

84 results from Opt-Cases 1 and 2 support this conclusion. Minimization of TFmax very nearly minimizes TW4 and vice versa. Not surprisingly, the designs for the two cases are also quite similar. Both designs are shear coaxial ( = 0) with the hydrogen flow area, HA, and the oxidizer post tip thickness, OPTT, at their maximum values. Table 4.7 shows that OA 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 moderate-to-long combustion lengths, Xcc, are consistent with the relatively slow mixing expected from a shear coaxial element. Optimization Study0 0.2 0.4 0.6 0.8 1 1234 Opt-CaseObjective TFmax TW4 TTmax Xcc Optimization Study0 0.2 0.4 0.6 0.8 1 1234 Opt-CaseDesign Variable %& '& OPTT 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 Opt-Cases 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

PAGE 97

85 designs are impinging-like with the hydrogen flow angle, at or near its maximum value and HA and OPTT at their minimum values in the chosen design space. Again, Table 4.7 shows that OA 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 HA and OPTT are consistent between the two pairs of dependent variables but OA 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 OA. A similar plot for TTmax and Xcc is shown in Figure 4.11B. The maximum injector face temperature, TFmax, is least sensitive to OA. 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 multi-objective optimization. The trends for TW4 and Xcc are similar to each other, with OA zero for both cases. The trend for TTmax relative to OA is opposite. Table 4.7 and Figure 4.10 show the resulting design is different for each of the four single-objective 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.

PAGE 98

86 Response of TFmax and TW4 0 0.04 0.08 0.12 0.16 0.2 00.20.40.60.81 '&Respons e TFmax TW4 Response of TFmax and TW4 0 0.04 0.08 0.12 0.16 0.2 00.20.40.60.81 '&Respons e TFmax TW4 Response Functions0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 00.20.40.60.81 '&Response TTmaxXcc D Response Functions0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 00.20.40.60.81 '&Response TTmaxXcc D Response Functions0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 00.20.40.60.81 '&Response TTmaxXcc D A B Figure 4.11: Variation of objectives with respect to OA for fixed HA 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 =0, HA=1 and OPTT=1. B) TTmax and Xcc for =1, HA=0 and OPTT=0. Multi-objective optimization To concurrently evaluate component life/survivability and performance considerations, a multi-objective optimization study is carried out. Starting with performance, Xcc, 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 Opt-Case 4, with Xcc minimized is repeated as the starting point. The design for this case is an impinging element ( = 0.917) with minimum flow areas and the thinnest oxidizer post tip. The consequences of this design are a minimum Xcc, a low TTmax and very high values of TFmax and TW4. When TTmax is minimized along with Xcc in Opt-Case 5, the optimum design does not change significantly from that in Opt-Case 4. The hydrogen flow angle increases

PAGE 99

87 form 0.917 to 1.0. The values of HA, OA and OPTT remain unchanged from OptCase 4. The combustion length, Xcc, is unchanged as compared to Opt-Case 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 OA was 1.0 for the single-objective 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). OptCase HA OA OPTT TFmax TW4 TTmax Xcc 4 0.917 0 0 0 0.987 (0) 0.926 (0) 0.182 (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) 0.928 (0) 0.151 (1) 0.00097 (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 (1) 0.224 (0) 0.155 (1) 0.279 (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 (1) 0.224 (1) 0.155 (1) 0.279 (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.11B. Although the individual TTmax and Xcc drive the design towards the opposite ends of the range of OA, the composite desirability function, D, has marginal variation with respect to OA. The optimizer picks the value of OA that gives the maximum value of the desirability function, D. The optimum design at this point is affected mostly by the

PAGE 100

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 Study0 0.2 0.4 0.6 0.8 1 4567 Opt-CaseObjective TFmax TW4 TTmax Xcc Optimization Study0 0.2 0.4 0.6 0.8 1 4567 Opt-CaseDesign Variable %& '& OPTT 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, Xcc). Normalized values shown). A) Objectives. B) Design variables. Opt-Case 6 minimizes TFmax simultaneously with TTmax and Xcc. In terms of the design variables, HA has now shifted from 0 to 1, while the other three remain at the same values as in Opt-Case 5. This shift in HA is consistent with the single-objective minimization of TFmax shown as Opt-Case 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 Xcc from 0 to 0.279 is a negative impact of the design change. Figure 4.13 sheds more light on the situation.

PAGE 101

89 Figure 4.13A plots the individual objectives as a function of HA. Figure 4.13B is similar, with the individual objectives plotted as a function of OA. These figures clearly illustrate that both TFmax and TW4 are considerably more sensitive to HA than OA. Accordingly, the optimizer uses the hydrogen flow area to most efficiently regulate TFmax. Figure 4.13A also shows the positive slope of Xcc relative to HA, which is responsible for the fairly large performance drop seen in this design. Response Functions 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 '&Respons e TTmax Xcc D TFmax TW4 Response Functions 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 '&Respons e TTmax Xcc D TFmax TW4 Response Functions0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 %&Response TTmax Xcc D TFmax TW4 Response Functions0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 %&Response TTmax Xcc D TFmax TW4 (a) (b) Figure 4.13: Variation of objectives with respect to HA or OA with other design variables fixed. (normalized values shown and D is the desirability function. All objectives are equally weighed in the composite function). A) TFmax, TW4, TTmax and Xcc, with respect to HA for =1, OA=0 and OPTT=0. B) TFmax, TW4, TTmax and Xcc, with respect to OA for =1, HA=1 and OPTT=0 Opt-Case 7 adds the final objective, TW4. Table 4.8 shows that addition of TW4 has no effect on the design and consequently no effect on any of the objectives. The chamber wall temperature, TW4, has previously been shown to be linked with TFmax, so this result is not surprising. The injector design resulting from inclusion of all the objectives in the optimization process (Opt-Case 7) is different from Opt-Case 4 where only the performance indicator,

PAGE 102

90 Xcc was optimized. In the design space evaluated for this effort, HA has a very strong effect on the system, especially for TFmax and TW4. Reference to Figure 4.13 shows that at the zero value for both HA and OA, Xcc is still decreasing. Performance could probably be further improved (i.e., Xcc shortened) by enlarging the design space to include smaller values of HA and OA. In the process of designing a real injector, the large drop in performance affected by the consideration of TFmax and TW4 would probably require a compromise design with certain variables weighted more heavily than others. To further probe the interplay between injector performance and component life/survivability, a composite optimization study is conducted with varying weights on the individual objective functions embodied by the composite desirability function. This is a functionality that the new design optimization methodology must possess. The designer must be able to weight, or favor, one or more dependent variables to have maximum flexibility in addressing design requirements. The use of geometric mean and desirability functions is an approach which addresses such a requirement. Additionally MOGA have also been shown to work well for such problems. The results of the MOGA will be presented in the later sections of this chapter. The results of the current study are shown in Table 4.9 and Figure 4.14. Reference to Table 4.9 (the numbers in parenthesis indicate the weight in the composite desirability function for that variable) shows that the designs suggested by Opt-Case 8 and Opt-Case 10 are similar. In Opt-Case 8 the life/survivability influencing objectives TFmax, TW4, and TTmax are weighted by a 2:1 ratio over the performance indicator (Xcc). Opt-Case 10 is the opposite, with performance having a 2:1 weighting over the life/survivability objectives. Reference to Tables 4.8 and

PAGE 103

91 4.9 shows that the design variables and objectives for these two cases are almost identical to those in Opt-Case 7, where all four objectives are weighted equally. Figure 4.13B shows the composite desirability function D, for Opt-Case 7, to be reasonably flat as a function of OA. Thus, small weightings have little or no effect on the injector design. Table 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 shown). OptCase HA OA OPTT TFmax TW4 TTmax Xcc 8 1 1 0.192 0 0.346 (1) 0.210 (1) 0.163 (1) 0.334 (0.5) CFD 0.343 0.200 0.185 0.326 Error (%) 0.04 0.12 0.91 0.22 9 0 1 0.497 0 0.0451 (5) 0.0823 (5) 0.471 (5) 0.627 (0.1) CFD 0.0401 0.0767 0.404 0.613 Error (%) 0.08 0.07 2.53 0.32 10 1 1 0 0 0.376 (0.5) 0.224 (0.5) 0.155 (0.5) 0.279 (1) CFD 0.369 0.214 0.252 0.264 Error (%) 0.10 0.12 3.93 0.38 11 0.612 0 0 0 0.919 (0.1) 0.898 (0.1) 0.282 (0.1) 0.0132 (5) CFD 0.911 0.890 0.361 0.0070 Error (%) 0.10 0.09 3.10 0.16 The other two cases shown in Table 4.9, Opt-Case 9 and Opt-Case 11, have relative weightings of 50:1 and 1:50, respectively for life/survivability (TFmax, TW4, and TTmax) and performance indicators (Xcc). For Opt-Case 9, the design tends toward the cases where TFmax and TW4 were individually minimized. The exception is that OPTT is equal to zero, which is required to minimize TTmax. For Opt-Case 11, the resulting design tends toward Opt-Case 4 where Xcc is minimized. Even with this strong weighting, the effect of

PAGE 104

92 including TFmax and TW4 is felt with the value of equal to 0.612. When Xcc is minimized by itself, is equal to 0.917. Optimization Study0 0.2 0.4 0.6 0.8 1 891011 Opt-CaseObjective TFmax TW4 TTmax Xcc Optimization Study 0 0.2 0.4 0.6 0.8 1 891011 Opt-CaseDesign Variable %& '& OPTT A B Figure 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 shown). A) Objectives. B) Design variables. Figure 4.15 shows the joint desirability function, D, for Opt-Case 9 (with a weighting of 50:1 in favor of life/survivability over performance), to be fairly flat. Accordingly, small weightings have little effect on the design. Significant weightings must be applied to push the design very far toward either life/survivability or performance. Opt-Case 11 has good performance and a low injector tip temperature. However, the chamber wall and injector face temperatures are very high. In an actual design, if this level of performance was required, active cooling could be required for the chamber wall and injector face. Additional CFD solutions were executed to confirm the optimum designs obtained from RSA for all 11 Opt-Cases. These results are shown in Tables 4.7, 4.8 and 4.9. The

PAGE 105

93 objective error, as defined by Eq. (4.5), is also shown for each case. The error is on the order of four percent or less in the case of TTmax. The error for the other three objectives is less than one percent. The larger discrepancy in TTmax is likely due to the steady state assumption made in the CFD analyses for this effort. The injector problem is actually unsteady. Response Functions 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 '&Respons e TTmax Xcc D TFmax TW4 Response Functions 0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 '&Respons e TTmax Xcc D TFmax TW4 Figure 4.15: Variation of TFmax, TW4, TTmax and Xcc, with respect to OA for =0, HA=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). Preliminary Observations Based on the preliminary optimization studies specific observations can be made. Single-objective optimization. The optimizer reliably locates single objective minimum values. Minimization of the individual objectives yields a different injector design. This establishes the fact that designing an efficient long life/survivability injector requires concurrent evaluation of all the relevant design variables. For the current study,

PAGE 106

94 1. The four objectives, based on the design obtained during their individual minimization, largely fall into two groups: (TFmax, TW4) and (TTmax, Xcc). However, there are differences between them. In each group, the responses of the two design objectives to OA are different. The observation also indicates that the three life/survivability-related objectives require compromises between design variables. 2. Minimizing TFmax and TW4 leads to a design with equal to zero (shear coaxial element), maximum fuel flow area and thickest post tip. This design also yields moderate to poor performance due to the slow mixing across the shear layer. 3. Minimizing TTmax and Xcc results in an impinging-like design with equal to one. It also has the minimum fuel flow area and the thinnest post tip thickness. This design performs well, but has very high wall and injector face temperatures. Multi-objective optimization. The injector design, when multiple objectives are included, is different from any of the single objective minimizations. For example, although the individual objectives, TTmax and Xcc drive the design towards the opposite ends of the range of OA, the composite desirability function accounting for all their effects exhibits only marginal variation with respect to it. The optimum design is affected mostly by other design variables. One can clearly see the benefit of using a global optimization procedure, which helps identify and interpret the trend of the composite function over the design space, especially when the individual objectives have opposing patterns. For the current study the following specific observations can be made. Equal weights 1. The design with all 4 objectives included (Opt-Case 7) is still an impinging element with low to moderate temperatures, but marginal performance. 2. Figure 4.13 indicates that lowering either propellant flow area beyond the current limit would probably decrease Xcc and thus increase performance. Lowering the oxidizer area would probably have less adverse effect on TFmax and TW4.

PAGE 107

95 Unequal weights 1. With modest weights on either the performance (Xcc) or life/survivability (TFmax, TW4 and TTmax) the design does not change appreciably. 2. Opt-Case 11 (with a large emphasis on performance) gives very good performance, modest TTmax, but very high values of TTmax and TW4. If enlarging the design space does not help, active cooling may be required. Sensitivity and Trade-Off Analysis The analysis is divided into two parts. Firstly, the global analysis is carried out where the entire design space is explored to estimate the sensitivity of the objectives to the design variables. Following this the POF is explored to understand the trade-offs between the objectives and also the different design trends. Global analysis The global analyses addresses, both the sensitivity of objective variability to design variables and the interaction between objectives over the complete design space. The global sensitivity indices are computed using the variance-based non-parametric approach proposed by Sobol (1993) and described in chapter 3. The code written using Matlab (2002) was validated using a well-known benchmark problem (Mckay, 1997). Table 4.10 summarizes the results of this study and lists the essential and non-essential design variables with respect to individual objective variability. A design variable is considered essential if it is responsible for at least 5% of the objective variability. Figure 4.16 shows the percentage of main factor ( Si) contribution of different design variables to individual objective variabilities. The variability of TFmax is largely influenced by HA and moderately by (Figure 4.16A). The effect of the other design variables is marginal suggesting that they are non-essential and in principle could be fixed. The variability of TW4 is considerably influenced only by HA (Figure 4.16B).

PAGE 108

96 TTmax is influenced considerably by OPTT and marginally by (Figure 4.16C). For Xcc, HA, OA and have considerable influence (Figure 4.16D). The total sensitivity indices ( Si total) were computed and compared with the main factors ( Si) and it was found that the contributions of the cross-interactions among the design variables to the objectives variability were negligible. Table 4.10: List of essential () and non-essential () 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. Design Variables Mean Error (%) Objective HA OA OPTT TFmax 5.2 TW4 11.5 Life/Survivability indicators TTmax 6.7 Performance indicator Xcc 6.1 A B Figure 4.16: Main factor (Si,) influence on objective variability. A) TFmax. B) TW4. C) TTmax. D) Xcc.

PAGE 109

97 C D Figure 4.16: Continued. This global sensitivity analysis shows that different design variables can have varied effect on individual objectives. Such a study can help the designer fix some of the design variables while analyzing the design. For example, in the current injector design study, for objective TFmax, the design variables OA and OPTT can be fixed at their mean value (0.5) as this do not result in significant differences in the prediction. The mean error (difference between predictions) throughout the design space between modified RSAs (obtained by fixing the non-essential design variables at their mean values), and the RSAs listed in Eqs.(4.1)-(4.4), (referred to as original RSA in the rest of the text) are given in Table 4.10. The error for TW4 is about 12% suggesting that the cutoff for the non-essential variables may have to be lowered to capture additional features of the original model. The mean error in the modified RSA of the remaining objectives is about 6%, suggesting that they capture the original RSA reasonably well. This information can be used for dimensionality reduction and therefore, to ease the search for

PAGE 110

98 optimum designs. Similar concept has been addressed in the past using a different approach. For example, Knill et al. (1999) have used linear aerodynamic to identify the important terms in the polynomial-based RSA which were then used for creating the surrogate model from Euler analyses. This reduced considerably the number of points needed for the Euler design of experiments. For the current work since the number of design variables is small the remaining studies are carried out without fixing the nonessential variables. From an engineer’s viewpoint and interest, this study provides an insight into the physics of a design problem by highlighting the design features that govern the individual objectives. Linear surrogate models for the four objectives are constructed as a function of the four design variables. The corresponding coefficients are shown in Table 4.11. The magnitude of the coefficients agrees well with the essential and non-essential nature of the design variables for each objective. Additionally, the R2 values for the linear RSA are compared with those of the original RSAs (referred as R2 nonlinear in Table 4.11). Comparing R2 linear with R2 nonlinear shows that most of the variability is accounted for by the linear model and the additional terms in the nonlinear model give marginal improvement. Table 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. HA OA OPTT Intercept R2 linear R2 nonlinear TFmax 0.237 -0.634 -0.033 -0.066 0.733 0.990 0.999 TW4 0.0719 -0.775 0.114 -0.052 0.826 0.986 0.999 TTmax -0.222 0.108 -0.063 0.662 0.367 0.918 0.994 Xcc -0.234 -0.402 0.433 0.108 0.138 0.958 0.997

PAGE 111

99 A correlation analysis was then carried out by Goel et al. (2004) to observe the interactions between objectives. A total of 14641 (114) design points were generated over the complete design space by varying one variable at a time by a constant value. The objectives were calculated using the original RSAs. The correlation matrix Cdes and corresponding p-values are computed using Matlab. The correlation matrix, Cdes shows that there is a strong correlation between objectives TFmax and TW4, as the corresponding coefficient is very close to 1. P-values and 95% confidence intervals for the correlation coefficients also establish the statistical significance of the results. Low P-values (<< 0.05) confirms the significance of the correlation results. This finding is in agreement with the observations made in the preliminary optimization study. The combustion chamber wall temperature, TW4, is excluded and the optimization problem is formulated with the remaining three objectives. Cdes = maxcc4max max cc 4 max TF X TW TT TF 1.00 -0.773 0.947 -0.272 X -0.773 1.00 -0.583 0.263 TW 0.947 -0.583 1.00 -0.193 TT -0.272 0.263 -0.193 1.00!" #$ #$ #$ #$ #$ #$ #$ %& Pareto front analyses The trade-offs between objectives and sensitivity analyses are carried out at the POF. The Pareto fronts of (TFmax, Xcc) and (TFmax, TTmax) are of interest as these objectives conflict one another. Goel et al. (2004) have generated the POFs with the aid of NSGA-II. Figure 4.17A shows the relation between TFmax and Xcc. It can be seen that the POF in this case is linear over a large region. A small increase in the value of TFmax ((10%) reduces the combustion length, Xcc, by nearly 50%. Figure 4.17B shows the

PAGE 112

100 relation between TFmax and TTmax. It is obvious that the POF is non-convex. It is also seen that a small drop in the value of the face temperature ((10%) can reduce the tip temperature TTmax by nearly 60%. Hence at a small cost of TFmax both Xcc and TTmax improves considerably. TFmax vs TTmax0 0.2 0.4 0.6 0.8 1 00.10.20.30.4 TFmaxTTmax TFmax vs Xcc0 0.2 0.4 0.6 0.8 1 00.20.40.60.81 TFmaxXcc A B Figure 4.17: Two-objective Pareto optimal front. A) TFmaxvs Xcc. B) TFmax vs TTmax. Following the trade-off studies, the three-objective (TFmax, Xcc and TTmax) optimization problem is solved using the NSGA-II algorithm (Goel et al., 2004). The solutions obtained are further refined using the –constraint local search strategy coded in Matlab. In this strategy, first, one of the objectives TFmax, is treated as the objective function and rest of the objectives Xcc and TTmax are treated as equality constraints. The constraint value is set equal to the corresponding objective value as found by NSGA-II simulation. The corresponding design variable vector is used as initial guess. This procedure is repeated for all individuals in the population. This gives a set of Pareto optimal solutions, referred to as Set A in this study. Next, Xcc is chosen as the objective and the other objective functions TTmax and TFmax are treated as constraints. For this

PAGE 113

101 problem, the obtained Pareto optimal set is referred to as Set B. Similarly, the Pareto optimal set, Set C, is obtained with TTmax as objective and TFmax and Xcc as constraints. Now all the Pareto optimal solutions obtained so far; i.e., Sets A, B, C and the original NSGA-II set are combined. To find the true Pareto optimal front, a non-domination check is carried out on this set of 400 solutions. This yields 254 non-dominated solutions. After removing the duplicates, there are 249 solutions, which are on the POF. These solutions are shown in Figure 4.18. This solution set is the global Pareto optimal solution set. The optimal solutions listed in Tables 4.7-4.9 are also plotted in Figure 4.18 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 -0.2 0 0.2 0.4 0.6 0.8 1 Xcc TFmax TTmax 1 2 9 8 6, 7, 10 11 4 5 3 Figure 4.18: Pareto optimal solution set (+) and the multi-objective optima listed in Tables 4.7-4.9 ( o ). The optimal solutions shown in Tables 4.7-4.9 lie along the POF obtained using MOGA and the local search algorithm. Although identical solutions are not identified by the different multi-objective optimization approaches, the trend of the Pareto front is adequately captured. This shows the effectiveness of the different approaches.

PAGE 114

102 Goel et al. (2004) has then used the hierarchical clustering algorithm (Jain and Dubes, 1988) to divide the obtained POF into 9 clusters for sensitivity and trade-off analyses. Values of the design variables and objectives for these designs are shown in Table 4.12. Graphically, the solutions are shown on the POF in Figure 4.19. It is clear from Figure 4.19 that solutions are selected uniformly over the design space. In Figure 4.20, the box plots of the design variables for clusters 1, 3, 6 and 9 are shown. The box plots for the objectives are shown in Figure 4.21. These plots highlight the variability of the design variables and objectives in each cluster. Table 4.12: Objective function and design variables of nine (9) representative designs from the Pareto optimal solution set. Cluster HA OA OPTT TFmax Xcc TTmax 1 0.000 1.000 0.842 0.712 0.023 1.090 0.880 2 0.000 1.000 0.356 0.587 0.028 0.749 0.0890 3 0.000 1.000 0.442 0.014 0.054 0.750 0.452 4 0.094 1.000 0.000 0.015 0.126 0.453 0.466 5 0.668 1.000 0.732 0.000 0.259 0.681 0.229 6 0.600 0.670 0.000 0.000 0.489 0.264 0.226 7 0.295 0.108 0.000 0.354 0.719 0.129 0.641 8 0.314 0.066 0.000 0.055 0.776 0.097 0.357 9 1.000 0.014 0.680 0.000 0.935 0.138 -0.043 For cluster 1 it is seen that the value of is fixed at 0 (shear coaxial injector) (Figure 4.20A) and HA is fixed at 1 (Figure 4.20B). This suggests that in this cluster, the designs are sensitive to and HA both of which reach their extreme values. The remaining two design variables OA (Figure 4.20C) and OPTT (Figure 4.20D), vary over a range. It is observed that TFmax is minimized (Figure 4.21A) where as Xcc and TTmax lie near their maximum and have little variability. Hence the designs in cluster 1 tend to minimize TFmax and represent shear coaxial injector designs. The design variables OA and OPTT do not influence TFmax but affect the remaining objectives, Xcc and TTmax. Partial correlation coefficients are estimated to obtain the relationship between the

PAGE 115

103 variances of these design variables and objectives (Table 4.13). It is noticed that as OA increases, Xcc increases (Rcorr = 1.00) and TTmax decreases (Rcorr = -0.638), thereby requiring a compromise in the design. As OPTT decreases both Xcc and TTmax decrease (Rcorr = 1.00 for both). Xcc TFmax TTmax 3 6 9 Figure 4.19: Pareto optimal solution set and nine (9) representative solutions from the same set. The circles identify the representative of a specific cluster. A B Figure 4.20: Box plots for the design variables in clusters 1, 3, 6 and 9. A) B) HA. C) OA. D) OPTT.

PAGE 116

104 C D Figure 4.20: Continued. Similar observations can be made for clusters 3, 6 and 9. Figures 4.21A-4.21C show that with increase in TFmax, from cluster 1 to cluster 9, Xcc and TTmax decrease (based on the median of the box plots). This highlights the trade-off between the objectives. Cluster 9 provides information about the opposing trend as to what was observed in cluster 1. As objectives Xcc and TTmax are minimized (Figures 4.21B and 4.21C) an impinging injector design is obtained ( ~ 1, Figure 4.20A) with TFmax exhibiting high values (Figure 4.21A). The HA is near minimum (Figure 4.20B) contrary to the design in cluster 1 where high HA minimized TFmax. The OA has considerable variation which suggests that the variability of objectives, Xcc and TTmax are not largely affected by this design variable. Figure 4.20D shows that Xcc and TTmax are minimized for the minimum value of OPTT. Table 4.13 gives the partial correlation coefficients for the set of design variables and the objectives in each cluster which shows considerable variation. The partial correlation coefficients for the combinations left out are effectively zero. These measures

PAGE 117

105 can be used to tune the design variables in a chosen cluster so as to improve on the objectives as per design requirements. A B C Figure 4.21: Box plots for the objectives in clusters 1, 3, 6 and 9. A) TFmax. B) Xcc. C) TTmax. In this study an elaborate optimization study with different compromised design goals has been presented along with both a global and POF sensitivity analyses. The conclusion drawn from these studies will be presented in chapter 5.

PAGE 118

106 Table 4.13: Partial correlation coefficients (Rcorr) of design variables vs. objectives for different clusters along the POF Cluster Design variable TFmax Xcc TTmax 1 OA 1.000 -0.638 1 OPTT 1.000 0.991 3 OA 1.000 6 0.982 -0.735 -0.983 6 HA -0.999 0.994 -0.729 9 0.877 -0.203 -0.769 9 HA -0.992 0.983 0.816 9 OA -0.977 0.997 -0.940

PAGE 119

107 CHAPTER 5 SUMMARY, CONCLUSION AND FUTURE WORK In this chapter the work done is this dissertation is summarized and conclusions drawn. Firstly, the summary and conclusion of the work related to Navier-Stokes (NS) CFD code verification is presented. Following this the design optimization study is summarized and conclusions drawn. Finally possible future research topics are briefly mentioned. Navier-Stokes Code Verification In this study, least square extrapolation (LSE) method has been mainly implemented on Navier-Stokes computations solved in the finite volume (FV) formulation. Two coarse grid solutions are extrapolated onto a fine grid using constant weighting parameters estimated by minimizing the residuals of the NS equations on the fine grid in the least square sense. The lid-driven cavity flow with two different boundary conditions (constant and variable lid velocities) are used to study the robustness and efficiency of LSE in situations of singularities, non-linearity and coupled equations. The initial test for LSE is done using a 2-D turning point problem is tested. Following this, only pressure is extrapolated for the lid-driven cavity flows. Finally LSE is implemented on the complete NS equations to extrapolate the velocity components and pressure simultaneously. The key observations pertaining to the LSE approach can be summarized as follows:

PAGE 120

108 It has been noticed with the aid of a 2-D turning point problem and a laminar liddriven cavity flow that this method can yield solution improvement for linear as well as nonlinear PDEs. Additionally it has been shown with this method that the extrapolating weights based on minimization of the residual of the governing equation performs better than the weights used by RE based on assumed order of convergence. It has been shown with the aid of two different Reynolds number (Re = 5.33 and 1000) that LSE works well over a wide range provided the base grids captures the flow features adequately. Minimizing the residuals does not always lead to minimizing the extrapolated solution error. In particular, it has been shown that flow with singularities would introduce large variation in the eigenvalues of the set of residual equation that is used for least square approach. This would influence the accuracy of the solution even though least square might effectively reduce the L2 norm of the residual. In practice, neglecting the regions of singularity can improve the performance. This suggests the requirement of further study to address issues of singularities in different flows. The LSE is seen to work well for the pressure-velocity coupled system. The Picardlike iteration adapted for the nonlinear momentum equation converges well for the shown case. It is also seen that the individual extrapolation of pressure or the coupled pressurevelocity extrapolation does not affect the accuracy of the extrapolated pressure field drastically. Overall issues related to singularities in the lid-driven cavity flow were noticed. This resulted in an elaborate study to understand the problems related to it and its influence on LSE. In this study only a constant weighting parameter is used. The preliminary implementation of spatially dependent weighting function exhibits inconsistent performance, which requires further investigation. Design Optimization Study In this study, the integration of CFD and optimization tools to explore the design of a single element rocket injector is presented. The design variables indicating the thermal environment and performance of the element are identified and the design points selected

PAGE 121

109 using a DOE technique. CFD solutions are obtained for each design and appropriate orders of polynomial-based RSAs are generated for the individual objectives. Based on the performance of the RSAs a grid refinement study is carried out and a new, more efficient grid is generated. Utilizing the results obtained from CFD computations on the fine grid, a preliminary multi-objective optimization study is carried out using a composite function based on geometric mean approach and different design trends, correlations and trade-offs observed. Following the preliminary study, RSA, global sensitivity indices and genetic algorithms have been used to conduct an elaborate sensitivity and trade-off analyses of the injector design. Global analyses have been conducted to observe the influence of design variables on the objective variability and estimate the correlation between objectives. A Pareto optimal front (POF), generated by Goel et al. (2004) using a RSAmulti-objective genetic algorithm (MOGA) coupled approach, is used to conduct tradeoff studies between objectives. A three-objective POF has been divided into 9 clusters and box plots have been used to observe the variability of the designs in each cluster. Additionally, partial correlation coefficients have been calculated to derive a relationship between the objectives and design variables in each cluster. Such analyses provide the designer with enough information to filter out designs based on his specific needs. Some observations based on the tools used, their performance, and the different design aspects of the injector are highlighted below. The degree of spread of a dependent variable over a parametric space is an important criterion for validation experiment design. Evaluation of large parametric spaces provides insight into potential validation measurements. Correlations of the variables highlight the trends and provide insight into physical processes that regulate the objectives.

PAGE 122

110 The study confirms that the details of the injector design govern the performance. Directing the fuel flow towards the element axis and thinning the oxidizer post lead to better performance. The same details are shown to have a large impact on the injector and combustion chamber environmental variables as well. Using CFD to evaluate and include more design variables and objectives in the conceptual design phase has potential to create more robust designs. This increase in the dimensionality of the design problem also creates the requirement for an optimization tool to guide the designer through often opposing trends to an acceptable design compromise. Based on the preliminary optimization study it was found that the proposed CFDRSM approach worked well. It was fairly efficient in terms of the number of CFD solutions required and provided sensible results in view of the physics of the problem. For both singleand multi-objective optimization efforts, in many cases, the resultant designs are at an extreme of the design space. This indicates the design space could be enlarged subject to practical design considerations. Again, a merit of the RSM is that the solutions obtained can be used repeatedly when revising the design scope. Furthermore, different scenarios, with either singleor multiobjective optimization, can make use of the information supplied by the response surface repeatedly. The global sensitivity indices (main factor and total sensitivity indices) help identify which design variable is essential for objective variability. Additionally, point out the effect of cross interactions among the design variables on the objectives variability. The variability of TFmax depends on and HA, TW4 on HA, TTmax on and OPTT and Xcc on HA and OA. The rest of the design variables for each objective can, in principle, be fixed, as their contribution to the objective variability is marginal. It is noticed that the mean errors between the modified RSAs and the original RSAs are about 6-12%. It is observed that TFmax and TW4 are strongly correlated and hence TW4 is excluded from the multi-objective optimization study. The conflicting nature of the objectives (TFmax, TTmax) and (TFmax, Xcc) is observed by examining the corresponding POFs. For a small cost in TFmax(~10%) both Xcc and TTmax (~50%) can be decreased considerably. The POF obtained for the 3-objective study (TFmax, TTmax and Xcc) gives an idea of the various compromises among the different objectives. A hierarchical clustering approach helps divide the POF into regions of similar objective goals.

PAGE 123

111 Box plots and partial correlation coefficients can assist the design selection along the POF. Box plots identify the variability of design variables and objectives in each cluster. Partial correlation coefficient identifies linear trends between the design variables and objectives variabilities. To minimize TFmax it is essential to have a shear coaxial injector ( = 0) and maximum change in baseline H2 flow area. To minimize Xcc and TTmax, an impinging injector ( > 0) is required with minimum change in baseline H2 flow area. This study has offered in depth information about the different design aspects pertaining to the design variables and objectives. The importance of different design variables on individual objectives has been identified and the interaction between the objectives noticed. The POF of the 3-objective study gives different choices for the designer to look into but the key observation is that for a marginal cost in the injector face temperature, TFmax, the performance, Xcc and O2 post tip temperature, TTmax can be considerably reduced. The MDO tools, namely DOE and RSM, have been shown to work well along with CFD. Vaidyanathan et al. (2003b) have shown that these tools can be use to address model fidelity issues in the context of a turbulent cavitating flows. In this work they have used these tools to identify the optimum values of the empirical parameters present in the turbulence and the cavitation models. This highlights the broad scope of the tools presented in this dissertation. Future work Future research directions are briefly presented in this section. First issues related to LSE are addressed and possible strategies suggested. Next, ideas for effective use of design tools used in this dissertation are presented.

PAGE 124

112 Navier-Stokes Code Verification In this dissertation, the weight function in LSE has not been allowed to vary spatially. For the problems tested and the techniques implemented, spatially dependent weight function was found to deteriorate the extrapolation and hence only a constant weight was used for all the study presented. One of the key issues during the extrapolation is the interpolation of coarse grid solutions onto the fine grid. This interpolation introduces high spurious wave number terms which combined with the spatially dependent weight function deteriorates LSE. For flows like lid-driven cavity flows, such high wave number components are amplified during residual estimation. The LSE approach is implemented on these polluted residuals and hence does not guarantee the minimization of the solution error. More details regarding this issue can be found in the work of Garbey and Shyy (2004). They have proposed using a low mode approximation of the residual for LSE implementation. This smoothing of the residual is shown to be effective in improving the performance of LSE. One of the other issues of LSE implementation is that the coarse grids should capture the flow feature adequately for efficient extrapolation. At the same time these coarse grids should not be so fine as to increase the computational cost. As the grid resolution is increased the solution reaches an asymptotic trend; i.e., there is no noticeable change in the flow feature with further refinement of the grid. Unless the coarse grids are properly chosen, the performance of LSE will not be optimal. Judicious selection of the base grids is an issue that deserves to be investigated. Design Optimization In this dissertation the RSAs were first generated and then used for the sensitivity and trade-off analyses. Hence although the global sensitivity analyses provided incite into

PAGE 125

113 the influence of design variables on the individual objective variability, it was not used for dimensionality reduction. Additionally the injector design problem has small number of design variables and hence there is not a serious requirement for such a reduction. For future design studies the following direction is suggested to tackle the curse of dimensionality. Whenever possible the preliminary design study can be conducted using a reduced order model, including semi-empirical and experimental data. This reduced model can be coupled with DOE, RSM and global sensitivity tools to identify the influence of design variables on objective variability. This information can then be used to reduce the dimensionality of the design space. Higher fidelity models can then be used to analyze the designs in this reduced design space. Additionally, the reduced design space will help in reducing the number of search directions during the optimization study.

PAGE 126

114 LIST OF REFERENCES Ainsworth M, Oden J T, 1992, “A Procedure for A Posteriori Error Estimation for h-p Finite Element Methods,” Computer Methods in Applied Mechanics and Engineering, 101, pp 73-96. Ainsworth M, Oden J T, 1993a, “ A Posteriori Error Estimators for Second Order Elliptic sSstems. Part 1: Theoretical Foundations and A Posteriori Error Analysis,” Computers Math. Applications, 25 (2), pp 101-113. Ainsworth M, Oden J T, 1993b, “ A Posteriori Error Estimators for Second Order Elliptic Systems. Part 2: An Optimal Order Process for Calculating Self-Equilibrating Fluxes,” Computers Math. Applications, 26 (9), pp 75-87. Ainsworth M, Oden J T, 2000, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York. Allen D M, August 1971, “Mean Square Error of Prediction as a Criterion for Selecting Variables,” Technometrics, 13 (3), pp 469-475. Allen D M, February 1974, “The Relationship between Variable Selection and Data Augmentation and a Method for Prediction,” Technometrics, 16(1), pp 125-127. Archer G E B, Saltelli A, Sobol I M, 1997, “Sensitivity Measures, ANOVA-like Techniques and the USA of Bootstrap,” Journal of Statistical Computation and Simulation, Vol. 58, pp 99-120. Babuska I, Rheinboldt W C, 1978, “ A Posteriori Error Estimates for the Finite Element Method,” International Journal of Numerical Methods in Engineering, 12, pp 15971615. Barth T-J, January 1993, “Recent Developments in High Order K-exact Reconstruction on Unstructured Meshes,” AIAA Paper. Calhoon D F, Ito J I, Kors D L, July 1973, “Investigation of Gaseous propellant Combustion and Associated Injector-Chamber Design Guidelines,” Aerojet liquid rocket company, NASA Cr-121234, Contract NAS3-13379. Chen Y S, January 1989, “Compressible and Incompressible Flow Computation with a Pressure-Based Method,” AIAA89-0286, AIAA 28th Aerospace Sciences Meeting.

PAGE 127

115 Chen Y S, Farmer R C, June 1991, “CFD Analysis of Baffle Flame Stabilization,” AIAA 91-1967, AIAA 27th Joint Propulsion Conference. Cheng G C, Farmer R C, January 2002, “CFD Spray Combustion Model for Liquid Rocket Engine Injector Analyses,” AIAA Paper 2002-0785, 40th AIAA Aerospace Sciences Meeting & Exhibit. Dean E B, 1996, “Taguchi Methods from the Perspective of Competitive Advantage,” http://akao.larc.nasa.gov/dec/tm.html. Deb K, 2001, Multi-Objective Optimization Using Evolutionary Algorithms, Wiley, Chichester UK. Deb K, Goel T, 2000, “Multi-objective Evolutionary Algorithms for Engineering Shape Design,” Evolutionary optimization (Eds. Sarker R, Mohammadin M, Yao X), Kluwer, pp 147 – 176. Deb K, Agrawal S, Pratap A, Meyarivan T, 2000, “A Fast and Elitist Multi-objective Genetic Algorithm for Multi-objective Optimization: NSGA-II,” Proceedings of the parallel problem solving from Nature VI conference, Springer-Verlag, Paris, France, pp 849 – 858. Deb K, Goel T, 2001, “A Hybrid Multi-Objective Evolutionary Approach to Engineering Shape Design,” Proceedings of Evolutionary Multi-criterion Optimization Conference, Springer-Verlag, Zurich, pp 385 – 399. Demkowicz L, Oden J T, Strouboulis T, 1984 “Adaptive Finite Elements for Flow Problems with Moving Boundaries. Part 1: Variational Principles and A Posteriori Error Estimates,” Comput. Methods Appl. Mech. Engg., 46, pp 217-251. Dickerson R, Tate K, Nurick W, January 1968, “Correlation of Spray Injector Parameters with Rocket Engine Performance,” AFRPL-TR-68-11, U. S. Air Force Rocket Propulsion Laboratory. Ferziger J H, Peric M, 1999, Computational Methods for Fluid Dynamics, Springer, 2nd edition, New York. Garbey M, Shyy W, 2003, “A Least Square Extrapolation Method for Improving Solution Accuracy of PDE Computations,” Journal of Computational Physics, 186, pp 1-23. Garbey M, Shyy W, 2004, “A Least Square Extrapolation Method for the A Posteriori Error Estimate of the Incompressible Navier Stokes Problem,” Department of Computer Science, University of Houston, TX (unpublished). Gill D S, Nurick W H, March 1976, “Liquid Rocket Engine Injectors,” NASA SP-8089, NASA Lewis Research Center.

PAGE 128

116 Goel T, 2001, "Optimal Shape Design of Mechanical Components Using Single and Multi-Objective Genetic Algorithms", Master's Thesis, Dept. of Mechanical Engg., Indian Institute of Technology, Kanpur, India. Goel T, Deb T, November 2002, "Hybrid Methods for Multi-Objective Evolutionary Algorithms", SEAL'02, Proceedings of the 4th Asia-Pacific Conference on Simulated Evolution And Learning Computational Intelligence for the E-Age, Singapore, pp 188-192. Goel T, Vaidyanathan R, Haftka R T, Queipo N, Shyy W, Tucker K, 2004, “Response Surface Approximation of Global Pareto Optimal Front in Multi-Objective Optimization of Simulation-Based Models,” 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Albany NY. Gu X, Dchock H J, Shih T I-P, Hernandez E C, CHu D, Keller P S, Sun R L, 2001, “Grid-Quality Measures for Structured and Unstructured Meshes,” AIAA 20010652. Haftka R T, Scott E P, Cruz J R, 1998, "Optimization and Experiments: A Survey," Applied Mechanics Reviews, 51(7), pp 435-448. Ilinca C, Zhang X D, Trepanier J-Y, Camarero R, 2000, “A Comparison of Three Error Estimation Techniques for Finite-Volume Solutions of Compressible Flows,” Comput. Methods Appl. Mech. Engg., 189, pp 1277-1294. Ishibuchi H, Yoshida T, 2002, “Hybrid Evolutionary Multi-Objective Optimization Algorithms”, Soft Computing Systems: Design, Management and Applications (Frontiers in Artificial Intelligence and Applications), 87, pp 163 – 172. Jain A K, Dubes R C, 1988, Algorithms for Clustering Data, Prentice Hall College, New Jersey. Jasak H, Gosman, A D, 2001, “Residual Error Estimate for the Finite Volume Method,” International Journal for Numerical methods for Fluids, 39, pp 1-19. JMPTM, The Statistical Discovery SoftwareTM, Version 5, Copyright ) 1989-2002, SAS Institute Inc., Cary, NC, USA. Knill D L, Giunta A A, Baker C A, Grossman B, Mason W H, Haftka R T, Watson L T, 1999, “Response Surface Methods Combining Linear and Euler Aerodynamics for Supersonic Transport Design,” Journal of Aircraft, 36(1), pp 75-86. Lasdon L S, Waren A, Jain A, Ratner M., March 1978, “Design and Testing of a Generalized Reduced Gradient Code for Nonlinear Programming,” ACM Transactions on Mathematical Software, 4(1), pp 34-50. MATLAB, The Language of Technical computing, Version 6.5 Release 13. 19842002, The MathWorks, Inc.

PAGE 129

117 McKay M D, March 1997, “Nonparametric Variant based Methods of Assessing Uncertainty Importance,” Reliability and System Safety, 57, pp 267-279. Miettinen K M, 1999, Nonlinear Multiobjective Optimization. Kluwer, Boston. Mullin M, Sukthankar R, 2000, "Complete Cross-Validation for Nearest Neighbor Classifiers," 17th International Conference on Machine Learning (ICML), Stanford University, Palo Alto, California. Myers R H, Montgomery D C, 1995, Response Surface Methodology – Process and Product Optimization Using Designed Experiment, John Wiley & Sons, New York. Nurick J H, November 1990, “DROPMIX A PC Based Program for Rocket Engine Injector Design,” JANNAF Propulsion Conference, Cheyenne, Wyoming. Oberkampf W L, Trucano T G, 2002, “Verification and Validation in Computational Fluid Dynamics,” Progress in Aerospace Sciences, 38, pp 209-272. Oden J T, 2001, “ A Posteriori Error Estimation,” Verification and Validation in Computational Solid Mechanics, Edited by Hans Mair, ASME/USACM Standards. Oden J T, Wu W, Ainsworth M, 1994, “An A Posteriori Error Estimate for Finite Element Approximations of the Navier-Stokes Equations,” Computer Methods in Applied Mechanics and Engineering, 111, pp 185-202. Ollivier-gooch C F, 1996, “A New Class of ENO Schemes based on Unlimited DataDependent Least-squares Reconstruction,” AIAA 34th Aerospace Sciences Meeting and Exhibit, Reno, NV, US, AIAA-96-0887. Owen A, 1992, "Orthogonal Arrays for: Computer Experiments, Integration and Visualization," Statistica Sinica, 2(2), pp 439-452. Patankar S V, 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, Taylor & Francis Group, New York. Papila N U, 2001, Neural Network and Polynomial-Based Response Surface Techniques or Supersonic Turbine Design Optimization, PhD Dissertation, Department of Mechanical and Aerospace Engineering, University of Florida. Pavli A L, September 1979, “Design and Evaluation of High Performance Rocket Engine Injectors for Use with Hydrocarbon Fuels,” NASA TM 79319. Pelletier D, Turgeon E, Lacassse D, Borggaard J, 2001, “Adaptivity, Sensitivity and Uncertainty Towards Sandards in CFD,” AIAA 2001-0192. Pieper J L, October 1991, “Oxygen/Hydrocarbon Injector Characterization,” PL-TR 913029, Philips Laboratory, Propulsion Directorate, Air Force Systems Command, Edwards Air Force Base, CA.

PAGE 130

118 Roache P J, 1994, “Perspective: A Method for Uniform Reporting of Grid Refinement Studies,” Journal of Fluids Engineering, 116, pp 405-413. Roache P J, 1997, “Quantification of Uncertainty in Computational Fluid Dynamics,” Annual Review of Fluid Mechanics, 29, pp 123-160. Roache P J, 1998, Verifiction and Validation in Computational Science and Engineering. Hermosa Publishers, Albuquerque, NM. Rupe J H, July 1965, “An Experimental Correlation of the Nonreactive Properties of Injection Schemes and Combustion Effects in a Liquid Rocket Engine,” NASA TR 32-255. Shyy W, 1994, Computational Modeling for Fluid Flow and Interfacial Transport, Elsevier Science Publishers, Amsterdam, Netherlands. Shyy W, Thakur S S, Ouyang H, Liu J, Blosch E, 1997, Computational Techniques for Complex Transport Phenomena, Cambridge University Press, New York. Shyy W, Garbey M, Appukuttan A, Wu J, 2002, “Evaluation of Richardson Extrapolation in Computational Fluid Dynamics,” Numerical Heat Transfer, Part b, 41, pp 1-27. Sobol I M, 1993, “Sensitivity Estimates for Nonlinear Mathematical Models,” Mathematical Modeling and Computational Experiment, New York, NY, John Wiley & Sons, 1, pp 407-414. Solver, 2000, Microsoft Corporation, 1985-1999, Microsoft Excel. Thakur S, Wright J, Shyy W, 2002, “STREAM: A Computational Fluid Dynamics and Heat Transfer Navier-Stokes Solver. Theory and Applications,” Streamline Numerics, Inc., and Computational Thermo-Fluids Laboratory, Department of Mechanical and Aerospace Engineering Technical Report, Gainesville, Florida. Tucker P K, Shyy W, Sloan J G, July 1998, “An Integrated Design/Optimization Methodology for Rocket Engine Injectors,” 34th AIAA / ASME / SAE / ASEE Joint Propulsion conference and Exhibit. Tucker P K, Shyy W, Vaidyanathan R, 2000, “An Optimization Based Approach to Injector Element Design,” AIAA-2000-3220. Turgeon E, Pelletier D, Borggaard J, January 2001, “Sensitivity and Uncertainty Analysis for Variable Property Flows,” AIAA-2001-0139, 39th AIAA Aerospace Science Meeting and Exhibit, Reno, NV. Unal R, Dean E B, 1991, "Taguchi Approach to Design Optimization for Quality and Cost: An Overview," Proceedings of the International Society of Parametric Analysts, 13th Annual Conference, ISPA, New Orleans, Louisiana.

PAGE 131

119 Unal R, Dean E B, 1995, "Design for Cost and Quality: The Robust Design Approach," http://mijuno.larc.nasa.gov/pap/robdes/robdes.html. Unal R, Braun R D, Moore AA, Lepsch R A, 1997, “Response Surface Model Building Using Orthogonal Arrays for Computer Experiments,” 19th Annual International Conference of the International Society of Parametric Analysis, New Orleans, Louisiana. Unal R, Lepsch R A, McMillin M L, 1998, “Response Surface Model Building and Multidisciplinary Optimization using D-Optimal Designs,” 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Paper No. 98-4759. Vaidyanathan R, Papila N, Shyy W, Tucker P K, Griffin L W, Haftka R T, Fitz-Coy N, September 2000, “Neural Network and Response Surface Methodology for Rocket Engine component Optimization,” 8th AIAA/NASA/USAF/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Long Beach, CA. Paper No. AIAA 2000-4880. Vaidyanathan R, Tucker P K, Papila N, Shyy W, January 2003a, “CFD-Based Optimization for a Single Element Rocket Injector,” 41st Aerospace Sciences Meeting & Exhibit, Paper No. 2003-0296 (in review for Journal of Fluids Engineering). Vaidyanathan R, Senocak I, Wu J, Shyy W, May 2003b, “Sensitivity Evaluation of a Transport-Based Turbulent Cavitation Model,” Journal of Fluids Engineering, 125(3), pp 447-458. Vaidyanathan R, Shyy W, Garbey M, Haftka R T, January 2004a, “CFD Code Verification Using Least Square Extrapolation Method,” 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, Paper No. AIAA-2004-0739. Vaidyanathan R, Goel T, Shyy W, Haftka R T, Queipo N, Tucker P K, July 2004b, “Global Sensitivity and Trade-Off Analyses for Multi-Objective Liquid Rocket Injector Design,” 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Fort Lauderdale, FL, Paper No. AIAA-2004-4007. Van Doormal J P, Raithby G D, 1984, “Enhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows,” Numerical Heat Transfer, 7, pp 147-163. Wang T S, Chen Y S, July 1990, “A United Navier-Stokes Flowfield and performance Analysis of Liquid Rocket Engines,” AIAA 90-2494, AIAA 26th Joint Propulsion Conference. Zhang X D, Trepanier J-Y, Camarero R, 1997, “An A Posteriori Error Estimation Method Based on Error Equations,” Proceedings of the AIAA 13th Computational Fluid Dynamics Conference, Snowmass, CO, US, AIAA-97-1889.

PAGE 132

120 Zienkiewicz O C, Zhu J Z, 1987, “A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis,” International Journal for Numerical Methods in Engineering, 24, pp 337-357. Zienkiewicz O C, Zhu J Z, 1992a, “The Superconvergent Patch Recovery and A Posteriori Error Estimates. Part 1: The Recovery Technique,” International Journal for Numerical Methods in Engineering, 33, pp 1331-1364. Zienkiewicz O C, Zhu J Z, 1992b, “The Superconvergent Patch Recovery and A Posteriori Error Estimates. Part 2: Error Estimates and Adaptivity,” International Journal for Numerical Methods in Engineering, 33, pp 1365-1382.

PAGE 133

121 BIOGRAPHICAL SKETCH Rajkumar Vaidyanathan was born on April 27, 1976, in Chennai, India. He spent the major part of his life and schooling in Kolkata, India. He obtained his bachelor’s degree in naval architecture from Indian Institute of Technology, Chennai, in 1998 and master’s degree in aerospace engineering at the University of Florida in 2000. He has, since then, been pursuing his doctoral studies under the guidance of Dr. Wei Shyy. His areas of interest include multidisciplinary optimization techniques and computational fluid dynamics code verification.


Permanent Link: http://ufdc.ufl.edu/UFE0005362/00001

Material Information

Title: Investigation of Navier-Stokes Code Verification and Design Optimization
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0005362:00001

Permanent Link: http://ufdc.ufl.edu/UFE0005362/00001

Material Information

Title: Investigation of Navier-Stokes Code Verification and Design Optimization
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0005362:00001


This item has the following downloads:


Full Text












INVESTIGATION OF NAVIER-STOKES 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

Fitz-Coy 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 NAVIER-STOKES 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 pressure-velocity................................ ... ..................24
Results and Discussion.................................................25
Two-Dimensional Turning Point Problem ............................................. 26
Lam inar Lid-Driven Cavity Flow ........................................ ............... 29
LSE for pressure only .............................. ........................... 33
LSE for coupled pressure-velocity 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 ulti-Objective 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
Single-objective optimization ............................................................82
M ulti-objective optim ization...................................... ............... 86
Prelim inary O observations ........................................ .......................... 93
Sensitivity and Trade-Off Analysis ..................................... .....................95
Global analysis ...............................................95
Pareto front analyses ..................................................... 99

5 SUMMARY, CONCLUSION AND FUTURE WORK ......... ................................107

N avier-Stokes Code V erification ........................................ ....................... 107
D esign Optim ization Study ............................................................................. 108
Future w ork ............................................. ............. ............... .111
Navier-Stokes 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 (non-normalized values shown). .................................. .................74

4.5: Performance of full and reduced quadratic RSAs for the four objective functions
(non-norm 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 non-essential (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, x-axis is the number of grid points
in each direction for the fine grid on which extrapolation is done...........................28

2.5: Two dimensional lid-driven cavity flow (dimensions 1 x 1). ...................................30

2.6: u-velocity contour for grid 121 x 121 and Re = 5.33. ..............................................31

2.7: u-velocity 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 lid-velocity. ........................................ .......................... 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 lid-velocity. ............................................35

2.12: Comparison of interpolated and extrapolated pressure with CFD solutions for Re =
1000 and variable lid-velocity ......... .............................................. ............... 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 lid-velocity (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 lid-velocity (reduced domain). ..................38









2.16: Comparison of interpolated and extrapolated pressure with CFD solutions for Re =
5.33 with constant lid-velocity. ............................. ...... ........................... 39

2.17: Comparison of interpolated and extrapolated pressure with CFD solutions for Re =
1000 with constant lid-velocity. ................................. ............... 40

2.18: Extrapolated pressure and velocity compared with CFD solution on grid 121 x 121
(M,), Re = 5.33 and variable lid-velocity. .........................................42

2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 x 171
(M,), Re = 1000 and variable lid-velocity.....................................43

3.1: Schematic of the procedure for global design optimization (Papila, 2001).............46

3.2: Full factorial 3-level 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: Near-injector 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: Two-objective Pareto optimal front ............................ ..... ........... .... 100

4.18: Pareto optimal solution set (+) and the multi-objective 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 NAVIER-STOKES 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 Navier-Stokes 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 CFD-based

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 model-based optimization

and sensitivity evaluation.

For Navier-Stokes (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 CFD-based 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 k-eturbulence closure. With the aid of these surrogate models,

sensitivity and trade-off 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 multi-objective optimization

study is carried out using a geometric mean approach. Following this, sensitivity analyses

with the aid of variance-based non-parametric approach and partial correlation

coefficients are conducted using data available from surrogate models of the objectives

and the multi-objective 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 Navier-Stokes 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 multi-disciplinary optimization (MDO) which

in turn promotes CFD-based design optimization studies.

Computational fluid dynamics-based 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 Navier-Stokes (NS) equations

and other related models introduces challenges that include physical modeling

uncertainty, geometric complexities, non-uniform and non-orthogonal 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

round-off 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 well-resolved 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 multi-objective problems. Multi-objective 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 non-dominated

solutions in the Pareto-optimal set is termed the Pareto-optimal front (POF). This set

provides the designer with a clear idea of the trade-offs 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 multi-objective

evolutionary algorithms (MOEA) have been proposed in literature (Deb, 2001). Insight

into the trade-offs between different objectives can be obtained from the POF.

In this dissertation different aspects of NS code verification and CFD-based 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 lid-driven 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 CFD-MDO 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 trade-offs 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.

* CFD-based 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 CFD-based 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 multi-element 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 trade-off analyses for the design are then carried out using the data

obtained from surrogate models of the design objectives, and the POF (multi-objective

optima) generated by Goel et al. (2004) with the aid of a multi-objective genetic

algorithm (MOGA) and a local search method (e -constraint strategy). The regions of the

POF that represent different trade-offs 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
variance-based non-parametric 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 two-dimensional turning point problem and a laminar

lid-driven 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
NAVIER-STOKES 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

Navier-Stokes (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 two-part 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; Ollivier-gooch, 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 flow-field output, say a velocity component or pressure, of a CFD

simulation on a grid size, h, and assuming second-order 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 higher-order terms leads to


U(x)= I4u x; h-u(x;h)] +O(h3) (2.3)


This is a second-order RE which results in a gain of order one. Similarly if the

solution is assumed to have first-order convergence apriori and O(h) term is eliminated,

it results in a first-order RE with a gain of order one.

U(x) = 2ux; h-u(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)
(4-h2)

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,+ 1-a 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((j-1) xr) (2.8)
J=1

with x,j = 0,... ,Mreals. The weighting function, ais dependent on the spatial co-

ordinate, x for a one-dimensional problem. For two dimensional problem it is a function

of spatial co-ordinates, 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, +(l-a)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 one-dimensional problem using two Fourier modes, Eq.

(2.12) that has to be minimized can be rewritten as










Lh [a +(1-ao)2 ] +Lh[a cos (xr) +(l-a cos (xr))2]-fh (2.13)

This approach can be generalized for nonlinear PDE problems via a Newton-like 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 n-times 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 non-exact 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 I||a1l = 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 non-exact 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 non-exact 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 (non-exact 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 in-house code STREAM which is

based on the pressure-based 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 second-order 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 v-component 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 non-linearity in the momentum equations

(a's have velocity components in them which make the equation non-linear) and the

coupling of pressure-velocity 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

Navier-Stokes problem is tackled in two steps. In the first step, only the pressure is

extrapolated and in the second step the pressure-velocity 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 Navier-Stokes 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 + (1-a)2

* InputpLsE, f2, v2 into Eq. (2.31).

* Find that minimizes the residual of Eq. (2.31).

* Calculate PLE.

LSE for coupled pressure-velocity

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 Picard-like

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, + (l-a) 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 =Y-7; + (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 + (1-f3)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, two-dimensional turning point problem (Figure 2.2) and

* Laminar, lid-driven cavity flow (Figure 2.5).

1. variable lid-velocity

2. constant lid-velocity











Two-Dimensional 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 -y-y


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(n-y) 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 two-dimensional. A sample solution is shown in Figure 2.3.

In the finite volume implementation, central difference is used for the diffusion term and









first-order 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 first-order 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 ofNavier-Stokes 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 second-order, 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, x-axis 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 under-resolved 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, lid-driven cavity flow is used and the

details are given in the following section.

Laminar Lid-Driven Cavity Flow

Garbey and Shyy (2003) have addressed the lid-driven cavity flow problem in a

vorticity-streamfunction formulation with a finite difference implementation. In the

current study the complete NS equations are solved in 2D for the lid-driven cavity flow

and the pressure along with u and v-components 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 lid-velocity uhd in the positive x-direction

(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 u-component of the velocity at the lid. The v-component 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 u-velocity 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 lid-driven cavity flow (dimensions 1 x 1).

The Reynolds number for the flow is defined based on the mean lid-velocity. The

u-velocity contours are shown in Figures 2.6 and 2.7 for Reynolds numbers 5.33 and

1000, respectively. Some difference is noticed in the u-velocity contours at the upper

corner regions of the cavity for both the Reynolds numbers. The negative variable uhd

results in a u-velocity 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: u-velocity contour for grid 121 x 121 and Re = 5.33. A) ud = -16x2(l-x)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: u-velocity 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(l-x)


-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(1-x)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 Pressure-Velocity (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 lid-velocity

is presented and in part (b) the case with constant lid-velocity 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 lid-velocity (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


3-4 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 lid-velocity. 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 lid-velocity (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 lid-velocity. 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 lid-velocity 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 lid-velocity. 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 lid-velocity (Re = 5.33). At this point it would be interesting to take a

look at a case where the lid-velocity 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 = 0-1 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 lid-velocity (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

S-6.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 of-4.0. This suggests

that a sub-domain 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 lid-velocity (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
U-2.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 lid-velocity. 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 lid-velocity (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.2-3.2
E -3.4 EMORE2 EM 111
EM0121 I-3.4 O--- -- EMn111
S-3.6 3 Mn

-3. EM E -3. 6 LSE

-4-1 -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 lid-velocity. 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 pressure-velocity computations

The coupled algorithm is implement on a lid-driven cavity flow with Re = 5.33 and

1000 and the variable lid-velocity 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 3-4

iterations of this algorithm results in converged values of a yand f with about 3-4 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, u-velocity and v-velocity 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 u-velocity 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
u-velocity


80 90 100 110
Grid Points (N)


Emn norm error vs N
Re = 5.33
v-velocity
-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 lid-velocity. A) Pressure. B) u-velocity. C)
v-velocity.


-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
0g-3.0
o
-3.1


-3.2


-3.3


Emn vs N
Re = 1000
u-velocity


130 140 150 160 170 180 130 140 150 160

Grid Points (N) Grid Points (N)


170 180


Emn vs N
Re = 1000
v-velocity


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 lid-velocity. A) Pressure. B) u-velocity.

C) v-velocity.
















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 multi-disciplinary 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 single-objective optimization study

which is dependent on two design variables. This is a general case and can be extended

for multi-objective 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 polynomial-based 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 second-order 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 k-1 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)
n-n


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 rms-error 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 cross-validation 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 k-fold cross-

validation (Mullin and Sukthankar, 2000). In k-fold cross-validation, the data points are

divided in k equally sized subsets. One of the subset is kept for testing and the remaining

k-1 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 k-1 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

rms-error. 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 rms-error 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 rms-eror, 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, (n-n ) n- (3.1
R=I-1- =- (I-R (3.10)
SS, /(n 1) n-nP


For a good fit, the value of R2 should be closer to one.

The polynomial-based 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

polynomial-based 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 3-level 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...q-1 }. 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 nr-row-by-t-column sub-matrix, 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 n-dimensional 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 variance-based non-parametric 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 multi-objective

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 variance-based

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 gradient-based 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)
B-A

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)
SC-E

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
s-8

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



Multi-Objective Genetic Algorithm

One of the other optimization approaches is to use Multi-objective 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

multi-objective 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

hyper-surfaces.

One of the MOGA that has been shown to be effective in finding the Pareto-

optimal solutions is the elitist non-dominated sorting genetic algorithm (NSGA-II) (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 non-domination 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 y-axis 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 "H-spread". The inner fence is located 1 step beyond the hinges, which is

equal to 1.5 times the H-spread. 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

empirical-based (Rupe, 1965; Dickerson et al., 1968; Pavli, 1979; Nurick, 1990; Pieper,










1991). Extensive sub- and full-scale hot-fire 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 one-dimensional 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 cold-flow and hot-fire 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 half-angle, 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 three-dimensional geometry of multi-element 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 CFD-based

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 trade-off 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 multi-objective optima (Pareto optimal front,

POF) generated by Goel et al. (2004) with the aid of multi-objective genetic algorithm

(MOGA) and a local search method. The regions of the POF that represent different

trade-offs 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 trade-off 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 variance-based

non-parametric 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 head-on interaction between the oxidizer and fuel (Gill

and Nurick, 1976). The second type of injector consists of non-impinging elements where

the propellant streams flow in parallel, typically in coaxial fashion (Figure 4.1B). Here,

mixing is accomplished through a shear-mixing 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 F-O-F

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 non-impinging 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) F-O-
F impinging element. B) Coaxial element.

One can combine the above-mentioned 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 F-O-F

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 Baseline-40% 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 pressure-based, finite difference, Navier-Stokes solver, FDNS500-CVS (Chen,

1989; Wang and Chen, 1990; Chen and Farmer, 1991) is used in this study. The Navier-

Stokes equations, the two-equation 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 7-species and 9-reaction 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 cross-validation 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 polynomial-based 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 non-normalized 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 (non-normalized 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: Near-injector 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 injector-generated 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 axi-symmetric

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 j-line 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 (non-normalized 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).


TFmax-Quadratic RS


1


0.8


0.6


0.4


0.2


0 0.2 0.4 0.6 0.8
RS


TTmax-Quadratic RS


0 0.2 0.4 0.6
RS


TW4-Quadratic RS












RS-CFD
Optimum


0 0.2 0.4 0.6 0.8 1
RS


Xcc-Quadratic 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. RS-CFD 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 4-design 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
(non-normalized 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, multi-objective optimizations are

performed with equal weights. Finally, the multi-objective 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. Non-normalized values of the

objectives refer to the scaled values and the normalized values refer to the re-scaled

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.

Single-objective 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 Opt-Cases 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 moderate-to-long 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
Opt-Case Opt-Case

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 Opt-Cases 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 impinging-like 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

multi-objective 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

single-objective 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.

Multi-objective optimization

To concurrently evaluate component life/survivability and performance

considerations, a multi-objective 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 Opt-Case 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 Opt-Case 5, the optimum design does

not change significantly from that in Opt-Case 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 Opt-Case 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 single-objective

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 --- a-------|1 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
Opt-Case Opt-Case

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.

Opt-Case 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 Opt-Case 5. This shift in AHA is consistent with the single-objective

minimization of TFmax shown as Opt-Case 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.