Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0009260/00001
## Material Information- Title:
- Numerical Modeling of Transpiration Cooled Rocket Injectors
- Creator:
- Tully, Landon Rothwell
- Place of Publication:
- Gainesville, Fla.
- Publisher:
- University of Florida
- Copyright Date:
- 2005
- Language:
- English
- Physical Description:
- online resource
## Subjects- Subjects / Keywords:
- Cooling ( jstor )
Fluids ( jstor ) Geometry ( jstor ) Inlets ( jstor ) Momentum ( jstor ) Porous materials ( jstor ) Porous plates ( jstor ) Pressure reduction ( jstor ) Simulations ( jstor ) Velocity ( jstor ) - Genre:
- Academic theses ( fast )
Academic theses ( lcgft )
## Notes- Abstract:
- 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 ad-hoc 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. ( ,,,,,, )
- Abstract:
- 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.
- General Note:
- Includes vita.
- General Note:
- Includes bibliographical references.
- General Note:
- Title from title page of source document.
- General Note:
- Index terms: Darcy; Forchheimer ; injector ; LPRE; porous ; rocket.
- Statement of Responsibility:
- by Landon Rothwell Tully
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Embargo Date:
- 5/1/2005
- Resource Identifier:
- 436098657 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads:
tully_l ( .pdf )
tully_l_Page_61.txt tully_l_Page_17.txt tully_l_Page_40.txt tully_l_Page_07.txt tully_l_Page_93.txt tully_l_Page_51.txt tully_l_Page_73.txt tully_l_Page_86.txt tully_l_Page_34.txt tully_l_Page_69.txt tully_l_Page_27.txt tully_l_Page_92.txt tully_l_Page_37.txt tully_l_Page_79.txt tully_l_Page_01.txt tully_l_Page_66.txt tully_l_Page_25.txt tully_l_Page_36.txt tully_l_Page_72.txt tully_l_Page_49.txt tully_l_Page_75.txt tully_l_Page_38.txt tully_l_Page_28.txt tully_l_Page_62.txt tully_l_Page_52.txt tully_l_Page_94.txt tully_l_Page_43.txt tully_l_Page_45.txt tully_l_Page_50.txt tully_l_Page_24.txt tully_l_Page_47.txt tully_l_Page_85.txt tully_l_Page_02.txt tully_l_Page_41.txt tully_l_Page_14.txt tully_l_Page_67.txt tully_l_Page_06.txt tully_l_Page_65.txt tully_l_Page_70.txt tully_l_Page_88.txt tully_l_Page_22.txt tully_l_Page_81.txt tully_l_Page_64.txt tully_l_Page_19.txt tully_l_Page_12.txt tully_l_Page_04.txt tully_l_Page_63.txt tully_l_Page_08.txt tully_l_Page_35.txt tully_l_Page_23.txt tully_l_Page_59.txt tully_l_Page_13.txt tully_l_Page_44.txt tully_l_Page_77.txt tully_l_Page_39.txt tully_l_Page_20.txt tully_l_Page_89.txt tully_l_Page_87.txt tully_l_Page_26.txt tully_l_Page_05.txt tully_l_Page_82.txt tully_l_Page_57.txt tully_l_Page_21.txt tully_l_Page_48.txt tully_l_Page_53.txt tully_l_Page_76.txt tully_l_Page_80.txt tully_l_Page_91.txt tully_l_Page_29.txt tully_l_Page_11.txt tully_l_Page_54.txt tully_l_Page_83.txt tully_l_Page_16.txt tully_l_Page_32.txt tully_l_Page_31.txt tully_l_Page_10.txt tully_l_Page_74.txt tully_l_Page_33.txt tully_l_Page_15.txt tully_l_Page_55.txt tully_l_Page_78.txt tully_l_Page_60.txt tully_l_Page_56.txt tully_l_Page_03.txt tully_l_Page_68.txt tully_l_Page_09.txt tully_l_Page_42.txt tully_l_Page_58.txt tully_l_Page_84.txt tully_l_Page_46.txt tully_l_Page_18.txt tully_l_Page_90.txt tully_l_Page_71.txt tully_l_pdf.txt tully_l_Page_30.txt |

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 2-1. Properties of common porous materials. ............................................................16 2-2. Heat transfer coefficients with respective fluid to solid specific surface area. .........26 5-1. Values of permeability and Forchheimer's coefficient. .........................................55 5-2. Fluid and porous solid properties, along with inlet velocities examined. ................56 5-3. Local thermal non-equilibrium models, w. respective coefficients.........................64 LIST OF FIGURES Figure p 1-1. Pratt and W hitney's R L 10 LPR E ..................................................... ..................... 3 1-2. Cutaway of a regeneratively cooled tubular thrust chamber .............. ............ 5 1-3. N on-im pinging injector elem ents .............. ......... ..................................................6 1-4. U nlike-im pinging injector elem ents ........................................ ........................ 7 1-5. Like-impinging injector elements ................. ..................... ................. 8 1-6. Steady state cooling m ethods .................................. ...................................... 9 2-1. Flow schem atic for D arcy's law .......................................... .......................... 12 2-2. Illustration of capillary model for a bundle of identical capillaries ........................17 2-3. Transition from the Darcy regime to the Forchheimer regime in unidirectional flow through an isothermal saturated porous medium ............................................. 20 2-4. A schematic of a representative elementary volume and the position vectors used..24 3-1. Generic grid setup for staggered arrangement of variables .................................28 3-2. G rid setup, w ith notation. ........................................ ............................................29 3-3. Outline of SIMPLE algorithm ............................. ............... 32 3-4. Geometry of the problem under investigation with labels corresponding to input variables in the code. ......................... ...... ...................... ........... ............. 37 3-5. Sample grid displaying the important characteristics of the grid generation process.38 4-1. Benchmark Case 1, lid driven cavity flow, problem setup....................................40 4-2. Benchmark Case 1, discretization scheme comparison on a 20 x 20 grid for R e = 1 0 0 0 ......................................................................... 4 1 4-3. Benchmark Case 1, discretization scheme comparison on a 120 x 120 grid for R e = 1 0 0 0 ......................................................................... 4 2 4-4. Benchmark Case 2, pressure driven flow in a pipe, problem setup.........................43 4-5. Benchmark Case 2, results comparison........................................... ...............46 4-6. Benchmark Case3, forced convection in a duct partially filled with a porous m material, problem setup .................. .......................... ...................... 47 4-7. Benchmark Case 3, comparison of velocity profiles to the analytical results of P oulikakos [20]. .................................................... ................. 50 4-8. Benchmark Case 3, grid with refined region at the porous/open channel interface..50 4-9. Benchmark Case 3, results of local grid refinement, with 80 nodes in the radial direction ......... ...... ............ ..................................... ...........................5 1 4-10. Benchmark Case 3, validation of solution to the energy equation........................52 5-1. Domain setup for both the experimental and numerical investigations of flow through a drilled orifice plate ................................. ............... ............... 54 5-2. Grid setup corresponding to the geometry seen in Figure 5-1 ................................56 5-3. Pressure drop across the orifice plate ............................................. ............... 57 5-4. Pressure drop across the orifice plate. ............................................ ............... 59 5-5. Velocity profiles 1 inch downstream of the orifice plate, Case C4 .........................60 5-6. N on-isotherm al porous plate setup .................................................. ..... .......... 61 5-7. Temperature profiles on the upstream and downstream faces of the porous slug.....62 5-8. Tem perature profiles, and AT ...................................................................... 63 5-9. Temperature profiles for LTNE models on downstream face of the porous slug .....65 5-10. Main injector assembly of the Space Shuttle Main Engine (SSME)....................66 5-11. Domain setup for both the numerical investigations of flow through a drilled orifice plate with a dynamically influence center jet .....................................67 5-12. Temperature profiles on the upstream and downstream injector faces ................68 5-13. Tem perature Profiles and AT...................................... ............... ............... 69 5-14. Percent mass flow through the center jet for flow speeds, V1-V9 .....................70 5-15. Results from 2UDS simulations with the centerjet ............................................. 70 6-1. U-velocity 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 x-direction or axial direction y y-direction or radial direction z z-direction 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 ad-hoc 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 hot-firing 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 jet-assisted take off (JATO) engines. Several of these were successfully flown in the F-84 fighter bomber, the PB2Y-3 flying boat, and the B-29, B-45, and B-47 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 upper-stage 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 stainless-steel 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 1-1. 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 [4-7]. 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 -` Cuh-n> stilinaf binds i tjector plate U=rated 180t iWt ThMat .^a^ QwnW Figure 1-2. 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: [4-7]. sitratiIle r ,-D"rm lu lredt w fueled Farid Id' - Fuel cooled tubulal mil Tlrn-larwrl 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: non-impinging, unlike-impinging, and like impinging. Non-impinging elements inject the fuel and oxidizer separately into the combustion chamber and include: showerhead, fan former, and coaxial elements as seen in Figure 1-3. 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 1-3. Non-impinging injector elements, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. Unlike-impinging 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 1-4. Injector face Injector face FIE IFuel Oxidizer Manifold Oxidizer Manifold Fuel Manifold Fuel Manifold (A) Unlike Doublet (B) Unlike Triplet Figure 1-4. Unlike-impinging injector elements, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. Like-impinging elements are very similar to the unlike-impinging 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 1-5. 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 1-5. Like-impinging 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 1-6. Figure 1-6. 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: [4-7]. 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 well-known 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 2-1) 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 cross-sectional area of the Figure 2-1. 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 (centi-Poise) 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 K-13 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 2-1. Table 2-1. Properties of common porous materials. Porous Material Porosity Permeability (m2) Foam Metal (also made of other materials) 0.98 Fiberglas 0.88-0.93 2.4E-11 5.1E-11 Berl saddles 0.68-0.83 1.3 E-07 3.9 E-07 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.56-0.59 9.5 E-14 1.2 E -13 Granular crushed rock 0.44 0.45 Soil 0.43-0.54 2.9E-13 1.4E-11 Sand 0.37-0.50 2.0E-11 1.8E-10 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.12-0.34 4.8 E -15 2.2 E -13 Sandstone (oil sand) 0.08-0.38 5.0E -16 3.0E -12 Limestone, dolomite 0.04-0.10 2.0 E -15 4.5 E -14 Coal 0.02-0.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 Navier-Stokes equations for internal flows. One model is based on laminar flow through a bundle of circular tubes of equal length and diameter D (Figure 2-2), in which the corresponding porosity F, is nid p2 /4 for n tubes per unit area. Capillaries / Conduits (DI I SOO o -0 Figure 2-2. Illustration of capillary model for a bundle of identical capillaries. The permeability can be derived beginning with the well known Hagen-Poiseuille 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 cross-sectional 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 high-speed 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 = u-D -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 2-3. 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]. N-H-^r-T-- .....t.-H.-- --h+.. fK=------ ,ReK fK 1 +0.55 S.1 ---__ ReK I A I II I [ I [ll- __---.-t--1-[-44 1---.--. .- 0.1 0.1 1 10 100 vK1/2 RCK=- Figure 2-3. 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 2-3, 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 1-5.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 Carman-Kozeny 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 Navier-Stokes equation. The first attempt at this was by Brinkman, who included a term analogous to the Laplacian term of the Navier-Stokes 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 Navier-Stokes 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 2-4. 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 2-2. 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 2-2. 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(1-l 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 Navier-Stokes 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 pressure-correction method for solving the Navier-Stokes equations on a Cartesian or axi-symmetric 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 3-1. Internal Node Location (Center of CV) Figure 3-1. 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 3-1. 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 3-2 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 (i-1 ,j) (i,j) (i+1 ,j) Yi-1 sw s se S y o o o (i,j-1) L Yi-2 Xi-2 Xi-1 Xi Xi+1 S x Figure 3-2. 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 v-momentum 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 u-momentum equations only and can be easily extended for use in the v-momentum 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 )1-1 Fe = u +h[eU -ue ) (3.7) where the explicit terms are denoted by the superscript "m-l", 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 lower-upper (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 Semi-Implicit Method for Pressure-Linked Equations. This procedure was originally put forth by Patankar and Spalding [17]. As a reference, and introduction, the reader is referred to Figure 3-3, the details of which will be discussed subsequently. Yes Set calculated T as Solve Discretized Energy "guessed" field Equation -No Convergence Yes Stop Figure 3-3. 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 3-1. 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 (P-Pp) (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 3-4. 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 sub-region. First, the number of control volumes in the x-direction is needed. Once the number of CV's for each sub-region (entry, porous, and exit) has been entered, the y grid lines are generated at each x-location. 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 y-direction. A sample grid can be seen in the Figure 3-5. 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 x-direction (m.) Figure 3-5. 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 well-known analytical solution for fully developed pipe flow. This test case was used to validate the axi-symmetric 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 well-known 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 4-1. Moving Lid UL Figure 4-1. 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 U-H pfULH The following Figures 4-2 and 4-3, 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 4-2 while the results for a 120 x 120 grid can be seen in Figure 4-3. 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 u-velocity (m/sec) 1st Order Upwind Central Difference 2nd Order Upwind o Ghia etal. (1982) x(m) Figure 4-2. Benchmark Case 1, discretization scheme comparison on a 20 x 20 grid for Re = 1000. A) U-Velocity through y-geometric. B) V-Velocity 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 u-velocity (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 4-3. Benchmark Case 1, discretization scheme comparison on a 120 x 120 grid for Re = 1000. A) U-Velocity through y-geometric. B) V-Velocity 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 4-4. This problem has several well-known 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 hydro-dynamically 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 4-4 .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 x-direction, because the flow was given ample room to become fully developed. Consequently, only the number of nodes in the y-direction will be noted. The results with 22 nodes in the y- direction can be seen in Figure 4-5, 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.2E-5, 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 u-velocity (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 4-5. 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 4-6. SP,:,rous Media Un SOpen Channel Figure 4-6. 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 sDa-Y2 ) 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 x-axis. 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 4-7. As Figure 4-7 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 4-8. 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 4-9. 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 u-velocity (Dimensionless) Figure 4-7. 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 x-direction (m.) Figure 4-8. 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 4-10A. Respectively, the results of the code under investigation can be seen in Figure 4-10B. 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 u-velocity (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 4-9. 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=0-5- s=0.5 J ./, -~-- 10" 10-5 10-4 10-3 Darcy Number (Da) Figure 4-10. 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. 10-2 10-1 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 well-known 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 5-1. 54 R 1.036" 2.63144 cm 000000000000 0000000000000000 ooo00000000000 000000 00000000000000000 Hole-Pattern / 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 / -- i--w^ ^ ^ ^^ ^ ^ ^^ ^ ^ Inlet Porous 7- 0.0263144 m, Velocity p i n n t p Plate -- - ---- Pi P2 Figure 5-1. 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 5-1 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 5-1. 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 5-1. Values of permeability and Forchheimer's coefficient. Permeability K (m2) Forchheimer Coefficient CF Case Model Value Model Value ed 2 d Cl P 1.16324E-9 m2 0.55 1-5.5 0.520801 32 De d2E 1.8 C2 p 5.87529E-12 m2 18 2.44905 180(1- c)2 1801/2 3 C3 Experimental 7.71735E-10 m2 Experimental and 0.160154 K from C3 C4 P 1.16324E-9 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 x-nodes, 40 in the porous region, and 120 y- nodes, resulting in the mesh seen in Figure 5-2. 0.025 0.02 0.015 0.005 0.oo------------ ------------- 0 0 0.01 0.02 0.03 0.04 0.05 x-direction (m) Figure 5-2. Grid setup corresponding to the geometry seen in Figure 5-1. 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 5-1 with fluid properties and boundary conditions given in Table 5-2. The results can be seen in Figures 5-3 and 5-4. 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 5-3 for three of the four above models; additionally the experimental results are provided. Table 5-2. 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.8048e-5 kg/m-s V2 7.766 Specific Heat (Cp) 1006.2 J/kg-K V3 10.543 Thermal Conductivity (k) 0.025913 W/m-K 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 5-3. Pressure drop across the orifice plate. Cases Cl, C3, and C4 are labeled in Table 5-1. As can be seen in Figure 5-3, the results of Case Cl provided pressure drops much greater than the real values. Additionally, by inspection of the coefficients in Table 5-1, 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 6-9%. 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.74035E-9 m2 and CF = 0.487469) were extracted and the corresponding simulations were run. The results compared to Case 3 can be seen in Figure 5-4. 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 y-axis at zero. This shift was ignored when K and CF were extracted. Furthermore, when the curve fit was forced to intersect the y-axis 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 5-4. 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.74035E-9 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 5-5. 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 5-5. 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 non-equilibrium (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 5-6. 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/m-K 500 K0 k500 K= 110 W/m-K k-n = 94.1417 W/m-K P1 P2 Figure 5-6. Non-isothermal porous plate setup. Domain of interest and corresponding boundary conditions are as seen, with inlet velocities ranging from V1-V9. 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 5-6, kff becomes, 94.1417 W/m-K. The results of this model for the same inlet velocities previously discussed (V1-V9), can be seen in Figures 5-7 and 5-8. Additionally it should be noted that the boundary conditions seen in Figure 5-6 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 5-7, 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 5-7. 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 5-8A 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 5-8B. 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 5-8. Temperature profiles, and AT. A.) Temperature profiles on the downstream face for all inlet velocities tested, V1-V9. B.) AT between the upstream and downstream faces for all inlet velocities tested, V1-V9. 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 5-3, where the value of asf associated with that particular model has been determined. Table 5-3. Local thermal non-equilibrium 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 dp--e + 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 5-1 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 5-3 is the best match. However, for the results discussed herein, a constant value of asf, as in equation 5-2, 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 5-9. 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 5-9. 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 hydro-dynamically verified, a more accurate model of an SSME injector can be explored. To begin Figure 5-10 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 5-11. 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 5-10. 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 non-isothermal porous plate, simulations were run for a range of flow speeds, VI V9. The results given here, Figures 5-12 through 5-14, 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/m-K 500 Ko 0.0263144 m ks =110 W/m-K keff = 94.1417 W/m-K 0.00635 m P1 P2 Figure 5-11. 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 5-1. As seen in Figure 5-12, a similar trend as for the porous plate is observed. Likewise, Figure 5-13A and 5-13B follow the similar patterns to those previously discussed. One additional comment should be made on the percent mass flow through the centerjet, Figure 5-14. 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 5-12. 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 5-15. 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 5-13. Temperature Profiles and AT. A.) Temperature profiles on the downstream face for all inlet velocities tested, V1-V9. B.) AT between the upstream and downstream faces for all inlet velocities tested, V1-V9. 0.0263 60 50 E S40 E 30 II I- < 20 0.0263 8.5 0 Figure 5-14. 500- 450- 400- E (, I- 5 10 15 20 25 Darcian Velocity Percent mass flow through the center jet for flow speeds, V1-V9. 0 0.005 0.01 0.015 0.02 0.025 r (m) Figure 5-15. 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 ^ m-1 m-1 u =7-AQe 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 6-1. 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 C-2 c 0.9 0.8 0.7 0.6 0 0.01 0.02 Porous 0.04 0.05 0.06 x (m.) Figure 6-1. U-velocity 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 ad-hoc 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 so-called 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 full-fledged 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 fine-tuned. 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 non-equilibrium 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 SP-125, 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 14-28. [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 1475-1478. [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, Springer-Verlag, 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 Three-Dimensional Parabolic Flows," International Journal ofHeat and Mass Transfer, Vol. 15, 1972, pp 1787-1805. [18] Thakur, S., Wright, J., and Shyy, W., STREAM, A Computational Fluid Dynamics and Heat Transfer Navier-Stokes 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., "High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method," Journal of Compuational Physics, Vol. 48, 1982, pp. 387-411. [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. |

Full Text |

PAGE 1 NUMERICAL MODELING OF TRANSPIRA TION COOLED ROCKET INJECTORS By LANDON ROTHWELL TULLY A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2005 PAGE 2 Copyright 2005 by Landon Rothwell Tully PAGE 3 In loving memory of Randolph R. Tully, Sr. PAGE 4 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 NASAs 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 NASAs 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 iv PAGE 5 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. v PAGE 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................iv LIST OF TABLES ...........................................................................................................viii LIST OF FIGURES ...........................................................................................................ix NOMENCLATURE .........................................................................................................xii ABSTRACT .......................................................................................................................xv CHAPTER 1 INTRODUCTION LIQUID PROPELLANT ROCKET ENGINES..........................1 1.1 Liquid Propellant Rocket History..........................................................................1 1.2 Liquid Propellant Rocket Engine Components.....................................................4 1.2.1 Injectors.......................................................................................................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 Modeling Efforts.......................................................................................12 2.2 Permeability and Porosity....................................................................................14 2.2.1 Porosity......................................................................................................14 2.2.2 Permeability...............................................................................................14 2.2.3 Permeability Models..................................................................................16 2.2.3.1 Capillary models.............................................................................17 2.2.3.2 Drag models....................................................................................18 2.2.3.3 Experimental methods.....................................................................18 2.3 High Reynolds Number Flow..............................................................................19 2.3.1 Forchheimer models..................................................................................21 2.3.2 Experimental methods...............................................................................22 2.4 Semi Heuristic Equations....................................................................................22 2.4.1 Semiheuristic Continuity and Momentum Equations................................23 2.4.2 Energy Equation........................................................................................25 vi PAGE 7 3 FINITE VOLUME METHOD, SIMPLE ALGORITHM, & GRID GENERATION.27 3.1 Algorithm Overview............................................................................................27 3.2 Discretization of the Governing Equations..........................................................28 3.3 SIMPLE Method..................................................................................................32 3.4 Discretization Schemes........................................................................................35 3.4.1 First Order Upwind Scheme......................................................................35 3.4.2 Central Difference Scheme........................................................................35 3.4.3 Second Order Upwind Scheme.................................................................36 3.5 Grid Generation...................................................................................................36 4 BENCHMARK TEST CASES / CODE VALIDATION............................................39 4.1 Benchmark Case 1, Lid Driven Cavity................................................................39 4.2 Benchmark Case 2, Developing Pipe Flow.........................................................43 4.3 Benchmark Case 3, Forced Convection in a Duct Partially Filled w ith a Porous Material..................................................................................................................47 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 Recommendations for Future Work....................................................................73 7 CONCLUSIONS..........................................................................................................74 LIST OF REFERENCES...................................................................................................76 BIOGRAPHICAL SKETCH.............................................................................................78 vii PAGE 8 LIST OF TABLES Table page 2-1. Properties of common porous materials....................................................................16 2-2. Heat transfer coefficients with respective fluid to solid specific surface area..........26 5-1. Values of permeability and Forchheimers coefficient.............................................55 5-2. Fluid and porous solid properties, along with inlet velocities examined..................56 5-3. Local thermal non-equilibrium models, w. respective coefficients...........................64 viii PAGE 9 LIST OF FIGURES Figure page 1-1. Pratt and Whitneys RL10 LPRE................................................................................3 1-2. Cutaway of a regeneratively cooled tubular thrust chamber.......................................5 1-3. Non-impinging injector elements................................................................................6 1-4. Unlike-impinging injector elements............................................................................7 1-5. Like-impinging injector elements................................................................................8 1-6. Steady state cooling methods......................................................................................9 2-1. Flow schematic for Darcys law................................................................................12 2-2. Illustration of capillary model for a bundle of identical capillaries..........................17 2-3. Transition from the Darcy regime to the Forchheimer regime in unidirectional flow through an isothermal saturated porous medium.....................................................20 2-4. A schematic of a representative elementary volume and the position vectors used..24 3-1. Generic grid setup for staggered arrangement of variables.......................................28 3-2. Grid setup, with notation...........................................................................................29 3-3. Outline of SIMPLE algorithm...................................................................................32 3-4. Geometry of the problem under investigation with labels corresponding to input variables in the code.................................................................................................37 3-5. Sample grid displaying the important characteristics of the grid generation process.38 4-1. Benchmark Case 1, lid driven cavity flow, problem setup........................................40 4-2. Benchmark Case 1, discretization scheme comparison on a 20 x 20 grid for Re = 1000.................................................................................................................41 4-3. Benchmark Case 1, discretization scheme comparison on a 120 x 120 grid for Re = 1000.................................................................................................................42 ix PAGE 10 4-4. Benchmark Case 2, pressure driven flow in a pipe, problem setup...........................43 4-5. Benchmark Case 2, results comparison.....................................................................46 4-6. Benchmark Case3, forced convection in a duct partially filled with a porous material, problem setup............................................................................................47 4-7. Benchmark Case 3, comparison of velocity profiles to the analytical results of Poulikakos [20]........................................................................................................50 4-8. Benchmark Case 3, grid with refined region at the porous/open channel interface..50 4-9. Benchmark Case 3, results of local grid refinement, with 80 nodes in the radial direction....................................................................................................................51 4-10. Benchmark Case 3, validation of solution to the energy equation..........................52 5-1. Domain setup for both the experimental and numerical investigations of flow through a drilled orifice plate...................................................................................54 5-2. Grid setup corresponding to the geometry seen in Figure 5-1..................................56 5-3. Pressure drop across the orifice plate........................................................................57 5-4. Pressure drop across the orifice plate........................................................................59 5-5. Velocity profiles 1 inch downstream of the orifice plate, Case C4...........................60 5-6. Non-isothermal porous plate setup............................................................................61 5-7. Temperature profiles on the upstream and downstream faces of the porous slug.....62 5-8. Temperature profiles, and T....................................................................................63 5-9. Temperature profiles for LTNE models on downstream face of the porous slug.....65 5-10. Main injector assembly of the Space Shuttle Main Engine (SSME).......................66 5-11. Domain setup for both the numerical investigations of flow through a drilled orifice plate with a dynamically influence center jet...............................................67 5-12. Temperature profiles on the upstream and downstream injector faces .................68 5-13. Temperature Profiles and T...................................................................................69 5-14. Percent mass flow through the center jet for flow speeds, V1-V9..........................70 5-15. Results from 2UDS simulations with the center jet................................................70 x PAGE 11 6-1. U-velocity at the centerline........................................................................................72 xi PAGE 12 NOMENCLATURE A cross sectional area a coefficient in the discretization equation C F Forchheimer coefficient c P specific heat Da Darcy number D e defined by equation 2.16 d p particle diameter F convective flux g acceleration due to gravity H height in Benchmark Case 1 h height of channel h sf convective heat transfer coefficient I1 and I2 interpolation factors for 2UDS I o ,I 1 modified Bessel functions of the first kind K, K permeability scalar and tensor K 0 K 1 modified Bessel functions of the second kind k thermal conductivity L length of porous slug, or length in Benchmark Case 1 m mass flow rate n number of capillaries P 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 R o radius Re d Reynolds number based on capillary diameter Re H Reynolds number based on cavity height, Benchmark Case 1 KRe Reynolds number based on permeability S momentum source term S surface area in discretized equations s radius at porous interface T temperature T m mean temperature T s surface temperature, wall temperature q s heat flux at surface xii PAGE 13 u x or axial velocity component u' x or axial velocity correction to conserve continuity u D x darcian velocity component u D darcian velocity vector u L lid velocity in Benchmark Case 1 u m mean velocity V volume v y or radial velocity component v' y or radial velocity correction to conserve continuity v D y darcian velocity component w width of channel w D z darcian velocity component x x-direction or axial direction y y-direction or radial direction z z-direction or distance from datum level Greek Symbols volume of CV dynamic viscosity effective dynamic viscosity, here = density porosity any scalar of vector quantity interpolation factor kinematic viscosity hydrodynamic boundary layer thickness T 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 xiii PAGE 14 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 xiv PAGE 15 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 ad-hoc 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 authors knowledge, all prior studies have been heuristic. The focus of this thesis is on xv PAGE 16 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. xvi PAGE 17 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 hot-firing 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 jet-assisted take off (JATO) engines. Several of these were successfully flown in the F-84 fighter bomber, the PB2Y-3 flying boat, and the B-29, B-45, and B-47 1 PAGE 18 2 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/LH 2 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&Ws most successful LPRE [1]. Deemed by P&W as the most reliable, safe and high performing upper-stage 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 stainless-steel 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. PAGE 19 3 Figure 1-1. Pratt and Whitneys 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. PAGE 20 4 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 [4-7]. 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 injectors success. To date, most designs have been empirical. The primary focus of this thesis is the thermal PAGE 21 5 Figure 1-2. 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: [4-7]. PAGE 22 6 Throughout the history of LPREs, a number of injector elements have been used. These are commonly split into three categories: non-impinging, unlike-impinging, and like impinging. Non-impinging elements inject the fuel and oxidizer separately into the combustion chamber and include: showerhead, fan former, and coaxial elements as seen in Figure 1-3. The coaxial elements are similar to those mentioned previously for P&Ws RL10 LPRE. Fuel Ox. Fuel Ox. Liquid Oxygen Gaseous Hydrogen Spacer Outer Sleeve Inner Sleeve Injector face Injector face Oxidizer Manifold Fuel Manifold(A) Shower Head(B) Spray or Fan Former(C) Hollow Post and Sleeve Element (Coaxial) adopted from [Sutton & Biblarz] Injector face Figure 1-3. Non-impinging injector elements, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. PAGE 23 7 Unlike-impinging 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 1-4. Fuel Ox. Fuel Ox. Fuel Ox. Fuel Ox. Injector face(A) Unlike Doublet(B) Unlike Triplet Oxidizer Manifold Fuel Manifold Oxidizer Manifold Fuel Manifold Injector face Figure 1-4. Unlike-impinging injector elements, recreated from Huzel [4], Sutton [5], Huzel [6], and Brown [7]. Like-impinging elements are very similar to the unlike-impinging 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 1-5. One common aspect to all types of injector elements is the injector face. In the case of the SSME and the P&Ws 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. PAGE 24 8 Fuel Ox. Fuel Ox. Injector faceLike Doublet (Self Impinging) Oxidizer Manifold Fuel Manifold Figure 1-5. Like-impinging 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 PAGE 25 9 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 1-6. Figure 1-6. Steady state cooling methods. Left: regenerative cooling jacket. Right: regenerative cooling jacket with nozzle extension for radiation cooling. PAGE 26 10 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: [4-7]. 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&Ws RL10 the injector face is made of a sintered stainless steel material, Rigimesh. The ultimate goal of this study is to obtain a method of acquiring the appropriate material properties to accurately model flow through the Rigimesh. To do so, the methods must first be applied to a porous material with well-known 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 PAGE 27 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. PAGE 28 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 Darcys investigations of hydrology (Figure 2-1) 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 Darcys Law, LKAQf P (2.1) 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 cross-sectional area of the UDUD P1P2 Porous Bed L Figure 2-1. Flow schematic for Darcys law, equation 2.1, recreated from Kaviany [9]. 12 PAGE 29 13 channel normal to the mean flow direction, L is the length of the porous slug in the flow direction, and lastly, P is the piezometric pressure which is defined as: gzp P (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 p P Furthermore, it should also be noted that Darcys law is commonly generalized in the following form: DfpuK (2.3) where K is a second order tensor which, again, is discussed in detail in the following section. Additionally, u D is the Darcian velocity vector or filter velocity vector. The Darcy velocity was originally defined by, A muD (2.4) thus, it is the average velocity of the flow at a given cross section [10]. Since Henry Darcys 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, PAGE 30 14 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, 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, VVf (2.5) where and the subscripts f and s represent the fluid and solid phases respectively. sfVVV 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 PAGE 31 15 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 length 2 or the unit of the darcy, which was created for practicality [10]. In the words of Dullien, a porous material has permeability equal to 1 darcy if a pressure difference of 1 atm will produce a flow rate of 1 cm 3 /sec of a fluid with 1 cP (centi-Poise) viscosity through a cube having sides 1 cm in length [10:78], hence, 21222310*987.0987.01 11 sec1 1mmcmatmcmcPcmdarcy In the case of anisotropic porous media, the permeability is represented by a second order tensor that has nine components as follows, 333231232221131211KKKKKKKKKK (2.6) 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., 332313232212131211KKKKKKKKKK (2.7) 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]. Darcys law then becomes, PAGE 32 16 DxfuKdxdp (2.8a) DyfKdydpv (2.8b) DzfwKdzdp (2.8c) where K x K y and K z 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 2-1. Table 2-1. Properties of common porous materials. Porous Material Porosity Permeability (m 2 ) Foam Metal ( also made of other materials ) 0.98 Fiber g las 0.88 -0.93 2.4E-11 5.1E-11 Berl saddles 0.68 -0.83 1.3 E -07 3.9E-07 Wire crim p s 0.68 -0.76 Black slate p owde r 0.57 -0.66 4.9 E -14 1.2E-13 Raschi g rin g s 0.56 -0.65 Leathe r 0.56 -0.59 9.5 E -14 1.2E-13 Granular crushed roc k 0.44 0.45 Soil 0.43 -0.54 2.9 E -13 1.4E-11 San d 0.37 -0.50 2.0 E -11 1.8E-10 Silica p owde r 0.37 -0.49 1.3 E -14 5.1E-14 S p herical p ackin g s well shaken 0.36 -0.43 Ci g arette filters 0.17 -0.49 1.1 E -09 Bric k 0.12 -0.34 4.8 E -15 2.2E-13 Sandstone ( oil sand ) 0.08 -0.38 5.0 E -16 3.0E-12 Limestone dolomite 0.04 -0.10 2.0 E -15 4.5E-14 Coal 0.02 -0.12 Concrete ( ordinar y 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, PAGE 33 17 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 Rigimesh 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 Navier-Stokes equations for internal flows. One model is based on laminar flow through a bundle of circular tubes of equal length and diameter D (Figure 2-2), in which the corresponding porosity is 42pdn for n tubes per unit area. Capillaries / Conduits dp dp Figure 2-2. Illustration of capillary model for a bundle of identical capillaries. The permeability can be derived beginning with the well known Hagen-Poiseuille equation, dxdpdufpm322 (2.9) PAGE 34 18 where u m is the mean velocity in the pore. Subsequently, the mean pore velocity can be related to the Darcian velocity by dxdpdnuufpmD1284 (2.10) Then, when Equation 2.10 is compared to Equation 2.3, the permeability is found to be 322pdK (2.11) Similar models have been proposed for ducts of varying cross-sectional 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 Darcys 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 PAGE 35 19 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 Darcys 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 Darcys law and further investigation may be necessary [10]. 2.3 High Reynolds Number Flow Darcys law, as discussed above, is only applicable to low speed flows, i.e. u D 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, )1(ORefpDfddu (2.12) where d p 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 fDfKKuRe (2.13) Based on this definition, Darcys law is valid for KRe up to approximately 0.2. For the purposes of this report, unless otherwise specified, Re d will be used. The limited applicability of Darcys law is based on the fact that it only accounts for the viscous resistance of the flow. For high-speed 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 Forchheimers equation, PAGE 36 20 DDfFDfCpuuKuK (2.14) where the C F is the so called Forchheimer coefficient. The transition from the Darcy (viscous) regime to the Forchheimer (inertial) regime based on KRe is illustrated in Figure 2-3. 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]. Figure 2-3. 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 2-3, the divergence from linearity begins at KRe 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 (u D 2 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 PAGE 37 21 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 C F 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, epFDdC5.5155.0 (2.15) where, hwwhDe2 (2.16) with w and h being the width and height of the channel, respectively. In the case of a cylinder D e is simply the diameter. Other models, beyond Forchheimers extension, have also been proposed. One such example is an extension of the Carman-Kozeny theory, for details see Kaviany [9] and Nield [14], where the fluid and matrix properties are related to the pressure drop, by 23322118.11180DfpDpfududdxdp (2.17) Once again, d p 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 Erguns PAGE 38 22 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 2321180pdK (2.18) 23211808.1 FC (2.19) 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 Navier-Stokes equation. The first attempt at this was by PAGE 39 23 Brinkman, who included a term analogous to the Laplacian term of the Navier-Stokes equation, DDuuK2 fp (2.20) where is an effective viscosity. Some authors have set equal to f while others have defined it using the Einstein formula, given by 15.21f (2.21) However, the governing equations used herein are based on the premise that equals f Other authors have attempted to derive a set of governing equations based on the local volume averaging of the Navier-Stokes 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 VkkdVV1 (2.22) 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 VkkkkdVV1 (2.23) This leads to the continuity and momentum equations for Newtonian fluid in steady, incompressible flow in porous mediums, PAGE 40 24 Figure 2-4. A schematic of a representative elementary volume and the position vectors used. The fluid phase is shown as continuous, scanned from Kaviany [9]. 0fu (2.24) Sfffffffpuuu2 (2.25) Here, 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 fffxFxffxuuKCKuS (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, 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. PAGE 41 25 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 TkTceffffpu (2.27) where, fseffkkk 1 (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 fseffkkk 11 (2.29) 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 ffsssfsffffeffffffPTTahTkTCu (2.30) and for the solid phase, as ffsssfsfssseffTTahTk0 (2.31) Here the effective thermal conductivities for the solid and fluid phases respectively, can be written as PAGE 42 26 ffeffkk (2.32) sseffkk 1 (2.33) where a sf is the specific surface area of the interface and h sf 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 2-2. 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 2-2. Heat transfer coefficients with respective fluid to solid specific surface area. Model h sf a sf Notes 1 pfdk6.031RePr1.12 pd 16 2 59.033.0RePr064.1pfdk pd21346.20 Re > 350 3 1323110RePr2555.0spfpkdkd pd 16 Models taken from Alazmi [13]. PAGE 43 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 Navier-Stokes 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 pressure-correction method for solving the Navier-Stokes equations on a Cartesian or axi-symmetric 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 Stones Method is used for solving the resulting system of linear equations. 27 PAGE 44 28 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 3-1. Grid Line Internal Node Location(Center of CV) Scalar Boundary Node Scalar Internal Node u, v Boundary Node, Respectivelyu, v Internal Node, RespectivelyScalar CVu Momentum CVv Momentum CVControl Volumes Nodes Figure 3-1. Generic grid setup for staggered arrangement of variables; denote scalar variables, denote u velocities, and 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 3-1. 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 3-2 for an understanding of the notation used herein. The notation discussed is based on the PAGE 45 29 W E P S Nnsweseswnwnexixi-1xi+1xi-2yi-1yi-2yiyi+1(i,j)(i+1,j)(i-1,j)(i,j-1)(i,j+1) x y yx Figure 3-2. 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 v-momentum 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, x and y are the same for all CVs. 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: WESNnbQaaQuauatotalnb*nbnb*PputotalnbnbunbPup,,, Wherevv vvv** (3.1) PAGE 46 30 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 u-momentum equations only and can be easily extended for use in the v-momentum equations by simply changing the superscripts. For a first order upwind scheme, the coefficients of equation 3.1 are: PNnnunuNyySma0, (3.2a) SPssusuSyySma 0, (3.2b) PEeeueuExxSma0, (3.2c) WpwwuwuWxxSma0, (3.2d) (3.2e) nbuporousunbuPWESNnbQaa,,, Where Here the symbol, has been used to denote the greater of A and B, the S terms represent the surface areas of the control volumes, and is the volume of an individual CV. The source term due to the porous media are represented in equation 3.2e as In final discretized form the viscous and inertial sources become, BA, uporousQ uKCKQxFfxfuporous (3.3) Finally, the source term Q total of equation 3.1 contains the pressure and the convective flux resulting from the deferred correction approach, uconvectiveupressureutotalQQQ (3.4) Where the pressure and convective flux components are, respectively, 1mwweeupressureSpSpQ (3.5) PAGE 47 31 12,mUDSCDSuconvectiveUDSuconvectiveuconvectiveFFQ (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 (3.7) 12,mUDSeUDSCDSeeUDSeeceuumumF where the explicit terms are denoted by the superscript m-1, 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 equation 3.6, and the implicit terms in the coefficients 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. uconvectiveQ unba Finally the discretized equations can be solved. Many methods are available, here the strongly implicit procedure (SIP) or Stones Method is used. It is a form of incomplete lower-upper (LU) decomposition. More detail can be found in Ferziger [15]. PAGE 48 32 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 Semi-Implicit Method for Pressure-Linked Equations. This procedure was originally put forth by Patankar and Spalding [17]. As a reference, and introduction, the reader is referred to Figure 3-3, the details of which will be discussed subsequently. Start Solve Discretized Momentum Equations Solve Pressure Correction Equation Correct Pressure and Velcoties Convergence Solve Discretized Energy Equation Stop Yes Set calculated p, u, v as guessed fields No Guess Velocity, Pressure, and Temperature Fields Convergence Yes Set calculated T as guessed field No Figure 3-3. Outline of SIMPLE algorithm. PAGE 49 33 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, ppp* (3.8a) uuu* (3.8b) vvv* (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, nnssfeewwfnbnbpnbPpPSvSvSuSupApa****'' (3.9) The coefficients of equation 3.9 are, nvPfpNaSa2 (3.10a) svPfpSaSa2 (3.10d) euPfpEaSa2 (3.10c) wuPfpWaSa2 (3.10d) WESNnbaanbpnbpP,,, Where (3.10e) Additionally, the following variables can be written as ePuPeEuPeuPaaa1 (3.11) PAGE 50 34 where e is an interpolation factor defined as PEPeexxxx (3.12) 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 CVs as seen in Figure 3-1. 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 Stones Method to solve equation 3.9. After the pressure corrections have been determined the interface velocities are updated by '''PEuPeeppaSu (3.13) or, ''*PEuPeeeppaSuu (3.14) Then the pressures are corrected by '*PPPppp (3.15) PAGE 51 35 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 at an interface is equal to the value of its upstream or upwind neighbor. Mathematically the convective flux can be written as 0,0,eEePeeFFF (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, EPeeEPPEPeeeeFxxxxFF (3.17) PAGE 52 36 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 0,0,21eEEEEEeWWPeeFIFIF (3.18) where, WPWexxxxI1 (3.19a) EEEeEExxxxI 2 (3.19b) 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. PAGE 53 37 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. Total Radius of Pipe Radius of Center Jet: y Sub Domain 1 Inlet Velocity y Sub Domain 2:Rtotal RCenter Jet xstart:x location of the inlet (typically 0.0) pstart:x location of the beginning of the porous region pend:x location of the end of the porous region xend:x location of the exit Figure 3-4. 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 sub-region. First, the number of control volumes in the x-direction is needed. Once the number of CVs for each sub-region (entry, porous, and exit) has been entered, the y grid lines are generated at each x-location. 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 y-direction. A sample grid can be seen in the Figure 3-5. PAGE 54 38 0 0.01 0.02 0.03 0.04 0.05 0 0.005 0.01 0.015 0.02 0.025 Grid SetupRadius (m.)x-direction (m.) Figure 3-5. Sample grid displaying the important characteristics of the grid generation process. PAGE 55 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 well-known analytical solution for fully developed pipe flow. This test case was used to validate the axi-symmetric 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 well-known 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 4-1. 39 PAGE 56 40 UL H L Moving Lid Figure 4-1. 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 HUHULfLHRe The following Figures 4-2 and 4-3, offered as a measure to the 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 4-2 while the results for a 120 x 120 grid can be seen in Figure 4-3. On a course grid the 2UDS solution was best; however, as the grid was refined, both the CDS and 2UDS schemes provided accurate solutions. PAGE 57 41 -0.2 0 0.2 0.4 0.6 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u-velocity (m/sec)y (m) 1st Order UpwindCentral Difference2nd Order UpwindGhia et al. (1982) A) 0 0.2 0.4 0.6 0.8 1 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 x (m)v-velocity (m/sec) 1st Order UpwindCentral Difference2nd Order UpwindGhia et al. (1982) B) Figure 4-2. Benchmark Case 1, discretization scheme comparison on a 20 x 20 grid for Re = 1000. A) U-Velocity through y-geometric. B) V-Velocity through x-geometric center. PAGE 58 42 -0.2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u-velocity (m./sec.)y (m.) 1st Order UpwindCentral Difference2nd Order UpwindGhia et al. (1982) A) 0 0.2 0.4 0.6 0.8 1 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 x (m.)v-velocity (m./sec.) 1st Order UpwindCentral Difference2nd Order UpwindGhia et al. (1982) B) Figure 4-3. Benchmark Case 1, discretization scheme comparison on a 120 x 120 grid for Re = 1000. A) U-Velocity through y-geometric. B) V-Velocity through x-geometric center. PAGE 59 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 4-4. This problem has several well-known 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 hydro-dynamically and thermally. The resulting profiles in the fully developed region are compared to their analytical solutions. dd Hydrodynamic Entrance Region Fully Developed Region DUinlet = Um dTdT Thermal Entrance Region FD Region Tinlet Constant Wall Temperature, Twall > Tinlet Ro ru Figure 4-4. Benchmark Case 2, pressure driven flow in a pipe, problem setup. As shown above both the hydrodynamic and thermal boundary layers are developing simultaneously. PAGE 60 44 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: 212omRruru (4.1) where u m is the mean velocity of the flow, which for incompressible flows is defined by dAu A uAm1 (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 fDkhDNu (4.3) where D is the diameter of pipe, k f is the thermal conductivity of the fluid, and h is the heat transfer coefficient defined by mssTThq (4.4) In equation 4.4, is the heat flux at the wall or surface, sq oRrfsrTkq (4.5) while T s is the wall temperature, and T m 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 PAGE 61 45 oRommuTrdrRuT022 (4.6) Finally, combining equations 4.3 through 4.6 leads to the following definition of the Nusselt number, msRrDTTDrTNuo (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 x-direction, because the flow was given ample room to become fully developed. Consequently, only the number of nodes in the y-direction will be noted. The results with 22 nodes in the y-direction can be seen in Figure 4-5, 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.2E-5, or 0.002 % 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. PAGE 62 46 0 0.05 0.1 0.15 0.2 0.25 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 u-velocity (m/sec)r (m) -0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 % Error (Analytical Numerical)/Analyticalr (m) % Error AnalyticalNumerical A) 10-2 10-1 3 4 5 6 7 8 9 10 15 Nu1/Gz = (x/D)/(ReD*Pr)B) Figure 4-5. 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. PAGE 63 47 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 4-6. rxsRo Porous Media Open Channel Um Figure 4-6. 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, 2oRKDa (4.8) 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. PAGE 64 48 Fluid Region, 0 r s Ayu42 (4.9) Porous Region, s r 1, DayDaCKDaBIuoo2121 (4.10) where, 212112112121212112DaKsDaIsDaKDaIDaKsDasDaDaKBooo (4.11) 42 221121212112121121sDasDaKsDaKsDasDaKsDaKsDaIsDaIBAooo (4.12) 211212112112sDaKsDasDaKsDaIBC (4.13) Here, u is the dimensionless velocity defined as **2*dxdPRuufo (4.14) 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, PAGE 65 49 Fluid Region, 0 r s Ayu42 (4.15) Porous Region, s r 1, DayDaCKyDaBIuoo2121 (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 21DaBIo to give 21yDaBIo 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 x-axis. 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 4-7. As Figure 4-7 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 4-8. 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 4-9. Not only did the error decrease, but the time to run the simulation, diminished. PAGE 66 50 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u-velocity (Dimensionless)y (m.) -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 % Error (Annalytical Numerical)/Analytical % Error 2nd Order UpwindAnalytical Figure 4-7. Benchmark Case 3, comparison of velocity profiles to the analytical results of Poulikakos [20]. 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Radius (m.)x-direction (m.) Figure 4-8. 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. PAGE 67 51 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 4-10A. Respectively, the results of the code under investigation can be seen in Figure 4-10B. By inspection, the results of the current code are within reason. No noticeable error can be seen, deeming the solution correct. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u-velocity (Dimensionless)y (m.) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 % Error (Annalytical Numerical)/Analytical 2nd Order UpwindAnalytical % Error Figure 4-9. Benchmark Case 3, results of local grid refinement, with 80 nodes in the radial direction. PAGE 68 52 10-5 10-4 10-3 10-2 10-1 0 1 2 3 4 5 6 NuDarcy Number (Da) s=1.0s=0.5s=0.7s=0.9B.) Figure 4-10. 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. PAGE 69 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, Rigimesh. However, due to the large number of unknowns of the Rigimesh, such as material properties and pore geometry, a set of experiments must first be performed on a porous medium with a well-known geometry. This will provide the proper methods of determining the material constants, and should lay the foundation for future studies of the Rigimesh. 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 Forchheimers 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 Rigimesh. To begin, the numerical and experimental setup studied can be seen in Figure 5-1. 53 PAGE 70 54 R 1.036" / 2.63144 cmHole Pattern / REVa = 0.033" / 0.08382 cmb = 0.066" / 0.16764 cm b b 0.020" / 0.0508 cm R 0.010" / 0.0254 cm a a PorousPlate 0.0263144 m 0.0254 m 0.03175 m 0.05715 m Inlet Velocity P1P2 Figure 5-1. 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 5-1 yields a porosity of 1442.016764.00508.042Volume TotalVolume Void22 (5.1) 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 5-1. 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 PAGE 71 55 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 Darcys 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 C F for Case 4; however, it was based on the permeability from the capillary model. Table 5-1. Values of permeability and Forchheimers coefficient. Permeability K (m 2 ) Forchheimer Coefficient C F Case Model Value Model Value C1 322pd 1.16324E-9 m 2 epDd5.5155.0 0.520801 C2 2321180pd 5.87529E-12 m 2 23211808.1 2.44905 C3 Experimental 7.71735E-10 m 2 Experimental and K from C3 0.160154 C4 322pd 1.16324E-9 m 2 Experimental and K from C4 0.205578 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 x-nodes, 40 in the porous region, and 120 y-nodes, resulting in the mesh seen in Figure 5-2. PAGE 72 56 0 0.01 0.02 0.03 0.04 0.05 0 0.005 0.01 0.015 0.02 0.025 r (m)x-direction (m) Figure 5-2. Grid setup corresponding to the geometry seen in Figure 5-1. 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 5-1 with fluid properties and boundary conditions given in Table 5-2. The results can be seen in Figures 5-3 and 5-4. 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 5-3 for three of the four above models; additionally the experimental results are provided. Table 5-2. 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 () 1.1875 kg/m 3 V1 3.58 Dynamic Viscosity () 1.8048e-5 kg/m-s V2 7.766 Specific Heat (c p ) 1006.2 J/kg-K V3 10.543 Thermal Conductivity (k) 0.025913 W/m-K V4 13.140 V5 16.301 V6 18.122 V7 20.0795 V8 23.306 V9 25.754 PAGE 73 57 5 10 15 20 25 0 1 2 3 4 5 6 7 8 9x 104 Darcian Velocity (m/sec)Pressure Drop (Pa) Case C1Case C3Case C4Experimental (Carroll and Omar) Figure 5-3. Pressure drop across the orifice plate. Cases C1, C3, and C4 are labeled in Table 5-1. As can be seen in Figure 5-3, the results of Case C1 provided pressure drops much greater than the real values. Additionally, by inspection of the coefficients in Table 5-1, 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 Rigimesh the pore geometry is unknown, PAGE 74 58 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 6-9%. 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 Forchheimers coefficient (K = 5.74035E-9 m 2 and C F = 0.487469) were extracted and the corresponding simulations were run. The results compared to Case 3 can be seen in Figure 5-4. 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, did not intersect the y-axis at zero. This shift was ignored when K and C 89.364363.20516.482xxy F were extracted. Furthermore, when the curve fit was forced to intersect the y-axis at zero, negative coefficients resulted; therefore, was not investigated. Regardless, by varying the values of permeability and Forchheimers coefficient, the error was reduced. Moreover, Forchheimers constant my only be piecewise constant depending on flow regime. To PAGE 75 59 better understand this further investigation is necessary. As stated before, as more data is collected for investigations herein, the results should improve. 5 10 15 20 25 0.5 1 1.5 2 2.5 3x 104 Darcian Velocity (m/sec)Pressure Drop (Pa) Case 3K & CF from Curve FitExperimental (Carroll and Omar) Figure 5-4. Pressure drop across the orifice plate. Case C3 and results from K and C F based on a curve fit of the experimental data. The curve fit resulted in the following values, K = 5.74035E-9 m 2 and C F = 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 5-5. 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 Rigimesh, 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 PAGE 76 60 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. 0 5 10 15 20 25 30 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 Axial Velocity (m/sec)r (m) NumericalExperimental V1 V9 Figure 5-5. 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 non-equilibrium (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 PAGE 77 61 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 5-6. 0.0263144 m 0.0254 m 0.03175 m 0.05715 m Hot Air 500 K Cooled Wall, Constant Temp, 273.15 K kf = 0.025913 W/m-Kks = 110 W/m-Kkeff = 94.1417 W/m-K TDownstream P1P2 TUpstream Figure 5-6. Non-isothermal porous plate setup. Domain of interest and corresponding boundary conditions are as seen, with inlet velocities ranging from V1-V9. The effective conductivity was found from equation 2.28 in Chapter 2. It is also provided here for convenience, fseffkkk 1 (2.28) Using the solid and fluid phase conductivities provided in Figure 5-6, k eff becomes, 94.1417 W/m-K. The results of this model for the same inlet velocities previously discussed (V1-V9), can be seen in Figures 5-7 and 5-8. Additionally it should be noted that the boundary conditions seen in Figure 5-6 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. PAGE 78 62 Here the temperature profiles on the upstream and downstream faces of the injector are examined. As seen in Figure 5-7, the cooled wall boundary condition has a greater affect on temperature near the centerline at lower flow speeds, V1 as opposed to V9. 0 0.005 0.01 0.015 0.02 0.0263 273.15 300 320 340 360 380 400 420 440 460 480 Temperature (K)r (m) Upstream FaceDownstream Face V9 V1 Figure 5-7. Temperature profiles on the upstream and downstream faces of the porous slug. Two test cases are shown here, V1 corresponds to an inlet velocity of 3.58 m/sec and V9 to 25.754 m/sec. Similarly, Figure 5-8A 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 5-8B. 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. PAGE 79 63 0 0.005 0.01 0.015 0.02 0.025 280 300 320 340 360 380 400 420 440 460 480 500 Temperature (K)r (m) Increasing Velocity V1 V9 A) 0 0.005 0.01 0.015 0.02 0.0263 0 10 20 30 40 50 60 70 r (m)T = Tupstream Tdownstream (K) V1 V1 V9 V9 B) Increasing Velocity Figure 5-8. Temperature profiles, and T. A.) Temperature profiles on the downstream face for all inlet velocities tested, V1-V9. B.) T between the upstream and downstream faces for all inlet velocities tested, V1-V9. PAGE 80 64 Next, the effect of LTNE will be examined. The model previously discussed in Chapter 2 can be seen again in Table 5-3, where the value of a sf associated with that particular model has been determined. Table 5-3. Local thermal non-equilibrium models, w. respective coefficients h sf a sf Model D escri p tion D escri p tion Value ( m -1 ) H0 Thermal Equilibrium NA NA H1 pfdk6.031RePr1.12 pd 16 10107.9 H2 59.033.0RePr064.1pfdk pd21346.20 712.9 H3 1323110RePr2555.0spfpkdkd pd 16 10107.9 Because of the known geometry of the porous plate, the actual value of a sf as opposed to those based on the models, can be determined analytically as follows: 12222 2.132721222*222mdbdLdbLdVAappppssfsf (5.2) where b is defined in Figure 5-1 and L is the thickness of the porous region. Alternatively a sf could be written, 122bdVAapssfsf (5.3) Based on this value of a sf Model H2 seen in Table 5-3 is the best match. However, for the results discussed herein, a constant value of a sf as in equation 5-2, will be used. For simplicity, only the results from two flow speeds, V1 and V9, on the downstream PAGE 81 65 face of the injector are given here. These can be seen in Figure 5-9. 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. 0 0.005 0.01 0.015 0.02 0.025 280 300 320 340 360 380 400 420 440 460 480 500 Temperature (K)r (m) Fluid Phase, H1Solid Phase, H1Fluid Phase, H2Solid Phase, H2Fluid Phase, H3Solid Phase, H3LTE All From V9 All From V1 Figure 5-9. Temperature profiles for LTNE models on downstream face of the porous slug. Two test cases are shown here, V1 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 hydro-dynamically verified, a more accurate model of an SSME injector can be explored. To begin Figure 5-10 shows the assembly of injectors PAGE 82 66 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 5-11. Figure 5-10. 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 non-isothermal porous plate, simulations were run for a range of flow speeds, V1 V9. The results given here, Figures 5-12 through 5-14, 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. PAGE 83 67 R 1.036" / 2.63144 cm R 0.20 / 0.635 cm 0.00635 m PorousPlate 0.0263144 m 0.0254 m 0.03175 m 0.05715 m Hot Air 500 K Cooled Wall, Constant Temp, 273.15 K TDownstream TUpstream P1P2 kf = 0.025913 W/m-Kks = 110 W/m-Kkeff = 94.1417 W/m-KNote: Each individual hole seen here are not representative of each injector element. The coaxial injector is represented by the larger diameter center jet. The small holes act as pores in the porous face and are not to be confused with an injector element. Center Jet Figure 5-11. 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 5-1. As seen in Figure 5-12, a similar trend as for the porous plate is observed. Likewise, Figure 5-13A and 5-13B follow the similar patterns to those previously discussed. One additional comment should be made on the percent mass flow through the center jet, Figure 5-14. As the flow speed increased, the percent mass flow through the center jet 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 PAGE 84 68 0 0.0032 0.01 0.015 0.02 0.0263 273.15 300 320 340 360 380 400 420 440 460 480 500 Temperature (K)r (m) Upstream FaceDownstream Face Porous Region Jet V9 V1 Figure 5-12. Temperature profiles on the upstream and downstream injector faces. Two test cases are shown here, V1 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 5-15. 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. PAGE 85 69 0 0.0032 0.01 0.015 0.02 0.0263 273.15 300 320 340 360 380 400 420 440 460 480 500 Temperature (K)r (m) Porous Region Jet Increasing Velocity V9 V1 A) 0 0.0032 0.01 0.015 0.02 0.0263 0 10 20 30 40 50 60 r (m)T = Tupstream Tdownstream (K) Increasing Velocity V1 V9 V9 V1 Porous Region Jet B) Figure 5-13. Temperature Profiles and T. A.) Temperature profiles on the downstream face for all inlet velocities tested, V1-V9. B.) T between the upstream and downstream faces for all inlet velocities tested, V1-V9. PAGE 86 70 0 5 10 15 20 25 30 8.5 9 9.5 10 10.5 11 Darcian Velocity% Mass Flow through Center Jet Figure 5-14. Percent mass flow through the center jet for flow speeds, V1-V9. 0 0.005 0.01 0.015 0.02 0.025 300 350 400 450 500 Temperature (K) Upstream FaceDownstream Face r (m) V9 V1 Figure 5-15. Results from 2UDS simulations with the center jet. The unrealistic spike in temperature in the center jet region can be seen at flow speed V1. PAGE 87 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 71 PAGE 88 72 11**1memeeuPeeexpxpauu (6.1) 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 6-1. Refining the mesh locally near the porous interface resulted in decreased amplitudes of these oscillations; however, this increased the simulation time drastically. 0 0.01 0.02 Porous 0.04 0.05 0.06 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 x (m.)u velocity (m/sec.) Figure 6-1. U-velocity 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. PAGE 89 73 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. PAGE 90 CHAPTER 7 CONCLUSIONS The use of porous models to replace ad-hoc 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 so-called 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 Rigimesh 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 full-fledged 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 fine-tuned. 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. 74 PAGE 91 75 Of utmost importance to the development of LPREs 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 non-equilibrium 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. PAGE 92 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/rl10_high.jpg [4] Huzel, D. K., and Huang, D. H., Modern 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, 7 th ed., Wiley, New York, 2000. [6] Huzel, D. K., and Huang, D. H., Design of Liquid Propellant Rocket Engines, 2 nd ed., NASA SP-125, 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, Masters Thesis, University of Florida. [9] Kaviany, M., Principles of Heat Transfer in Porous Media, 2 nd Edition, Springer-Verlag, New York, 1995. [10] Dullien, F. A. L., Porous Media Fluid Transport and Pore 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 14-28. [12] Guin, J. A., Kessler, D. P., Greenkorn, R. A., The Permeability Tensor for Anisotropic Nonuniform Porous Media, Chemical Engineering Science, Vol. 26, 1971, pp 1475-1478. 76 PAGE 93 77 [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, 2 nd Edition, Springer-Verlag, New York, 1999 [15] Ferziger, J. H., and Peric, M., Computational Methods for Fluid Dynamics, 3 rd Edition, Springer-Verlag, 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 Three-Dimensional Parabolic Flows, International Journal of Heat and Mass Transfer, Vol. 15, 1972, pp 1787-1805. [18] Thakur, S., Wright, J., and Shyy, W., STREAM, A Computational Fluid Dynamics and Heat Transfer Navier-Stokes 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., High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method, Journal of Compuational Physics, Vol. 48, 1982, pp. 387-411. [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, 5 th Edition, Wiley, New York, 2002. PAGE 94 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. 78 |