UFDC Home  myUFDC Home  Help 



Full Text  
NUMERICAL MODELING OF TRANSPIRATION COOLED ROCKET INJECTORS By LANDON ROTHWELL TULLY A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2005 Copyright 2005 by Landon Rothwell Tully In loving memory of Randolph R. Tully, Sr. ACKNOWLEDGMENTS This research was made possible by the contributions of numerous professional colleagues, friends, and family. Their combined expertise, encouragement, and support have kept me on track, engaged, and committed to high quality research. I am forever thankful for all their efforts and for their belief in my work. Dr. Jacob N. Chung was unfailing in his wisdom and guidance, without which this research would not have been possible. Dr. Bruce F. Carroll and doctoral candidate Ahmed F. Omar contributed their support and parallel experimental efforts which were vital to the outcome of this research. Dr. Siddharth S. Thakur had no direct ties to this work but he graciously gave up his time to offer his suggestions when I asked. His teachings and help were invaluable. Dr. Wei Shyy agreed to invest considerable time and effort to be on my supervisory committee and is responsible for the continued success and quality of the Mechanical and Aerospace Engineering Department at the University of Florida. This study would not have been possible without the contributions of NASA's Institute for Future Space Transport, the University Research Engineering and Technology Institute, and the Constellation University Institute Program. Additionally the guidance of Paul Kevin Tucker from NASA's Marshall Flight Space Center was imperative in the progress of this research. My parents, brother, and extended family were always on hand with support and encouragement. My friends have given me a lifetime of memories and helped make my time at the University of Florida an enjoyable one. Ryan Avenall has worked through the courses with me and helped critique my research. Finally, I am grateful to Marci Gluckson for her constant encouragement, reassurance, and inspiration in my work. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ......... .................................................................................... iv L IS T O F T A B L E S ........ .. ................................................................. .... .. .. v iii LIST OF FIGURES ......... ........................................... ............ ix N O M E N C L A T U R E ......................................................................................................... x ii A B S T R A C T .........x.................................... ....................... ................. xv CHAPTER 1 INTRODUCTION LIQUID PROPELLANT ROCKET ENGINES..........................1 1.1 L iquid Propellant R ocket H history .........................................................................1 1.2 Liquid Propellant Rocket Engine Components ............................................. 4 1 .2 .1 Injecto rs .................. ... ............ ................... ........................ .. 4 1.2.2 Combustion Chamber and Nozzle............................................................8 1.3 Objectives of the Study..................................................................... 10 2 HISTORY OF POROUS MODELING .............................................. ...............12 2 .1 Initial M odeling E efforts .......................................................................... ... ... 12 2.2 Perm ability and Porosity........................................................ ............. 14 2 .2 .1 P oro sity ................................................................................. 14 2.2.2 P erm ability ......... .......................................................... .. .... .... ..... 14 2.2.3 Perm ability M odels.......................................... ............................ 16 2.2.3.1 C apillary m odels ........................................ ......... ............... 17 2.2.3.2 D rag m odels ............................ ...... .. .......... .. .......... ....18 2.2.3.3 Experim mental m ethods................................. ....................... 18 2.3 H igh R eynolds N um ber Flow .................................... ........................... ......... 19 2.3.1 Forchheim er m odels ............................................................................ 21 2.3.2 Experim ental m ethods ........................................ .......................... 22 2.4 Sem i H euristic E quations ................... ... ..................... ............... .... 22 2.4.1 Semiheuristic Continuity and Momentum Equations.............................. 23 2 .4 .2 E n ergy E qu ation ............................................................. .....................2 5 3 FINITE VOLUME METHOD, SIMPLE ALGORITHM, & GRID GENERATION.27 3.1 A lgorithm O overview ................................................................ ................... 27 3.2 Discretization of the Governing Equations............................. ..................28 3.3 SIM PLE M ethod ......... .................. ......... ... ..................... 32 3.4 D iscretization Schem es......... .................. .................................. ............... 35 3.4.1 First Order Upwind Scheme ............................ ..... .. ............... 35 3.4.2 Central D difference Schem e..................................... ....................... 35 3.4.3 Second Order Upwind Schem e ...................................... ............... 36 3 .5 G rid G generation .............................. .......................... .. ........ .... ..... ...... 36 4 BENCHMARK TEST CASES / CODE VALIDATION ................. .... ............39 4.1 B enchm ark Case 1, Lid D riven Cavity ...................... ........................ ..........39 4.2 Benchmark Case 2, Developing Pipe Flow ............ ................. ..................... 43 4.3 Benchmark Case 3, Forced Convection in a Duct Partially Filled with a Porous M material ...................... .. .. ......... .. .. ................................................. 4 7 5 TRANSPIRATION COOLED INJECTOR..................... ...... ............... 53 5.1 Drilled Orifice Plate Isothermal Simulations .............. .................................. 53 5.2 Drilled Orifice Plate Non Isothermal Simulations..........................................60 5.3 Drilled Orifice Plate with Center Jet Dynamic Interaction ...........................65 6 OBSERVATION AND RECCOMENDATIONS...............................71 6.1 Stability Observation using a Central Difference Scheme ................................71 6.2 Stability Observations on a Collocated Grid ..................... .................. .......... 71 6.3 R ecom m endations for Future W ork ........................................ .....................73 7 CON CLU SION S....... ............................................................. ...... .. 74 LIST OF REFEREN CES ............................................................................. 76 B IO G R A PH IC A L SK E TCH ..................................................................... ..................78 LIST OF TABLES Table p 21. Properties of common porous materials. ............................................................16 22. Heat transfer coefficients with respective fluid to solid specific surface area. .........26 51. Values of permeability and Forchheimer's coefficient. .........................................55 52. Fluid and porous solid properties, along with inlet velocities examined. ................56 53. Local thermal nonequilibrium models, w. respective coefficients.........................64 LIST OF FIGURES Figure p 11. Pratt and W hitney's R L 10 LPR E ..................................................... ..................... 3 12. Cutaway of a regeneratively cooled tubular thrust chamber .............. ............ 5 13. N onim pinging injector elem ents .............. ......... ..................................................6 14. U nlikeim pinging injector elem ents ........................................ ........................ 7 15. Likeimpinging injector elements ................. ..................... ................. 8 16. Steady state cooling m ethods .................................. ...................................... 9 21. Flow schem atic for D arcy's law .......................................... .......................... 12 22. Illustration of capillary model for a bundle of identical capillaries ........................17 23. Transition from the Darcy regime to the Forchheimer regime in unidirectional flow through an isothermal saturated porous medium ............................................. 20 24. A schematic of a representative elementary volume and the position vectors used..24 31. Generic grid setup for staggered arrangement of variables .................................28 32. G rid setup, w ith notation. ........................................ ............................................29 33. Outline of SIMPLE algorithm ............................. ............... 32 34. Geometry of the problem under investigation with labels corresponding to input variables in the code. ......................... ...... ...................... ........... ............. 37 35. Sample grid displaying the important characteristics of the grid generation process.38 41. Benchmark Case 1, lid driven cavity flow, problem setup....................................40 42. Benchmark Case 1, discretization scheme comparison on a 20 x 20 grid for R e = 1 0 0 0 ......................................................................... 4 1 43. Benchmark Case 1, discretization scheme comparison on a 120 x 120 grid for R e = 1 0 0 0 ......................................................................... 4 2 44. Benchmark Case 2, pressure driven flow in a pipe, problem setup.........................43 45. Benchmark Case 2, results comparison........................................... ...............46 46. Benchmark Case3, forced convection in a duct partially filled with a porous m material, problem setup .................. .......................... ...................... 47 47. Benchmark Case 3, comparison of velocity profiles to the analytical results of P oulikakos [20]. .................................................... ................. 50 48. Benchmark Case 3, grid with refined region at the porous/open channel interface..50 49. Benchmark Case 3, results of local grid refinement, with 80 nodes in the radial direction ......... ...... ............ ..................................... ...........................5 1 410. Benchmark Case 3, validation of solution to the energy equation........................52 51. Domain setup for both the experimental and numerical investigations of flow through a drilled orifice plate ................................. ............... ............... 54 52. Grid setup corresponding to the geometry seen in Figure 51 ................................56 53. Pressure drop across the orifice plate ............................................. ............... 57 54. Pressure drop across the orifice plate. ............................................ ............... 59 55. Velocity profiles 1 inch downstream of the orifice plate, Case C4 .........................60 56. N onisotherm al porous plate setup .................................................. ..... .......... 61 57. Temperature profiles on the upstream and downstream faces of the porous slug.....62 58. Tem perature profiles, and AT ...................................................................... 63 59. Temperature profiles for LTNE models on downstream face of the porous slug .....65 510. Main injector assembly of the Space Shuttle Main Engine (SSME)....................66 511. Domain setup for both the numerical investigations of flow through a drilled orifice plate with a dynamically influence center jet .....................................67 512. Temperature profiles on the upstream and downstream injector faces ................68 513. Tem perature Profiles and AT...................................... ............... ............... 69 514. Percent mass flow through the center jet for flow speeds, V1V9 .....................70 515. Results from 2UDS simulations with the centerjet ............................................. 70 61. Uvelocity at the centerline................................................ ............................ 72 NOMENCLATURE A cross sectional area a coefficient in the discretization equation CF Forchheimer coefficient cp specific heat Da Darcy number De defined by equation 2.16 dp particle diameter F convective flux g acceleration due to gravity H height in Benchmark Case 1 h height of channel hsf convective heat transfer coefficient II and 12 interpolation factors for 2UDS 1o,I1 modified Bessel functions of the first kind K, K permeability scalar and tensor Ko, K1 modified Bessel functions of the second kind k thermal conductivity L length of porous slug, or length in Benchmark Case 1 rh mass flow rate n number of capillaries Y piezometric pressure P piezometric pressure p pressure p' pressure correction to conserve continuity Pr Prandtl number Q volumetric flow rate Q source term in discretized equations r radial direction Ro radius Red Reynolds number based on capillary diameter ReH Reynolds number based on cavity height, Benchmark Case 1 Re[ Reynolds number based on permeability S momentum source term S surface area in discretized equations s radius at porous interface T temperature Tm mean temperature Ts surface temperature, wall temperature qs" heat flux at surface u x or axial velocity component u' x or axial velocity correction to conserve continuity UD x darcian velocity component UD darcian velocity vector UL lid velocity in Benchmark Case 1 Um mean velocity V volume v y or radial velocity component v' y or radial velocity correction to conserve continuity VD y darcian velocity component w width of channel WD z darcian velocity component x xdirection or axial direction y ydirection or radial direction z zdirection or distance from datum level Greek Symbols AQ volume of CV [t dynamic viscosity [t' effective dynamic viscosity, here Ct' = t p density ; porosity Y any scalar of vector quantity k interpolation factor 1 kinematic viscosity 6 hydrodynamic boundary layer thickness 6T thermal boundary layer thickness Subscripts and Superscripts * values may not conserve continuity * dimensional values in Benchmark Case 3 2UDS second order upwind differencing scheme CDS central differencing scheme D Darcian or Filter E east node EE east node, two away from P e east interface eff effective f fluid phase k generic phase, solid or fluid m iteration number N north node NN north node, two away from P n north interface nb neighboring nodes S south node SS south node, two away from P s south interface s solid phase UDS first order upwind differencing scheme W west node WW west node, two away from P w west interface x x component of a vector y y component of a vector z z component of a vector Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science NUMERICAL MODELING OF TRANSPIRATION COOLED ROCKET INJECTORS By Landon Rothwell Tully May 2005 Chair: Jacob N. Chung Major Department: Mechanical and Aerospace Engineering As society has a growing demand for space travel, the need to advance current technology is evident. With the increasing demand comes a heightened need for a cost effective and reliable means of transport. In order to achieve this, a better understanding and new methods of controlling the temperature of certain components of liquid propellant rocket engines (LPREs) is necessary. Currently a number of component models are based on adhoc assumptions of the fluid flow characteristics. One such example is the injector face. Injector plates for LPREs are often made of porous materials to aid in cooling. One such material is Rigimesh, made from sintered stainless steel. One process of cooling the injector face, termed transpiration cooling, occurs when a small amount of fuel is bled through the Rigimesh, downstream of which combustion takes place. Previous studies have shown the fuel flow rate through the porous media to be proportional to the pressure drop across the plate and the local temperature of the metallic material. However, to this author's knowledge, all prior studies have been heuristic. The focus of this thesis is on numerically modeling the flow through porous materials in an attempt to gain an efficient method of determining the flow field in LPRE injectors and temperatures across the injector face plate. CHAPTER 1 INTRODUCTION LIQUID PROPELLANT ROCKET ENGINES The idea of liquid propellant rocket engines (LPREs) was conceived over 100 years ago. Since then, many designs have been built and tested in the United States alone; and their significance lies in the fact that they are the heart of modern day space exploration [1]. Before presenting the work of this thesis, a brief introduction to the history and design of LPREs is discussed. 1.1 Liquid Propellant Rocket History Over the past century LPREs have become a well developed and established technology; as a consequence they are now the main source of propulsion for space launch vehicles. Robert H. Goddard has been attributed as the first developer of a LPRE. In 1921 he constructed the first LPRE, this lead to the first static hotfiring test in 1923 and later the first flight on March 16, 1926. Goddard went on to achieve many other technological breakthroughs in the development of LPREs; however, much of his work went unrecognized within the industry as a whole. Goddard published very little during his lifetime and was hesitant to release his work for fear that others would utilize unproven concepts. Much of his work has become known after the fact, primarily through the publication of various works by his wife [1]. Most of the designs today are due in part to the work of Aerojet, Rocketdyne, and Pratt and Whitney, to name a few. Aerojet, which originated in 1942, was the developer of many jetassisted take off (JATO) engines. Several of these were successfully flown in the F84 fighter bomber, the PB2Y3 flying boat, and the B29, B45, and B47 bombers. Furthermore, Aerojet also developed the booster LPRE for the Bomarc Area Defense system and is most widely known for their development of the Titan LPREs. Rocketdyne, which started in 1945, is now owned by the Boeing Company, and has been the largest LPRE company in the United States. As stated by Sutton, "up to June 1, 2001 Rocketdyne engines had boosted 1516 vehicles" [1:997] This remarkable accomplishment includes the of the Saturn I and the Saturn V from the United States Apollo program, as well as the historic LOX/LH2 Space Shuttle Main Engine (SSME) developed in 1972 [1]. Pratt and Whitney (P&W), a United Technologies Company, came onto the scene in 1957. One of their most notable accomplishments was the development of the RL10 LPRE, Figure 1.1, in 1959. Since then it has become P&W's most successful LPRE [1]. Deemed by P&W as the "most reliable, safe and high performing upperstage engine in the world" [2:1], the RL10 has been the workhorse of the industry. The injectors of the RL10 are hollow post and sleeve elements, essentially two concentric pipes (see Section 1.2.1 of this thesis). Moreover, P&W engineered a porous stainlesssteel material called Rigimesh, which was used to cool the injector face via transpiration cooling [1]. This method of cooling the Rigimesh is the primary focus of this report and will later be discussed in detail. Figure 11. Pratt and Whitney's RL10 LPRE, from Pratt and Whitney [3]. TheRL10 produces a thrust from 16,500 lb to 22,300 lb and weighs 310 lb to 370 lb depending on the model. The fuel and oxidizer are liquid hydrogen and liquid oxygen, respectively. 1.2 Liquid Propellant Rocket Engine Components The primary function of an LPRE is to generate thrust through combustion. High temperature and high pressure gasses are produced in the combustion chamber (Figure 1 2) and ejected at very high velocities through the nozzle. In the many mechanisms of a LPRE, the thrust chamber is the main component responsible for the generation of the thrust. Within the thrust chamber, the liquid propellants are injected, the droplets are atomized then vaporized, combusted, accelerated to sonic velocities at the throat, and then supersonic velocities within the diverging section of the nozzle. A number of components are responsible for the above mentioned events, for example, the injector is responsible for the introduction and atomization of the propellants into the combustion chamber. As expected, the combustion chamber vaporizes the propellants, combusts them, where they are accelerated through the nozzle and finally ejected [47]. Brief discussions of the components of a thrust chamber are presented here, with emphasis on the particular components applicable to this report. 1.2.1 Injectors The primary function of the injectors is to introduce and gauge the flow of the oxidizer and fuel to the combustion chamber. Over the decades of successful LPRE development many styles and designs of injectors have been offered, some more successful than others. Several factors, such as, engine performance, combustion stability, reliability, structural integrity, the weight and cost of the injectors, thermal protection, and hydraulic characteristics all contribute to an injector's success. To date, most designs have been empirical. The primary focus of this thesis is the thermal Onxyg Fuel int' ,F MFeniifo C' Seakig ringmlj pyroldikic ivpilera EUkd leadl / CoinbutWoWi sectn (pCfoWIcKig) 7/ cikautivc ianl (Iik*uD Lc'fir 1 Tirwi a, ,ril tuiil ` Cuhn> stilinaf binds i tjector plate U=rated 180t iWt ThMat .^a^ QwnW Figure 12. Cutaway of a regeneratively cooled tubular thrust chamber, from reference Sutton [1] and Sutton [5]. Originally used in the Thor missile. The nozzle diameter is about 15 inches and originally produce a thrust of 120,00 lbf, which was later upgraded to 165,00 lbf. modeling of the injector face by numerical methods. For a more detailed discussion regarding the other design parameters the reader is referred to the following references: [47]. sitratiIle r ,D"rm lu lredt w fueled Farid Id'  Fuel cooled tubulal mil Tlrnlarwrl s lbrfing bid monkl Fuel rIfIig Fnifmfllji Throughout the history of LPREs, a number of injector elements have been used. These are commonly split into three categories: nonimpinging, unlikeimpinging, and like impinging. Nonimpinging elements inject the fuel and oxidizer separately into the combustion chamber and include: showerhead, fan former, and coaxial elements as seen in Figure 13. The coaxial elements are similar to those mentioned previously for P&W's RL10 LPRE. Oxidizer Manifold Fuel Manifold FL (A) Shower Head Gaseous Hydrogen Liquid Oxygen Injector face (B) Spray or "Fan Former" ,'.uler Sli. Inner Sleeve Injector face (C) Hollow Post and Sleeve Element (Coaxial) adopted from [Sutton & Biblarz] Figure 13. Nonimpinging injector elements, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. Unlikeimpinging elements direct a stream of fuel and a stream of oxidizer to the same point, where they impinge upon each other and mix. These include unlike doublets and unlike triplets, as seen in Figure 14. Injector face Injector face FIE IFuel Oxidizer Manifold Oxidizer Manifold Fuel Manifold Fuel Manifold (A) Unlike Doublet (B) Unlike Triplet Figure 14. Unlikeimpinging injector elements, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. Likeimpinging elements are very similar to the unlikeimpinging elements in design. However, rather than directing fuel and oxidizer at one another; fuel is injected at fuel and oxidizer at oxidizer. This produces a fan shaped spray of the separate propellants, which helps in atomization of the propellants. As expected the types include, like doublets (commonly called self impinging) and like triplets. Like doublets are the most common and can be seen in Figure 15. One common aspect to all types of injector elements is the injector face. In the case of the SSME and the P&W's RL10 the injector face is made of a porous material, through which fuel is bled to aid in cooling the face. This process has been termed transpiration cooling; and the numerical modeling of this process is the primary focus herein. Later, a simplified version of a single injector element will be thoroughly discussed and analyzed. Injector face Oxidizer Manifold Fuel Manifold LI _ Like Doublet (Self Impinging) Figure 15. Likeimpinging injector elements, also known as, like doublet or self impinging, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. 1.2.2 Combustion Chamber and Nozzle The combustion chamber is the portion of the thrust chamber where the propellants are mixed and burned. The nozzle is typically designed as an integral part of the combustion chamber, and has a converging diverging shape to accelerate the propellants to supersonic velocitites. Overall, the design of the combustion chamber and nozzle is very dynamic. The volume and shape of these chambers must provide complete combustion of the propellants, provide adequate cooling, as well as accelerate the flow to the necessary speeds. Cooling of the thrust chamber is extremely important. Without proper cooling the walls of the combustion chamber and nozzle could become too hot, causing the materials to weaken, possibly to the point of failure. Fundamentally, all internal faces of the LPRE are exposed to hot gases, specifically the injector face and walls of the combustion chamber and nozzle. A typical point of concern is the wall temperatures at the throat, where the heat transfer rate is greatest. However, the temperatures at all locations within the LPRE must be controlled. Ordinarily there are two methods of cooling the combustion chamber and nozzle; the steady state method and heat sink cooling. One method of steady state cooling is to line the thrust chambers are with a cooling jacket. The throat region, previously mentioned, requires the greatest care because of the extreme heat transfer rate hence cooling jackets are designed to provide the highest coolant velocity at the throat. Additionally, the coolant, usually fuel, is run through the cooling jacket and then fed into the injectors; this process is called regenerative cooling. It recycles the heat absorbed by coolant rather than letting it go to waste. Another type of steady state cooling is radiation cooling. Here an additional section is added to the nozzle exit where it becomes extremely hot and radiates the extra heat into space. An example of the two steady state cooling process, cooling jackets and radiation cooling, can be seen in Figure 16. Figure 16. Steady state cooling methods. Left: regenerative cooling jacket. Right: regenerative cooling jacket with nozzle extension for radiation cooling. The primary reason for mentioning these cooling methods is an indirect relation to the studies herein. Some researchers, such as Ryan Avenall [8] have been studying the use of metallic foams lining the cooling jacket, particularly in the throat region, to increase the heat transfer rate as in a heat exchanger. The governing equations and developed code in this report may be of future use in the modeling of metallic foam lined cooling jackets, injectors, or fluid flow through any porous material. Alternatively, cooling with transient heat transfer has been used. In some short duration engines the chambers are made thick enough to absorb the heat of the combustion process. Other methods include ablative cooling and heat sink cooling. These methods are only briefly mentioned for totality; a more detailed discussion can be found in references: [47]. 1.3 Objectives of the Study The primary objective of this study is to numerically investigate the conjugate heat transfer and flow characteristics of transpiration cooled injectors. In the cases of the SSME and P&W's RL10 the injector face is made of a sintered stainless steel material, RigimeshTM. The ultimate goal of this study is to obtain a method of acquiring the appropriate material properties to accurately model flow through the RigimeshTM. To do so, the methods must first be applied to a porous material with wellknown material characteristics and a known physical geometry. Here, flow through an orifice plate made of Beryllium copper will be examined numerically, and compared to the respective physical experiments performed by Ahmed F. Omar, which are to be published in the near future in his dissertation. The results of this study will aid in future simulation and design of rocket combustors. Additionally, by gaining a better understanding of flow 11 through porous materials, the computational models put forth herein can be extended to handle cooling jackets lined with metallic foams. This too, will aid in the future development of LPRE combustion chambers. CHAPTER 2 HISTORY OF POROUS MODELING This chapter presents a brief history of modeling flow through porous media. First, a presentation of the models and a discussion of important terms are provided, then the details and modern equations for modeling are offered. Additionally, a brief review of previous modeling efforts by some other authors is presented. 2.1 Initial Modeling Efforts Since 1856 people have been investigating methods of modeling flow through porous materials. In fact, it was Henry Darcy's investigations of hydrology (Figure 21) that revealed a linear relationship between the rate of flow through a porous bed and the pressure drop. This linear relationship is now referred to as Darcy's Law, KA A9 Q : (2.1) p/ L In equation 2.1, K is the specific permeability which will be discussed in detail in the following section, Q is the volumetric flow rate, A is the crosssectional area of the Figure 21. Flow schematic for Darcy's law, equation 2.1, recreated from Kaviany [9]. channel normal to the mean flow direction, L is the length of the porous slug in the flow direction, and lastly, Y is the piezometric pressure which is defined as: S= p + pgz. (2.2) In equation 2.2, z is the distance measured upwards (against the gravitational direction) from some datum level [10]. For the purposes of the work described herein the gravitational and buoyancy effects are ignored, hence equation 2.2 reduces to J = p. Furthermore, it should also be noted that Darcy's law is commonly generalized in the following form: Vp = u (2.3) K where K is a second order tensor which, again, is discussed in detail in the following section. Additionally, UD is the Darcian velocity vector orfilter velocity vector. The Darcy velocity was originally defined by, UD = (2.4) A thus, it is the average velocity of the flow at a given cross section [10]. Since Henry Darcy's initial efforts, the modeling of porous materials has gained much attention. His model is constantly being verified and extended to handle a wide range of flows. From geological aspects such as ground water modeling, to high speed flows through metal foams used in heat exchangers, the limits and use for porous models are endless. As with all engineering applications, efficient and affordable experiments are not always attainable. Therefore, scientists and engineers revert to these models to gain an understanding of systems before they are used for real life applications. As in many cases, such as in this thesis, the applicability of these models is still being verified, therefore herein, the results of both numerical and physical experiments are discussed simultaneously. Furthermore, it is important to note that the physical experiments, and their resultant data used herein, were not performed by the author. A majority of the data from physical experiments referenced within thesis was provided by the work of, Dr. Bruce F. Carroll and Ahmed F. Omar, or in some cases a specific reference may be mentioned. All of the experiments performed by this author are numerical, and the results of physical experiments may be referenced for comparison and validation of the computational results. 2.2 Permeability and Porosity Before any modeling efforts can be presented, a more detailed discussion of the permeability and porosity of the porous structure is needed. Both properties are based on the geometry and structure of the pores in the medium and are necessary to modeling the flow through such media. 2.2.1 Porosity To commence, the porosity, F, of any material is defined as the volume fraction of the porous sample occupied by voids. In a saturated fluid all of the volume of the voids is occupied by the fluid, hence the porosity is defined as, E= (2.5) V where V = Vf + VJ, and the subscripts f and s represent the fluid and solid phases respectively. 2.2.2 Permeability Next, the general term, permeability, refers to the conductivity of a fluid through the porous medium. This, however, is not very useful because it may vary with fluid properties. As an alternative, the specific permeability is referred to as the conductivity of the fluid through the porous material, independent of the fluid properties. For simplicity, for single phase flow the specific permeability will be referred to as permeability. Based on this definition, the permeability is a material property based on the geometry of the pore and is generally a second order tensor. It has units of length2, or the unit of the darcy, which was created for practicality [10]. In the words ofDullien, "a porous material has permeability equal to 1 darcy if a pressure difference of 1 atm will produce a flow rate of 1 cm3/sec of a fluid with 1 cP (centiPoise) viscosity through a cube having sides 1 cm in length" [10:78], hence, 1(cm3 sec)1(cP) 1 darcy = m c)l(t = 0.987 n2 = 0.987* 10 12 2 1(cm2 1(atm/cm) In the case of anisotropic porous media, the permeability is represented by a second order tensor that has nine components as follows, K,, K,, K,3 K = K2 K22 K3. (2.6) K,3 K32 K33 A great deal of effort has been put forth by a number of authors, Whitaker [11] and Guin [12] for example, to prove that K is symmetric under the assumption that the anisotropic porous media is orthotropic (has three mutually orthogonal principal axes), ie., K,, K,, K, K= K, K2, K3 (2.7) K K,, K22 K13 K23 K33 As previously stated, the above tensor is symmetric, moreover a diagonal matrix is formed when the coordinate axes are aligned with the principal axes of the medium [10]. Darcy's law then becomes, u (2.8a) dx K V (2.8b) dy KY dp if WD (2.8c) dz K where Kx, Ky, and Kz have been termed the directional permeabilities. To gain a better feel for the permeability and porosity several values of common materials are given in Table 21. Table 21. Properties of common porous materials. Porous Material Porosity Permeability (m2) Foam Metal (also made of other materials) 0.98 Fiberglas 0.880.93 2.4E11 5.1E11 Berl saddles 0.680.83 1.3 E07 3.9 E07 Wire crimps 0.68 0.76 Black slate powder 0.57 0.66 4.9 E 14 1.2 E 13 Raschig rings 0.56 0.65 Leather 0.560.59 9.5 E14 1.2 E 13 Granular crushed rock 0.44 0.45 Soil 0.430.54 2.9E13 1.4E11 Sand 0.370.50 2.0E11 1.8E10 Silica powder 0.37 0.49 1.3 E 14 5.1 E 14 Spherical packing, well shaken 0.36 0.43 Cigarette filters 0.17 0.49 1.1 E 09 Brick 0.120.34 4.8 E 15 2.2 E 13 Sandstone (oil sand) 0.080.38 5.0E 16 3.0E 12 Limestone, dolomite 0.040.10 2.0 E 15 4.5 E 14 Coal 0.020.12 Concrete (ordinary mixes) 0.02 0.07 Data taken from Kaviany [9]. 2.2.3 Permeability Models The value of the permeability is often determined empirically, however, several models have been proposed based on the geometry of the porous medium. In general, the above mentioned models can be split into two categories, capillary models and drag models. In capillary models, the flow is considered in complex channels or conduits, where, as in drag models, the flow is considered over objects. In this thesis the capillary models are of greater interest due to their applicability to the previously mentioned RigimeshTM and the orifice plates used in the physical and computational experiments performed herein. However, for diligence both are mentioned and several models will be described. 2.2.3.1 Capillary models First, capillary models are based on known solutions of the NavierStokes equations for internal flows. One model is based on laminar flow through a bundle of circular tubes of equal length and diameter D (Figure 22), in which the corresponding porosity F, is nid p2 /4 for n tubes per unit area. Capillaries / Conduits (DI I SOO o 0 Figure 22. Illustration of capillary model for a bundle of identical capillaries. The permeability can be derived beginning with the well known HagenPoiseuille equation, d2 dp u d p (2.9) 32/, fdx where Um is the mean velocity in the pore. Subsequently, the mean pore velocity can be related to the Darcian velocity by n7d4 dp UD = sum (2.10) S 128/@ 8 dx Then, when Equation 2.10 is compared to Equation 2.3, the permeability is found to be ed2 K= (2.11) 32 Similar models have been proposed for ducts of varying crosssectional area and a network of conduits. For more information on capillary models, the reader is referred to Kaviany [9] and Dullien [10]. 2.2.3.2 Drag models In drag models, flow is considered around objects, spheres, cylinders, etc. rather than through a network of capillaries as previously discussed. Here the total resistance to the flow is compared to Darcy's law as a means of obtaining the permeability of a particular material. One example is for flow over cylinders. Several authors have investigated this numerically and have extracted the permeability based on a curve fit; one such example is cited in Kaviany [9]. Other authors have investigated flow over spheres, different arrangements of cylinders, and other various shapes. Generally, a great deal of work has been put forth into the determination of the permeability based on these drag models, however, for brevity, these models are merely mentioned. For additional details on this subject, the reader is referred to Kaviany [9] and Dullien [10]. 2.2.3.3 Experimental methods Experimental rigs used to determine the permeability of a material are commonly called permeameters. Many styles of these devices have been designed, but in general the methodology remains the same. Measure the pressure drop for a number of flow rates, of which Reynolds number is on the order of unity; then fit a straight line to the data points. If the line passes through the origin at zero, then Darcy's law is being obeyed; however, because of the error in the experimental efforts this may not always be the case. If a straight line cannot be fit to the data, then the system is not obeying Darcy's law and further investigation may be necessary [10]. 2.3 High Reynolds Number Flow Darcy's law, as discussed above, is only applicable to low speed flows, i.e. uD is sufficiently small. Typically this means that the Reynolds number, based on the average pore diameter, has an order of magnitude near unity or less, puDd Red, 0(1) (2.12) Pf where dp is the average particle diameter. In this case, the diameter of the capillary is used. However, tt should be noted that other authors have defined the Reynolds number as Re = (2.13) 'If Based on this definition, Darcy's law is valid for Re / up to approximately 0.2. For the purposes of this report, unless otherwise specified, Red will be used. The limited applicability of Darcy's law is based on the fact that it only accounts for the viscous resistance of the flow. For highspeed flows the inertial contributions to the flow resistance become noticeable. A number of authors have attempted to account for the inertial effects, and a method that has gained wide spread use has been termed Forchheimer's equation, Vp = uD pF uD (2.14) where the CF is the so called Forchheimer coefficient. The transition from the Darcy (viscous) regime to the Forchheimer (inertial) regime based on Re K is illustrated in Figure 23. Additionally, a number of different equations have been proposed to handle the inertial effect. For a detailed discussion on the effect the different Forchheimer terms, the reader is referred to Alazmi [13]. NH^rT .....t.H. h+.. fK= ,ReK fK 1 +0.55 S.1 __ ReK I A I II I [ I [ll __.t1[44 1.. . 0.1 0.1 1 10 100 vK1/2 RCK= Figure 23. Transition from the Darcy regime to the Forchheimer regime in unidirectional flow through an isothermal saturated porous medium, scanned from Nield [14]. As seen in Figure 23, the divergence from linearity begins at Re K of approximately 0.2, with a greater influence in the range of 1 10. Additionally, at high Reynolds numbers, the inertial resistance dominates because at high Darcian velocities the inertial term (uD2 term) is much larger than the viscous term, as can be inferred from equation 2.14, [14]. The high velocity regime is of utmost importance for the current research; hence, both the viscous (Darcian) and the inertial (Forchheimer) coefficients will be considered in the subsequent calculations. 2.3.1 Forchheimer models As for permeability, models are necessary for determining the inertial coefficient. These models are not as common as those for permeability; however, two models will be discussed here. Originally it was thought that the value of CF was constant and approximately 0.55. However, later studies have shown that it varies with different porous mediums. One study for flow over spheres showed that, CF= 0.55 15.5 P (2.15) De where, 2wh D = 2 (2.16) w +h' with w and h being the width and height of the channel, respectively. In the case of a cylinder De is simply the diameter. Other models, beyond Forchheimer's extension, have also been proposed. One such example is an extension of the CarmanKozeny theory, for details see Kaviany [9] and Nield [14], where the fluid and matrix properties are related to the pressure drop, by dp 180(1 )2 1 . 1 2(2.17) dx= f d 2+3 .8 (2.17)dPf dx dS e dp P P Once again, dp is the particle diameter or appropriate characteristic length scale of the porous medium. Also, the numerical constants, 180 and 1.8, in equation 2.17 have been written by Nield as 150 and 1.75, respectively. The relationship then becomes Ergun's correlation. However equation 2.17 as it stands will be used herein. Additionally, the permeability and Forchheimer coefficient can be extracted from equation 2.17 as d E3 K P (2.18) 180(1 e)2 1.8 CF = (2.19) 1801/2 E3/2 It should be noted that this is an ad hoc procedure for purposes of numerical solutions only. 2.3.2 Experimental methods Similar methods to those performed for determining the permeability experimentally can be applied here. In the experimental procedures applicable to this report the pressure upstream and the pressure downstream of the plate were varied in order to apply an appropriate mass flow rate. The appropriate mass flow rate means the flow rate is large enough for inertial effects to dominate. Then, with the permeability known from the previous low speed experiments, the inertial or Forchheimer coefficient can be extracted based on equation 2.11. Note that this does not include the pressure loss due to any boundary layer effect felt in the real situation. 2.4 Semi Heuristic Equations In order to model flow through porous media in conjunction with flow through plain media, i.e. an open channel, a single set of governing equations is desired. As stated by Kaviany, this would be "too complicated to be of practical use" [9:66]. Instead, many authors have attempted to develop an equation for modeling flow through porous media comparable to the NavierStokes equation. The first attempt at this was by Brinkman, who included a term analogous to the Laplacian term of the NavierStokes equation, Vp= UD + / D (2.20) K where p' is an effective viscosity. Some authors have set p' equal to p/u, while others have defined it using the Einstein formula, given by / = ,f[1+ 2.5(1 )]. (2.21) However, the governing equations used herein are based on the premise that p' equals "If. Other authors have attempted to derive a set of governing equations based on the local volume averaging of the NavierStokes equation. This however leaves a large number of unknowns that require experimental verification [9]. A number of different equations have been proposed; however, one typically accepted form based on semiheuristic techniques is as follows. 2.4.1 Semiheuristic Continuity and Momentum Equations First, the volume average of any scalar of vector quantity can be defined as (WVk = skdV. (2.22) V V Here V is the volume of a representative elementary volume, (REV) as seen in Figure 2 4. Additionally the intrinsic volume average is defined as (WVk)k = I ikdV. (2.23) VkV This leads to the continuity and momentum equations for Newtonian fluid in steady, incompressible flow in porous mediums, interfacial Solid Volume V g \ area) Fluid Volume V1 Contained in A  Figure 24. A schematic of a representative elementary volume and the position vectors used. The fluid phase is shown as continuous, scanned from Kaviany [9]. V.())= 0 (2.24) uPf (U). Vuj p= V f)p + Pf' V2 uf)+S. (2.25) 6 6 Here, c is the previously defined porosity (equation 2.5) and S is a momentum source term representing the flow resistance containing two terms. For the x direction S is represented by S fKx u f (2.26) In equation 2.26 the first and second terms are respectively, the previously discussed Darcy and Forchheimer terms. In the plain media or open channel, the porosity, 6, is unity and the source term is zero. Equation 2.25 then reduces to the well known Navier Stokes equation for the incompressible constant property flow of a Newtonian fluid. 2.4.2 Energy Equation Finally, the steady energy equation with no generation, and assuming thermal equilibrium between the solid and fluid phases of the porous media can be written as (c,()u). 1 T) = V [ ke, V(T, (2.27) where, ke, = (I )k, + skf. (2.28) Additionally, Nield [14], states that the effective thermal conductivity, equation 2.28, is applicable when the conduction between the solid and fluid phases occurs in parallel. If, on the other hand, the conduction takes place in series, the effective thermal conductivity should be written as ( + ) (2.29) kl k kk This effectively gives a bound for the actual effective thermal conductivity, parallel being the upper bound and series being the lower bound. Subsequently, the energy equation can be extended to allow for heat transfer between the solid and fluid phases. For the fluid phase the energy equation can be written as (POC)f uf) )f eff f ]+ fa+h ) Tf ), (2.30) and for the solid phase, as 0o= .[ks (T,)"]hfaT,) Tf ). (2.31) Here the effective thermal conductivities for the solid and fluid phases respectively, can be written as kfe = ck, (2.32) k,e = (1 c)k,, (2.33) where asf is the specific surface area of the interface and hsf is the heat transfer coefficient between the solid and fluid phase. As expected, the heat transfer coefficient must be based on another model. Alazmi and Vafai give three different models seen in Table 22. The effect of these models, along with the effect of the Forchheimer and Brinkman terms are discussed in Alazmi [13], their respective paper. However, for the purposes of this thesis the applicability of all three models to the porous media of interest will be further explored. Table 22. Heat transfer coefficients with respective fluid to solid specific surface area. Model hf asf Notes k(2 +1.1Pr1/3 Re06) 6(1 ) dp dp 4k,4f 033 Re059 20.346(1l E)E2 2 1.064 Pr33 Re059 Re > 350 ydp d dE d 6(1 g) 3 +  0.2555Pr1/3 Re2/3 kf 10k dp Models taken from Alazmi [13]. CHAPTER 3 FINITE VOLUME METHOD, SIMPLE ALGORITHM, & GRID GENERATION As previously discussed, the use of porous injector plates is common in rocket engines. However, full fledged experiments, in order to gain temperature profiles and pressure drops through the porous media, are not always practical. Additionally, for many practical applications, the NavierStokes equations, with or without the addition of the momentum source term for porous media, can not be solved analytically. A computer simulation on the other hand, can be fast, easy, and most important, affordable. Throughout this chapter the methods used to discretize the partial differential equations presented in Chapter 2, into a set of algebraic equations and then solve those equations are discussed. 3.1 Algorithm Overview The algorithm and notation within is based primarily on the methods presented in Chapter 7 of Ferziger [15], as well as Patankar [16]. The code employs the SIMPLE based pressurecorrection method for solving the NavierStokes equations on a Cartesian or axisymmetric staggered grid. The discretization is based on the finite volume technique, where the diffusive terms are approximated using a central difference scheme (CDS) and the convective terms are approximated by one of three schemes, first order upwind scheme (UDS), second order upwind scheme (2UDS), or CDS. A deferred correction approach of Stone's Method is used for solving the resulting system of linear equations. 3.2 Discretization of the Governing Equations The finite volume technique is based on dividing the domain of interest into a finite number of control volumes (CVs), or cells. The governing equations are then discretized on a structured grid which is created by defining each grid line in space. As a result each scalar control volume has been defined by its' grid lines, see Figure 31. Internal Node Location (Center of CV) Figure 31. Generic grid setup for staggered arrangement of variables; denote scalar variables,  denote u velocities, and T denote v velocities. The locations of scalar variables such as pressure and temperature are on each node, while the location of the u and v velocities are staggered with respect to the pressure as seen in Figure 31. The location of each scalar node is determined by finding the geometrical center of each CV and the locations of the u and v nodes are on their respective grid lines. The reason for selecting the staggered grid was to prevent impractical oscillations in the solution, in particular when a porous source term is added; this will be discussed later in Chapter 6. Once a computational domain has been created the discretized form of the governing equations can be determined. To begin, the reader is referred to Figure 32 for an understanding of the notation used herein. The notation discussed is based on the Nodes * Scalar Boundary Node o Scalar Internal Node S u, v Boundary Node, Respectively Su, v Internal Node, Respectively Control Volumes SScalar CV W u Momentum CV Sv Momentum CV Ax y,.   Yi+1 N o 0 0 Ay (ij+1) Yi nw n ne W P E O w 0 e 0 (i1 ,j) (i,j) (i+1 ,j) Yi1 sw s se S y o o o (i,j1) L Yi2 Xi2 Xi1 Xi Xi+1 S x Figure 32. Grid setup, with notation. respective control volume being examined. For example, when looking at the u momentum equation the notation will be based on that CV, likewise for the vmomentum equation and for scalar variables. In general, the equations herein well reference the control volume being examined by a superscript. To begin, the node of interest is labeled P; its neighbors are labeled N, S, E, and W, for North, South, East, and West, respectively. For instructive purposes a uniform discretization will be assumed, Ax and Ay are the same for all CV's. However, this is not necessary, nor is it required in the code presented. Through the methods presented by Ferziger [15] and Patankar [16] we arrive at the following discretized form of the momentum equations: apup anb nb total nb Where nb = N, S, E, W (3.1) v vZa'v* +Q't v n = av nb tta nb where the superscripts refer to the fact that the velocities may not conserve continuity, and must later be corrected. Additionally, the coefficients of equation 3.1 can be determined as follows. Here they are presented for the umomentum equations only and can be easily extended for use in the vmomentum equations by simply changing the superscripts. For a first order upwind scheme, the coefficients of equation 3.1 are: a" = r ,0 + S (3.2a) YN  a" =[ .',0+ P'S (3.2b) a = I O,0+ I eS(3.2c) XE X a" =[ho + ,0 S" (3.2d) XP XW a, = a Qpoous A Where, nb N, S,E, W. (3.2e) nb Here the symbol, [[A, B]], has been used to denote the greater of A and B, the S terms represent the surface areas of the control volumes, and AQ is the volume of an individual CV. The source term due to the porous media are represented in equation 3.2e as Q" rou In final discretized form the viscous and inertial sources become, SP f F Pf CF ) Q =Kuu (3.3) porous Kx Finally, the source term Qtotal of equation 3.1 contains the pressure and the convective flux resulting from the deferred correction approach, Qotai = Q pressure + Qo,.ectve. (3.4) Where the pressure and convective flux components are, respectively, Q =(PeSe pw.S\)m' (3.5) Qcuonvechve (Fonvecve U (Fconvecve CDS2US 1 (3.6) The deferred correction approach is used to obtain a high order of accuracy, while only adding a small amount of effort compared to low order schemes. This can be done by treating the high order terms explicitly and the low order terms both implicitly and explicitly. Essentially, the high order terms are moved to the right hand side and lower order terms are placed on both sides of the equation. The terms on the right hand side are determined explicitly, from the previous iteration, while the terms on the left hand side are found implicitly. As the solution converges the lower order terms (which are found on both sides of the equation) drop out and the final solution is that of the higher order scheme [15]. The method used above was to treat the convective flux as Fc iLUUS C (CDS2DS 2US UDS )11 Fe = u +h[eU ue ) (3.7) where the explicit terms are denoted by the superscript "ml", meaning values found from the previous iteration. Additionally, the superscripts UDS, CDS, and 2UDS refer to a first order upwind difference scheme, a central difference scheme, and a second order upwind difference scheme respectively. These are discussed further in Section 3.4 of this thesis. Furthermore, the explicit terms can clearly be seen in Qonvece, equation 3.6, and the implicit terms in the coefficients a", where nb = N,S,E, W. Here the solution effort and progress should be similar to that of an upwind difference scheme, while the final solution should have the accuracy of the higher order scheme. Finally the discretized equations can be solved. Many methods are available, here the strongly implicit procedure (SIP) or Stone's Method is used. It is a form of incomplete lowerupper (LU) decomposition. More detail can be found in Ferziger [15]. 3.3 SIMPLE Method As a result of the above approximations being based on values from previous iterations, continuity may not be satisfied. Thus, we must correct the resulting velocity and pressure fields to account for the change in mass flux. The procedure used has been given the name SIMPLE, which stands for SemiImplicit Method for PressureLinked Equations. This procedure was originally put forth by Patankar and Spalding [17]. As a reference, and introduction, the reader is referred to Figure 33, the details of which will be discussed subsequently. Yes Set calculated T as Solve Discretized Energy "guessed" field Equation No Convergence Yes Stop Figure 33. Outline of SIMPLE algorithm. To begin, an initial guess of the velocity, pressure, and consequent interface velocities and mass flow rates, must be made. Once these fields have been guessed, the solution of the imperfect velocity fields may be obtained by the methods outlined in Section 3.2. Next these fields must be corrected, to preserve continuity, as follows, p = P +p' (3.8a) S=u + u' (3.8b) v = + v' (3.8c) with p', u', and v' representing the corrections to their respective variables. First the pressure correction must be determined. Using continuity one can arrive at equation 3.9 to determine the pressure correction. For brevity the derivation of this correction is not presented here, however, it can be found in either Ferziger [15] or Patankar [16]. The resulting pressure correction is, aPpP = APn +p f [SS u*Sj]+Pf [vS vS ]. (3.9) nb The coefficients of equation 3.9 are, aj = Pf a (3.10a) a= Pf 2 (3.10d) a e SP S2 aw =f p (3.10d) p w aP = aP Where, nb= N, S, E, W. (3.1Oe) nb Additionally, the following variables can be written as (au)e = (au,)Ee +(a, (1 e ) (3.11) where ke is an interpolation factor defined as A, = xy xP (3.12) XE XP Furthermore, in order to determine the pressure correction the interface velocities on a pressure (scalar) CV, must be evaluated. This leads to the importance of the use of a staggered grid. Here, the interface velocities for a scalar CV have already been defined as the nodal velocities for the u and v CV's as seen in Figure 31. Alternatively, had a collocated grid been used, these velocities would be unknown and would have to be interpolated from neighboring values. This can be achieved by applying the so called Rhie and Chow momentum interpolation technique, however, the steep gradient across the interface of an open channel and a porous zone leads to some problems. This will be discussed in more detail in Chapter 6. With the interface velocities known, the pressure corrections may be found by applying Stone's Method to solve equation 3.9. After the pressure corrections have been determined the interface velocities are updated by =a (PPp) (3.13) ap or, e =U pE p) (3.14) Then the pressures are corrected by PP = pP +( p) (3.15) Finally, check for convergence; if the sum of the absolute values of the residuals is less than the designated tolerance, the solver stops. If not, the method repeats using the velocity and pressure fields obtained above as the initial guess for the next iteration. 3.4 Discretization Schemes For purposes of robustness, there separate discretization schemes for the convective terms are offered in the current code. All three schemes can be directly implemented in equations 3.6 and 3.7 by selection of the proper scheme denoted by the superscripts UDS, CDS, and 2UDS. Here only a brief description of the schemes is provided, for more information the reader is referred to Ferziger [15], Patankar [16], and Thakur [18]. 3.4.1 First Order Upwind Scheme The first order upwind scheme assumes the value of 4 at an interface is equal to the value of its upstream or upwind neighbor. Mathematically the convective flux can be written as Feq = op [[Fe,0]] E [[ Fe,0]]. (3.16) The results of this scheme can be seen in equation 3.2. The first order upwind scheme is the basis upon which the deferred correction approach, previously discussed, is based. The first order upwind scheme is treated implicitly and the higher order schemes (central difference and second order upwind) are treated explicitly. In order to properly apply this approach the higher order schemes must be discussed. 3.4.2 Central Difference Scheme In the central differencing approach the interface values are approximated by a simple linear interpolation between the neighboring nodes, FAe = Fe xe (p )= Fe (p 0). (3.17) XE X The central difference scheme has been applied to all of the diffusion terms as discussed in Sections 3.2 and 3.3. It can also be applied to the convective terms through the aforementioned deferred correction approach without any modification to the algorithm. 3.4.3 Second Order Upwind Scheme The second order upwind scheme, on the other hand, is based on a linear extrapolation of the two upstream neighbors. This results in Fe 0e = 0P Ow)i + W e,0]] {(E EE),2 + EE[[ F,0]] (3.18) where, I,= e (3.19a) XP X I2 X E (3.19b) EE E As with the CDS the 2UDS can be applied through the use of the deferred correction approach. 3.5 Grid Generation In order to perform the procedures outline in the previous sections, a computational grid must first be created. This grid may be created by a number of means; the specific methods are of no particular importance. However, obtaining a proper grid with refined cell sizes in certain regions is crucial. In order to obtain a proper mesh, the grid generation methods used within are discussed in order to familiarize the reader with the methods used herein. First, the dimensions of the geometry, based on the setup and description shown in Figure 3.4, must be entered. Also, it should be noted that the grid generation described here is only applicable to the results seen in Chapter 5 of this report. The grids used in the test cases seen in Chapter 4 were created specifically for those applications and need not be discussed here. Additionally, is should be noted that all inputs are assumed to be SI, hence all lengths should be entered in meters. start: x location of the beginning of the porous region start: xend: x location of the pend: x location of inlet (typically 0.0) x location of the end the exit of the porous region y Sub Domain 2: InletRtotal Rcenter Jet V elocity Total Radius of Pipe Rtotal Renteret Radius of Center Jet: y Sub Domain 1 Figure 34. Geometry of the problem under investigation with labels corresponding to input variables in the code. Once the dimensions have been entered correctly the user must then enter the number of control volumes in each subregion. First, the number of control volumes in the xdirection is needed. Once the number of CV's for each subregion (entry, porous, and exit) has been entered, the y grid lines are generated at each xlocation. The code starts by determining the location of the grid lines in the porous region assuming a uniform discretization. Once these grid lines have been defined the neighboring cells of the open channel to the porous zone (in both entry and exit regions) are defined as the same size as those in the porous region. Next, expansion factors are determined based on the number of cells and the length of the entry and exit regions, respectively. Then the grid lines are defined. A similar procedure is performed in the ydirection. A sample grid can be seen in the Figure 35. Grid Setup 0.025 0.02 E 0.015 (U) C 0.01 0.005 0 0 0.01 0.02 0.03 0.04 0.05 xdirection (m.) Figure 35. Sample grid displaying the important characteristics of the grid generation process. CHAPTER 4 BENCHMARK TEST CASES / CODE VALIDATION As a means of verifying the code developed herein, results of several well known test cases are presented and compared to their respective current simulations. Here three benchmark cases have been selected; each to individually validate a specific aspect of the code. First, the results of the standard lid driven cavity are compared to those of Ghia [19]. The second test case is the wellknown analytical solution for fully developed pipe flow. This test case was used to validate the axisymmetric grid in addition to the solution of the standard energy equation. In this particular case, both the fully developed velocity profile and Nusselt number are compared to their respective analytical solutions. Finally, a less wellknown test case, based on the analytical solution for "Forced Convection in a Duct Partially Filled With a Porous Material," [20], is presented. 4.1 Benchmark Case 1, Lid Driven Cavity The first benchmark case mentioned above is the popular lid driven cavity flow of Ghia [19]. This problem has been well documented and a number of accurate solutions are available within the literature. The results used here, for comparison, are those provided by Ghia [19]. Before proceeding, the problem description is as follows: flow inside a square cavity of length and height both 1 unit, is driven by motion of the top boundary of the domain as seen in Figure 41. Moving Lid UL Figure 41. Benchmark Case 1, lid driven cavity flow, problem setup. For this case the lid velocity was set to 1 unit/sec. The results of Ghia [19] are categorized according to Reynolds number. Here the Reynolds Number is defined as UrH PfUH ReH UH pfULH The following Figures 42 and 43, offered as a measure to the V /P performance of the code under development, illustrate the u and v velocity profiles through the geometric centers in the x and y directions, respectively, for the three discretization schemes, first order upwind (UDS), central difference (CDS) and a second order upwind scheme (2UDS). The results on a 20 x 20 grid can be seen in Figure 42 while the results for a 120 x 120 grid can be seen in Figure 43. On a course grid the 2UDS solution was best; however, as the grid was refined, both the CDS and 2UDS schemes provided accurate solutions. A) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.2 0.4 0.6 0.8 uvelocity (m/sec) 1st Order Upwind Central Difference 2nd Order Upwind o Ghia etal. (1982) x(m) Figure 42. Benchmark Case 1, discretization scheme comparison on a 20 x 20 grid for Re = 1000. A) UVelocity through ygeometric. B) VVelocity through x geometric center. 1st Order Upwind  Central Difference S 2nd Order Upwind o Ghia etal. (1982) 0.1 M 0 E 5 0.1 >, > _ 42 A) 1 0.9 0.8 0.7 0.6 E 0.5 0.4 0.3 1st Order Upwind 0.2  Central Difference S 2nd Order Upwind 0.1 o Ghia etal. (1982) 01 0.2 0 0.2 0.4 0.6 0.8 1 uvelocity (m./sec.) B) 1 st Order Upwind 0.3 + Central Difference 2nd Order Upwind 0.2 o Ghia etal. (1982) 0.1 O 0.1 0 > 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 x (m.) Figure 43. Benchmark Case 1, discretization scheme comparison on a 120 x 120 grid for Re = 1000. A) UVelocity through ygeometric. B) VVelocity through x geometric center. 43 4.2 Benchmark Case 2, Developing Pipe Flow The second benchmark case is pressure driven flow in a circular duct, described graphically in Figure 44. This problem has several wellknown analytical solutions for particular aspects of the flow field. For example, the fully developed profile is a well known function of radius and mean velocity. However, analytical solutions to the entire flow field, including the developing flow regime, are not as straightforward, particularly when both the velocity and thermal boundary layers are developing simultaneously. For this reason, the geometry was created with substantial room to allow the flow to reach its fully developed conditions, both hydrodynamically and thermally. The resulting profiles in the fully developed region are compared to their analytical solutions. Hydrodynamic Entrance Region Fully Developed Region Uinlet = Urm Tilet Figure 44 .6 .......) R. Thermal Entrance Region FD Region .. .. .. .. .. .. .. .. .. .. . ......... ............... 7 ConstantNNWall Temperature, Twall > Tinlet Benchmark Case 2, pressure driven flow in a pipe, problem setup. As shown above both the hydrodynamic and thermal boundary layers are developing simultaneously. ......II; The analytical solution of the velocity profile for fully developed laminar pipe flow can be found in most fluids or convective heat transfer texts, for example Incropera [21], and is as follows: u(r)= 2um 1 r (4.1) where um is the mean velocity of the flow, which for incompressible flows is defined by u =IudA. (4.2) Next, the dimensionless temperature gradient at the wall, the Nusselt number, is a well known constant when the flow is thermally fully developed. Therefore, it is used to validate the solution of the energy equation. The Nusselt number is defined as hD NuD = (4.3) kf where D is the diameter of pipe, kf is the thermal conductivity of the fluid, and h is the heat transfer coefficient defined by q" = h(T T). (4.4) In equation 4.4, q" is the heat flux at the wall or surface, OT q = k, O (4.5) Sr f = R , while Ts is the wall temperature, and Tm is the mean or bulk temperature. For this case the wall temperature is known, as it is a prescribed boundary condition, however, the mean temperature must be computed. Similar to the mean velocity, the mean temperature for incompressible flows can be determined by T 2 JuTrdr. (4.6) 0 0 umR2o (46) Finally, combining equations 4.3 through 4.6 leads to the following definition of the Nusselt number, Tr D O r=Ro NuD (= T .) (4.7) It is important to note that the reason for deriving the above form of the Nusselt number is its ease of calculation in the numerical solution. As it stands, the dimensional temperature gradient at the wall and mean temperature can be easily computed resulting in a simple computation of the Nusselt number. Here the results shown are independent of the grid size in the xdirection, because the flow was given ample room to become fully developed. Consequently, only the number of nodes in the ydirection will be noted. The results with 22 nodes in the y direction can be seen in Figure 45, for a Reynolds number of 20. As can be seen, the results are very accurate, with a maximum absolute error in the velocity profile of approximately 4.2E5, or 0.00200. Also, it should be noted that the Nusselt number converged to an acceptable value, Nu = 3.66, which is accurate within 2 decimal places of the analytical solution. % Error (Analytical Numerical)/Analytical 0.01 0.008 0.006 0.004 0.002 SAnalytical Numerical uvelocity (m/sec) E 0.05 0.04 0.03 0.02 0.01 0 B) 15 10 9 8 7 z 6 5 4 3 Figure 45. Benchmark Case 2, results comparison. A) Fully developed velocity profile in comparison with the analytical solution and % error. B) Nusselt number vs. dimensionless axial distance. 0 0.002 S % Error 0.09 0.08 0.06 0.05 0.04 0.03 0.02 0.01 102 10 1/Gz = (x/D)/(ReD*Pr) n\ \*, 4.3 Benchmark Case 3, Forced Convection in a Duct Partially Filled With a Porous Material The third and final test case is based on the analytical solution provided by Poulikakos [20]. In this thesis flow is assumed fully developed both hydraulically and thermally for the geometry depicted in Figure 46. SP,:,rous Media Un SOpen Channel Figure 46. Benchmark Case3, forced convection in a duct partially filled with a porous material, problem setup, recreated from Poulikakos [20]. The results provided are velocity profiles as a function of the Darcy number, K Da= (4.8) R2o and porous substrate thickness, s. In this particular case the porous medium is assumed to be isotropic; consequently, the permeability is the same in all directions as denoted by the lack of subscript on K in equation 4.8. Additionally, several plots of the Nusselt number, also as a function of Da and s are provided. However, no closed form solution was presented, only a representative figure. The resulting solution of Poulikakos and Kazmierczak for the fully developed velocity profile is as follows. Fluid Region, 0 O r < s 2 u = + A (4.9) 4 Porous Region, s 5 r < 1, u =Blo Da2 + CK yDa Da (4.10) where, A JKsDa lj  B IDaK Da a + YK, Da (4.1 B= 2(4.11) I0 Da Y )K sDa +I sDaY2 KO Da S((4.12) sD2 K sDa SDa+2 2 sP, KI sDaY2 ) 4 II sDa,Y2) s Y u= (4.14) R dP. /Jf dx, with the denoting the standard dimensional terms. Here, several discrepancies were noticed. The results of equations 4.9 through 4.13 do not match the results shown in Figure 2 (b) of the report. For this reason the equations were scrutinized and properly written as, Fluid Region, 0 O r < s u= +A (4.15) Porous Region, s 5 r < 1, u =IBIo yDa + CKo yDa2 Da (4.16) Equations 4.15 and 4.16 were not derived from first principles by the author; they were simply inferred from the Figures provided by Poulikakos and Kazmierczak. The corrections were obtained by matching the velocity at the porous interface, r = s, of equations 4.9 and 4.10. Here the addition of a y in the term BI Da 2 to give BIo (yDa 1/2 provided matching velocities at the interface. Also, the negative sign seen in both equations 4.15 and 4.16 was added to correct for the direction of the velocity profiles along the positive xaxis. The results of a single test case using a second order upwind scheme, with s = 0.8, and 240 evenly spaced grid points in the radial direction can be seen in Figure 47. As Figure 47 shows, the numerical results of the code under development are within reason. The maximum error of approximately 0.12% occurred at the interface between the solid and porous region, which lead to a local refinement of the grid at the interface. The second grid investigated consisted of only 80 nodes in the radial direction, with local refinement at the interface as seen in Figure 48. The results of the local refinement are excellent. With only a third of the total grid points used in the previous simulation, the maximum error was decreased by slightly more than one half, from 0.12% to 0.056%, as seen in Figure 49. Not only did the error decrease, but the time to run the simulation, diminished. % Error (Annalytical Numerical)/Analytical 0.02 0 0.02 0.04 0.06 0.08 0.1 0.1 % Error 2nd Order Upwind Analytical 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 uvelocity (Dimensionless) Figure 47. Benchmark Case 3, comparison of velocity profiles to the analytical results of Poulikakos [20]. 1    0.9 0.8 0.7 0.6 S0.5 0.7  rQ 0.4 0.3 0.2 0.1 0 0 10 20 30 40 50 60 70 80 90 100 xdirection (m.) Figure 48. Benchmark Case 3, grid with refined region at the porous/open channel interface. Here only a total of 80 nodes were used in the radial direction. Finally, the results of the solution to the energy equation are presented. Here Poulikakos provided no analytical solution, however, a solution was found numerically. The results of the above mentioned numerical solution can be seen in Figure 410A. Respectively, the results of the code under investigation can be seen in Figure 410B. By inspection, the results of the current code are within reason. No noticeable error can be seen, deeming the solution correct. % Error (Annalytical Numerical)/Analytical 0 0.01 0.02 0.03 0.04 0.05 0.06 % Error 2nd Order Upwind Analytical 0 0 0.02 0.04 0.06 0.08 0.1 0.12 uvelocity (Dimensionless) T1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.14 0.16 0.18 Figure 49. Benchmark Case 3, results of local grid refinement, with 80 nodes in the radial direction. i6 A.. 2  _. .!! , S n _ '' .*'.'' (b) Iio' 10" .0' 10 Da B.) 6 5 /s=1.0 s=0.9  s=0.7 s=05 s=0.5 J ./, ~ 10" 105 104 103 Darcy Number (Da) Figure 410. Benchmark Case 3, validation of solution to the energy equation. A.) Numerical results of Poulikakos for a constant wall temperature boundary condition, scanned from Poulikakos [20]. B.) Results from current code under investigation. 102 101 CHAPTER 5 TRANSPIRATION COOLED INJECTOR In Chapter 1 of this thesis, a number of different injector styles was presented. Here a simplified version of the hollow post and sleeve, or, coaxial element is investigated. A single stream of fluid is injected through a large orifice, while a smaller amount of the same fluid is bled through the injector face. The main focus of this report is the flow through the porous injector face and the respective temperatures. The investigations herein are based on a drilled orifice plate acting as the porous medium. Future investigations should include the actual material used in the SSME, RigimeshTM. However, due to the large number of unknowns of the RigimeshTM, such as material properties and pore geometry, a set of experiments must first be performed on a porous medium with a wellknown geometry. This will provide the proper methods of determining the material constants, and should lay the foundation for future studies of the RigimeshTM. 5.1 Drilled Orifice Plate Isothermal Simulations In order to accurately model flow through porous materials, the correct values of porosity, permeability (viscous coefficient), and the Forchheimer's coefficient (inertial coefficient) must be determined. Several methods were previously discussed, they are applied here. The results are compared to the corresponding experimental results, in an attempt to determine the appropriate method for extension to the RigimeshTM. To begin, the numerical and experimental setup studied can be seen in Figure 51. 54 R 1.036" 2.63144 cm 000000000000 0000000000000000 ooo00000000000 000000 00000000000000000 HolePattern / REV 00000000000 00000000 000000000000 00000 a R 0.010" / 0.0254 cm 000000000000 000000 0.066"0.16764cm 00000000000000000000 S000000000000000000 Inlet 1Porous 0000000000000000000000 0.020" / 0.0508 cm 0000000000000000000000 00000000000000000000000 b S0000000000000000000000 000000000000000000000 0000000000000000000 a8)20.033"0.08382cm 000000000000000000 b =a = 0=033"/4 08382 cm 00000000000000000 b = 0.066 / 0.16764cm 00000000000000 000000000000 00000000 0.05715 m < 0.03175 m 0.0254 m /  iw^ ^ ^ ^^ ^ ^ ^^ ^ ^ Inlet Porous 7 0.0263144 m, Velocity p i n n t p Plate    Pi P2 Figure 51. Domain setup for both the experimental and numerical investigations of flow through a drilled orifice plate. For a known geometry the porosity can be easily computed. Here, the representative elementary volume (REV) seen in Figure 51 yields a porosity of 2 '(0.0508)2 eVoid Volume 0.1442. (5.1) Total Volume 0.167642 With the porosity known, the models for permeability and the Forchheimer coefficient, discussed in Chapter 2, can now be applied. The results are given in Table 51. Furthermore, the constants for Case 3 were determined from the experimental results put forth by Dr. Bruce F. Carroll and Ahmed F. Omar. The permeability was extracted from equation 2.3, based on the pressure drop at low flow velocities, an inlet velocity of 3.58 m/s. This corresponds to Re on the order of 100, which is expected to introduce some error because this is much greater than unity, as required by Darcy's law. The Forchheimer coefficient was then extracted from equation 2.14 using the previously determined permeability and the pressure drop from the largest tested flow velocity of 25.754 m/s, corresponding to a Re of approximately 840. A similar procedure was applied to determine CF for Case 4; however, it was based on the permeability from the capillary model. Table 51. Values of permeability and Forchheimer's coefficient. Permeability K (m2) Forchheimer Coefficient CF Case Model Value Model Value ed 2 d Cl P 1.16324E9 m2 0.55 15.5 0.520801 32 De d2E 1.8 C2 p 5.87529E12 m2 18 2.44905 180(1 c)2 1801/2 3 C3 Experimental 7.71735E10 m2 Experimental and 0.160154 K from C3 C4 P 1.16324E9 m2 Experimental and 0.205578 32 K from C4 In order to accurately compare the results of the four test cases, the simulations must be run such that there is no dependence on the grid. To facilitate a grid independent solution, the number of nodes in the domain of interest was increased until the absolute difference between the current results and those computed on the previous grid was negligible. This was achieved, with 120 xnodes, 40 in the porous region, and 120 y nodes, resulting in the mesh seen in Figure 52. 0.025 0.02 0.015 0.005 0.oo  0 0 0.01 0.02 0.03 0.04 0.05 xdirection (m) Figure 52. Grid setup corresponding to the geometry seen in Figure 51. The grid is refined at the open channel / porous interface and the near wall regions. Subsequently, simulations were run for the four cases of Table 51 with fluid properties and boundary conditions given in Table 52. The results can be seen in Figures 53 and 54. All results seen herein are based on a second order upwind scheme, 2UDS, unless otherwise specified. Pressure drops for a range of velocities, and the velocity profiles downstream of the orifice plate have been compared to the experimental results. The pressure drop is plotted against the Darcian Velocity in Figure 53 for three of the four above models; additionally the experimental results are provided. Table 52. Fluid and porous solid properties, along with inlet velocities examined. Fluid and Material Properties Inlet Velocities Fluid, Air @ 24.24 C Case Velocity (m/sec) Density (p) 1.1875 kg/m3 VI 3.58 Dynamic Viscosity ([t) 1.8048e5 kg/ms V2 7.766 Specific Heat (Cp) 1006.2 J/kgK V3 10.543 Thermal Conductivity (k) 0.025913 W/mK V4 13.140 V5 16.301 V6 18.122 V7 20.0795 V8 23.306 V9 25.754 4 x 10 9 Case C1 8  Case C3 SCase C4 Experimental (Carroll and Omar) 7 6 t5 5 10 15 20 25 Darcian Velocity (m/sec) Figure 53. Pressure drop across the orifice plate. Cases Cl, C3, and C4 are labeled in Table 51. As can be seen in Figure 53, the results of Case Cl provided pressure drops much greater than the real values. Additionally, by inspection of the coefficients in Table 51, the coefficients of Case C2 would provide pressure drops much greater than the actual; therefore, the simulations were not run for this case. Cases C3 and C4 provided comparable results; however, the coefficients of Case C4 were slightly better. Additionally, the experimental methods used to determine the permeability and Fochheimer coefficient put forth in Case C3 were slightly flawed. The Reynolds number, used to extract the permeability was much larger than unity, and is the reason for the greater discrepancy. In the case of the RigimeshTM the pore geometry is unknown, making experimental determination of the material constants necessary. To obtain accurate values of these coefficients experiments must be run at lower speeds. As stated above, Case C3 provided the best results and is therefore further scrutinized. The greatest absolute error occurred for flow speed V2 and was approximately 1700 Pa, corresponding to 88%. A similar absolute error occurred for flow speeds V7 and V8; however, because of the drastic pressure drop the percent error was in the range of 69%. The attractiveness of there results lie in the fact that absolute error was a maximum of approximately 1700 Pa, or 0.25 psi. In particular at the high flow speeds the absolute error remained about the same. It is the belief of this author that further experimental investigation and the methods of Case 4 would provide better values of the coefficients and hence, reduced errors. In hopes of reducing the previously discussed error, a quadratic curve fit to the experimental data was applied. From this, the values of permeability and Forchheimer's coefficient (K = 5.74035E9 m2 and CF = 0.487469) were extracted and the corresponding simulations were run. The results compared to Case 3 can be seen in Figure 54. Here error was reduced to a maximum value of 65% occurring for V2. The percent error of all other flow speeds decreased as well. Some of the discrepancy seen here is based on the fact that the curve fit, y = 48.516x2 + 20.363x 364.89, did not intersect the yaxis at zero. This shift was ignored when K and CF were extracted. Furthermore, when the curve fit was forced to intersect the yaxis at zero, negative coefficients resulted; therefore, was not investigated. Regardless, by varying the values of permeability and Forchheimer's coefficient, the error was reduced. Moreover, Forchheimer's constant my only be piecewise constant depending on flow regime. To better understand this further investigation is necessary. As stated before, as more data is collected for investigations herein, the results should improve. 4 x 10 Case 3 3 K & CF from Curve Fit Experimental (Carroll and Omar) 2.5 S2 Q 2 1.5 09 0.5 5 10 15 20 25 Darcian Velocity (m/sec) Figure 54. Pressure drop across the orifice plate. Case C3 and results from K and CF based on a curve fit of the experimental data. The curve fit resulted in the following values, K = 5.74035E9 m2 and CF = 0.487469. In addition to the pressure drop across the plate, velocity profiles one inch downstream of the plate were compared to the experimental results, Figure 55. Here much more discrepancy in the trends can be seen. This is due largely in part to manufacturing defects in the physical porous plate; the hole pattern was not perfect (Ahmed F. Omar, personal communication, Summer 2004). It is expected that the results of the RigimeshTM, which is a more natural porous material, will be a better match for the results of the porous models used herein. Additionally, by the nature of the hole pattern in porous plate, the orifices did not extend completely to the boundary. Essentially, a small layer near the wall was not acting as a porous zone. This may have caused some flow separation downstream of the plate, resulting in the errors seen here. Finally, because of the high Reynolds numbers examined in this work, it is expected that the addition of a turbulence model would help account for the downstream velocity profile. Numerical 0.03 Experimental 0.02 0.01 V1 S0 \ 0.01 0.02 * V9 0.03 0 5 10 15 20 25 30 Axial Velocity (m/sec) Figure 55. Velocity profiles 1 inch downstream of the orifice plate, Case C4. 5.2 Drilled Orifice Plate Non Isothermal Simulations With the inertial and viscous coefficients known, attention can be turned to determining the proper form of the energy equation, local thermal equilibrium (LTE), or local thermal nonequilibrium (LTNE) between the solid and fluid phases. Furthermore, if an LTNE model is appropriate, the proper value of the heat transfer coefficient must be determined. Generally, when the difference between the thermal conductivities of the solid and fluid phases is large, the assumption of thermal equilibrium is insufficient. Another possibility would be to examine the Biot number, which compares the internal thermal resistance of the solid to boundary layer thermal resistance. In the investigations herein, the solid phase conductivity is about 5000 times that of the fluid phase, hence the LTNE effects are expected to be significant. To begin, simulations corresponding to LTE were investigated for the domain seen in Figure 56. 0.05715 m 0.03175 m Cooled Wall, Constant 0.0254 m Temp, 273.15 K T,,, T,,,,, ^^^^^^^M ,,,,. Hot Air H0i 0.0263144 m k, = 0.025913 W/mK 500 K0 k500 K= 110 W/mK kn = 94.1417 W/mK P1 P2 Figure 56. Nonisothermal porous plate setup. Domain of interest and corresponding boundary conditions are as seen, with inlet velocities ranging from V1V9. The effective conductivity was found from equation 2.28 in Chapter 2. It is also provided here for convenience, ke = ( e)k, + Esk. (2.28) Using the solid and fluid phase conductivities provided in Figure 56, kff becomes, 94.1417 W/mK. The results of this model for the same inlet velocities previously discussed (V1V9), can be seen in Figures 57 and 58. Additionally it should be noted that the boundary conditions seen in Figure 56 included an inlet fluid temperature of 500 K, a constant wall temperature of 273.15 K in the porous region, and adiabatic walls in the open channel. Here the temperature profiles on the upstream and downstream faces of the injector are examined. As seen in Figure 57, the cooled wall boundary condition has a greater affect on temperature near the centerline at lower flow speeds, VI as opposed to V9. I , 480  460 v 9 440 420 400 V1 380 1 S360 1 340 320  Upstream Face Downstream Face 300 273.15 0 0.005 0.01 0.015 0.02 0.0263 r (m) Figure 57. Temperature profiles on the upstream and downstream faces of the porous slug. Two test cases are shown here, VI corresponds to an inlet velocity of 3.58 m/sec and V9 to 25.754 m/sec. Similarly, Figure 58A shows the effect of the boundary condition on the temperature profiles of the downstream face for a range of flow velocities. Additionally, the effect of the flow speed on the temperature difference of the upstream and downstream face can be seen in Figure 58B. At high flow speeds the temperature gradient in the axial direction was much greater in the near wall regions as noted by the large temperature difference. However, at the centerline, the axial temperature gradient approaches zero as the flow speed increases. A) 500 Increasing Velocity 0 , 400 S380 360 E  360 0.005 0.015 0.02 0.025 r(m) Increasing Velocity \// Vi 0 V0 V9 0 0.005 0.01 r(m) 0.015 r (m) 0.02 0.0263 Figure 58. Temperature profiles, and AT. A.) Temperature profiles on the downstream face for all inlet velocities tested, V1V9. B.) AT between the upstream and downstream faces for all inlet velocities tested, V1V9. E S40 F E 30 20 II a 20 Next, the effect of LTNE will be examined. The model previously discussed in Chapter 2 can be seen again in Table 53, where the value of asf associated with that particular model has been determined. Table 53. Local thermal nonequilibrium models, w. respective coefficients hsf asf Model Description Description Value (m1) HO Thermal Equilibrium NA NA k(2+ 1.lPr1/3 Re06) 6(1 s) H1 10107.9 d d H2 1.064 Pr33 Re059 20346( 712.9 dp dp H3 dpe + 2 s ) 10107.9 0.2555Prl/3 Re2/3 k 10k, dp Because of the known geometry of the porous plate, the actual value of asf, as opposed to those based on the models, can be determined analytically as follows: a 2(2 2)*L 2 1327.2 (m ), (5.2) sfb 2 2[(di /2)2 b 2 1/2 2 =' where b is defined in Figure 51 and L is the thickness of the porous region. Alternatively asf could be written, A4/ 2na a== p (5.3) SV, b 2(1 ) Based on this value of asf Model H2 seen in Table 53 is the best match. However, for the results discussed herein, a constant value of asf, as in equation 52, will be used. For simplicity, only the results from two flow speeds, VI and V9, on the downstream face of the injector are given here. These can be seen in Figure 59. From this, a basic understanding of the trends can be identified. As expected; the solid and fluid phases were at significantly different temperatures. In all cases the temperature of the solid phases was lower than that when LTE was assumed. Furthermore, Model H3 had the larger heat transfer coefficient; therefore, the fluid and solid phases were close to being in thermal equilibrium. Finally, as experimental data is collected more conclusions will be drawn regarding the conditions for the presence of LTE. 500 4 8 0 h __ 460 All From V9 All From V1 420 360 Fluid Phase, H1 . *'. \ 360 ** Solid Phase, H1 :. ' Fluid Phase, H2 : 340 Solid Phase, H2 a Fluid Phase, H3 S 4 320 Solid Phase, H3 LTE 300 280 0 0.005 0.01 0.015 0.02 0.025 r(m) Figure 59. Temperature profiles for LTNE models on downstream face of the porous slug. Two test cases are shown here, VI corresponds to an inlet velocity of 3.58 m/sec and V9 to 25.754 m/sec. 5.3 Drilled Orifice Plate with Center Jet Dynamic Interaction With the porous models now hydrodynamically verified, a more accurate model of an SSME injector can be explored. To begin Figure 510 shows the assembly of injectors for the SSME. The model used herein is to examine a single element of this assembly. Additionally, the real injector introduces both fuel and oxygen into the combustion chamber as separate streams that mix in the combustion chamber. However, for the investigation here a single fluid, air, will be examined as seen in Figure 511. MAIN INJECTOR ASSEMBLY Spark Fluted oxidizer posts where igniter  hot hydrogen evaporates the oxygen Fuel intet from hot gas manifold Oxygen leCold hydrogen cavity Oxygen inlel manifolds Five compartment baffle with 75 cooled Thrust infection posts load transmitting Primary injection cone plate (transpiration cooled) with 525 main injection elements ignition flame tube Oxygen / from main' oxygen valve Figure 510. Main injector assembly of the Space Shuttle Main Engine (SSME), scanned from Sutton [5]. Image shows baffle with five outer compartments. As with the isothermal porous plate and the nonisothermal porous plate, simulations were run for a range of flow speeds, VI V9. The results given here, Figures 512 through 514, include the temperature profiles on the upstream and downstream faces, the temperature difference across the porous plate, and the percent mass flux through the center jet. However, it should be noted, that the results here are based on a first order upwind scheme, UDS, rather than a 2UDS. The reason for this will be discussed later. 00000 R 1.036" /2.63144 cm oooooo00000000000 0000000000000 00000000000000 00 0000o0o0o0oooo 0o o 00R 0.20" / 0.635 cm 0000000000000 00000 0000000000000 000000 0 0 0 0 0 0 0 00 0 0 00 0 Note: Each individual hole seen here are 000400 00 o0 000oo not representative of each injector 000000000 0 0 0 0 0 0 0 000 0 o element. The coaxial injector is ooooooooo 0 o0000000 represented by the larger diameter "center 00000000000000 0 0 000 jet." The small holes act as pores in the 0oooooooo 0 0o o o 0 obo 0 0 porous face and are not to be confused 000000000000000000 0 00000000000000 o0o00 with an injector element. 000000000000000000 00000000000000000 000000000000000/ Center Jet 000000000000 00000000 0.05715 m S0.03175 m Cooled Wall, Constant 0.0254 m Temp, 273.15 Ko TUpstream TDownstream Hot Air  0Pate kf = 0.025913 W/mK 500 Ko 0.0263144 m ks =110 W/mK keff = 94.1417 W/mK 0.00635 m P1 P2 Figure 511. Domain setup for both the numerical investigations of flow through a drilled orifice plate with a dynamically influence center jet. The drilled hole pattern representing the porous injector face is the same as previously investigated, see Figure 51. As seen in Figure 512, a similar trend as for the porous plate is observed. Likewise, Figure 513A and 513B follow the similar patterns to those previously discussed. One additional comment should be made on the percent mass flow through the centerjet, Figure 514. As the flow speed increased, the percent mass flow through the centerjet decreased. Initially this was somewhat counter intuitive. The percent mass flow through the center jet was expected to rise as the flow speed increased due to the additional resistance to the flow through the porous region. However, the extra effort to 400 ^ 380 V1 E 360 340 Jet Porous Region 320   Upstream Face 300 Downstream Face 273.15 0 0.0032 0.01 0.015 0.02 0.0263 r(m) Figure 512. Temperature profiles on the upstream and downstream injector faces. Two test cases are shown here, VI corresponds to an inlet velocity of 3.58 m/sec and V9 to 25.754 m/sec. drive the flow through the center jet at high speeds was more than that caused by the porous region, making the percent mass flow through the center jet decrease. With the results now presented, attention can be turned to the reason for using a UDS for this particular problem. Initially, this problem was examined using a 2UDS; however, slight overshoots in the fluid temperature in the center jet region were observed, as seen in Figure 515. The UDS, on the other hand, produced results that did not include the temperature overshoot. Currently, the cause of this overshoot cannot be offered; however, more investigation into its origin is underway. A) 500 Increasing Velocity S400 S380 Q360 E  360 340 Jet Porous Region 320 300 h 273.15 0 0.0032 0.015 0.02 r(m) Increasing Velocity v / ,,C, / Jet  Porous Region 0 0.0032 0.01 0.015 0.02 r(m) Figure 513. Temperature Profiles and AT. A.) Temperature profiles on the downstream face for all inlet velocities tested, V1V9. B.) AT between the upstream and downstream faces for all inlet velocities tested, V1V9. 0.0263 60 50 E S40 E 30 II I < 20 0.0263 8.5 0 Figure 514. 500 450 400 E (, I 5 10 15 20 25 Darcian Velocity Percent mass flow through the center jet for flow speeds, V1V9. 0 0.005 0.01 0.015 0.02 0.025 r (m) Figure 515. Results from 2UDS simulations with the centerjet. The unrealistic spike in temperature in the center jet region can be seen at flow speed V1. CHAPTER 6 OBSERVATION AND RECCOMENDATIONS Here a number of recommendations for future studies of transpiration cooled injectors are proposed. Additionally, several findings not previously discussed are presented. The recommendations and lessons learned are by no means complete; however, several observations became apparent to this author in the course of this work and are presented here in hopes of furthering this research. 6.1 Stability Observation using a Central Difference Scheme It is well known that the central differencing scheme can lead to physically impossible solutions at high Peclet Numbers. However, in the cases examined here, the central difference scheme did not converge. The added numerical dissipation of the first and second order upwind schemes helped to keep these schemes stable, while the lack of dissipation in the CDS kept the respective simulations from converging (Dr. Siddharth Thakur, personal communication, November 2004). 6.2 Stability Observations on a Collocated Grid Previously, in Chapter 3, the differences between a staggered and a collocated grid were discussed. Most importantly, a momentum interpolation technique must be implemented when using a collocated grid, whereas, on a staggered grid, the staggering of the variables provides this effect. Prior to the use of a staggered grid, as in this thesis, the same simulations were run on a collocated grid, using the Rhie and Chow momentum interpolation technique, where interface velocity is found from 1 >a ^ m1 m1 u =7AQe 1 [p Kp r (6.1) e a"e x e Jx In equation 6.1 the overbar represents an interpolated value. The results on a collocated grid included spurious oscillations in the regions just upstream and downstream of the porous plate, see Figure 61. Refining the mesh locally near the porous interface resulted in decreased amplitudes of these oscillations; however, this increased the simulation time drastically. 1.5 1.4 1.3 1.2 C2 c 0.9 0.8 0.7 0.6 0 0.01 0.02 Porous 0.04 0.05 0.06 x (m.) Figure 61. Uvelocity at the centerline. The oscillations seen are on a collocated grid, a side effect of the Rhie Chow momentum interpolation technique. The cause of these oscillations is inherent to the Rhie Chow momentum interpolation technique. This method of momentum interpolation essentially couples the velocity and pressure gradient driving the flow by determining the interface velocity. However, when porous source terms are included, as those used in this thesis, the flow is driven by both the pressure gradient the momentum source. In order to model this on a collocated grid, an appropriate momentum interpolation technique to couple all three terms is necessary (Dr. Siddharth Thakur, personal communication, November 2004). 6.3 Recommendations for Future Work The first modification or addition to the current code would be the inclusion of a turbulence model. The existence and modeling of turbulent flows within porous media is still a subject of debate, however, of more importance here would be the turbulent flow downstream of the porous plate. The main effect of a turbulent model would be the change in temperature on the injector face caused by the recirculation of the hot flow near the face especially for the cases with a center jet. The mixing of the low speed fuel bled through the porous material and the high velocity center jet causes eddies and vortices to form, which could impart substantial differences on the temperature profile across the injector face. Subsequently, the exploration of a momentum interpolation technique for use on a collocated grid, capable of handling the added momentum source would be beneficial. Not only from an academic standpoint, but also because the added applicability of multi grid acceleration common to collocated grid. This of course leads to the addition of multi grid capability to the current code. This would decrease simulation time drastically, and allow for refined meshes to be considered as well the capability to efficiently add new models, such as the turbulent models previously discussed. CHAPTER 7 CONCLUSIONS The use of porous models to replace adhoc boundary conditions in previous computational models for liquid propellant rocket engine (LPRE) injectors has been investigated. Several methods of determining porous material constants such as permeability and the socalled Forchheimer coefficient were analyzed. In cases examined within this thesis, the geometry of the pores was known allowing for the exploration of several analytical models. The combination of experimental methods and analytical models resulted in coefficients that provided superlative results. However, as more experimental data is collected, it is the belief of this author that the experimental methods will provide more accurate coefficients. Furthermore, based on the unknown geometry of the pores in most materials, in particular the RigimeshTM used in LPRE injectors, experimental determination is believed to be necessary. However, the experimental determination of the constants can be done on a much smaller scale than any fullfledged experiment. Additionally, although a maximum error of 88% was seen in the pressure drop across the injector, it is the belief of this author, that with further investigation and additional experimental data, the pressure and velocity curves can be finetuned. Additionally, it is important to note that these large percentage errors were seen at low flow speeds. Throughout the range of flow speeds examined, the absolute error was in the range of 400 to 1800 Pa, or roughly 0.06 to 0.26 psi, corresponding to approximately 2% error at the high flow speeds. 75 Of utmost importance to the development of LPRE's are temperature profiles across the face of the injector assembly. Due to time constraints, no comparison to experimental data was provided within this thesis; however, a number of curves based on local thermal equilibrium and local thermal nonequilibrium models have been presented. It is the hope of this author that these curves and developed code can be used to validate these models as experimental results are found, and aid in the development of LPRE injectors. LIST OF REFERENCES [1] Sutton, G. P., "History of Liquid Propellant Rocket Engines in the United States," Journal of Propulsion and Power, Vol. 19, No. 6, 2003, pp. 978 1007. [2] Louden, P., Pratt & Whitney Space Propulsion RL10. Retrieved November 2, 2004, from http://www.pw.utc.com/presskit/factsheets/space_2003_status rl10.pdf [3] Pratt & Whitney, RL10 High Resolution Image. Retrieved November 2, 2004, from http://www.pw.utc.com/presskit/images/rl 10 high.j pg [4] Huzel, D. K., and Huang, D. H., "Moder Engineering for Design of Liquid Propellant Rocket Engines," Vol. 147, Progress in Astronautics and Aeronautics, AIAA, Washington DC, 1992. [5] Sutton, G.P., and Biblarz, O., Rocket Propulsion Elements, 7th ed., Wiley, New York, 2000. [6] Huzel, D. K., and Huang, D. H., Design ofLiquid Propellant Rocket Engines, 2nd ed., NASA SP125, 1971. [7] Brown, C. D., Spacecraft Propulsion, AIAA Education Series, AIAA, Washington DC, 1996. [8] Avenall, R. J., 2004, "The Use of Metallic Foams for Heat Transfer Enhancement in the Cooling Jacket of a Rocket Propulsion Element," Master's Thesis, University of Florida. [9] Kaviany, M., Principles of Heat Transfer in Porous Media, 2nd Edition, Springer Verlag, New York, 1995. [10] Dullien, F. A. L., Porous Media Fluid Transport andPore Structure, Academic Press, Inc., New York, 1979. [11] Whitaker, S., "Advances in Theory of Fluid Motion in Porous Media," Industrial & Engineering Chemistry, Vol. 61, No. 12, 1969, pp 1428. [12] Guin, J. A., Kessler, D. P., Greenkor, R. A., "The Permeability Tensor for Anisotropic Nonuniform Porous Media," Chemical Engineering Science, Vol. 26, 1971, pp 14751478. [13] Alazmi, B., and Vafai, K., "Analysis of Variants Within the Porous Media Transport Models," Journal of Heat Transfer, Vol. 122, 2000, pp. 303 326. [14] Nield, D. A., and Bejan, A., Convection in Porous Media, 2nd Edition, Springer Verlag, New York, 1999 [15] Ferziger, J. H., and Peric, M., Computational Methodsfor Fluid Dynamics, 3rd Edition, SpringerVerlag, Berlin, Germany, 2002. [16] Patankar, S. V., Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC, 1980. [17] Patankar, S. V., and Spalding, D. B., "A Calculation Procedure for Heat, Mass and Momentum Transfer in ThreeDimensional Parabolic Flows," International Journal ofHeat and Mass Transfer, Vol. 15, 1972, pp 17871805. [18] Thakur, S., Wright, J., and Shyy, W., STREAM, A Computational Fluid Dynamics and Heat Transfer NavierStokes Solver Theory and Applications, Version 4.5.2, Retrieved November 2, 2004, from http://aemes.mae.ufl.edu/cfdweb/cgi bin/main.cgi?index=0&altmenu=4 [19] Ghia, U., Ghia, K. N., and Shin, C. T., "HighRe Solutions for Incompressible Flow Using the NavierStokes Equations and a Multigrid Method," Journal of Compuational Physics, Vol. 48, 1982, pp. 387411. [20] Poulikakos, D., and Kazmierczak, M., "Forced Convection in a Duct Partially Filled With a Porous Material," Journal of Heat Transfer, Vol. 109, 1987, pp. 653 662. [21] Incropera, F. P., and DeWitt, D. P., Fundamentals of Heat and Mass Transfer, 5th Edition, Wiley, New York, 2002. BIOGRAPHICAL SKETCH Landon Rothwell Tully was born August 29, 1980, in Abington, Pennsylvania. He graduated from Cypress Lake High School, Fort Myers, Florida, in 1998. He attended the University of Florida and received a Bachelor of Science, with honors, in mechanical engineering in the fall of 2002. Since then he has been pursuing a Master of Science degree in mechanical engineering while working as a graduate research assistant. 