UFDC Home  myUFDC Home  Help 



Full Text  
TWODIMENSIONAL MODELING OF A CHEMICALLY REACTING, BOUNDARY LAYER FLOW IN A CATALYTIC REACTOR By PATRICK D. GRIFFIN 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 2006 Copyright 2006 by Patrick D. Griffin ACKNOWLEDGMENTS Siemens and the National Aeronautics and Space Administration supported this research. I thank Dr. David Mikolaitis, Dr. David Hahn, and Dr. Corin Segal for their assistance. I also thank my parents for their continued support and involvement in and out of my educational development. TABLE OF CONTENTS A C K N O W L E D G M E N T S ................................................................................................. iii LIST OF TABLES .............. ............................................... ........ vi L IST O F F IG U R E S .... ...... ................................................ .. .. ..... .............. vii A B S T R A C T .......................................... .................................................. v iii CHAPTER 1 IN T R O D U C T IO N ............................................................................. .............. ... 2 REDUCTION OF CONSERVATION EQUATIONS ...........................................7 A p p ly in g A ssu m option s ...................................................................... .....................8 C ontinuity E quation ............................................ ............................. ...........9 Species C continuity E quations.................................... ......................... .. ......... 10 M om entum E qu ation s ........................................... ........................................ 12 Energy Equation .................. .................................... ................. 15 O rder of M agnitude A analysis ......................................................... ............... 21 C ontinuity E quation ........... ................................................ ....... .... .... .... ... 22 Species C continuity Equations.................................... ......................... .. ......... 23 A xial M om entum Equation ........................................ ........................... 25 V vertical M om entum Equation ................................... ............................. ....... 27 Energy Equation .................. .................................... ................. 30 U nit A naly sis ........................................... ........................... 34 C ontinuity E quation ........... ................................................ ....... .... .... .... ... 34 Species C continuity E quations................................... .......................... ... ........ 35 M om entum Equation ............................................................. ............... 36 Energy Equation .................. .................................... ................. 37 Sum m ary of Governing Equations ........................................ ........................ 40 3 PROGRAM METHODOLOGY ........................................ .......................... 41 D iscretization ..................................... ................................ ......... 43 P aram eters and C conditions .............................................................. .....................45 Input and O utput Files ........................................ ................... ..... .... 47 Initial C conditions of a Stage ............................................... ............................ 49 S ta g e O n e ........................................................................................................ 5 0 B lasius Solution................................................ 51 Subsequent Stages ............................. .................... .. ........ .. .............53 Solving G governing Equations........................................... ............... .. ............. 53 M om entum Equation ............................................................. ..................57 C ontinuity E quation.......... ..... .................................... .... .... ...... .. 59 Species C continuity Equations.................................... ............................. ....... 61 E energy E qu action ................................................................ ............. .... 63 Species/Energy System of Equations ............................................... ...............67 4 T E S T IN G .......................................................................................7 1 C a se O n e ................................................................................................................ 7 2 R results of C ase O ne ................................................................ .. .......... ......73 C a se T w o ............................................................................. 7 5 Results of Case Tw o ......................... ..... ......... ........... ....... 76 C ase T h ree ............................................................................ 7 9 R results of C ase T three ....................................................................... ... .... ...80 C a se F o u r ...................................... ..................................................... 8 3 R esu lts of C ase F ou r ................................................................... ........ .......... 83 5 PROGRAM LIMITATIONS AND IMPROVEMENTS ............................................90 6 C O N C L U S IO N ................................................................................................. 9 1 L IST O F R E F E R E N C E S ...................................... .................................... ....................93 B IO G R A PH IC A L SK E T C H ...................................................................... ..................95 v LIST OF TABLES Table page 21 Equations m odeling the flow ......................................................... ............... 40 22 Units of the governing equations. ........................................ ....................... 40 41 Param eters and conditions of case one.......................................... ............... 73 42 Param eters and conditions of case two.......................................... ............... 76 43 Parameters and conditions of case three........................................ ............... 79 44 Parameters and conditions of case four.................. ............................................. 83 LIST OF FIGURES Figure page 21 D im ensionless variables. ................................................ ............................... 21 31 Flow chart for single stage modeling. ........................................... ............... 55 32 Flux components in the species/energy system......................................................68 33 Source components in the species/energy system.................................................69 34 Boundary conditions of the species/energy system............... ........ ............... 70 41 Axial velocity profiles of case one. ........................................ ....................... 74 42 Pressure plot of case one. ............................................... .............................. 75 43 Axial velocity profiles of case tw o................................................ ........ ....... 77 44 Pressure plot of case tw o. ............................................... ............................... 78 45 Reduction in methane concentrations of case two. .............................................78 46 Axial velocity profiles of case three.......................... ...................... .................. 81 47 Pressure plot of case three. .............................................. .............................. 82 48 Reduction in methane concentrations of case three. ..........................................82 49 Axial velocity profiles of case four. ........................................ ...................... 84 410 Pressure plot of case four. ............................................................. .....................85 411 Temperature profiles of case four. ........................................ ....................... 86 412 Methane concentrations of case four.....................................................................86 413 Hydrogen concentrations of case four. A) Mass fractions of atomic hydrogen. B) Mass fractions of diatomic hydrogen. ...................................... ...............87 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 TWODIMENSIONAL MODELING OF A CHEMICALLY REACTING, BOUNDARY LAYER FLOW IN A CATALYTIC REACTOR By Patrick D. Griffin August 2006 Chair: David Mikolaitis Major Department: Mechanical and Aerospace Engineering Problems associated with fossil fuels are increasing interest in alternative forms of energy production. Hydrogen is quickly becoming a popular option, but the efficient, affordable production of hydrogen is needed for it to become a viable source of energy. Catalytic reformation of hydrocarbons and alcohols appears to be a promising means of hydrogen production, but little is known about the surface chemistry. Research on heterogeneous catalyst and their reaction mechanisms is growing. A greater understanding of the surface chemistry could yield cheaper, more effective catalysts. The evolving chemistry of the surface catalyst is in need of a flexible software program to test new surface mechanisms. A program is developed to model chemically reacting flow through a catalytic reactor. The reactor is represented in twodimensional Cartesian coordinates with negligible body forces acting on the fluid. The flow is characterized as a steady, low Mach number, boundary layer flow of a Newtonian fluid. Basic principles of mass, species mass, momentum, and energy conservation are expressed mathematically and simplified. These principles are transformed into the equations controlling the behavior of the fluid and its motion through a process of applying assumptions, an order magnitude analysis, and a unit analysis. A code is written to numerically solve the resulting system of coupled governing equations. The methodology of constructing the program is decomposed into developing an orthogonal computational mesh, quantitatively defining the reactor and flow, locating chemical data and solutions, establishing initial boundary conditions, and solving the governing equations. The program is used to model four different flows: one with no chemistry, the second with only gas chemistry, and the third and fourth with gas and surface chemistry. Calculated solutions from each case are examined to confirm that the software produces reasonable results and is operational. The software is found to predict the point of ignition when the initial temperature is great enough to cause catalytic combustion. CHAPTER 1 INTRODUCTION The world is becoming increasingly aware of its dependence on fossil fuels. This fuel is meeting over eightyfive percent of our country's energy demands, which includes everything from electricity to transportation [1]. The power of fossil fuels lies in the atomic bonds of the hydrocarbons that make up these fuels. Energy is released by breaking these bonds in the process of combustion. The burning of fossil fuels also releases harmful byproducts that include: carbon monoxide, carbon dioxide, and nitrogen oxides. The carbon released into the atmosphere is originally trapped underneath the earth's surface, leading to an overall increase of carbon oxides in the atmosphere. Some believe these byproducts are leading to weather changes and health problems around the world. Energy extraction from fossil fuels is a relatively easy process and the fuel is readily available in deposits beneath the earth's surface. For these reasons, fossil fuels have become the main source of the world's energy production. The finite source is nonrenewable and will eventually run out. Decreasing supplies will lead to a rise in fuel cost and alternative forms of energy will become cheaper than fossil fuels. Economics involved with the decrease in fuel supplies will dictate that the world turn to alternative forms of energy. Whether for ecological or economical reasons, the world will need to find alternative forms of energy. Some look to the most abundant element in the universe, hydrogen. Hydrogen is a clean, renewable source of energy that can be used in combustion engines and fuel cells. Fuel cells are very efficient at producing electricity from hydrogen with the byproduct being water. A major obstacle in this alternative fuel is the affordable production of the energy carrier. Hydrogen rarely stands alone in its pure form. Most of the earth's hydrogen is bonded to oxygen and carbon, in the form of water, alcohols, and hydrocarbons. Water is an extremely stable molecule and takes a great deal of energy to extract hydrogen atoms. This energy must come from renewable sources if we wish to address the problems associated with fossil fuels. Hydrogen extraction from alcohols and hydrocarbons is much easier. However, fossil fuels are currently the main source of hydrocarbons. About ninetyfive percent of the hydrogen supply comes from the catalytic steam reforming of natural gas according to the US Department of Energy [2]. Natural gas is a relatively clean fossil fuel consisting mostly of methane. But natural gas is still a finite resource that will eventually run out. A very promising renewable source of hydrogen comes from ethanol. Ethanol is an alcohol that can be derived from biomass such as corn. Fuels produced from biomass release carbon into the atmosphere that is originally in the atmosphere leading to zero netproduction of carbon oxides [3]. There are many promising energy alternatives to fossil fuels. However, fossil fuels are so entrenched in our way of life, economically and politically, that few expect a quick transition away from fossil fuels. Most believe that hydrogen production will initially come from fossil fuels, with a gradual transition to renewable sources of hydrogen production. The world's attraction to the hydrogen economy is leading to an increased interest in heterogeneous catalyst for converting hydrocarbons and alcohols into the energy carrier. Surface catalysts are useful in increasing the reaction rates in combustors and reformers. Catalytic combustors burn the fuel over a catalyst. This bums fuel at a lower temperature, which decreases the amount of nitrogen oxides produced in the exhaust [4]. Catalytic reformers transform complex hydrocarbons and alcohols into hydrogen by stripping the fuel of their hydrogen atoms. In either case, the fuel molecule is adsorbed by the catalytic surface. The molecule forms a bond with the surface, usually through an oxygen or carbon atom. This weakens the adjacent bonds between the oxygen or carbon atom and the hydrogen atoms. The hydrogen atoms now begin to break off the molecule. The product molecule will detach from the surface once it is finished reacting with the catalyst. This leaves the surface free to adsorb a new reactant molecule. The catalyst provides reaction pathways with lower activation energies. In effect, the catalyst lowers the energy needed to break a molecule apart [5]. The efficient, affordable production of hydrogen is needed for this alternative fuel to become a viable source of energy. The efficiency of a metal to catalyze a given molecule is defined by how well the catalyst adsorbs the reactants and desorbs the products. Silver, for example, isn't a good catalyst because it doesn't form strong enough attachments with reactant molecules. Tungsten, on the other hand, isn't a good catalyst because it adsorbs too strongly. Metals like platinum and nickel make good catalysts because they adsorb strongly enough to hold and activate the reactants, but not so strongly that the products can't break away. [6] The efficiency of the catalyst has no affect on the metal's price. The price is dependant on the demand and rarity of the metal. As mentioned above, platinum and nickel are two common metals used in catalyst. Platinum cost approximately $1000 per ounce, where nickel cost around $0.4 per ounce [7]. With such a large disparity between efficiency and price, a greater understanding of the surface chemistry could lead to cheaper, more effective catalyst. Catalytic reformation of hydrocarbons and alcohols appears to be a promising means of hydrogen production, but little is know about the surface reactions. Surface catalysts are not fully understood because the chemistry around the surface is difficult to measure, especially in normal operating conditions. In the past, catalysts were treated as a black box. The black box representation of the catalyst usually consists of one global surface reaction or a small series of reduced mechanisms. Modifications to the black box can be made until the model accurately reproduces the experimental data. While this method is adequate for engineering applications, it does not accurately represent the chemistry involved [4]. Many studies have recently taken place in attempts to understand the surface reaction mechanisms of the heterogeneous catalyst [811]. The studies are mostly concerned with determining the reaction pathways and the stepbystep chemical degradation process of molecules. This is leading to new chemical reactions being added to the surface chemistry. The evolving chemistry of the surface catalyst is in need of a flexible software program to test the new mechanisms being added. A program adaptable to the changing surface chemistry is developed in this study. The program models a twodimensional, chemically reacting flow though a catalytic reactor. The fluid motion is characterized as a steady, low Mach number, boundary layer flow. The catalytic reactor consists of a heterogeneous catalyst covering the inside surface of a pipe or channel. The fluid motion is modeled as a flow through two flat plates with a pressure gradient. The two flat plates are modeled as catalytic surfaces and are identical. Basic principles of mass, species mass, momentum, and energy conservation are employed to generate the model. These principles are expressed mathematically and simplified for this specific problem. The process of reducing the principles into the governing equations consists of applying assumptions, an order magnitude analysis, and a unit analysis. A software code is written to numerically solve the resulting governing equations. The calculated solutions thermodynamically and kinetically define the fluid and its motion. The code is written in MATLAB, a programming language created by MathWorks that is used in many chemical flow simulations. MATLAB provides several builtin capabilities that make the software well suited for this problem. It's compatibility with Cantera being one such capability. Cantera is a free software package developed by Professor David Goodwin at the California Institute of Technology to solve problems concerning chemical reactions [12]. The program utilizes Cantera software to manage the chemistry. The methodology of the program's development includes the creation of an orthogonal computational mesh to resolve the equations. Then the establishment of parameters and conditions that quantify the reactor and fluid flow is performed. The location of input and output data is defined and initial boundary conditions are set. Finally the governing equations are solved and the flow in the catalytic reactor is modeled. The program is tested with four different cases: one with no chemistry, another with only gas chemistry, and two with gas and surface chemistry. Each case is modeled and the results are examined to confirm that the software produces reasonable results and is operational. In the future, calculated solutions can be compared to experimental measurements. New surface mechanisms can be tested with the program resulting from this study. Improved chemical kinetic data are updated in Cantera. The new chemistry is processed by Cantera and incorporated into the program. Any change in the chemistry being used to model the flow is done so inside the separate software of Cantera and not the main program. Because the data is stored separately, the program is able to remain flexible with the type of catalyst and fuel being used. This also allows the type of reaction pathways to change as our understanding of catalyst grows without altering the code. As a result the program is adaptable to the varying surface reaction pathways. Comparing the twodimensional model to experimental data provides a means of validating the accuracy of the new chemistry. With a better understanding, heterogeneous catalyst might be the key to the clean, renewable source of energy. CHAPTER 2 REDUCTION OF CONSERVATION EQUATIONS The chemically reacting flow through the reactor is modeled by numerically solving the governing equations. Four of the equations are derived on the principles of mass, momentum, and energy conservation [13]. The velocity field, pressure, and temperature in the reactor are determined with these four equations. Knowing two independent thermodynamic properties would adequately model the flow if it were not for the chemistry taking place. The catalytic surface is expected to induce chemical activity changing the fluid composition. A set of equations is needed to determine this changing chemical composition. The species continuity equations satisfy this need and are called upon to calculate the composition of the flow. One equation is needed to determine the mass fraction of a single atom or molecule. As a result, the number of equations inside this set is equal to the number of species used to model the flow, denoted as N. The mass, momentum, species and energy equations along with an equation of state are all that is needed to determine flow properties through out the reactor. These governing equations are coupled to one another in several different ways. All of them contain properties dependent on the flow variables. For example, the momentum, species, and energy equations contain transport properties such as viscosity, diffusion coefficients, and thermal conductivity. Most of these properties are dependent on the pressure, temperature, and composition of the flow. Properties dependent on the flow variables indirectly couple the equations to one another. The equations are also directly coupled to one another. Vertical gradients of the species mass fraction not only appear in the set of N species continuity equations, but also the final form of the energy equation. In addition to this, the velocity and velocity gradients can be found in all of the equations. This makes for a group of highly coupled equations that control the behavior of the flow. Equations of mass, momentum, species, and energy conservation are broken down and modified to reflect this specific model while Cantera processes the equation of state for an ideal gas. Conservation equations are transformed into the governing equations by applying assumptions characterizing the reactor and flow. Governing equations are individually examined in an order magnitude analysis after the assumptions are made. Dominant terms in a given equation are found by comparing their magnitude to the magnitude of other terms in the equation. Neglecting the weak terms and retaining the strong terms further reduce the equations. A unit analysis or unit check is applied to the resulting system of equations to ensure the validity of the equations. The process also establishes the units of each variable, property, and solution. Applying Assumptions Turns goes through a similar process of simplifying the conservation equations for a steady onedimensional flow [14]. Instead of one dimension, the computational space of the catalytic reactor is modeled in twodimensional orthogonal space. These two dimensions are the rectangular coordinates x and y, which represent the axial direction and vertical direction respectfully. Once the governing equations are reduced to their twodimensional form, they are simplified by making assumptions about the fluid and its motion. The chemically changing fluid is always considered a Newtonian fluid, which carries many assumptions with it. Most importantly of which is that the shear stress is linearly proportional to the rate of deformation. Another important assumption concerns the fluid's motion. The flow is modeled as a steadystate flow, meaning all fluid properties are independent of time. As a result, a partial derivative of any quantity with respect to time is zero. More assumptions are made in order to reduce the governing equations and are discussed during that process below. For the most part, the analysis mirrors that of a boundary layer flow. However, the catalytic surface creates density variations in the flow and compressibility must not be ignored. Continuity Equation The analysis begins with the reduction of the continuity or mass equation. The continuity equation is a mathematical representation of the conservation of mass that states that mass cannot be created or destroyed. In a Eulerian method of description, the conservation of mass is described as the time rate of change of mass in a control volume being equal to the net flux of mass through the control surface. Equation 21 is the vector form of the continuity equation. + V (p)= 0 (Equation 21) at The steady flow assumption leads to the partial derivative with respect to time being zero. The first term in Equation 21 is dropped as a result, leaving only the mass flux in vector notation. The catalytic reactor is being modeled in a twodimensional space. Therefore, the mass flux is written out into its twodimensional form with u and v representing the x and y component of the velocity, respectfully [15]. ( + )=p 0 (Equation 22) Dx Qy Further simplification is restricted due to the fact that density (defined asp) variations occur in the flow. Equation 22 represents the continuity equation for the flow field being modeled. This equation proves to be very important in the reduction of the other governing equations. However, it is not the form used by the program. The computer code uses the mass equation to determine the vertical velocity and some mathematical manipulation is needed before the equation reaches its final form below. Qv 8(pu) Sp v (pu) p (Equation 23) ay x y Species Continuity Equations Much like the continuity equation, the species continuity equation requires that the rate of gain of a single species mass in a control volume equals the net flux of the species mass in through the control surface. Dissimilarity in the two equations arises due to the chemistry. Instead of equaling zero, it equals the net chemical production of that species in the control volume. The continuity equation of species i is shown as Equation 24 and the set consists of one equation for each species. The time rate of change of the species mass is zero because the flow is steady. The species mass flux is expanded out into its twodimensional Cartesian coordinate form and the result is Equation 25. a(pY ) = m +V = i"'" "(Equation 24) S+ = "' (Equation 25) 9x 9y On the right hand side of the species continuity equation lies the net chemical production of species i in the control volume (ih"). This is determined using Cantera, which gives the chemical production in moles. Therefore, the species mass chemical production is replaced with the molar chemical production times the molecular weight of the species. The species mass entering the control volume, know as the species mass flux, transpire as a result of two modes, bulk flow and diffusion [14]. S+ '= AbMW (Equation 26) mhx = pYu + miD (Equation 27a) m" = Pv + m1yD (Equation 27b) The first term in Equations 27a and 27b is the mass flux due to the bulk flow. It is equal to the product of the density, species mass fraction, and fluid velocity component corresponding to the direction of the mass flux. The second term is the mass flux due to diffusion. The species mass flux can now be placed in the reduced species continuity equation and the two modes separated from each other. (pu)+ (pv)+ iD + (lf yD ) = AczMW (Equation 28) The chain rule is applied to the bulk flow terms and the process is shown in Equation 29. This leaves the continuity equation being multiplied by the mass fraction plus two mass fraction gradient terms being multiplied by the density and velocity. The continuity equation is equal to zero via Equation 22. After dropping the mass equation and replacing the two bulk flow terms with Equation 29, the species continuity equation reduces to Equation 210. a a 8a(pu) a(pv) (qai (pYu)+(pYv)= Y+ + pv (Equation 29) ax ay Ox ay ax Qy Coninuity Equanon pa +a (I,,f ) ( ,D' ) D= bMW (Equation 210) Dx ) + Qy C)C C Diffusion is a result of concentration gradients, temperature gradients, pressure gradients, and uneven body forces. Ordinary diffusion from concentration gradients is the only mode of diffusion considered in this model. The species mass diffusion is approximated using a mixtureaveraged diffusion coefficient [14]. The mass diffusion terms are replaced with Equations 211 a and 21 lb inside the governing equation. Equation 212 is the species continuity equation after all the assumptions are applied. Further simplification is possible with an order magnitude analysis. xDff = pD,  (Equation 211 a) rhyf,Df = PD, (Equation 21 lb) ay pu + ~+ D +a P = aMW, (Equation 212) ox dy dx^ ax } y oy Momentum Equations The momentum or NavierStokes equation is analogous to Newton's law of momentum conservation. The momentum equation states that the rate of change of linear momentum per unit volume equals the net momentum flux through that volume plus the sum of forces acting on the volume. This is mathematically written in vector form as Equation 213. Forces acting on the control volume are broken up into the surface forces and body forces. Surface forces are defined as the divergence of the stress tensor. The stress tensor is shown below in its rectangular, twodimension form as Equation 214. i + VV+V. ( = V + pFB (Equation 213) at 0 = j (Equation 214) 1^I Ua~ T\ I~ ^p) The stress tensor is made up of the shear stress acting on the surface plus the pressure acting normal to the volume. Again, the flow is considered steady with negligible body forces. Therefore, the time derivative and body force terms are dropped. The vector form of the momentum equation is separated into its xcomponent and y component equations. Equation 215a represents the twodimensional Cartesian momentum equations in the xdirection, while Equation 215b is in the ydirection. 8(puu,) 8(puv) )d 8a +uu + + (Equation 215a) ax ay ax ay 8(pvu) 8(pwv) o 9o + ) + (Equation 215b) The xmomentum equation is used as one of the governing equations in the computer program. Equation 215b, on the other hand, is not explicitly used in the computer code. It is used to gain some insight into the behavior of the pressure. Both equations simplify in a similar manner so both are reduced collectively. In this process, the xcomponent equation is given first followed by the ycomponent equation. The chain rule is applied to the momentum convection on the left hand side of the two Navier Stokes equations. F(pu) 9(pv) au au Qao ao S+ + pu+ pv + (Equation 216a) ax ay ] x Qy ax y Conhnurty Equahon 0(pu) 8(pv) av av BC7 ao, v +pu(+pvj + (Equation 216b) Conhnurty Equahon This results in the continuity equation being multiplied by the xvelocity in Equation 216a and by theyvelocity in Equation 216b. The continuity equation is equal to zero and the first term in these equations is dropped accordingly. The two momentum equations are reduced to Equations 217a and 217b. pu + pv + (Equation 217a) ax ay Ox ay 8v 8v Oa 8a pu+ pv + (Equation 217b) The flow through the reactor consists of a Newtonian fluid. The stress acting on a Newtonian fluid has no preferred direction, meaning the stress tensor matrix is symmetric and its components are defined below. By definition the shear stress of a Newtonian fluid is proportional to the rate of deformation. Shear stresses of a Newtonian fluid are written below as functions of the velocity gradients [15]. ,x = rT p (Equation 218a) Oy = Ty p (Equation 218b) O, = y = ry = ry (Equation 218c) = u 2 Pa u +av c 3 OaxOy av 2(au~av 1 , = 2u 2P + y3 Ox O IX= U+ YOX (Equation 219a) (Equation 219b) (Equation 219c) Replacing the stress components with their definitions above and pressure gradient, the two momentum equations become, separating the au au ap pu+ pv + ax ay 8x a r[2p u + i (Equation 220a) 8 8u 8F 2(au 8v 8 8u a v 2pj p + + p + 8x 8x 8x 3 8x 8y 8y By 8x av 8v ap pu+ pv  + 8x ay ay Sau v + v] a 2(u+v] (Equation 220b) S+ +[ 2au 0) a+ 8x L ay ax ) y 8y Ly 3 a8x ay Equations 220a and 220b represent the momentum equation in the xdirection and ydirection respectfully. Both are reduced to their final form via an order of magnitude comparison. Energy Equation The Energy Equation requires that the rate of change per unit volume is equal to the net energy flux into that volume due to convection, heat, and work [13]. 8 V2 V2 p pe+ +V pV e+V =Vq+ +pVF (Equation 221) at 2 )] 2 From left to right in Equation 221, the first term is the time rate of change of energy, which is zero because of the steady flow assumption. The second term represents the flux of energy due to convection and equals the heat transferred into the control volume plus the work done by the surface forces and the body forces. The work done by the surface forces is determined using the stress tensor of Equation 214. The body force is assumed to be negligible; therefore, the work done by the body force is neglected. Equation 221 is written out in its twodimensional Cartesian form with the assumed simplifications. 8 [ V2 8 V2 Spu e+ + pv e+ = LLv (Equation 222) Oq< Oq y a u + 9 + (x a Vu" + x 8) The stress components are replaced with the shear stress and pressure. The pressure is separated from the shear stress terms, leaving two pressure work terms at the end of the energy equation. pu e+ + v e+ = c 2 Sy 2 J [ 2 [) { y 2] (Equation 223) S + U(r +VT) ( + Ur, +VT') () (Vp) Placing the x and y pressure work terms on the far right of Equation 223 into the corresponding x and y energy convection terms on the left, Equation 223 becomes Equation 224. The internal energy and pressure is replaced by the enthalpy, defined in Equation 225 as the internal energy plus the product of the pressure and the specific volume. The energy transfer due to the shear stress work is replaced by a variable called (r work) to save space. Simplification of this term is possible via an order of magnitude analysis of the governing equation, but first the energy convection and heat flux terms are reduced. As of now the energy equation can be written out as Equation 226. pI e+p2+ + pI e++ = workor) (Equation 224) aL \ p 2} Oy \ p 2} a ' h=e+p (Equation 225) P Lpu[ 2 a\+ + pv 2 h+  +(z work) (Equation 226) 8x 2 y 2 O8x 8y  First the two energy convection terms on the left hand side of Equation 226 are simplified. These expressions consist of an enthalpy flux and kinetic energy flux, both due to bulk flow. The two convection terms are separated into enthalpy convection and kinetic energy convection. Performing the chain rule on the two kinetic energy convection terms produces four separate terms. Equation 228 illustrates the process. Two of these terms are combined to form the kinetic energy multiplying the continuity equation, which equals zero. After dropping this term, the kinetic energy flux is replaced with the last two terms of the equation above and the energy equation now takes the form of Equation 229. 8(puh) 8(pvh) a rv2+ 8 VV2 + + pu\ + pv = S ay ax Ly2 ay J (Equation 227) qx aqy + (+r work) x ay a [Pu V2 % a [ V2K ] pu + pv = 8x 2 J y 2 v2 r(pu) (pv) +L 2 +Pv V2 (Equation 228) V2 \8 pu) 8pv)\ pu 8 ^ ^ pv 8 ^ + + V + V 2 [ x 8y 2 ax 2 ay ConhnuityEquaton a(puh) + (pvh) pu 8 (v+ pv = (V2) 8x 8y 2 8x 2 8y x y 2 x 2(Equation 229) aqx aq + ( work) ax ay Moving over to the right hand side of the energy equation, the heat flux terms are now simplified. The heat flux is determined using Fourier's Law of Heat conduction plus the flux of enthalpy [14]. The enthalpy flux here is due only to diffusion. The flux of enthalpy from the bulk flow is already accounted for in the convection term. The vector equation of the heat flux is broken up into the two Cartesian coordinate components, x and y. Equations 27a and 27b are used to replace the species diffusion mass flux inside the sum of the heat flux. N q = kVT + ( D, ) (Equation 230) i=1 NT aTN N q = k + ( T ) k + ( h) (puYh (Equation 231a) x D1i x ,= ,=i T N N N y = k + (1"yDh) =k + (,yh ) (pvy~) (Equation 23 1b) o z=1 y 1=1 1=1 Taking the partial derivative of the two equations above produces Equation 232a and 232b. The partial derivative is not affected by the species sum and therefore can be moved inside the sum. Similarly, the mass flux due to the bulk flow is not affected by the sum and can be moved outside of the sum. The sum located inside the partial derivative of the last term contains the product of the species mass fraction and species enthalpy. Since the species enthalpy is given on a mass basis, the sum is equal to the specific enthalpy of the flow. The last term can be rewritten as the partial derivative of the density, velocity component, and enthalpy product. __, aT_ N(a >h aF aqx a L k dT]+," h 1 pu (Y ) (Equation 232a) Ox Ox x x LX) x. ' aq= a aT ] a k+ 8q_," yh, pv (Yh,) (Equation 232b) 9y 8y 9y O y ak + ,xh [+ puh] (Equation 233a) ax ax L x]x axL Ax 8qy [k aT] + a, hi:hj, [pvh] (Equation 233b) ay a ay 8yy ) ay The two heat flux terms of the energy equation are replaced with Equations 233a and 233b above. Once this is complete the energy equation takes the form of Equation 234 shown below. ( p 8, pu(V2 PV ,(V2) /x ./2 8x 2 By [k 2]+ [k~l Fm h [mFh,,] + (Equation 234) ax 8x ay ay ,8ax1x j jy j uh]+ vh] +(r work) The enthalpy convection cancels with the modified enthalpy diffusion of the heat flux. The chain rule is performed on the expressions inside the two sums. The process, shown in Equations 235a and 235b, leaves the species enthalpy times the partial derivative of the species mass flux plus the species mass flux times the enthalpy gradient. Equations 27a and 27b are used again to replace this species mass flux. This procedure takes the original partial derivatives and splits it into three terms each. a[,xa] =hO ax Ox+ =a h ,x + puY +hxiff" (Equation 235a) )x ax ax axx Lk, h] +'h + a =h h + Pv +, h Df,y h (Equation 235b) y Oy ,"dy  Y Replacing the two partial derivatives with their expanded expressions above, the process of simplifying the sums can begin. For the last two expressions in Equations 2 36a and 236b, the gradient of the species enthalpy is equal to the product of the species specific heat and the species temperature gradient. Every species comprising the fluid at a given point in the flow is assumed to have the same temperature, which means the temperature gradient can be moved outside of the sums. This process is done for all four of the terms containing enthalpy gradients. Oa1 h h 1 O aha OX( hj= puY a AxDif (Equation 236a] ].,, h m+ pyYh l + (Equation 236b) \ ax 9 ax OpuY = pu (Yc =p (Equation 237a) pv Y = pv Y(c =pV aTc (Equation 237b) Sx =p = D c (Equation 238a) f OX j (1 Oax NP OX X 9a (rnDif yDf (" )p p(D, c (Equation 238b) In Equations 237a and 237b, the product of the species mass fractions and multiplied by the mass flux, which equals the flow density times the proper velocity component. The diffusion mass flux inside the sum of the last two expressions is approximated using the mixtureaveraged diffusion equation, Equation 211. After replacing the four terms with four equations above, Equations 236a and 236b are added together and rearranged before being placed into the energy equation. ,a aT D h, h = Pc + pVC + S1 c (Equation 239) h am + a j pJ D, K c +m L A 1 1=1ax Sy a Lx ( x o y y Species _Continuty The two dimensional gradient of the species mass flux is replaced with the species chemical production via the species continuity equation, Equation 26. The energy chemical production via the species continuity equation, Equation 26. The energy I equation is reduced to the Equation 240. Further simplification is performed with an order magnitude analysis in the next section. pu V2 ) pv a(v2) T a T_ 2 8x 2 By 8xL 8xx By By puc T pvc (hMW)+ (Equation 240) ax 8y 'I S8x 8x y y Order of Magnitude Analysis An order of magnitude comparison between terms in a given equation determines which terms must be reserved and which terms can be neglected. Governing equations that are modified for this specific model are simplified further by eliminating the insignificant terms. It is necessary to nondimensionalize the equation prior to comparing terms. Variables are nondimensionalized with the uniform properties of the flow entering the reactor. Most of the properties are chosen such that the resulting magnitudes are on the order of one. The dimensionless variables and their magnitude are shown below in Figure 21. p* 0,o (1) T*= = 0(1) 0(1) k* k 0(1) u*= 0(1) p* P =0(1) y*=Y= (1) D,= = D = 0(1) U pmU H D v* H c v*= =0(?) P = 0(1) L] H =0(S) c = =0(1) U Lo L CP Figure 21. Dimensionless variables. All but two of the dimensionless parameters have a magnitude on the order of one. The unknown magnitude of the vertical velocity is found with the continuity equation. The characteristic distance in the axial direction, L, is much greater than the characteristic distance in the vertical direction, H. A dimensionless parameter with a very small magnitude, denoted by 0(6), is produced when the characteristic height is divided by the characteristic length. Continuity Equation Equation 22 is the twodimensional continuity equation that is reduced based on the steady flow assumption. Flow properties are replaced with their appropriate dimensionless variables. After some algebraic rearranging, the mass equation is rewritten in its dimensionless form. a( + =( 0 (Equation 22) 8(p*u*) L 8(p*v*) + = 0 (Equation 241) ax* H Qy* All of the known dimensionless variables have an order magnitude of one. It has already been noted that the characteristic length is much larger than the characteristic height. This produces a relatively small value that divides the vertical mass flux term. In order to balance the mass equation, the dimensionless yvelocity must have the same order magnitude as the division of the height by the length. O(1) 1 v* 0(+ = 0 (Equation 242) 0(1) 0(3)0(1) v* = 0(3) (Equation 243) While the continuity equation remains unchanged, the comparison of terms reveals that the vertical velocity of the flow is small compared to axial velocity. This is a common result in boundary layer flow analysis. Growth of the boundary layer is dictated by the viscosity, or momentum transfer, and does not affect the entire flow until farther downstream. The flow does not consist entirely of a boundary layer flow. However, a vertical velocity does not exist at the entrance of the reactor, on the surface, or at the centerline of the pipe or channel. The vertical velocity remains much smaller than the axial velocity through out the reactor because of these boundary conditions. Species Continuity Equations The species continuity equation is reduced based on the assumptions of a steady, twodimensional flow, with ordinary diffusion being the only mode of diffusion. This equation is shown below. pu +pv ~+ pD y PDm = MW (Equation 212) The dimensionless form of the species mass equation is obtained by replacing the flow variables with their proper dimensionless counterpart. The species mass fraction is exempt from this part of the process because it is already a dimensionless quantity that varies between zero and one. The unknown magnitude of the mass fraction does not pose a problem since it is found in every term on the left hand side of the equation. As a result it affects the magnitude of each term equally. After some algebraic manipulation, the left hand side is rewritten in its dimensionless form as Equation 244. The right side of the species equation, the species chemical production, is not compared to the rest of the equation. Neglecting this term would result in the modeling of a nonreacting flow. Therefore, the convection and diffusion terms are the only terms considered. S a8Y L( _Y p*u* +L p*v* \+ 8x*) v (Equation 244) p*D 2 a *Kp* UL Lx dx*) H2 ay* D" Dy *J Note that the massaveraged diffusion coefficients are nondimensionalized by an arbitrary value. This value is chosen such that the quantity of the dimensionless property is roughly one. This ensures that the dimensionless diffusion coefficients have an order magnitude of one, but the size of the value relative to the product of the incoming velocity and characteristic length is unknown. The species diffusion terms cannot be compared to the species bulk flow terms as a result. However, the comparison between the diffusion terms inside the brackets is still possible. O(1)+ 0 )+D (1)+f (Equation 245) 0(8) UL 0(S) Both of the bulk flow terms have a magnitude on the order of one. The first term inside the brackets, corresponding to diffusion in the axial direction, also has a magnitude of one. The second term corresponds to the diffusion in the vertical direction and has a magnitude much greater than one. The order magnitude comparison of the species continuity equation shows that the xcomponent of the species diffusion is much smaller than the vertical diffusion and can be neglected. Information about the size of the characteristic diffusion coefficient relative to the product of the characteristic velocity and length is needed to determine the parameter multiplying the vertical diffusion inside the brackets. The parameter must be very small, on the order of 0(6)2, in order for the vertical diffusion to be of the same magnitude as the two bulk flow terms. This means that the species mass transfer from the bulk flow is much greater than the species mass transfer due to diffusion. Though this is most likely the case for the flow being modeled, further restricting the flow to this assumption does not simplify the equation. The species continuity equation for the flow through the reactor is now reduced down to Equation 2 46 after dropping the axial diffusion term. pa + a + pD  = cbMW, (Equation 246) tx ay ay ay) Axial Momentum Equation Analysis of the axial momentum equation is performed in the same manner as the other equations. Flow properties are nondimensionalized by their characteristic variables. The comparison begins with the momentum equation governing a steady, two dimensional flow of a Newtonian fluid in the xdirection. Characteristic variables are rearranged, and the dimensionless form of the xmomentum equation is shown as Equation 247. The viscous term inside the brackets is compared separately from the momentum flux and pressure gradient terms due to the length of the expression. The comparison of the dimensionless momentum flux and pressure gradient terms is now possible. Excluding the vertical velocity, all of the dimensionless variables have a magnitude on the order of one. The division of the characteristic length by the characteristic height produces a relatively large value. This value is multiplied by the dimensionless vertical velocity, which is a small quantity. The overall effect produces a momentum flux and pressure gradient terms that all have the same order magnitude of one. As a result, none of these terms is less important than the other and none of the three can be ignored. pu+ pv a + 8x 8y &x \ 8 S u 8 2 u 8v 8 u v ra 2 a avx (Equation 220a) 8x x L x x 38x 8y 8y y y 8x au* L au* p*u* + p*v* Ox* H ay* (1)+ ( (1) =0(1) 0(3) ap* L = + 2rT term} + P.U x LPU + {2\r term ) pm2U ~mx (Equation 247) (Equation 248) Viscous terms inside the brackets are transformed into dimensionless variables and compared to each other. Equation 249 represents the dimensionless form of the viscous term. The characteristic properties are reorganized and the viscous term is now multiplied by the inverse of the Reynolds number. The Reynolds number is a common dimensionless parameter used to compare inertial forces to viscous forces. The Reynolds number in Equation 250 is based on the length of the reactor and therefore is a comparison of these two forces in the axial direction. L ( L r term}\= pmU2  L a u u 2 u av + u +v 2 2 + + + 9 L ]_u* a 9 (2Su* +L av* + 2*  + + P a* a* 3* H y*) P, UL a *( L2 u* L av*) *H *\ a+* 9y* H2 9y* H 9x*) (Equation 249) (Equation 250) Magnitudes of each term that comprise the viscous momentum transfer expression can now be compared to one another. Every expression inside the brackets is of the order of one, except for a single term. This term is underlined twice in Equation 251 and has a magnitude much greater than one. The result is a significant reduction of the viscous term. With the exception of the highlighted term, every expression is neglected and the complex viscous expression is simplified to just one term. Information about the magnitude of the Reynolds number is needed to compare the viscous term with the rest of the momentum equation. The inverse of the Reynolds number must have a magnitude of 0(6)2 for the remaining viscous term to be of a similar size as the momentum flux and pressure gradient. A large Reynolds number assumption forces the viscous term to balance with the other terms in the equation. It also forces the inertial forces of the flow to be more significant than the viscous forces. This is a reasonable assumption because the Reynolds number is based on the axial direction, where inertial forces are expected to be greater than the viscous forces [16]. 1 O 0(3) 1 O (3) 1 (1) O(1)+ +0 (Equation 251) Re O (3)1+ () 0 + () The order of magnitude comparison of the NavierStokes equation in the axial direction produces a couple of useful results. The complex momentum transfer due to viscosity is simplified to a single term. This reduces the xmomentum equation to its final form used in the computer code. In addition, the Reynolds number must be large for the viscous momentum transfer to be of a size comparable to the rest of the momentum equation. The large Reynolds number result is used later in the magnitude comparison of the ymomentum equation. pu + p +/ (Equation 252) 1 Re = (Equation 253) 0(3)2 Vertical Momentum Equation The momentum equation in the vertical direction is not used directly in the computer program. It is used to gain some insight into the behavior of the pressure through out the flow. The process of comparing terms in theymomentum equation is the same as the process for the xmomentum equation. The NavierStokes equation for a steady, twodimensional flow of a Newtonian fluid is nondimensionalized by the characteristic scales. The resulting dimensionless equation takes the form of Equation 2 54 after some algebraic manipulation. Again, the viscous term is broken down separately because of the length of the expression. 8v 8v 8p PU + PV =+ O8x y Dy ( [ ] _ [ 22u + v ) (Equation 220b) a u 8v} a 8v \ \l )+ + 2/I/  _ + O8x \y 8x y 8y 8y 3\ 8x 8y 8v* L 8v* L 8p* L p*u* +p*v* + (Equation 254) 8x* H Dy* H Dy* p.U2 term Comparison of the momentum flux and pressure terms is postponed until the magnitude of the viscous momentum transfer is known. After the flow properties are converted to their dimensionless form, the viscous term becomes Equation 255. Similar to the viscous term in the axial NavierStokes equation, the parameter multiplying the viscous term is the Reynolds number. Analysis of the axial momentum equation determined that the Reynolds number is on the order of 1/0(6)2. Knowing the magnitude of the Reynolds number allows the viscous term to be compared to the momentum flux and pressure gradient terms. But first an order magnitude comparison must be performed on all the terms inside the brackets. 08 ~ J L 0u 8v* ] L2 8 a 2,u * + +H2 2/* SL ; "a H + :J]kA; "t L *" I (Equation 255) PUL L 8 [ 2Su* L v*) IH O  +H H 8y 38x H Dy* 1 + 0(5)]+ 0 )[O(1)+o0)1} (Equation 256) Re O (5) O ([5) 0(5) 0(5) As a reminder, all of the dimensionless properties, excluding the vertical velocity, have a magnitude of one. The dimensionless vertical velocity and the division of the characteristic height by the length have a magnitude much less than one. The result is the magnitude of all but a single term inside the brackets of Equation 255 reducing to an order of 1/0(6). The single remaining term, underlined twice in Equation 256, has a magnitude of 0(6). This is much smaller than the other terms and can be neglected. As a result, the dimensionless viscous term in the vertical momentum equation has an order magnitude of 1/0(6). The process ofnondimensionalizing the equation is complete. An analysis of the vertical momentum equation is now possible with the knowledge of the viscous term's magnitude. ,v* L v* L 9p* 1 p*u* + + v* O + r term}) (Equation 257) Ox* H 9y* H 9y* Re 0( ()+ 0(3)= +0()20 )} (Equation 258) S0(,) 0(,) 0(,) The magnitude of every term, with the exception of one, is a very small quantity, 0(6). The exception is underscored in Equation 258 and is relatively large compared to the rest of the equation. The highlighted value corresponds to the pressure gradient in the ydirection. The pressure gradient is considerably larger than the other terms and dominates this equation. Neglecting all of the irrelevant terms, the dimensionless pressure gradient is the only term remaining. This results in the pressure gradient in the ydirection being essentially zero. O 0 (Equation 259) 8y The order of magnitude comparison of the vertical momentum equation reveals that pressure through out the flow is a weak function of the vertical position. Pressure can be treated as strictly a function of the axial direction, x. Energy Equation The magnitude comparison of the energy equation begins with Equation 240. The equation is simplified based on the assumptions of a steady, twodimensional flow with negligible work done on the fluid by the body forces. Magnitudes of each similar expression are compared to one another separately. First, the two kinetic energy convection terms are compared to each other. Then the two heat conduction terms are compared and so on. The work done by the shear stress is evaluated last. pu (V) V pva(V2) 8 FkT] a FkT 2 8x 2 y 8xL 8xx y By pu8T pvcp (h, cMW) + (Equation 240) PDrcp + O] + workr) pDc x 8Ox Ox ay ay Analysis of the energy convection on the left hand side of Equation 240 is a qualitative process. The couple of convection terms are not nondimensionalized, but a simplification is possible with the understanding of the behavior of the velocity field. The energy being transferred consists of kinetic energy, which is proportional to the magnitude of the flow velocity squared. Balancing the continuity equation proved that the vertical component of the velocity is much smaller than the axial component. Squaring both components only make this difference more pronounce. V2 =2 V2 2 2 = U2 [ *2 *2] = U2 [O(1)2 +0()2] V2 = 2 (Equation 260) The contribution of the vertical component to the magnitude of the flow velocity is neglected. The kinetic energy is calculated using only the axial component of the flow velocity and the kinetic energy convection is rewritten as Equation 261. The two partial derivatives of the xvelocity squared are performed to obtain the final form of the kinetic energy transfer. puS C(2) pvS(U2) 2u 8u pUa(2+ = (u2) p + puv (Equation 261) 2 Ox 2 Oy Ox Sy Moving over to the right hand side of Equation 240, the two heat conduction terms are now compared in a much more quantitative procedure. Variables are replaced with their dimensionless counterparts in Equation 262 so a dominant term can be found. These terms cannot be compared to the rest of the energy equation because the magnitude of the parameter outside of the brackets of Equation 263 is unknown. Comparing the magnitudes of the two expressions reveals the dominant term. It is underlined twice and corresponds to the heat conduction in the vertical direction. Being much smaller than vertical conduction, the heat conduction in the axial direction is neglected in the final form of the energy equation. a k a k =O kk*T* 2 [k (Equation 262) k. T{ ( ( T1 2 ( 2T* L2 [ H2 y* kT O (1)+ 1 (Equation 263) L 0(,)2 The energy transfer due to the diffusion of enthalpy is the next term to be reduced. Every variable and property is replaced with its dimensionless representation in Equation 264. The order magnitude analysis reveals that the second term, underscored in Equation 265, is much larger than the other. The lesser of the two terms is the enthalpy diffusion in the axial direction and is neglected in the energy equation. "F(Y ~ T aY8Ty Z\PDLmPC a x OT]2 pD=1 c 8 8x x y oy) (Equation 264) Dp, c,T, c. [ Y OT +* L2 O OT* * L2 H 2 yx* y *H s" a' __H_ DpocPT No D^^cT 0(1) o(1)+ 1 v ~ ( ) (Equation 265) The last expression to be analyzed in the energy equation is the energy transfer due to work done by the shear stress. It is greatly simplified by comparing the magnitude of each term that comprises the work. Shear stress is defined as Equations 219ac for a Newtonian fluid. The energy transfer expression becomes Equation 267 after the shear stress definitions are substituted. The flow properties are replaced with their proper dimensionless variables and characteristic scales. ( _work) = (urxx + VT ) + (ux + VT ) ax *y 8F vu 2 ( QPu ov a(u + iv1+ Sx2up aup u+.+vI+ +x ax ax 3 ax V y ay ax Q v 2 u 8v u v 2vv + \ + up + ay ay 3 ax ay ay ax (Equation 266) (Equation 267) S* u* 2 u* u* L av*> 2u*u*u** + + U2P. a dx* 3 dx* H dy* L2 ax* *( *L du 8v* v* *\ + H Oy* 8x* ___ v* + 9 2 OH du 8v* 2v* *v** + + U2P, L2 8 y2* ~ 3 L 8x* y* L2 H2 Iy (yu H 8v*) 8y L dx* (Equation 268) The dimensionless form of the viscous work can now be used to determine which terms dominate the energy transfer expression. Most of the terms have an order magnitude of one. There are two terms that do not have this magnitude. One is very small with a magnitude on the order of 0(6)2, while the other has a large magnitude and is underlined twice. The highlighted term is substantially greater than the other terms inside the brackets and can be considered the only dominant term. Neglecting all the other weak terms, the energy transfer due to the shear stress work is reduced to Equation 270. ( o(1) o() 0o() o() o()2 u 0(2)2 0()2 0(8)21 0 ( 8)2 L0(8)2 0(5)2 0(g)2 0(g)2 0(g)2 (Equation 269) y upy (Equation 270) The order magnitude analysis of the energy equation has greatly simplified the governing equation. The kinetic energy convection, heat conduction, enthalpy diffusion, and shear stress work are analyzed individually. These four modes of energy transfer are not compared to one another because no assumptions are made about which mode is more important. Analysis of the velocity field has revealed that the axial component can be used to determine the magnitude of the velocity at any point of the flow. The heat conduction in the xdirection and many of the terms that comprise the shear stress work are neglected due to the analysis. As a result of all this simplification, the energy equation becomes Equation 271. 2 u au 8 T aT aT pu +puv k p puc P ox By By By ox By L (Equation 271) N N\ ]a YaT7 8 au S(hc)A'MW)+ pDc  + I upu zK1 1L ayay y) ay L ay j Unit Analysis Units are substituted into the governing equations to ensure each expression in an equation balances with the other expressions. Replacing variables and properties with their units reveals several important aspects of these quantities. The inspection validates that the equations were reduced without misplacing any variables or properties. It also locates properties that require a unit conversion and determines units of the calculated solutions. Note that Cantera calculates the properties with the International System (SI) of measurement [17]. In order to keep unit conversions to a minimum, the variables also use this system of measurement. Continuity Equation The main program calculates the vertical velocity component with the mass equation. Solving the mass equation for the yvelocity produces Equation 23. Each variable and property is replaced with its units. The density is given in kilograms per meter cubed by Cantera. Therefore, the velocity and differential distances are measured in meters per second and meters, respectfully. From Equation 272, it is clear that the mass equation balances with equivalent units on both side of the equation. av 8(pu) Sp pv (u) (Equation 23) y x y kg X 1 kg 6 1 6 kg 1 kg kg A^ __, (Equation 272) mn s mn s s m3 ) m s nm3 .s Unit analysis of the mass equation reveals that no unit conversion of the density is necessary and the two differential step sizes should be given in similar units. Units of the calculated vertical velocity depend on the units of the axial velocity, which is defined by the initial condition. Although units cancel each other out in the mass equation, the other governing equations prove that SI units should be used for the variables. Species Continuity Equations The reduced species continuity equation is algebraically reorganized in a form the program can solve. This form is discussed more in section Solving Governing Equations. Equation 246 is transformed into Equation 273. International System of measurement is used for the velocity and differential distances, and their units are meters per second and meters, respectfully. The density is still kilograms per meter cubed and the mass fraction is dimensionless. Cantera gives a mixtureaveraged diffusion coefficient in meters squared per second. The unit of the net production rate is kilomoles per second meter cubed, and the molecular weight is given in kilograms per kilomole. All of these units are placed into Equation 273 to complete the analysis. pu = pD,  pv + aAMW (Equation 273) ax 9y Qy j y kg 1 1 kg W/ 1 kg X 1 I6 kg S3 S 3 S (Equation 274) kg kg m3 S m3 S The unit check of the species continuity equation shows that converting the molecular production rate into a mass production rate is indeed necessary. No other unit conversion is required if SI units are used for the velocity components and differential step sizes. The process proves the equation is reduced correctly from a unit analysis point of view, and its solution is dimensionless. Momentum Equation The reduced momentum equation is reorganized into a form similar to Equation 2 73. Equation 275 is the form of the momentum equation solved by the program. Units of the velocity components, differential distances, and density remain unchanged. The pressure and dynamic viscosity is also given in SI units. The SI unit of measurement for pressure is the Pascal, which equals a kilogram per meter per second squared. Cantera reports the dynamic viscosity in units of Pascalsecond. A Pascalsecond is equivalent to a kilogram per metersecond. Variables and properties are replaced with these units and result in balanced Equation 276. pu / du pv (Equation 275) Ox oay ) Ox ay kg j 1 1 kg X 1 kg 1 kg / 1 na s s i mnm s / n s n na s s i /  m/2 S (Equation 276) kg kg 2 2 2 2 m *s ms *ss The analysis finds that none of the properties determined by Cantera necessitate a unit conversion. Units of the simplified momentum equation balance accurately and its solution, the axial velocity, is calculated in meters per second. Energy Equation Like the last two governing equations, the reduced energy equation is organized into the form of Equation 277. 8T 8 k T + u 8T pucP k p+ual vc + SY L y (Equation 277) Z OYT A' Ou Ou \pDc [ (hcMW)pu puv D' "Ic y 8y x 8y Units of the flow properties already discussed in the previous governing equations remain the same. Several new properties are encountered in the energy equation. The temperature, thermal conductivity, species enthalpy, and specific heat of the fluid and species i are used exclusively by this equation. The temperature is measured in degrees Kelvin, and the thermal conductivity is given in units of watts per meterKelvin. A watt per meterKelvin is equivalent to a Joule per meterKelvinsecond. Specific heat of the fluid can be determined on a mass basis in Cantera. The unit of the fluid's specific heat is Joules per kilogramKelvin. Each expression is analyzed individually with the units established. Moving left to right in Equation 277, the axial energy convection is analyzed first. Substituting units into the axial energy convection shows that the term has units of Joules per cubic meter second. Units of the other expressions must reduce to this unit to balance the energy equation. aT /J / J pucp  (Equation 278) SOx m s S ./ M 3 S The next expression evaluated is the energy conduction along with the viscous term. Units of the viscous term are equivalent to the energy conduction at a Joule per cubic meter second and Equation 279 balances with Equation 278. 8 F 8T 8u~ 1 J __ 1_\ J___ [L O +up O = J + = (Equation 279) oy Oy oy m m/ *sm s / / m S Units of the second expression reduce to the same units of the axial energy convection. The vertical energy convection also reduces to these units and is shown below in Equation 280. PVCT J J ay in3 s 3S 9./,M/ 4 3. (Equation 280) Cantera reports the species specific heat in a column vector that has been nondimensionalized by the universal gas constant. Multiplying the vector by the universal gas constant produces specific heats with the units of Joules per kilomole Kelvin. The sum of the energy diffusion expression includes the species specific heat. The term highlighted is added to convert the species specific heat from a molar basis to a mass basis. It is the inverse of the species molecular weight and must be added to the enthalpy diffusion expression for the units to conform to the rest of the energy equation. Equation 282 illustrates the modification. Z pDDC 3 I = (Equation 281) ay ay In S 3 S N NT PLDc [Y OT]pD , LLmOT (Equation 282) Much like the species specific heat, Cantera reports the enthalpy of each species in a dimensionless column vector. The vector is nondimensionalized by the universal gas constant and the temperature of the fluid. After multiplying the vector by these two properties, the resulting enthalpy has the units of Joules per kilomole. Unit analysis of the enthalpy production expression is done in Equation 283. Because the species enthalpy is reported on a molar basis, the chemical production rate does not need to be converted to a mass basis. The molecular weight term is underscored and is dropped such that the unit of this expression is consistent with the other terms of the energy equation. NJ 61 kg J kg 0(h)MW) = 3 M3 ,(=1 ht m3 s kmol m3 *s kmol N N I(haAMW4T ) :> Y(h0 ) (Equation 283) (Equation 284) The last expression in the unit check is the kinetic energy convection. Equation 2 85 shows that the expression needs no modification. pu2 puv   (Equation 285) )x C y m3 / f m3.  Analysis of the energy equation reveals that energy diffusion and enthalpy production expressions required modification. Molecular weights are added to the energy diffusion term and removed from the enthalpy production term. After these modifications, the energy equation becomes Equation 286 where the units of each expression are equal and the equation balances. Temperature being the dependent variable of the energy equation is calculated in degrees Kelvin. 8T 8 k T + u~1 8T pucp k +u I pvcp + 8x 8y 8y 8y 8 By OX YL OY Y (Equation 286) p D. c, LY,' T O a8u au P\D \ (hco pupuv I Iy 8y y 8x By Summary of Governing Equations In conclusion the governing equations are applied to the twodimensional modeling of the catalytic reactor. Each equation is reduced based on assumptions describing the fluid and its motion. Through an order magnitude analysis these equations are simplified further. The units of each term are verified and the equations are balanced. The resulting equations solved by the program are summarized in Table 21. Table 21. Equations modeling the flow. Equation Principle Equation number number Mass av a(pu) Sp p v 23 Conservation 2y ax y Species Mass a ai" SpeciespuMass pD MW i= 1,2,3,...N 273 Conservation Ox Sy[ 9 9y Momentum p u a u P u75 P pva 275 Conservation a c ( /T , COT __ FOT Ou PCT Energy Ox OY  + Y 1 28+ Conservation P. j i P, i Yr O~ T _I U 2 _u 7cu S DAm (ht, P 2 pu puv 1 MW 8y By 8x 8y The unit analysis also determines the units of properties and variables found in the governing equations. Units of each fluid property are compiled in Table 22. Table 22. Units of the governing equations. Property Variable Units Property Variable Units Differential Production km dx, dy m re, kmol/m3S step sizes d dy rates Ckmo Diffusion Thermal W/2 D,, m2/s k W/m K coefficients m conductivity Density kg/m3 Specific heats c J/kgK Enthalpies h, J/kmol Temperature T K Mass Velocity facts Y, dimensionless Veloci, m/s fractions components U, Molecular A k l Viscosity Pas weights PMW kg/kmol weights Pressure p Pa CHAPTER 3 PROGRAM METHODOLOGY Flow variables are calculated throughout the catalytic reactor via a stepbystep process of solving the governing equation. The process begins by creating a discrete mesh of points to numerically solve the equations. Solutions to the flow variables are recorded at these points. The next step involves establishing parameters and conditions of the reactor and fluid. These values characterize the reactor and initial conditions of the flow. Folders to import and export data must also be defined. Once the first three steps are completed, the code can begin to find the solutions. The program sets the initial conditions and solves the simplified governing equations in Table 21. The code is written in MATLAB and consists of a main program with three subprograms. One of the subprograms finds the initial velocity components. The main program creates the mesh, finds fluid properties, and sets conditions needed to solve the equations. Information from the main program is sent to the other two subprograms. Solutions are found by the subprograms and sent back to the main program where it is saved in the solution variables. MATLAB is chosen above other programming languages because of its builtin ability to handle vectors, vector operations, and partial differential equations. MATLAB incorporates several computational tools capable of solving partial differential equations. A function called pdepe is used to solve the momentum equation as a single equation. It is also used to solve the energy equation and species continuity equations as a set of coupled partial differential equations. This makes MATLAB well suited for modeling chemically reacting flows. Another useful property of MATLAB is its compatibility with Cantera. Cantera is a free software package developed by Professor David Goodwin at the California Institute of Technology to solve problems concerning chemical reactions. The main MATLAB program calls upon this software to determine the thermodynamic, transport and chemical kinetic properties of the flow and catalytic surface. Cantera is able to construct objects of different phases and tie the phases together through an interface. This allows for the chemical interaction between the gas and surface [17]. Several studies attempt to model catalytic combustion similar to this model. A study by the National Institute for Advanced Transportation Technology at the University of Idaho modifies an existing code. Lawrence Livermore National Laboratory provides the existing Hydrodynamics, Combustion, and Transport (HCT) code. The finite difference code, HCT, utilizes the same principles of conservation for its calculations. Dissimilarity occurs in the application of the governing equations to the onedimensional timedependent catalytic combustion, opposed to the twodimensional steadystate catalytic reactor modeled by this program. Still, the study offers some insight into the chemistry and equations involved with modeling a catalytic combustor [5]. In a second study, Chou et al. [4] uses CURRENT with CHEMKIN and SURFACE CHEMKIN software to model a twodimensional monolith catalytic combustor. CURRENT is a code developed by Winters et al. [18] for low Mach number chemically reacting flows. The study discusses the chemistry and boundary conditions of the model and compares the calculations to experimental data. This program uses a similar symmetric boundary condition at the centerline. Discretization The twodimensional computational space of the reactor is broken down and discretized before the equations can be numerically solved. The mesh is that of a planar geometry with the height determined by the reactor's radius. The upper boundary is moved to the centerline and the lower boundary is still the catalytic surface. This reduces the height of the computational space and in turn reduces computation time and memory used by the computer. Now is a good time to mention that the centerline is assumed to be a streamline and symmetric conditions are assumed to exist at this boundary due to the two identical plates modeling the surface of the pipe or channel. This assumption affects the boundary conditions discussed in the section Solving Governing Equations. The length of the reactor is broken down into stages, the first stage being the entrance. This is also intended to reduce the time needed to calculate a solution. It is expected that the flow changes relatively fast in the beginning of the reactor when the catalyst is first encountered. This corresponds to the first few stages of the computational space. To help resolve the solution in these stages a smaller differential step size in the xdirection is chosen. Once the properties reach a quasisteady state, the step size can be increased to help lower the computation time. An orthogonal mesh is created for every stage. Each stage has its axial direction discretized in a linear manner, where every point is an equal distance apart. The distance is set for a given stage but can change from stage to stage. This allows the user to adjust the axial step size of a stage if the program cannot converge on a solution. A possible source of this problem is a significant change of flow properties in the xdirection. Recall that the governing equations are simplified based on the assumption that the characteristic length is much larger than the characteristic height. In other words, the vertical gradients are much larger than the axial gradients. While this assumption is still applicable, there may be areas where a change in step size is needed, such as the reactor's entrance. The point separation in the ydirection, in contrast to the axial point placement, is the same throughout the reactor. Although the vertical point placement must remain the same for every stage, it is not restricted to only a linear displacement. The point displacement is set as a power of the point location. For example, setting the power to one would position the points linearly. Setting the power to two creates a quadratic point displacement, leading to more points near the surface. A larger power places more points near the surface. Varying the power allows the user to control the location of the points in the vertical direction. This aids the program in resolving the varying chemical composition near the surface. The catalytic surface serves as the main source of the chemical reaction in the flow. Therefore, it is expected that most of the chemical change will occur near the surface. More points are needed near the surface to determine the change in the chemical composition in the vicinity of the catalyst. A tight mesh near the surface also helps resolve the fluid velocity boundary layer. Velocity, temperature, and composition variables are not found for the entire stage at once. Instead the stage's mesh is broken up further into minimeshes. A minimesh contains all the vertical points for a group of three axial locations. Governing equations are solved one minimesh at a time due to the coupling of the equations. The pressure, temperature and mass fraction of a minimesh must be approximated prior to solving the equations. Jumping ahead might seem premature because the governing equations meant to calculate the variables have not been solved yet. However, fluid properties dependent on the solution are imbedded inside the equations. These properties must be established in order to solve the equations. One can now being to appreciate the complex coupling of the governing equations. Once the equations are solved, the program updates the variables and moves downstream to the next minimesh. Parameters and Conditions Parameters and conditions of the catalytic reactor and incoming flow are set inside the code of the main program. All of the values, composition being the only exception, must equal a real scalar. The computer code begins by setting parameters of the catalytic reactor, such as the radius, stage length, stage number, surface temperature, and the distance of the nonreactive surface. Dimensions of the computational space are constructed with the height and length of the stage. The stage number is simply the sequential numbering of each stage for which a solution is calculated. The value of this number determines the data used to set the incoming conditions and the output folder in which the export files are stored. This is discussed further in the sections Initial Conditions of a Stage and Input and Output Files. In the energy equation, the temperature at the wall or surface boundary is held constant at the value entered. The distance of the nonreactive surface refers to the entrance of the reactor where there is no catalyst on the surface. This is only important for the first stage and can be ignored for any other stage. The differential step sizes in the vertical and axial direction are also set at this point, along with the power used to discretize the vertical direction. These values are used to construct the twodimensional mesh described in the section Discretization. After the parameters of the reactor are entered, the conditions of the incoming flow are defined. The speed, temperature, pressure, and composition of the flow entering the reactor are established. The incoming flow is assumed to be a uniform flow where the velocity is purely in the axial direction. Initially there is no vertical component to the fluid's velocity and the speed of the xvelocity is the same at every point. Therefore, only a single quantity is needed to define the velocity vector entering the reactor. The flow's chemical composition is initially modeled as a wellmixed fluid. This simply means that the species mass fractions are also the same at every point entering the reactor. The composition is the only value entered as a string variable. This string contains the name and mass faction of the species present in the incoming flow. Cantera reads the string to set the composition of the gas. The program uses these values to set the initial conditions of the variables. The last parameter to set is the PC variable, also a scalar. This variable controls whether the program iterates on a solution and if so, how many times the iteration takes place. An inherent delay in the solution process exists because the governing equations are decoupled. The delay is exaggerated by properties that are dependent on the solution inside the equation. Some of these properties include density, viscosity, and diffusion coefficients. With no iteration (PC equal to one), the program numerically solves the governing equations for one minimesh. Then the program updates the variables and properties and moves one differential step downstream to solve the equations at the next three axial locations. The program is continually updating the properties prior to moving downstream; therefore the delay is expected to be small. To improve the calculation one may choose to iterate on a solution using a predictor/corrector type method. To iterate on a solution the PC variable is set to quantity greater than one. For example, the program iterates once on a calculation if the variable is equal to two. Iteration occurs by solving the equations and updating the properties with the known solution. This could be seen as a predictor step, now to correct the calculation. Instead of moving one differential step downstream the program recalculates the solution for the same three axial locations based on the updated properties. If the PC variable is three, the iteration occurs twice, and so on. Input and Output Files Input and Output file names are given prior to operating the program to direct import and export data. The input text file contains the chemical data Cantera require to model the gas and solid of the catalytic reactor. Naming the input file informs the program where the chemical data are located to import. Data determine properties found within the governing equations. Only the filename of the input file is needed if it is located in Cantera's current working directory. This directory is initially set as the data folder inside Cantera's main folder, which is installed with the free software. The pathname of the output folder provides the program the location of the export folders. Export folders must be created inside the output folder and given the name Stagel, Stage2, Stage3... etc. The solutions of a stage are recorded in the folder with the corresponding number. Therefore, the export folder of a stage must exist before seeking the solution of that stage. The entire pathname is stored in the string variable saveFile. Not only is this string variable used to export solutions of a stage; it is also used to import initial conditions for most of the stages. This is discussed further in the section Initial Conditions of a Stage. Considerable amounts of data are required to model the gas and solid of the catalytic reactor. Cantera accesses this data via the input text file specified. These files contain information on the chemical kinetics, thermodynamics, and transport properties of many different species. Data consistent with the modified Arrhenius function determines the chemical kinetic properties of the gas phase reactions. This data include activation energy, preexponential coefficients and temperature exponent. In addition to this, surface reactions apply reactive sticking probability. The thermodynamic properties are determined using a NASA polynomial parameterization or Shomate parameterization. Coefficients of either parameterization are incorporated in the data of the input file. Information needed to calculate transport properties based on either a multicomponent or mixtureaveraged transport model is also included. The multicomponent transport model provides a more accurate solution than the mixtureaveraged model. However, the multicomponent model requires more data and computation time than its counterpart [17]. The program saves several variables to the output folder for every stage. The value of each variable is saved as a double precision scalar, vector, or matrix in an ASCII file. The axial location, axial velocity, pressure, temperature, mass fraction of each species, pressure gradient, and vertical velocity are all stored in the export folder. The axial location is saved in order to keep track of which discretized points in the mesh the various solutions correspond. The xlocation is saved as a vector that begins at zero and ends at the length of the stage. The axial velocity is recorded at every point in the stage and the variable is saved as a matrix. It is possible to record the vertical velocity as a matrix in a similar manner with little addition to computation time. This is due to the fact that the variable is already determined to solve the governing equations. However, the y component of the velocity is so small when compared to the xcomponent that it is not recorded as part of the solution. This will help retain memory space for the other properties. Two independent thermodynamic properties are recorded to thermodynamically define the fluid. One of these properties is the pressure, which does not vary in the vertical direction. Being only onedimensional, the pressure at each axial location is recorded and the variable is saved as a vector. The other thermodynamic property is temperature and it remains a function of both dimensions, x and y. Temperature at every discretized point is calculated and saved as a matrix in the stage's export folder. The composition of the fluid is recorded as mass fractions of each species. Like the temperature, the mass fractions are a function of both dimensions. The mass fraction of each species is saved as a matrix into its own file. As a result, the number of mass fraction files saved in the export folder is equal to the number of species, N. All the properties needed to kinetically and thermodynamically define the flow are recorded. The only other variables saved are the pressure gradient and vertical velocity. The pressure gradient is recorded as a scalar and the vertical velocity is saved as a vector. Both variables correspond to the last axial position of the stage and are used as initial conditions of the next stage. Initial Conditions of a Stage The program can operate once all of the parameters, conditions, and file names are designated. Initial boundary conditions of the stage's first minimesh are established before solving the governing equations. Velocity components at all vertical points in the first axial location are required to define the momentum equation and its initial boundary condition. The pressure, temperature, and composition in the first minimesh must also be defined to estimate the properties inside the governing equations. The model requires only one pressure value per xlocation, because the pressure is independent of the vertical direction. The result is only three scalars being required to define the pressure in the mesh. Temperature and species mass fractions are twodimensional and must be set for every point in the stage's first minimesh. The process used to define these variables depends on the stage number. Conditions of the initial stage or stage one are based on the values discussed in the section Parameters and Conditions. Every other stage uses the solution of the previous stage to set these initial boundary conditions. The program can begin to solve the governing equations once the initial conditions are set. Cantera creates a gas object and surface object prior to defining initial conditions. The gas is adjusted to the pressure, temperature, and composition entering the reactor and the two objects are connected through an interface. The gas object is created at this time because Cantera provides a simple means to set the composition variable of the first stage. Only the composition's string variable is needed to establish the mass fraction of all the species initially present. Cantera can take the composition of the gas object and return the mass fraction of every species. This is much easier than searching for the species not present and setting their mass fraction to zero. Cantera also ensures that the sum of the mass fractions equals one. Stage One The reactor is characterized by the absence of a catalytic surface at its entrance. The catalyst does not begin until further downstream. This is where stage one begins and the initial boundary conditions of the first minimesh are determined. Minimal change in the conditions should occur over the nonreactive surface with the exception of the two velocity components. Therefore, the initial conditions of the temperature, composition, and pressure remain the flow conditions entering the reactor. Temperature and mass fractions at the first three axial locations are approximated by the values entered as initial conditions. The surface temperature is set to the value entered as a parameter. The initial pressure is equal to the pressure of the incoming flow, and the pressure at the next two differential steps is calculated with the pressure gradient. In contrast to the other variables, the surface affects the velocity vector. A boundary layer develops changing the profile of the axial velocity, which produces a vertical velocity. Blasius solution is used to model the boundary layer and determine the two velocity components. Blasius Solution Axial velocity is quantified by two values at a point in the beginning of the reactor. The singularity point is located on the front edge of the reactor, where the incoming flow first encounters the surface. A finite value is given to the uniform velocity entering the reactor. The velocity at this point must also equal zero due to the boundary conditions of the velocity. To overcome the singularity point, Blasius solution is used to calculate the two velocity components at the end of the nonreactive surface, where the first stage begins. H. Blasius is well known for obtaining an exact solution to a laminar boundary layer flow over a flat plate. Blasius is able to find a similarity solution to the continuity and momentum equations through proper scaling and nondimensionalization of the two equations. In his solution, the dimensionless stream function replaces the two velocity components as the dependent variable. The two coordinates, x and y, are also combined into one dimensionless independent variable. Blasius transforms the two partial differential equations into one ordinary differential equation. A power series expansion or numerical methods can then be used to solve the third order, nonlinear equation. The dimensionless stream function and its derivative are used to calculate the axial and vertical velocity [15]. Blasius solution describes a twodimensional, steady, incompressible flow with no pressure gradient. Recall that the assumption of constant density is not applicable to this model due to the chemistry involved. A pressure gradient equal to zero is also not accurate because the flow is assumed to have pressure changes in the axial direction. However, Blasius solution is used as a reasonable estimate to the velocity profile over the nonreactive surface. The change in density is caused mostly by the catalytic surface, and the catalyst is not present in the region that Blasius solution is employed. A change in density from species diffusing upstream is possible, but the effect should be negligible. The production of a new species with a diffusion velocity great enough to overcome the axial velocity is needed for this to occur. The behavior of the pressure gradient also permits the use of Blasius solution. The pressure slowly decreases as the flow moves downstream. This produces a decreasing pressure gradient that has an initial value of zero. The change in pressure is small and a zero pressure gradient should be a reasonable approximation at the entrance of the reactor. Fortunately a nonreactive surface is located in the region of the singularity point. Blasius solution can be used to generate the velocity profile at the end of the nonreactive surface. The main program calls on one of the subprograms, a function called Blasius. The function imports the axial differential step size, stage length, ycoordinate vector, location where the catalytic surface begins, viscosity, and initial speed of the flow. A shooting method determines the dimensionless stream function and its derivative. Because the equation developed by Blasius is dimensionless, the calculated values of the stream function are independent of the values imported. The axial velocity vector and vertical velocity vector are determined using the dimensionless variables with the imported variables. The axial velocity is treated as the initial condition of stage one. The vertical velocity is used in the momentum equation. Note that the function Blasius is only used in the first stage, where the nonreactive surface is located. Other stages use the solutions of the previous stage to set the velocity profile. Subsequent Stages If the stage number is greater than one, initial boundary conditions are taken from the export folder of the previous stage. The program locates and loads initial quantities with the saveFile variable. The last value of the preceding stage's pressure vector becomes the initial pressure of the current stage. The vertical velocity vector and pressure gradient (scalar) are loaded from the preceding stage's output files. Initial values of the current stage's axial velocity are set to equal the xvelocity at the end of the last stage. Temperature and composition variables are all that remain to import. Recall that the temperature and mass fraction must be set for all three axial locations. The first axial location is equal to the last axial location in the previous stage's matrix, similar to the axial velocity variable. For the next two differential steps downstream, the xlocation vector of the last stage is imported. This vector along with the temperature and mass fractions of the previous stage determine the variable's gradient at the end of the last stage. A secondorder backwarddifference formula is used to estimate this gradient [19]. Values of the next two axial locations in the variable of the current stage are linearly extrapolated using the gradient. With the initial conditions of the stage set, the program is prepared to solve the equations controlling the behavior of the flow. Solving Governing Equations Governing equations are solved for three axial locations at a time. Remember that the equations contain properties dependent on the solution. Solving one minimesh at a time allows the program to update properties inside the equations prior to calculating the solution one step downstream. For the same reason, the pressure, temperature, and composition at these three places must be approximated before solving the equations. For the stage's first minimesh, the method for these approximations is described in the section Initial Conditions of a Stage. Approximations of the other minimeshes are calculated with the solution of the previous minimesh. Note that two of the axial locations of the next minimesh are repeated locations of the previous minimesh because the program moves only one differential step. Two MATLAB functions, or subprograms, are written to solve the momentum equation and the species/energy equations separately. These two functions are named Momentum and SpeciesEnergy. The process of the solving the governing equations for a stage is illustrated in the flow chart of Figure 31. Once the initial approximations for a single minimesh are established, the program begins with the momentum equation in the xdirection. The main program calls on Cantera to find properties embedded in the momentum equation. Information is exported to the subprogram Momentum, which is then used to solve the equation. The subprogram sends the solution, the axial velocity, back to the main program. The mass or continuity equation calculates the vertical velocity and confirms that the solution of the momentum equation is accurate. The pressure gradient is updated and the momentum equation is solved again if the boundary condition of the vertical velocity at the centerline is not reacted. Iftheyvelocity equals zero, the species continuity and energy equations are solved together with the Species_Energy function. Again, the data needed to solve the system of equations are provided by Cantera and sent to the subprogram. This time the composition and temperature are sent back as the solution. Once the solutions are calculated and any iteration correction is performed, the program updates the variables, saves the data in the export folder, and moves one differential step downstream. The temperature and composition of the next axial location are predicted using a linear extrapolation from the previous two xlocations. The program loops back and solves the equations for the next minimesh. This continues until the end of the stage is reached or the program cannot resolve the solutions. IninitiaInitial Approximations (next minimesh) Continue until ) End of Stage Conditions to a minimesh end of Stage of a Stage I T e1 Momentum Ia , II c Solutions saved SI Modify to export folder (xvelocity) dp dxeotode CANTERA ass yvelocity 01 (yvelocity) Test Boundary Condition yvelocity = 0 [ Species/Energy (comp/temp) Figure 31. Flow chart for single stage modeling. The subprograms Momentum and Species_Energy use a partial differential equation (PDE) solver provided by MATLAB to numerically solve the equations. The solver numerically computes the momentum equation as a single equation. It is also used to solve the energy equation and species continuity equations as a set of coupled partial differential equations. The solver, named pdepe, calculates the solution of partial differential equations with the form shown below in Equation 31. The gradient of the dependent variable with respect to x is multiplied by a coupling term, c. This term along with the flux term, J and the source term, s, are functions of the two independent variables, the dependent variable and its vertical gradient. a z}az a Q at z 8 z' c x,y, x f x,y,z +s x,y, z, (Equation 31) The symmetry parameter, m, along with the coupling, flux, and source term define the PDE. The mesh spacing of the two independent variables, one initial condition, and two boundary conditions are all that remain to solve the equation. The mesh spacing is simply the minimesh discussed in the section Discretization. Function pdepe selects the xmesh dynamically to resolve the solution, but only reports the answer at the mesh points specified. Strictly speaking the initial condition is a boundary condition. It is the value of the dependent variable at the first of the three axial locations. The initial boundary condition of the dependent variable needs to be given as a function ofy. The other two boundaries are found at the catalytic surface and centerline. Both must fit the form shown in Equation 32. Boundary conditions are expressed in terms ofp, q, andf The flux termfis already defined in the PDE above, so only p and q are needed to establish the boundary conditions. p(x,y,z)+q(x,y)f x,y,z,2 =0 (Equation 32) Some fluid properties, such as density, are converted into functions ofy to conform to Equation 31 above. Most of these properties are functions of both dimensions. However, properties dependent only on the vertical direction should be an acceptable representation for several reasons. Fluid properties vary more in the vertical direction than the axial direction and are a stronger function ofy. In addition to this, the functions only need to represent the properties for the three xlocations of the minimesh. Initial boundary conditions also need to be turned into functions ofy. Builtin MATLAB functions spline and unmkpp generate the function representations. Properties are found at every vertical point in the middle of the three xlocations and saved in a vector. In the case of the initial condition, the vector contains the initial values of the dependent variable. The function spline uses the vector to create twenty separate piecewise polynomials of the form of the cubic spline. Function unmkpp extracts the four coefficients of the each polynomial and saves it into a fourbytwenty matrix for each representation. The matrices are exported to either the Momentum function or SpeciesEnergy subprogram to reconstruct the piecewise polynomial. The coefficients and a heavyside step function connect the piecewise polynomials inside the subprogram. The result is a smooth function representation of the initial conditions or fluid properties embedded in the governing equation. Momentum Equation Solving the momentum equation begins by guessing the pressure gradient. Pressure at the three axial locations is determined with the guessed pressure gradient. The pressure along with the approximated temperature and composition are used to determine the density and viscosity. These are the properties found in the simplified momentum equation. 9u 9 Su1 9p 9u pu / u pv (Equation 275) Ox dy y Ox dy Properties are determined at every vertical point in the middle of the three x locations and transformed into functions ofy for the coupling, flux, and source terms. Comparing the momentum equation to the form used by the pdepe function, it is evident that the axial velocity replaces the dependent variable, z. The symmetry parameter, m, is zero. The coupling, flux, and source terms equal Equation 33, Equation 34, Equation 3 5, respectfully. c x,y,uI = pu (Equation 33) fxy, \= =P (Equation 34) s xyu = pv (Equation 35) Oy ) x oy The coupling term, c, contains the density and axial velocity. This term is allowed to be a function of the dependent variable. As a result, only the density must be transformed into a function. The flux term,f equals the viscosity multiplied by the vertical gradient of the axial velocity. To fit Equation 31, the gradient will remain but the viscosity is represented by a function of the vertical direction. The last two terms in the momentum equation are combined into the source term. These two terms consist of the pressure gradient and the product of the density, yvelocity, and ygradient of the dependent variable. Recall that the pressure is only a function of the axial location. Therefore, the pressure gradient remains constant at a given xlocation and does not need to be transformed into a function ofy. The vertical gradient of the xvelocity is allowed inside the source term. The density and yvelocity product is the only element of the source term transformed into a function. The axial velocity's initial boundary condition is also transformed into a function ofy. Boundary conditions at the surface and centerline are all that remain to solve the momentum equation. Appling the noslip assumption, the axial velocity is zero on the catalytic surface. The centerline of the reactor is assumed to be a streamline with a vertical gradient of the xvelocity equal to zero. Equation 36a and 36b show these two conditions in a form recognized by the pdepe function. u= Ou+(0).f =0 (Equation 36a) = 0 (0)+(1) p =0 (Equation 36b) ay 9y Boundary conditions of the momentum equation are defined inside the subprogram Momentum. At the surface, ory equal to zero, p equals the dependent variable and q is zero. The centerline condition dictates that equals zero and q equals one. This sets the flux term to zero at the boundary. The flux term is the product of the viscosity and vertical gradient of the axial velocity. Since the viscosity is finite, the gradient must equal zero, which is the condition sought. The Momentum subprogram can now be used to solve the momentum equation. The function imports the discretized mesh and guessed pressure gradient. Coefficients of the initial boundary condition, coupling, flux, and source terms are also imported. These are the coefficients of the piecewise cubic polynomial. The secondorder, nonlinear PDE is solved and the axial velocity at the three axial locations is returned to the main program. Continuity Equation The axial velocity solution must be verified because the value of the pressure gradient is assumed. This value directly affects the momentum equation by being part of the source term. It also indirectly affects the solution by changing the properties dependent on the pressure. Equation 23 is the mass or continuity equation that calculates the vertical velocity. At the same time, the solution of the mass equation acts as a check to the momentum equation. First, the gradient of the density and xvelocity product is determined at every vertical point at the end of the minimesh. This partial derivative is calculated with a secondorder backwarddifference formula. The ygradient of the density is found with a secondorder central difference formula with varying spacing. Once these two gradients are found, the partial derivative of the vertical velocity is approximated with another secondorder central difference formula with varying spacing in Equation 37 [19]. Solving for the velocity at the next mesh point produces Equation 38. (v Q(pu) 8p  V) (Equation 23) ay 8x ay v l+(a2 1)V a2 1 P a(a +)(y y 1) where a (= yJ) y y ) a(pu) ap ax 8y (Equation 37) a(pu) _ap 8x yJ 8 a(a +1)(y y) (Equation 38) (Equation 38) I P(a2l )v pa2v1 P a(a+l)(yJ yJ _1) Equation 38 is used to find the vertical velocity component at the last of the three xlocations. This new vertical velocity becomes the variable used by the next minimesh downstream. The species and energy equations still use the original yvelocity for their calculations. Once the yvelocity is found at every vertical point, its value at the centerline is checked. Being a streamline boundary condition, there should be no flow across the boundary and the yvelocity should roughly equal zero. If the velocity does not meet this requirement, the pressure gradient is adjusted and the program loops back to the momentum equation. The amount of the adjustment is proportional to the size of the yvelocity at the centerline. A weighted correction modifies the pressure gradient. This continues until the centerline yvelocity is less than one tenthousandths. At which point the axial velocity of the minimesh is saved or spliced to the axial velocity variable of the entire stage. The program then moves on to the remaining two governing equations. Now is an excellent moment to discuss the reasoning behind breaking apart the momentum and mass equations from the other governing equations. It has already been shown that all of the equations are highly coupled and should be solved as such. However, that approach leads to a very problematic and timeconsuming calculation due to the unknown pressure gradient. Solving the entire group of equations until the correct pressure gradient is found would take a great amount of computing time. Decoupling the momentum and mass equations significantly reduces the time of the calculation. This method does not come without its disadvantages. Separating the governing equations creates a delay in the solution. This delay can be overcome with the iteration process already discussed in the section Parameters and Conditions. Species Continuity Equations The remaining equations are not decoupled, but instead are solved simultaneously by the function pdepe. The set of species equations are solved for the mass fraction of each atom or molecule. The number of equations in this set is equal to the total number of species in the model, defined as N. Equation 246 shows the simplified species equation. pU  pD \pv +a, iMW = 1,2,3,...N (Equation 246) ax oy oy oy Mass fraction of species i is the dependent variable of the species equations. Again, the symmetry parameter, m, is zero. Cantera determines the density, diffusion coefficients, net gas production rates, and molecular weights found in the equation. These quantities and the two velocity components are used to create the coupling, flux, and source terms shown in the three equations below. c x,y,Y, \= pu (Equation 3 f x, y, =pDm y (Equation 31 s x, y,Y Y = pV +aMW, (Equation 31 S 9y } 9y 9) 0) 1) The density and axial velocity make up the coupling term. Axial velocity is no longer the dependent variable as it is in the momentum equation. This means that the multiple of the density and axial velocity must be transformed into a function ofy. The flux term,f, can be found inside the parenthesis of Equation 246. It equals the density multiplied by the species mixtureaveraged diffusion coefficient and the vertical gradient of the mass fraction. The flux term is allowed to be a function of dependent variable's vertical gradient. Therefore, only the density and diffusion coefficient product is represented by a function. The source term becomes the combination of the last two terms of the species equation. This term equals the species mass production minus the product of the density, yvelocity, and vertical gradient of the dependent variable. The product of the density and yvelocity is transformed into a function representation, while the vertical gradient is left unaltered. This yvelocity is the original vector and not the velocity found from the mass equation. The species mass production and initial boundary condition are also changed to a function of the ydirection. Two boundary conditions of the species equations can be connected to the species mass flux. At the surface boundary is a heterogeneous catalyst where species can be created or destroyed. Assuming a steadystate model, the species flux can be equated to the production rate on the catalyst. Species mass flux into the surface equals the destruction rate and the flux away from the surface is the creation rate. Mathematically written in Equation 312 and reorganized into Equation 313a to fit the form defined by the PDE solver. Cantera determines these production rates for each species. A symmetric boundary condition is applied to the upper boundary. Resulting in the vertical gradient of a species mass fraction approximately equaling zero at the centerline. Equation 313a and 313b show these two conditions in a form recognized by the pdepe function. 8Y S,sufaceMW = pDm __ (Equation 312) (,suaceW )+ (1). pD, = 0 (Equation 313a) Y( )Y = 0 (0)+ (1) .D 0 (Equation 313b) Parameterp equals the mass production rate and q is one for the lower boundary at y equal zero. The upper boundary condition hasp equal to zero while q equals one. This sets the gradient equal to zero because neither the density nor the diffusion coefficient of the flux term equal zero. The system of partial differential equation is defined along with their initial and boundary conditions. Before the species equations are solved, the energy equation is added to the group. Energy Equation The energy equation contains kinetic energy terms defined by the velocity field. Kinetic energy terms are in the form ofx andy gradients of the axial velocity. The solution of the momentum equation is used for the xvelocity. Both partial derivatives are calculated at every vertical point in the minimesh. A simple secondorder central difference formula is applied to estimate the axial gradient. The vertical gradient is a little more complicated because the spacing in the ydirection may vary. A secondorder centraldifference formula that is modified for varying point spacing is used for the core of the calculations. Theygradient at the surface is found with either a firstorder or secondorder forwarddifference formula for equal spacing. If the y spacing is linear, then the secondorder formula is used. The firstorder equation is used if the spacing is nonlinear [19]. Although it is first order, the error should be small because the spacing near the surface is tight. The gradient at the centerline equals zero due to the boundary condition of the momentum equation. Equation 286 is the simplified energy equation with the two kinetic energy terms at the end of the equation. OT 8 F T + u OT puc k+up pc + X Y L (Equation 286) C LY IT a 2u Bu p1 D,  (ho pu puv =1 1 By 8y 8x By Equation 286 contains many thermodynamic and transport properties that need to establish. Cantera retrieves the properties at every vertical point in the middle of the three xlocations and saves them to vectors. The main program uses the vectors to create the function representations of the coupling, flux, and source terms. It is apparent that the temperature is now the dependent variable of the PDE. The symmetry parameter is zero and the coupling, flux, and source term are listed below. c x,y,T, = puc (Equation 314) I qy ) ( 9T' T ku f x, y,T, k+ (Equation 315) OT) OT c OY BT s x,y,T, = VC P +P D Cv+m  d S(Equation 316) Y (how) pu puv The coupling term is the multiple of the density, axial velocity, and specific heat. The entire expression is transformed into a function y. The flux term, found inside the brackets, is made up of two parts. The first equals the thermal conductivity times the temperature's axial gradient. Second is the combination of the viscosity, xvelocity and its ygradient. The two parts must be kept separate for the flux term to remain a function of the temperature gradient. Thermal conductivity is turn into one function, while the second part is turned into another. The remaining five terms are grouped into the source term. For the first representation, the product of the density, yvelocity, and specific heat are changed to a function of y, and the temperature gradient remains a variable. The second term consists of a complicated sum containing the species mass fraction gradient. This is where the coupling between the governing equations directly takes effect. Calculation of the sum is addressed in the section Species/Energy System. The third expression is the other sum in the equation. However, it is not nearly as difficult as the last because it does not contain any of the system's dependent variables. This sum is simply calculated in the main program and added with the last two kinetic energy terms. The last three terms in the energy equation are combined into a function representation. Initial boundary condition of the temperature is transformed into a function for the subprogram. Boundary conditions are established for the temperature at the surface and centerline. The temperature at the catalytic surface (y=0) is held constant. This lower boundary condition is shown in Equation 317a. The upper boundary condition is characterized by no heat flux. Temperature of the flow is uniform when it enters the reactor. At which point, the catalyst induces chemistry in the flow and heat production occurs at the surface. The temperature begins to increase at the surface and slowly expand up to the centerline. A thermal boundary layer is created and heat flux across the streamline is zero until the layer reaches the centerline. A long distance is needed for this to occur and the heat flux remains zero at the streamline for the short distance of the reactor. Equations 317a and 317b show the conditions in a form recognized by the pdepe function. T= ,ace (T ac)+(0). f = 0 (Equation 317a) OT OT Ou T= 0 (0)+(1). k +u/ = 0 (Equation 317b) In Equation 317a, p is the dependent variable minus the temperature at the surface and q equals zero. This produces the constant value at the surface. The upper condition is created with parameterp equal to zero and q equal to one. This sets the flux term, which consists of two parts, to zero at the boundary. The first part is the heat flux and the second contains the gradient of the axial velocity. A problem arises because only the heat flux should be zero. However, the second term vanishes at the centerline due to the boundary condition of the momentum equation. The result is the proper symmetry condition at the centerline boundary. Species/Energy System of Equations The species and energy equations are combined for the pdepe function to solve. The resulting system of (N+1) equations is shown below. The flux and source terms are split into three parts to accommodate the various expressions in each equation. [c] '= [fl]z+ [f2]+ [f3] [sl] Z+ [s2] + [s3] (Equation 318) ax Oy y y Jy ay where Z= (Y,Y2,...,YN, T) Components of the dependent variable vector, Z, consist of each species mass fractions and the temperature. Coefficients of the cubic spline generated for the coupling, flux, and source terms are grouped together. Each component of the vectors c,fl,J3, sl, and s3 represents the cubic spline function of that component. The other two expressions (f2 and s2) generate the sums involving the system's dependent variables. The coupling term of each equation is combined into one group, c. Separating this term is not necessary because it does not contain any mass fractions or temperature variables. pu pu [c]= (Equation 319) pu pUCp The flux term is broken up into three separate collections. The first group,fl, is the combination of the species equations' flux terms and the energy equation's thermal conductivity expression. These terms are multiplied by the axial gradient of the dependent variable. The Nh component offl is set to zero and replaced with a sum inJ2. Group j3 is a result of the additional flux term in the energy equation. The first N components of this group are zero because none of the species equations contain an additional term. pD, 0 0 PD2, 0 0 [f1]= [f2] dotFnl, I [f3]= 0 1 Ou k 0 _ Figure 32. Flux components in the species/energy system. The group,J2, is a result of "the fact that the sum of all the species diffusion fluxes must be zero" [13:227]. This fact is rearranged and shown as the sum in Equation 318. It suggested that this equation "be applied to the species in excess, which in many combustion systems is N2" [13:227]. Diatomic nitrogen is the last species listed in the input file used in the tests. This corresponds to the Nth component in the dependent variable vector. Note that the Nth component infl is zero and negative one inJ2. The result is this sum replacing the massaveraged diffusion flux in the last species equation. MATLAB's builtin dot product is utilized to create the sum. The last two values ofFnl are zero because the sum does not include the Nh specie (diatomic nitrogen) or temperature. SPDmj = dot Fn l, (Equation 320) where Fnl =(pDm, pDm, ..pD, ,m0, 0,O) The source term is also grouped into three expressions. Group sl is a combination of the terms multiplying the gradient of the dependent variable. The second group, s2 produces the complicated sum discussed in the section Energy Equation of this chapter. All the remaining terms are compiled into s3. PV 0 r n MWI pv 0 1 1 pv 02 2 sl= [s2]= dot Fn2 s3= pv 0 pvc, 1 ,) PU2 u puv Figure 33. Source components in the species/energy system. The sum in the energy equation, shown in Equation 321, is a function of the species mass fractions. The energy equation is not decoupled from the species equation and the mass fractions are dependent variables in the system. To reproduce this sum, all of the properties (this excludes the species mass fractions) are transformed into functions ofy for each species and stored in a vector. The dot product of this vector and the dependent variable vector gradient recreates the sum inside the subprogram Species_Energy. The sum is then multiplied by the temperature gradient. Note that the temperature gradient is not part of the sum, consequently the last value in Fn2 is zero. yPCprm = dotl Fn2, l (Equation 321) ( pcpCDlm pcp2D2m pcpNDNm where Fn2= PC ,PCp2D ..p N,0 The boundary conditions of both equations are also grouped together. Figure 34 shows the conditions for the system of equations. For the lower boundary condition, parameterp is broken up into two groups, because the dependent variable is part of the temperature's lower boundary condition. )l,surface ~ 0 1 0 1 2surfaceMW 0 1 0 1 + : Z+ f=0 + f=0 NsufaceMIW 0 1 0 1 Tuace 1 0 0 1 Figure 34. Boundary conditions of the species/energy system. Function Species_Energy imports information to solve this system of governing equations. The subprogram imports the discretized minimesh and all the polynomial coefficients. Parametersp and q of the system's boundary conditions, which are compiled in the main program, are also sent to the subprogram. The function then calculates the solution of each component in the dependent variable vector and returns it to the main program. A mass fraction less than 1E20 is treated as error and the value is set to zero. Solutions for the minimesh are spliced to their corresponding variables for the entire stage. The program moves one differential step downstream and loops back to the momentum equation. Recall that before solving the equations for a given minimesh the temperature and composition must be defined at all three axial locations. The program linearly extrapolates these values from the previous two xlocations. The same process solves the governing equations for the next minimesh. Its solution is spliced to the stage's solution variable and the program moves on. This continues until the end of the stage is reached or the differential xstep is not small enough to resolve the solution. CHAPTER 4 TESTING A process of running the code for several cases and examining the solutions is performed in order to test the program. Four cases are used to test the software. Beginning with no chemistry in the first case, chemistry is slowly introduced to the other cases. The second case involves only gas chemistry while the last two tests include both gas and surface chemistry. Gas and surface chemical reactions are modeled at five hundred and seven hundred degrees Kelvin. Slowly introducing chemistry to the model will aid in locating errors during the debugging process should any problem arise. Each case uses the input file named ptcombust.cti, which is provided by Cantera as part of the software package. This file contains data for the methane/oxygen surface mechanism on platinum developed by O. Deutschmann. The input file ptcombust.cti calls on the file gri30.cti, also part of the Cantera package, to manage the gas reactions. The file gri30.cti contains data for the optimized GRIMech mechanism and for this program calculates transport properties based on a mixtureaveraged transport model. Once the program finds the solution for a given case, the results are examined. No experimental data is available at this time to compare to the program's solutions. However the different tests can confirm that the software produces reasonable results and is operational. Several parameters and conditions that characterize the reactor and incoming flow are similar for the four cases. The reactor has a radius or thickness of two centimeters and a length of thirty centimeters. Therefore, the height of the mesh is two centimeters and the sum of stage lengths equals thirty centimeters. The distance of the nonreactive surface at the entrance of the reactor is given the variable name Lnocat and is equal to one centimeter. In case one, a catalytic surface is not present and the onecentimeter value only determines the location of the initial velocity condition found with Blasius solution. The differential step size in the vertical direction is set at fourhundredths of a centimeter. This mesh spacing in theydirection is not linear. In order to place more points near the surface, the power discussed in the section Discretization is set to four. The PC variable is equal to one so no iteration occurs. Temperature of the surface is set to four hundred degrees Kelvin for case one and two. The initial temperature of the flow is also four hundred degrees Kelvin for the first two cases. A mixture of air and methane at one atmosphere of pressure comprise the fluid entering the reactor of every test. Case One Case one models a chemically inactive gas passing through a reactor with no catalytic surface. The flow is essentially a nonreactive flow through a pipe or channel with a pressure gradient. With no chemical reactions taking place, density remains the same and the entire system of governing equations is altered. All of the equations could be simplified for an incompressible flow and the set of species equations could be removed all together. Although modifying the code in this way would defeat the purpose of the test. To test the program only fluid properties are changed, while the code remains unaltered. Turning the gas chemistry off is achieved by equating the species mass production rate to zero. This only affects the source term, s3 in Figure 33, in the species continuity equations. Creating a surface with no catalyst also exclusively affects the species continuity equations. Production rates at the surface are forced to zero changing the lower boundary conditions in Figure 34. One stage is used to model the flow of the first case. Recall that the computational space of the reactor can be broken up into stages. Being able to change the axial differential step size of each stage allows the program to resolve the changing composition. However there is no varying composition because the chemistry is removed in this test. The reactor does not need to be split into stages for this reason. The fluid mixture, given in mass fractions in Table 41, enters with a velocity of one meter per second. Other parameters and conditions of this test run are listed in Table 41. Table 41. Parameters and conditions of case one. Reactor Parameters Initial Flow Conditions Radius 0.02 m Velocity 1 m/s Lnocat 0.01 m Temperature 400 K Stage length 0.30 m Pressure 101325 Pa Surface temp 400 K CH4:0.004, dy 0.0004 m Composition 02:0.23, power 4 (mass fractions) N2:0.752, AR:0.014 No. of Stages 1 PC 1 dx (Stagel) 0.01 m Input file ptcombust.cti Results of Case One The solution should mirror that of a viscid twodimensional laminar flow through two flat plates with the pressure slowly decreasing. Figure 41 illustrates the axial velocity profile at four different locations. As expected the axial velocity solution resembles a boundary layer flow increasing from zero at the surface to the centerline velocity. The centerline velocity increases to compensate for the loss of mass flux near the surface. This is shown in Figure 41, where the centerline velocity increasing downstream as the boundary layer grows. Note the overshoot in the velocity profile at x equal zero. This is due to the program creating a function representation of the initial velocity condition. Other than the overshoot the velocity solution is a smooth continuous model of what is expected for flow over a flat plate with changing pressure. 0.020 i 0.018 SI 0.016 0.014 x=0 0.012 x=.1  x=0.2 ' y (m) 0.010 x=0.3  0.008 0.006 0.004 0.002 0.000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 xvelocity (m/s) Figure 41. Axial velocity profiles of case one. Increase in the centerline velocity should lead to a decrease in pressure. The pressure change of Figure 42 shows this to be the condition. The pressure slowly decreases downstream from its initial value of one atmosphere. The incompressibility of case one allows the use of Bernoulli's Equation to calculate the change in pressure. This provides an alternate means of finding the pressure with the axial velocity and ensures that the solutions of the program are consistent. The velocity at the streamline or centerline, where viscous effects are not present, is used in Bernoulli's Equation. The pressure difference calculated from both the program and Bernoulli's Equation is graphed in Figure 42. The change in pressure predicted by the program and Bernoulli's Equation is very similar and the behavior is typically found in the beginning stages of a pipe or channel flow. 0.00 0.02 0.04 0.06 p program P 0.08 Bernoulli (Pa) 0.10 0.12 0.14 0.16 0.18 0 0.05 0.1 0.15 0.2 0.25 0.3 x(m) Figure 42. Pressure plot of case one. Variations in the species mass fractions should not exist because all chemistry is neglected. Temperature should also have minor changes for the same reason in addition to the low Mach number of the flow. This is the case for the first test of the program. As expected, the calculated composition and temperature remain constant. The software produces the expected solutions for all variables in case one. Case Two A flow characterized by gas reactions and no surface reactions is modeled in case two. The flow is that of a chemically reacting fluid passing over a noncatalytic surface. The only difference between case two and case one is the presence of gas chemistry in the flow. Cantera determines the value of the species mass production rates. Unlike case one, these values are not forced to equal zero. Removing the surface chemistry is achieved by altering the lower boundary conditions of the species continuity equations. Surface production rates are set to zero just as they are in case one. Three stages are used to model the flow in case two. Because the temperature of the flow is relatively low, little change in the composition is expected. However, using more than one stage will test the process involved with multiple stages. This includes the saving and loading of variables and the smooth connection of the stages. The length of each stage is one centimeter. The same fluid composition enters the reactor, but the initial velocity is now half a meter per second. Parameters and conditions of this test run are listed in Table 42. Table 42. Parameters and conditions of case two. Reactor Parameters Initial Flow Conditions Radius 0.02 m Velocity 0.5 m/s Lnocat 0.01 m Temperature 400 K Stage length 0.30 m Pressure 101325 Pa Surface temp 400 K CH4:0.004, dy 0.0004 m H dy 0.0004 Composition 02:0.23, Power 4 power (mass fractions) N2:0.752, No. of Stages 3 AR0.014 AR:0.014 dx (Stagel) 0.01 m dx (Stage2) 0.01 m PC 1 dx (Stage3) 0.01 m Input file ptcombust.cti Results of Case Two The axial velocity calculated in the second test is graphed in Figure 43 for three x locations. Similar to the first test, the velocity is recognized as a typical pipe or channel flow solution. The initial velocity is half a meter per second and the centerline velocity increases from this value as the flow becomes fully developed. The presence of gas chemistry does not appear to affect the solution of the momentum equation. Variation in the composition is not anticipated and the velocity profile is comparable to that in case one. The reduction in the initial velocity does remove the overshoot found in Figure 41. 0.020 0.018 0.016 x=0 0.014 x=0.1 0.012 x=0.2 x=0.3 y (m) 0.010 0.008 .0 0.006 _,_ 0.004 0.002 0.000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 xvelocity (m/s) Figure 43. Axial velocity profiles of case two. Minor variations in the density do not change the velocity, meaning the pressure should also behave the same. Figure 44 illustrates that the pressure decreases downstream much like the pressure in case one. A difference in the program's solution and Bernoulli's solution is noticeable and there is almost a twentytwo percent difference between the two. The general behavior of the pressure is consistent with expectations; however, the software produces values that are not validated by Bernoulli's Equation. As expected, the temperature remains constant at four hundred degrees Kelvin. Some changes in the temperature do occur but are very small and can be considered numerical error. Some of the mass fractions also contain small fluctuations. Figure 45 shows the change in the mass fraction of the species methane. While some of the species mass fractions behave oddly, it is most likely a product of numerical error. The program has a secondorder accuracy and the largest step size equals one centimeter. The 78 resulting error has the size of one tenthousandths, which is greater than the error seen in Figure 45. 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0 0.05 0.1 0.15 0.2 0.25 x (m) Figure 44. Pressure plot of case two. 0.02 0.018 0.016 0.014 0.012 y (m) 0.01 0.008 0.006 0.004 0.002 0 3E12 2E12 1E12 0 1E12 reduction of CH4 mass fraction Figure 45. Reduction in methane concentrations of case two. 2E12 The software produces reasonable solutions for a chemically reacting flow without a catalyst. The low temperature in this test results in little gas reactions and the mass fractions remain nearly constant. A large incoming temperature will lead to combustion of the fuel/air mixture. Care is taken to avoid combustion because the governing equations are reduced based on the expectation that characteristic length scales in the axial direction are large. The smooth connection of multiple stages is also confirmed by the second test. This can be seen in Figure 44 where the pressure is a continuous function ofx. At this point the code calculates expected values for a flow with and without chemical gas reactions. Case Three A complete test of the software is performed in case three where gas and surface chemistry both exist. As it is originally intended, the model is that of a chemically reacting fluid flow over a catalytic surface. Cantera finds the gas and surface production rates used by the set of species continuity equations. Unlike the previous two tests, these values are not forced to equal zero. Table 43. Parameters and conditions of case three. Reactor Parameters Initial Flow Conditions Radius 0.02 m Velocity 0.5 m/s Lnocat 0.01 m Temperature 500 K Stage length 0.30 m Pressure 101325 Pa Surface temp 500 K CH4:0.004, dy 0.0004 m Composition 02:0.23, power 4 (mass fractions) N2:0.752, No. of Stages 2 AR:0.014 dx (Stagel) 0.001 m PC 1 dx (Stage2) 0.001 m Input file ptcombust.cti Minimal change in the composition is encountered due to the low temperature and only two stages are applied. The temperature of the gas and catalytic surface is increased to five hundred degrees Kelvin and the differential step size in the axial direction decreases to one millimeter. The other parameters and conditions are similar to the second test and all are found in Table 43. Results of Case Three At first the original software does not obtain a solution for the entire flow in case three and changes are made accordingly. Two centimeters into the reactor the program is unable to resolve the changing composition and the code prematurely terminates. This problem is found to be associated with the application of the dot products in the subprogram SpeciesEnergy. Once these dot products are removed, the program is able to solve the entire computational space. The dot product of Equation 320 is replaced with the mixtureaveraged diffusion coefficient. The Nh species equation is now similar to the rest of the species equations. The other sum, Equation 321, is no longer performed by the dot product but is found by adding the function representations of each species. These two dot products create the sums involving the gradients of the species mass fractions. The problem is not noticeable in case one because the change in the mass fractions is zero. This problem might be the source of the odd behavior seen in some of the mass fractions and pressure difference of the second test. The software is able to model the entire reactor after the corrections are made. The velocity profiles graphed in Figure 46 are typical of boundarylayer growth in the presence of a pressure gradient and are consistent with the models of the first two tests. Again, the velocity increases at the centerline and the boundary layer grows as the fluid moves downstream. 0.02 0.018 x=O , 0.016 x=0.05  0.014 . x=0.1 : _ 0.0 x=0.2 0.012 , x=0.3 " y (m) 0.01 r* 0.008 0.006 1 0.004  ''" 0.002 0 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 xvelocity (rms) Figure 46. Axial velocity profiles of case three. Figure 47 shows the pressure change calculated from the program and Bernoulli's Equation. Little difference is seen between the two solutions and both agree favorably. It is evident that after the corrections are made, the program produces reasonable values for the pressure in case three. The general behavior is also consistent with that of the other two tests. Temperature of the computational space remains constant at five hundred degrees Kelvin. Like case two the low temperature means gas reactions are at a minimum, but the presence of the catalytic surface generates chemical reactions. The reactions produce a slight increase in the temperature just above the surface but the change is minimal. The chemical decomposition of methane illustrated in Figure 48 also seems logical, but the values are small enough to be considered numerical error. The species mass fraction decreases from its initial value at the catalyst and the effect diffuses to the 82 centerline as the flow moves downstream. The decomposition is not sufficient to generate any significant chemical activity. 0.00  0.01 0.02 0.03 S0.04 (Pa) 0.05 0.06 0.07 0.00 0.05 0.10 0.15 0.20 0.25 x(m) Figure 47. Pressure plot of case three. 0.02 0.018 0.016 x=O x=0.1 0.014 x=0 .... x=0.2 0.012 x=0.3 y (m) 0.01  0.008  0.006  0.004 t """ 0.002 0 3.5E11 3E11 2.5E11 2E11 1.5E11 1E11 5E12 reduction of CH4 mass fraction Figure 48. Reduction in methane concentrations of case three. 0 5E12 The third case test reveals that the use of dot products in the subprogram, SpeciesEnergy, leads to resolution problems in the code. After removing the dot products, the software produces good results. However, case three does not produce a considerable amount of chemical activity and a higher temperature is use in case four. Case Four Similar to the third test, case four is another complete test of the software where gas and surface chemistry both exist. The incoming gas temperature and surface temperature is increased to seven hundred degrees Kelvin in an attempt to generate chemical reactions. Three stages are applied in an attempt to resolve the changing gas composition. Parameters and conditions are listed in Table 44. Table 44. Parameters and conditions of case four. Reactor Parameters Initial Flow Conditions Radius 0.02 m Velocity 0.5 m/s Lnocat 0.01 m Temperature 700 K Stage length 0.30 m Pressure 101325 Pa Surface temp 700 K CH4:0.004, CH4:0.004, dy 0.0004 m y 0.0004n Composition 02:0.23, Power 4 power (mass fractions) N2:0.752, No. of Stages 3 AR0.014 AR:0.014 dx (Stagel) 0.001 m dx (Stage2) 0.00001 m PC 1 dx (Stage3) 0.000001 m Input file ptcombust.cti Results of Case Four The program is not able to obtain a solution for the entire flow in case four. Nearly five centimeters into the catalytic reactor, rapid change in the fluid's composition is followed by a large increase in temperature. It appears that an initial temperature of seven hundred degrees Kelvin is sufficient to cause ignition of the air/fuel mixture over the catalytic surface. The software is unable to resolve the rapidly changing flow variables after this point. This is due to the fact that the code being tested is not designed to model a combustion process. Governing equations are reduced based on the assumption of a relatively large characteristic length. Large axial gradients involved with the ignition of the fuel will cause the code to terminate at the point of ignition. The velocity profiles, Figure 49, behave similarly to the other test and do not show any error prior to ignition. Inaccuracy in the axial velocity at the point of combustion is visible just above the surface in the boundary layer. The combustion of the fuel leads to a temperature increase in this same region. The large temperature change causes the density, found in the momentum equation, to change rapidly leading to error in the velocity solution. 0.02 0.018 0.016 x=0 0.014 x=0.025 ....... ignition 0.012 y (m) 0.01 0.008 0.006 " 0.004 . 0.002  . 0 _ . 0.00 0.10 0.20 0.30 0.40 0.50 0.60 xvelocity (m/s) Figure 49. Axial velocity profiles of case four. The change in pressure is graphed in Figure 410. This plot shows that the pressure decreases from its initial value of one atmosphere and appears more linear than the other pressure graphs. Although combustion occurs, significant change in density is not present until ignition and Bernoulli's Equation calculates a pressure difference nearly identical to the pressure change found by the program. Ignition is predicted just before five centimeters into the reactor where the pressure gradient becomes very steep. The pressure's behavior at this point is unexpected and is attributed to the resolution problems associated with combustion. 0 0.001 0.002 0.003 program P0.004 ,'_; Bernoulli (Pa) 0.005 Bern 0 .0 0 5^ ^  0.006 0.007 0.008 0.009 0 0.01 0.02 0.03 0.04 0.05 x (m) Figure 410. Pressure plot of case four. Temperature is graphed at four axial locations in Figure 411. The temperature continues to increase just above the catalytic surface as the flow moves downstream. The exothermic reactions induced by the catalyst lead to the temperature increase in the boundary layer. This variable becomes large and unstable just before the code terminates, which is visible in Figure 411. At the point of ignition, the temperature increases to over eight thousand degrees Kelvin. This value cannot be viewed as an accurate representation of the temperature. However, it appears that the ignition of the fuel is occurring just above the catalytic surface. 0.02 0.018 0.016 0.014 0.012 y (m) 0.01 0.008 0.006 0.004 0.002 0 650 700 750 800 850 900 Temp (K) 950 1000 1050 1100 Figure 411. 0.02 0.018 0.016 0.014 0.012 y (m) 0.01 0.008 0.006 0.004 0.002 0 Temperature profiles of case four. 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 CH4 mass fraction Figure 412. Methane concentrations of case four. The chemical decomposition of methane in case four, before resolution problems arise, behaves much like the mass fraction reduction graphed in Figure 48. The mass x=0 x=0.025   ignition i 87 fraction slowly decreases at the catalyst and the reduction effect diffuses up, away from the surface. Significant methane decomposition occurs right before combustion ends the program. This can be seen in the plot of the methane mass fraction in Figure 412. Many species are produced once significant amounts of methane are broken down. Two such species are atomic and diatomic hydrogen and their mass fractions are graphed in Figure 413. It seems that the greater initial temperature (700K) produces the desired effect of chemical activity. However, the large temperature also produces other species such as OH radicals, and temperature continually grows to the point of ignition. The code is not designed to process combustion and consequently ends at this point. As expected the catalyst aids in the production of hydrogen and the mass fraction of both the hydrogen atom and molecule increase at the surface. The "wiggle" found in the graph of Figure 413 at the last axial position is a result of the absolute value of an overshoot. Cantera cannot process negative mass fractions and any overshoot into the negative must be adjusted. 002 002 0018 0018 0016 0016 0 014 0014 0012 0012 y(m) 001 y(m) 001 0008 0008 0006 0006 0 004 0 004 0002 . 0002 001 0 001 002 003 004 005 006 007 008 002 002 006 01 014 018 022 H mass fraction A H2 mass fraction B Figure 413. Hydrogen concentrations of case four. A) Mass fractions of atomic hydrogen. B) Mass fractions of diatomic hydrogen. The fourth test predicts ignition of the fuel just above the catalytic surface nearly five centimeters into reactor. The greater initial temperature reveals that the program is unable to model catalytic combustion, but can forecast the point of ignition. At this point the software is unable to resolve the rapidly changing flow variables. However, solutions past this point are no longer physically realistic because the assumptions made to simplify the governing equations are not valid. Characteristic length scales in the axial direction become much shorter in the combustion process, which result in very large gradients. The program's resolution problems can be attributed to these large gradients found in some of the variables being determined. The incoming temperature of five hundred degrees Kelvin in case three is too low to produce any significant chemical activity. While the initial temperature of case four is too great and causes combustion. Two additional tests are preformed to better understand the temperature dependence of chemical activity in the reactor. Both tests are similar to case three and four, but use an initial temperature of five hundred fifty and six hundred degrees Kelvin respectfully. The solution of the case using a temperature of five hundred fifty is very similar to the solution of case three. The temperature and composition of the flow remain nearly constant. The other solution, using an initial temperature of six hundred degrees Kelvin, is similar to case four. The temperature continues to increases as the flow moves downstream until ignition is reacted. Comparable to case four, the composition begins to change at this point with the decomposition of methane and the production OH radicals and other species. The point of ignition is further downstream than case four due to the lower initial temperature. It is clear that the chemical activity is highly dependent on the initial temperature. An initial temperature in the range of five hundred fifty to six 89 hundred degrees Kelvin is the temperature needed to cause ignition in the reactor being modeled. CHAPTER 5 PROGRAM LIMITATIONS AND IMPROVEMENTS The program possesses several limiting characteristics when modeling a reacting flow. Calculated solutions are secondorder estimates due to the finite difference equations. Error from these estimates could propagate into the governing equations causing inaccuracies in the calculated solutions. The software uses a mixtureaveraged transport model in order to minimize the time needed to solve the system of equations. The temperature at the catalytic surface is held constant. Initial conditions of a mini mesh and fluid properties embedded in the equations are transformed into smooth cubic spline functions. Errors are undoubtedly produced in this process and rapid changes are not converted to smooth functions very well. Consequently, this program cannot model past the point of combustion and is only physically accurate for a relatively slow reformation process. To improve the program, the minimesh could be enlarged to include more than three axial locations and the use of higherorder finite difference equations would be possible. Increasing the size of the minimesh worsens the effect of the delay discussed in the section Parameters and Conditions and iteration would probably be needed. The code could also be modified to support a multicomponent transport model. Both changes would improve the accuracy of the solution but greatly increase the computation time. The lower boundary condition of the temperature could also be modified to represent a more realistic adiabatic surface or a surface with heat transfer. CHAPTER 6 CONCLUSION A program is created to validate new surface mechanisms of heterogeneous catalysts. The adaptable program models a chemically reacting flow over a catalytic surface. The catalytic reactor is represented in twodimensional Cartesian coordinate form with negligible body forces acting on the fluid. The flow is characterized as a steady, low Mach number, boundary layer flow of a Newtonian fluid. The principles of mass, species mass, momentum, and energy conservation are expressed mathematically and simplified into the governing equations. The model is constructed by numerically solving the system of coupled partial differential equations. The code, which consists of a main program with three subprograms, is written in MATLAB and uses Cantera to calculate chemical properties based on a mixtureaveraged transport model. Allowing Cantera to manage the chemistry independent of the main code allows the program to remain flexible with the varying reaction pathways. Four different cases are utilized to test the program. Calculated solutions from each case are examined to confirm that the software produces reasonable results and is operational. The software is found to predict the point of ignition in the fourth test where the initial temperature is great enough to cause catalytic combustion. Calculated values need to be compared to experimental data to truly determine the accuracy of the program. If the comparison between experimental data and the model reveals error in the program, improvements could be made to the code. Sacrificing computation time for accuracy might be necessary. Once the solutions of the program 