Citation
Numerical Modeling of Transpiration Cooled Rocket Injectors

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 )

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