UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 MODELING AND SIMULATION OF A FUEL CELL REFORMER FOR CONTROL APPLICATIONS By MOHUA NATH 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 2007 1 PAGE 2 2007 Mohua Nath 2 PAGE 3 To my parents for their love, sacrifice and steadfast support 3 PAGE 4 ACKNOWLEDGMENTS The completion of this work would not have been possible without the immense contribution of my committee, Dr. William Lear, Dr. Oscar Crisalle and Dr. James Fletcher. They have guided, advised and encouraged me with a lot of patience and supported me in every step during my masters program. I would like to take this opportunity to express my deepest gratitude to them. I would like to thank the most important people in my life, my parents and sister whose faith in me have brought me here. I would also like to mention my friends Alpana Agarwal, Jaya Das, Gaurav Malhotra, Nadeem Islam and Champak Das for being like a family to me and creating a home away from home. A special mention goes to Daniel Betts for sharing his knowledge of fuel cells with me and providing me with valuable suggestions. Last but not the least; I would like to thank Rana Dutta for motivating me to make this last leap possible. You are the light at the end of the tunnel. 4 PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES...........................................................................................................................7 LIST OF FIGURES.........................................................................................................................8 ABSTRACT.....................................................................................................................................9 CHAPTER 1 INTRODUCTION..................................................................................................................11 1.1 Fuel Cells......................................................................................................................11 1.1.1 Phosphoric Acid Fuel Cells...............................................................................12 1.1.2 Polymer Electrolyte Fuel Cells (PEMFC).........................................................13 1.1.3 Alkaline Fuel Cells (AFC)................................................................................14 1.1.4 Molten Carbonate Fuel Cells (MCFC)..............................................................15 1.1.5 Solid Oxide Fuel Cells (SOFC)........................................................................15 1.2 Hydrocarbons as Indirect Fuel......................................................................................16 1.3 Fuel Reforming.............................................................................................................17 1.3.1 Autothermal Reforming (ATR).........................................................................17 1.3.2 Partial Oxidation Reforming.............................................................................17 1.3.3 Steam Reforming..............................................................................................18 2 BACKGROUND AND LITERATURE REVIEW................................................................19 2.1 Reformer.......................................................................................................................19 2.2 Literature review...........................................................................................................24 2.2.1 Langmuir Hinshelwood Model.........................................................................25 2.2.2 Nakagaki Correlation........................................................................................27 3 MODELING OF A PACKED BED REFORMER................................................................28 3.1 Overview.......................................................................................................................28 3.2 Background of Thermal Model.....................................................................................28 3.2.1 Methanol Steam Fuel Cell Reformer................................................................28 3.2.2 General Description of the Reformer................................................................29 3.2.3 Partial Differential Equation.............................................................................30 3.3 Finite Difference Method..............................................................................................30 3.3.1 Discretization....................................................................................................30 3.3.2 Initial and Boundary Conditions.......................................................................35 3.3.2.1 Dirichlets Boundary Conditions........................................................35 3.3.2.2 Neumanns Boundary Conditions.......................................................35 3.3.2.3 Initial conditions.................................................................................36 5 PAGE 6 3.3.2.4 Treatment of an undetermined boundary condition............................36 3.3.2.5 Special Case of Fictitious Nodes at Neumanns Boundary Conditions...........................................................................................37 3.4 Model Validation..........................................................................................................39 3.5 NonDimensional Analysis...........................................................................................42 4 RESULTS AND DISCUSSION.............................................................................................46 5 CONCLUSIONS....................................................................................................................60 APPENDIX A MATLAB CODES.................................................................................................................63 B TEST BED BUS2 CONTROL LOGIC.................................................................................80 C DESCRIPTION OF CONTROL SCHEME...........................................................................91 LIST OF REFERENCES.............................................................................................................108 BIOGRAPHICAL SKETCH.......................................................................................................110 6 PAGE 7 LIST OF TABLES Table page 4.1 Constants used for finding solution to the reformer model. .............................................46 7 PAGE 8 LIST OF FIGURES Figure page 2.1 Generalized Fuel Cell Schematic.......................................................................................23 3.1 Generalized reformer schematic........................................................................................29 32 Discretization along length and radius of the reformer.....................................................32 33 Discretization along length and radius of the reformer and along time increments..........32 3.4 Fictitious nodes..................................................................................................................37 41 Temperature profile of reformer after 230 mins................................................................47 42 Temperature plots at midradius of reformer after 230 minutes........................................48 43 Temperature profile according to analytical solution........................................................49 44 Reformer temperature with zero boundary conditions (Numerical solution)....................50 45 Temperature profile of midradius according to Analytical Solution................................51 46 Transient analytical solution for axial temperature gradient.............................................52 47 Transient numerical solution for axial temperature gradient.............................................52 48 Analytical heat transfer solution for temperature as a function of axial position and time at radial distance = 0.75R...........................................................................................53 49 Numerical solution at different location along xaxis at r = 0.75 m. The temperature profiles are shown at different time instances....................................................................53 410 Non dimensional temperature profile using analytical method......................................54 411 Comparison of temperature profile of midradius according to Numerical Solution with and without heat generation.......................................................................................56 412 Temperature profile of reformer after 230 mins (with heat generation)............................57 413 Comparison of temperature profile of midradius according to numerical solution with and without convection..............................................................................................59 8 PAGE 9 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Master of Science MODELING AND SIMULATION OF A FUEL CELL REFORMER FOR CONTROL APPLICATIONS By Mohua Nath December 2007 Chair: William Lear Cochair: Oscar Crisalle Major: Mechanical Engineering The limited success in the hydrogen storage and distribution technology has driven the need for the development of an effective fuel processor. The dynamic performance of a reformer is of critical importance for the successful commercialization of hydrogen as fuel for stationary and transportation applications. The reforming technology is of particular interest to utilities that require a clean and efficient method of generating electricity from fuel cells. As an effort to achieve a better control of the fuel processor parameters, a dynamic model of a generalized reformer is built. The model successfully predicts the transient temperature gradient across a reformer catalyst bed and the reformate exit temperature and is capable of predicting the response of the reformer to disturbances and load fluctuations. The reaction and heat transfer in the catalyst packed bed was analyzed numerically using a generalized physical model. These results provide valuable insight into the transient response of a reformer. The numerical output was compared with analytical results which agreed well with each other. This confirmed the validity of the numerical method. Also as a part of this research, the control logic of a methanol powered fuel cell bus was studied. The overall control scheme shows that the catalyst bed temperature plays an important 9 PAGE 10 role determining the fuel flow into the reformer. The successful prediction of reformer parameters can thus be utilized to eventually design a reformer capable of quick starts and faster transient response. 10 PAGE 11 CHAPTER 1 INTRODUCTION The growing world economy and the limited resources of nonrenewable fuels emphasize the need for aggressive development of alternative fuels. Advanced power generation technology utilizing alternative fuels can become a factor in reducing emission of greenhouse gases, improving urban air quality and reducing dependency on foreign oil. Hydrogen can be used in fuel cell which is a promising technology, providing efficient and reliable source of energy for a wide range of applications [1,6]. Although hydrogen is abundantly found in nature, extraction of hydrogen from its compounds remains a challenge before it can be commercially viable as a fuel. Hydrogen can be used for multiple applications, ranging from power generation to transportation applications, or internal combustion engines to fuels cells. Reforming hydrogen from any hydrocarbon carrying fuel such as natural gas, biomass, coal or ammonia provides an attractive solution for hydrogen storage for portable applications [1,2]. The particular case of reforming to produce hydrogen for use in fuel cells for transportation applications is selected here for specification of geometric and thermodynamic parameters to be used in this study. 1.1 Fuel Cells Fuel cells are electrochemical devices that produce direct current electrical energy from chemical energy. Fuel and an oxidant are continually supplied to the fuel cell for the reaction to take place [2]. The main components of a fuel cell are an electrolyte, catalyst and a porous anode and cathode. In the presence of a catalyst, the fuel, particularly hydrogen, splits into a proton and an electron. The electrons flow through an external circuit generating electricity and thereafter combines with protons and oxidants to form byproducts at the cathode. Oxygen, acting as an 11 PAGE 12 oxidizing agent in this reaction, combines with the protons and electrons at the cathode to form a molecule of water. The overall reaction is as follows: OHOH22221 11 There are many types of fuels cells currently under investigation, including phosphoric acid fuel cells (PAFC), polymer electrolyte fuel cells (PEMFC), alkaline fuel cells (AFC), molten carbonate fuel cells (MCFC) and solid oxide fuel cells (SOFC). The classification of fuel cells is made based on the type of electrolyte. Among these, the PAFCs and PEMFCs hold the most potential for use as an alternative to internal combustion engines for transportation applications [1]. 1.1.1 Phosphoric Acid Fuel Cells Phosphoric acid fuel cells were the first type to be commercially investigated, other than for the U.S. space program. The electrolyte is phosphoric acid, usually contained in a silicon carbide matrix, and the electrodes made of Teflon bonded platinum or porous polytetradluoroethylene (PTFE) bonded carbon, which is a polymeric binder used to hold the carbon black particles together [1, 7, 5]. In the presence of a catalyst, usually platinum, the positively charged hydrogen ions migrate towards the cathode. Electrons generated at the anode travel through an external circuit towards the cathode, thus creating electric current. The hydrogen ions and electrons combine to form water which is a byproduct of this electrochemical process. The phosphoric acid fuel cells operate optimally between 150C to 220C, as at lower temperatures, carbon monoxide poisoning of the anode may occur. PAFC stacks need a heat sink during operation, usually a coolant, either liquid or gas, that passes through channels integral to the membrane electrode assembly [1, 2]. 12 PAGE 13 The chemical reactions are as below: 2Anode: 22HH e 221Catode: 222OHeH O 221Overall: 2 2 H OH OeO The PAFC stack operates at efficiencies between 37% and 40% and since heat is a byproduct of the electrochemical process, the overall efficiency of a combined process can reach 80%. PAFCs can tolerate a 0.5% concentration of carbon monoxide and are minimally affected by the presence of carbon dioxide. 1.1.2 Polymer Electrolyte Fuel Cells (PEMFC) Proton Exchange Membrane Fuel Cells are considered most suitable for transportation applications. A PEMFC is based on a solid polymer membrane which may be a thin plastic film of sulphonated fluropolymers that act as an electrolyte [1]. The two porous electrodes on either side of the membrane are made hydrophobic by coating with a compound like PTFE, thus helping the reactants to diffuse into a platinum layer that acts as a catalyst. The hydrogen ions diffuse through the membrane towards the cathode and the electrons flow through an external circuit thus producing electric current. Oxygen acts as an oxidizing agent and combines with the electrons and hydrogen ions and forms a byproduct of water. The chemical reactions are as below: 2Anode: 244HH 22Catode: 442OHeH 222Overall: 22 H OH O PEMFCs have higher volumetric energy density than other types of fuel cells, thus making them compact and suitable for vehicular applications. The optimal operating temperature is around 80C which allows quick startup and faster transient response. Recent developments 13 PAGE 14 however have elevated the operating temperature of PEMFCs beyond 150C to reduce the effect of CO poisoning, simplify water and thermal management and to recover high value heat. Other advantages are due to the fact that the electrolyte is solid therefore there is no spillage and corrosion thus contributing to its longer shelf life [5, 7]. One of the major disadvantages of the PEMFC is that the membrane is required to be continually hydrated to operate optimally, thus water management becomes a critical issue. PEMFCs use platinum as a catalyst, which has very low tolerance to CO poisoning. They can operate under a maximum of 10 ppm of CO, thus requiring a clean reformate gas to be used as a fuel. Also, increasing the fuel cell temperature beyond 100C can vaporize the water in the electrolyte which is essential for the conduction of ions, thus requiring tight control of fuel cell temperature and pressure [2, 7]. 1.1.3 Alkaline Fuel Cells (AFC) Alkaline fuel cells use a waterbased electrolyte solution of potassium hydroxide (KOH) with a concentration that can vary according to the fuel cell operating temperature, which could be from 65C to 220C [2, 7]. The hydroxyl ion (OH) acts as the charge carrier as in all other fuel cell types, However, the water formed at the anode travels towards the cathode and regenerates to hydroxyl ions. Therefore, in this type of fuel cell the byproduct is only heat. The chemical reaction is shown below: Anode: 2 H2 + 4 OH4 H2O + 4 eCathode: O2 + 2 H2O + 4 e4 OH222Overall: 22 H OH O The main disadvantage of AFCs is that hydrogen and oxygen have to be supplied to the cell with negligible concentrations of CO2, CO or CH4 as it can poison the electrolyte. This requirement makes it difficult to be used for transportation applications. The advantage is that 14 PAGE 15 they operate on comparatively lower temperatures does not require noble metals and have very high efficiency. 1.1.4 Molten Carbonate Fuel Cells (MCFC) These cells use a mixture of molten salt carbonates as an electrolyte which is contained in a porous, inert matrix made of ceramic. The mixture is usually of varied percentages of lithium carbonate and potassium carbonate. These fuel cells normally operate at a temperature of around 650C. The high operating temperature indicates that these fuel cells can operate directly on gaseous hydrocarbon fuels [1, 2]. At the high temperature of 650C, the alkali melt and become conductive to carbonate ions (CO32) which travel towards the anode. The ions flow from the cathode to the anode where they combine with hydrogen to produce water, carbon dioxide and electrons. Thus the byproducts are carbon dioxide, water and heat. The chemical reactions are shown below: Anode: CO32+ H2 H2O + CO2 + 2 eCathode: CO2 + 1/2O2 + 2eCO32 Overall: H2 + 1/2O2 + CO2 H2O + CO2 The main advantage is that at the higher temperature, fuel reforming can take place inside the fuel cell stack itself, thus eliminating the need of an external reformer. The disadvantage of these type of fuel cells is the time required to obtain the high temperature which greatly slows the startup process and makes the response sluggish [7]. 1.1.5 Solid Oxide Fuel Cells (SOFC) The Solid Oxide Fuel Cell (SOFC) operates under the highest temperature conditions, ranging from 600C to 1000C; thus it can operate with a number of different fuels. The electrolyte is a thin, solid ceramic material which conducts the charge carrying oxygen ions. The 15 PAGE 16 efficiency of the fuel cells are around 60% which is the highest among the fuel cell types. The byproduct of steam can be utilized for various purposes. The chemical reaction is as follows: Anode: 2 H2 + O2 2 H2O + 4 eCathode: O2 + 4 e2 O2Overall: 2 H2 + O2 2 H2O The main application of SOFCs is largescale industrial systems where the demand is higher power, and long startup times only minimally affect the performance or system requirements. Similar to MCFCs, the high temperature of the solid oxide fuel cells make them capable of operating on impure fuels and reforming occurs inside the fuel cell itself. SOFCs are being developed more than MCFCs because they have higher efficiency and ability to operate under higher temperatures [2, 7]. 1.2 Hydrocarbons as Indirect Fuel Hydrogen acts as an energy carrier; however it is not a source of energy. Thus there is a need for converting other sources of energy to hydrogen before it can be used directly in a fuel cell. Electrolysis of water is a well known method of hydrogen production, but it cannot be used for commercial production of hydrogen as it would be more economical to use the electricity directly as an energy source than to use it to produce an energy carrier in the form of hydrogen. The other method commonly used to generate hydrogen is by the reforming of hydrocarbons. Hydrocarbons in liquid form are easier to carry onboard and thus can be used for transportation applications. The most commonlyused process to produce hydrogen is the steam reforming reaction which is to react hydrocarbons with water at a high temperature. To improve the hydrogenproduction efficiency and remove impurities, a watergas shift reaction usually takes place after 16 PAGE 17 the steam reforming reaction. In this process, carbon monoxide produced from steam reforming reaction can be utilized to further break down water into hydrogen. 1.3 Fuel Reforming 1.3.1 Autothermal Reforming Among all the methods that produce hydrogen, autothermal reforming reaction is considered to be one of the most effective processes as it allows faster startup and response time. ATR uses liquid hydrocarbons as fuel that undergoes a reaction with steam or air in a single reactor. While operating in ideal condition, with the optimal amount of air, fuel and steam, the reactions efficiency can reach up to 93.9%. The ATR process is also capable of using hydrocarbons such as gasoline and diesel, which can make it more commercially viable. One of the more recent developments is the possible reduction of operating temperatures from 1,200C to 650900C by reducing oxygen to carbon ratios. The main advantage of reducing the operating temperature is that it allows for simpler reactor design, lowers cost in terms of material complexity and requires less fuel during startup conditions to preheat the reactor [8]. The basic reaction is as follows: 422122CHOHCO 4223CHHOHCO 1.3.2 Partial Oxidation Reforming Partial oxidation reaction is not as frequently used for commercial purposes as the other types. It is an exothermic reaction where a hydrocarbon reacts with controlled amount of oxygen. The advantages of this reaction are that it has a fast response time, high efficiency and does not require a catalyst. In addition, the byproduct of heat can be transferred via a heat exchanger for 17 PAGE 18 other applications. The disadvantages of this process are that it requires a high operating temperature and a high fuel /air ratio for the combustion. The overall reaction is as below: 422122CHOHCO 1.3.3 Steam Reforming The majority of hydrogen produced commercially is from steam reforming reaction. This reaction combines steam with hydrocarbon feedstock in a high temperature and pressure reactor. This is an endothermic reaction, thus requiring an external source of to maintain the temperature in the reactor. Generally, a nickel catalyst bed is used to speed up the reaction and increase the efficiency of the process. The advantages of this reaction are that it can achieve efficiency as high as 85% with heat recovery and can achieve the reaction efficiency of 80% [1]. The main drawbacks of this process are its plant size and slow startup. The overall reaction is as below: 4223CHHOHCO 18 PAGE 19 CHAPTER 2 BACKGROUND AND LITERATURE REVIEW 2.1 Reformer The development of a fuel reformer and its optimization may hasten the advent of widespread fuel cell deployment. The limited success so far, in the establishment of an infrastructure for hydrogen supply, can be overcome by the alternative solution of reforming hydrogen rich fuels onboard. Research on storage methods of hydrogen by physical or chemical means for fuel cell based vehicles has not yet provided a satisfactory practical solution, thus further supporting the need of reforming a hydrocarbon fuel such as methanol. The fuel used for reforming may vary according to the application. Whereas methanol, gasoline or ethanol may be a preferred fuel for transportation applications, natural gas or propane may have advantages for stationary applications. [12]. A reformer system typically consists of a premix tank, preheater, a reactor, a gas pretreatment unit and a burner. A mixture of liquid hydrocarbon and water is fed from the premix tank to the preheater where it is vaporized and superheated before it is sent to the reactor. The energy for the endothermic reaction is supplied by the catalytic burner which heats up the reactor to a preset temperature. The reformed gas is treated in the gas pretreatment unit to remove impurities like carbon monoxide and to further adapt it to meet the fuel cell requirements. The particular type of reforming taken as an example in this study is a catalytic steam process. It utilizes catalytic steam reforming to process the fuel mixture of methanol and water into a hydrogen rich gas. The power generation is accomplished by the oxidation of hydrogen in the fuel cell stack. The depleted fuel mixture is combusted in the steam reformer burner before is the products are finally liberated into the atmosphere. 19 PAGE 20 Several engineering issues hinder the commercialization of fuel cellreformer systems for various applications. One of the most important issues is that the steam reformer dynamic response is considerably slower than that of the fuel cell stack. The response of the reformer is limited by heat transfer rates between the burner gases and the catalyst bed. The response can be defined as a corresponding increase in hydrogen flow rate corresponding to a similar step change in the load at the fuel cell stack. A quick change in heat transfer rates is required to meet the variation of power levels. In addition, the temperature of the catalyst bed should remain constant at design temperature conditions. The increase of catalyst bed temperature from its design conditions can cause the permanent degradation of the catalyst. Thus an optimized reformer design and robust control scheme for the corresponding reformer is to be developed for a quick and efficient dynamic response of the system. The present work is to provide a foundation for such a development. Tighter control of the reformer is required in order for the fuel cell performance to follow the load changing conditions of the vehicle. Control of the reformer has to be achieved by keeping some variables as close to the set point as possible, mostly the catalyst bed temperature, while changing the reformate flow in accordance with the demand at the fuel cell stack which is the loadfollowing generation of electric power. This can be achieved by varying the temperature and flow of exit gases from the burner which is the source of heat for the reformer endothermic reaction. Figure 2.1 shows the general schematic of a fuel cell system with an onboard methanol reformer, where the overall output, in the form of electricity can be stored in batteries or can be used to drive the powertrain of a vehicle. 20 PAGE 21 The example system considered here for the sake of this study is a 50kW, 30ft Fuel Cell Power Test Bed Bus, often referred as TBB, at the Fuel Cell Research Laboratory at the University of Florida [9]. A 175 cell PACF stack is used in the TBB. Air at slightly above ambient pressure is supplied to the cathode to provide oxygen. Hydrogen rich gas, i.e. reformate, is generated by an onboard fuel processor which converts methanol and water to H2, CO2 and CO. During normal operating conditions, the fuel cell stack is provided excess reformate to ensure sufficient hydrogen is available to react at the electrodes. Approximately 80% of the hydrogen in the reformate is consumed and the remaining 20% is supplied to the reformer burner providing heat to continue with the endothermic reforming reaction. For a premix of methanol and water as reactants, the endothermic energy requirement is around 131 kJ/mole at 298C. Methanol is initially delivered by a pump to fuel the startup burner, which in turn provides heat to bring the fuel cell stack and reformer to an initial operating temperature. Once the setpoint temperature is reached in the reformer, the premix fuel is pumped into the reformer where the steam reforming reaction begins to take place. In this type of indirect methanol system, the fuel flow into the reformer is varied according to the load demands at the fuel cell stack output. In the first part of the reaction inside the reformer, the hydrocarbon, in this case methanol, undergoes a cracking reaction to decompose into carbon monoxide and hydrogen. This is an endothermic reaction where energy is continuously absorbed from the surroundings. As the reactants that enter the reformer are already preheated, there is an initial drop in the reactant temperatures, after which the heat required to maintain the reaction is provided solely by the hot gases burned at the catalytic burner. Carbon monoxide produced by the process of steam reforming may poison the noble metal catalysts in the fuel cell stacks. To reduce the concentration of carbon monoxide in the reformate, a water gas shift reaction is used to reduce 21 PAGE 22 the carbon monoxide concentration and additionally increase the hydrogen yield. The overall reaction is as follows: 232HCOOHCH (Cracking reaction) (21) 222HCOOHCO (Shift reaction) (22) This process produces hydrogen rich reformate which is sent directly to the anode of the fuel cell to be oxidized. The unused hydrogencontaining gas, called the flue gas, is returned via a feedback loop to a burner that burns the remaining hydrogen, providing the heat required for the endothermic reformer reaction. In case of increasing power demands from the fuel cell stack, the reformer must flow more reactants into the chamber and produce more hydrogen. This response is usually slow as the constraints are the convective heat transfer between the burner gases and reformer walls and conductive heat transfer from the walls to the catalyst inside the reformer [3]. Accompanied by an increase in fuel flow inside the reformer, and a subsequent demand for heat to the reformer walls, the fuel flow inside the burner is also increased. This helps to maintain the temperature inside the reformer, keeping the rate of reaction constant, as the premix fuel is maintained at uniform concentration. Since these processes have a large lag time, the result may be low quality reformate being directed to the fuel cell during transient conditions. On the other hand, during rampdown conditions of low power demand, the premix fuel flow is reduced to adjust with the low hydrogen requirement at the fuel cell stack. Subsequently, the fuel flow to the burner is also reduced as the hydrogen flue gas returning from the anode is adequate to provide heat to the reactor [9]. Due to the time delay in the reformer due to reaction kinetics, heat diffusion time and convection, the effect of the control action is not measurable for a period of time. Sometimes this 22 PAGE 23 long response time can exceed the quickchanging power demands at the fuel cell stack. This causes the feedback loop to be sluggish where feedback signal is crucial to automatic control of the reformer. The control action is thus inadequate since a change in the fuel flow will deliver the required hydrogen in the stack only after a certain time delay, thus already creating a hydrogen starvation at the stack or depleted hydrogen at the anode flue gas, again causing hydrogen starvation at the burner. Combustion Anode Flue Gas products Anode Reformate Air Air Cathode Reformer Steam Excess burner Reformer Heat exchanger plate Combustion products Vaporizer Heat Exchanger Neat Water methanol methanol Air Startup burner premix air and water Figure 2.1. Generalized Fuel Cell Schematic Modeling is an essential tool to understand the component level interactions of the reformer with the fuel cell stack and its implications on the overall system performance. A reasonable representation of the transient response will enable future development of design and control architectures for the reformer. The reformer model is to predict the concentration of species in the reformate and the reformer catalyst bed. The hydrogen output with respect to time will determine the transient response and time delay of the reformer. The hydrogen flue gas returning to burner also determines the methanol fuel flow rate to the burner. 23 PAGE 24 To meet the needs of future controloriented studies, another objective of modeling the reformer is to obtain a set of differential equations which will represent a state space equation. The state space equation will enable integration of the model into a simulation environment in order to numerically predict the behavior of the system under varying operating conditions. 2.2 Literature review According to Helms and Haley, a quick starting and fast transient response are the most important characteristics of a fuel cell power plant for transportation applications. By far the most popular reforming technology for onboard transportation application is catalytic steam reforming [13]. Geyer et al indicated that the methanol steam reforming technology is superior in its steam reforming technology. However, steam reforming transient response is slow in comparison to other components of a fuel cell system, thus limiting the overall effectivity of a plant in terms of dynamic response [17]. Thus it is crucial to increase the dynamic response of a fuel cell reformer. The first step towards improving the design and response characteristics of a reformer is to develop a model that can be utilized to study efficiency and thermodynamics of the system. The model should also be able to be used in conjunction with a controller design that can help achieve faster response to the system dynamics and disturbances. Kumar et al, Vanderborgh et al., and Geyer et al developed steady state models of the steam reformers, which do not predict the steam reformers transient behavior [16, 18]. Ohl et al developed a first principles dynamic model of a steam reformer for use in parametric design studies. The reformer was assumed to be a well mixed tank reactor. A series of rate reactions for each constituent of the reformate and fuel mixture consisted of the state variables of the system. The heat transfer equation consisted of the change in specific enthalpy due to the reaction kinetics within the reformer. These equations were then represented in a state 24 PAGE 25 equation form to describe a dynamic model of a steam reformer. Model results were then coupled with an optimization process to determine the design parameters of a steam reformer [8, 20, 21]. The study however, does not take into consideration the heat transfer between the walls of the reformer and any external heating agent such as burner hot gases or electric heaters. Convective heat transfer between the gases and the reformer wall temperatures are neglected. Ahmed et al developed a performance model of a reformer that predicted the performance and temperature of a reformate gas mixture. A complete conversion of fuel was assumed as there was a lack of chemical kinetics equation for the partial oxidation reaction for which the model was built. The result of the analysis showed that there existed a linear relation between the exit gas temperature and the inlet temperature of the fuel gas mixture. This model was however specific to only partial oxidation reaction and was not a general reformer model that could be applied for other types of fuel reforming [14]. Choi et al showed the results of kinetics of methanol decomposition along with methanol steam reforming and the water gas shift reaction. A non linear least squares optimization method was used to obtain expressions for rate of reactions [15]. Numerical analysis concluded separate rate of reactions for methanol steam reaction, water gas shift reaction and CO selective oxidation. The three reactors were then integrated and modeled in MATLAB. This study helped in observing the behavior of the reactor by changing its volume and temperature. Chois study however did not take into account the heat transfer in the reactor or the dynamic response of the reformer. 2.2.1 Langmuir Hinshelwood Model Ohl et al developed a dynamic model of the methanol reformer using the Langmuir Hinshelwood (LH) reaction rate for the methanol decomposition reaction. The use of this reaction rate is particularly advantageous because of the wide range of pressures that it covers. 25 PAGE 26 As the pressure increases the adsorption by the surface of catalyst bed and products slow down the reaction. This reaction rate takes all this into account. According to B. A Pepply, at constant temperature the rate of reaction increased with operating pressure although the final conversion decreased with pressure at thermodynamic equilibrium [17]. LH reaction kinetics is one of the simplest reaction mechanisms and it describes most catalytic surface reactions. The LH mechanism assumes that all reactants are adsorbed before the actual reaction can take place. Reactions occur between the adsorbed molecules following a fast diffusion process. The adsorbed molecules then are desorbed. The LH mechanism consists of many reaction steps taken together, ranging from 2 and may extend upto 30. Each step is assumed to be an elementary step, meaning that the reaction is supposed to occur exactly as it is written. Although the general success of the LH mechanism has been accepted, inaccuracies can apparently occur due to certain reactions that may be autocatalytic. Since the water gas shift reaction is itself an autocatalytic reaction, the appropriateness of the use of this rate equation for the water gas shift reaction is debatable. In addition to it, there are speculations that the water gas shift reaction may not occur in the presence of methanol on the catalyst bed. Since the LH mechanism assumes that one step is the rate limiting step while the others are at equilibrium, thus it reduces the rate of the overall reaction defined by the rate of the decomposition of methanol. Ohl expressed the methanol decomposition rate as follows: wwmmmmmcpbpbpbkmr1 (23) 26 PAGE 27 where r is the reaction rate for the decomposition of methanol, mc is the mass of the catalyst bed, is the effectiveness factor, km is the methanol decomposition rate, bm and bw are the methanol and water adsorption equilibrium constants respectively, pm and pw are the partial pressures of methanol and water, respectively. 2.2.2 Nakagaki Correlation Another way to express the rate of reaction has been developed by Nakagaki et al at Toshiba Power Systems. They carried out tests for evaluating the reaction rate of methanol decomposition with varying mass flow rates. Results showed that the diffusion resistance were more significant for lower mass fluxes thus the reaction rate varied at lower mass fluxes [11]. However it was found to be constant at higher mass fluxes, greater than 0.14kg/sm2 Also a correlation was found between pressure and the reaction rate. It was found that the reaction rate became lower at higher pressures, thus they expressed the rate of reaction as power law of the total pressure. Dependence of reaction rate to temperature was found to adhere to Arrheniuss law for temperatures, T<513 K. However for temperatures, T>513 K, conclusive results were not drawn. Based on these experimental observations, Nakagaki et al., derived the reaction rate of decomposition on Methanol, rm on Zn/Cu catalyst to be: (24) NRTEmlmOHCHeTPkr3/0)513/( where, ko = 1.35 X 106 mol / (gcat.s.atm) l = 0.13 E = 1 X 105 J/mol N = 1.3 M = 10 if T > 513 K else = 0 27 PAGE 28 CHAPTER 3 MODELING OF A PACKED BED REFORMER 3.1 Overview This chapter describes the basic principles behind the modeling of physical reformers. The theory governing the equations and the application of initial and boundary conditions are explored. In addition, the text explains the method of finite differences adopted to solve numerically the underlying partial differential equation. An independent analytical solution is presented to solve the heat transfer equation for a special case and thus validate the numerical model. The rate of reaction and the conversion of hydrocarbons to hydrogen depend on the temperature of the catalyst bed of the reformer. Consequently, the catalyst bed temperature and the outlet temperature of the product gases are the most critical process variables affecting the performance of the reformer [4]. The purpose of building the numerical model is to predict the temperature profile and reformate compositions as a function of time and space for a reactor of any cylindrical geometry under different operating conditions. This model should be useful for developing a controller that will manipulate the required process parameters of a methanol fuelcell reformer or a biomass gasifier reactor to optimize their performance. 3.2 Background of Thermal Model 3.2.1 Methanol Steam Fuel Cell Reformer The conversion of a methanol and steam mixture to hydrogen and other product gases takes place inside the reformer. This reformer is usually a cylindrical packedbed reactor with Cu/Zn catalyst pellets [10]. The fuel mixture is preheated in a superheater before being injected into the reformer with a feed pump. The heat for the endothermic reaction is supplied through the 28 PAGE 29 walls of the reformer by a catalytic burner. The rate of reaction or decomposition of methanol depends on the temperature of the reformer and the concentration of methanol. The motivation for developing a thermal model for a methanol fuelcell reformer is to provide a tool for designing a controller which will maintain the reformer at an optimal operational temperature. The thermal model consists of a heat transfer equation with conduction and convection terms. 3.2.2 General Description of the Reformer Due to the similarity between methanol reformers and other hydrocarbon reformers and the common goal of developing a controller, a single model of a general reformer is proposed. A cylindrical packedbed reformer where the reactions take place on the catalyst pellets is assumed [19]. Figure 3.1 shows the mass and energy flows, where r and x are respectively the radial and longitudinal axis of the reformer. The total length of the reformer in the axial direction is L and that in the radial direction is R. External heating is supplied through the walls of the reformer. The flow of gas mixture takes place in the axial direction and the radial mass flow is ignored. heat r x reactants reformate Figure 3.1. Generalized reformer schematic. 29 PAGE 30 3.2.3 Partial Differential Equation A cylindrical coordinate system is chosen in this problem to convert from threedimensions to twodimensions, as the heat conduction and convection is symmetrical about the r axis. The partial differential equation relevant to this model in the cylindrical coordinate system is x TCudVEd x Tk r Tk r T r ktTCggasgeeeee .2222 (32) where T = T(x,r,t) and where e is the effective mass density of the control volume, is the effective specific heat of the control volume, is the temperature inside the reformer, is the effective conductivity of the control volume, eC T ek r is the radial axis, is the axial axis, is the rate of heat generated in the reaction, Vis the total volume of the reformer, uis the velocity of gas, x .gEd gas is the density of the gas, and the void factor is the ratio of mass of gas to the total mass of solid inside the control volume[3]. For simplification of the numerical analysis, the convection term and heat generation term x TCudVEdggasg. gasgTuC x is initially omitted. 3.3 Finite Difference Method 3.3.1 Discretization The first order derivatives of a function f(x+h), in terms of the discrete differences can be expressed using Taylor series expansion as 30 PAGE 31 hffxfkk 1')( (34a) or, hffxfkk1')( (34b) Adding (34a) and (34b) yields a third possible approximation, namely kk f ffxh '()12 1 (34c) Equations (34a), (34b), and (34c) are finite difference approximations using backward, forward and central differentiation, respectively. Using the central difference scheme in the differential equation (32), the first order discretization can be written in the form TT (35a) Similarly (34a) second order discretization is of the form (35b) If the radius of the reformer is divided into M equal discrete elements of size r then represents a node on the radial axis, j r .Similarly if the length of the reformer is discretized into N elements of size x,en i represents a node on axial direction, x. Figure 32 shows the discretization of a crosssection of the upper half of the reformer. Only the upper half of the cylinder is considered as the heat transfer is symmetrical about the raxis of the reformer. Figure 33 introduces another axis with time,t, as an independent variable. Though while simulating, discretization in time is not carried out, we will further proceed to make a dynamic model of the reformer with the help of SIMULINK which utilizes an explicit method to find the transients of th r rTjiji 1,, 21,,2,22 2TT rTrTjijiji 31 PAGE 32 the temperature profile every instant of time, denoted as thn nt where is a finite time interval. t r Figure 32. Discretization along length and radius of the reformer. Figure 33. Discretization along length and radius of the reformer and along time increments. x (0,M) i T j+1 x i,T j i1, j T T i+1,j i,j1 T r (0,j) (0,0) (i,0) (N,0) n n+1 x r Ti j +1 1 j i T Ti+1 j Ti j Ti j +1 32 PAGE 33 Thus the governing PD(32) can be discretized in theand x r E directions into N and M number of nodes and written in the ODE form. For simplicity of numerical analysis, initially the heat generation term has been excluded. 1,1,,1,2 ,ij ,12,,1,2222nnnnnnnnijijijijijijijijeeeeffeffdTTTTTTTTkCkkdtrrrx (36) The above general equation form can be written as a set of ODEs by invoking the value and and the resulting set can be rearranged in terms of a banded pentadiagonal matrix of coefficient, as follows: T Ni...3,2,1,0 Mj...3,2,1,0 NMrlrtrlrtrlrtdtdT0 NMNMNMTTTTlrtrlrtsrlrtsrlrtsrlrtsrsrlrtsrlrsrldtdTdtdtdT1431761000000 (37a) where, = T similarly, lrtsrlrsrlrdTdtdT0 TsrlrdTdt22 T1 TMNNNTNMNNTTTTTTTTTTTTTTT]...........,,,,..,,,,[],.......,,....,,[,2,2,32,22,11,1,31,21,11321 33 PAGE 34 TMNNNTNMNNdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdTdtdT],...,,,..,,,,[],.......,,....,,[,2,2,32,11,1,31,21,11321 The nodes representing i = 0 and j = 0 have not been included in this matrix because i=0 is the boundary region at the entrance of the reformer. Since at this point, the boundary condition is assumed to be constant always, which is equal to the inlet flow gas temperature, it will have a fixed numerical value. Storing N number of numerical values for every computational iteration will take up a lot of computer storage space, thus the nodes are not taken into consideration in this matrix, which was built solely for iteration purpose. All the nodes representing i=0 is represented by one value in the computer program and is used every time the boundary nodes are required to calculate the temperature of its adjacent node. Similarly for j=0 which represents the radius of the cylindrical reformer, it is assumed that there is a no flux condition across the radius of the reformer, these radial nodes will have the same value as its adjacent nodes. In order to save computer space by storing duplicate values, these nodes are not included in the matrix g the iterations. The treatment of these boundary conditions will be dealt with more detail in a later section. Let the N x M pentadiagonal matrix in equation (37) be denoted as P. Then, while doin d dtPT (38) where, T NMNMTTTT,........,121T This equation can be integrated to yield tiiiid nnn 10 0TPTTT (39) 34 PAGE 35 where T0T is a vector of initial conditions TNMNM121Thus a new vector of temperature values as a function of time can be found using various algorithms and iteration methods. SIMULINK has been used to solve the above equations usian explicit RungeKutta method. 3.3.2 Initial and Bound ,........,00000T ng ary Conditions lar case of interest, so as to define appropriate problems with unique solutions. The initial condition gives the specific temperature distribution in the system at time zero, and the boundary conditions specify erature or the heat flow at the boundaries of the medium. 3.3.2.1 Dirichlets Boundary Conditions ons are fixed boundary conditions on the reformer wall and the inlet of the reformer, where the temperature is held constant by electric heaters supplying heat into the system through the walls or by the inflow of preheated gas and. This is expressed atical = at r = R TTTT Boundary conditions and initial conditions are prescribed for a particu the temp The Dirichlets Boundary Conditi mathemly as follows: Txrt(,,) wallT x L0 and = at x = 0, and3.3.2.2 Neumanns Boundary Conditions. Neumann Boundary Conditions are natural boundary conditions in the outlet where there is no heat flux across the radius because of symmetry and the exit of the reformer is assumed to be insulated so that zero heat flux can be assumed in that plane. This is expressed mathematically as follows: t0 Txrt(,,) inT rR0 t0 35 PAGE 36 rTxrt(,,) = 0, at r = 0, x L0 and t0 Txrt(,,) = 0, at x = L, where L is the length of the reformer, x rR 0 and t0 where and rxTxrtTxrtrx TT ghout the reformer. This is expressed mathematically as follows: 0TxrtxLR where is andary condition (,,)(,,) 3.3.2.3 Initial conditions The reformer is heated to a critical temperature by the startup burner. The initial condition is assumed to be uniform throu (,,), a t 0, 0initTtr initT constant value. 3.3.2.4 Treatment of an undetermined bou Let us consider the governing equation (32) again after neglecting the convection term, namely .22 22 g eeeeffeffCkktrrrxdV (310) For the firstor dEkTTTT der term in r the Neumannsun bodary conditions at r = 0 do not apply here because of a term 0/0 that becomes undefined. Applying LHospitals rule yields rrTrrrTr 1 (311) where the right hand side can be written as, 22r TTr (312) rrrReplacing this term in (310) yields, 22222 g eeCk effeffdETTTktrxdV (313) 36 PAGE 37 Furthermore, where r=0, the term g dE dVWhile bu is neglected. ilding a numerical simulation, this special case is accounted for separately to obtain the correct temperatures of the nodes at the axis of the reformer. .3.2.5 Special Case of Fictitious Nodes at Neumanns Boundarormer requires temperatures of a fictitious node, Since there is zero heat flux across the exit nodes of the reformer, we make use of this to assign a value to the fictitious nodes. 3y Conditions While implementing central difference scheme, the temperature of a node depends on the temperatures of its adjoining nodes. Thus the nodes at the exit edge of the ref nNjT,1 Figure 3.4. Fictitious nodes The noflux condition is used to find the value of the fictitious nodes. First the xdirection derivative is written as N jNjT TT x x ,,11 (314) which can be solved for to yield, NjT,1 nNjT,1 nNjT ,1 nNjT,1 nNjT ,1 37 PAGE 38 NjNj TTxT x ,,112 (315) Now at the zeroflux condition, T x 0 so it follows that 1,1,NjNjTT (316) Hence the temperature at fictitious node at ),1(jN can be replaced by the temperature at. A similar relationship can be developed at the other edge having fictitious nodes Finally after implementing all the boundary conditions, the governing pentadiagonal n the form, ),1(jNwhich lie at 0r matrix 36 can be rewritten i MNMNMTlrdt)1(1)1()1)(100 ldTT2102 rlr NTs1 lrrlr ltldTdt72002 MNNNTTTTTtrlrttrlrtrtstsltltsltssslslslsdtdTdTdtdTdtdTdtdtd4321)1((1600200200200000200200020000 (37b) 38 PAGE 39 In the above equation, the fictitious node N+1 and its multiples, which represent the exit end of the reformer has been included in the pentadiagonal matrix. Its value is calculated exactly as the value of N1 node and its multiples. The nodes along the radius of the reformer form the undetermined boundary conditions. Here, at the value of j=0, the equation (313) is implemented. Thus all along the radius, the node temperatures have different coefficients than that of the rest of the discretized reformer. As shown above, in the initial first N nodes of the reformer this boundary condition is imposed. This equation represents the governing matrix that was used to numerically model the generalized refiormer. 3.4 Model Validation To validate the accuracy of the numerical results, a comparison is made by solving the partial differential equation (32) using an analytical method. The physical model is simplified by neglecting the convection term to yield, 222211TTTrrrxt T (317) T= initTT and where initT is the initial temperature. where, rR0, x L0and t0, Equat ion (317) is subject to the boundary conditions, 0T, at 0rRxL (318a) 0T at 0, 0 x rR (318b) 0, at ,0T x Lr R x The problem (317) can be solved analytically by applying boundary conditions (318a) to (318c) to yield TxrtCRrx(,,)()( (318c) xe())22011 (319) mptmpmpmp 39 PAGE 40 where the value of constant ishe principle of orthogonality. Thus mpC determined by using t multiplying both sides by drrrRmR)(00 and dxxxLp0),( yields FdrdxxXrrRNNCpm RL1 rxpmmp)()()()(000 (320) Substituting (320) into (317) yields mp t RLmpeTxrtRrXxrFRrXxdrdx()(,,)()()()()220 (321) Now, it is necessary to find the eigenfunctions mpmpmprxNN()()01100 )(0rRm the norm)(mN and the e igenvalu em From Table 31, Case 3 in reference [ondition in question 10], it follows that for the boundary c )(0rRm )(0rJm = (322a) and )(2)(1202RJRNmm (322b) where m are the positive roots of 0)(0 RJm and J is the Bessel function of order 0. To find the eigenfunctions (xXp 0 ) the norm ) (pN and the eigenvalue p we refer to Table 22, Case 8 in [10] to find that ()sinpPXxx (323a) and LNp2)(1 (323b) 40 PAGE 41 where p are the positive roots of cos(Lp ) = 0. Substituting (322a) (323b) into (321) yields mptRLmpmpmpmprxeTxrtJrsinxrFJrsinxdrdNN x()()()()()()()22001100 (324) Separating the integrals produces the equation (,,) mptRLmpmpmpmprxeTxrtFJrsinxrJrdrsinxdxNN()(,,)()()()()()()22001100 (325) Integrating analytically the rightmost integral, and using the expressions (322a), (322b), yields (323a) and (323b) mptRmpmmpm pr e1TxrtFJrsinxrJrdrLRJR4 ()'(,,)()()()()220022110 (326) 0 mptRmpmmpmpre1TxrtFJrsinxrJrdrLRJR ()'(,,)()()()()22002211004 (327) where from [1], )(1)(10rrJdrrrJmm (328) Substituting (328) into (327), the generalized final analytical solution reduces to, mptmpm mpmmpLRJR()2111 4FeJrsinxJRTxrt()()()()(,,)2201 (329) 41 PAGE 42 3.5 NonDimensional Analysis A nondimensional an alysis is presented here to determine the heat transfer process in the reformer. The formulation of the problem in terms of nondimensional terms helps in studying the physics of the process and should provide a useful tool for data analysis and comparison. Consider the general partial differential equation 222211TTTtrrrx T In order to nond imensionalize along radius and axial direction, let Rrr/ (330a) L x30b) and x (3 0TTwall 0TTT (330c) Subst ituting these new variables into the governing equation yields t TTTxxxTxTTrrrTRrTTrTTTTwallwallwall wall 02 (331) or )()()())((00002 t TT x T L r T R r r T R 22222 022111 (332) aking 2 L tt T equation (332) can be expressed as t TTT xrrRLrTRL 222221 (333) 2 42 PAGE 43 Using the method of separation of variables )(),(),,(txrtxrT (334) Equation (332) can be represented as 22222)()(111tdtdtxrrzrz (335) where, .Then, taking the separated equations 2 )/(RLz 0)()(2ttdtd (336a) and 02222 x r r z r z (336b) xXr Again invoking separation of variables ),(Rxr )()( quation (335) becomes e 01 12zRdz 2222xdXdXrddRrrdR (337a) R and X Separating in 2221 x dX(337b) nd Xd a 222221rvdrdRrzrdRdzR (337c) ads to two equations that can be solved individually to yie leld 0222X x d Xd (,):sincosXxxandx (338a) 43 PAGE 44 0020202R r ddR r z r dRdz (338b) 0)()(2ttdtd (338c) Note: as equation has no dependence, it follows that, v=0 Thus the solutions of (338,ac) taken as, XrRepmt ),(x ),,(,02 e expressed as XrRCtxrT(),(),,(011 can b tpex2), (339) mmpmp Thus the final solution to equation (332) is rRrxXrRNNFetxrTmrxpmmppmt ,(),(),()()(),,(101000112Now to find the eigenfunctions from reference [1], rdxdxXp ),() (340) )()(00rJrRmm (341a) )( 21 )(JRNm 202Rm (341b) here ,m are the positive roots of 0),(0RJm w )(sin)(xxXpp (342a) LNp2)(1 (342b) where, p are the positive roots of 0)(cos Lp wh ere R and L are assumed to be 1. Substituting equations (341a) to (342b) in (340), 44 PAGE 45 rdxdxrJrxrJRJLRFetxr pmrxpmmpmt ),(sin),(),(sin),()(4),,(101000112022ion in nondimensional form, T (343) Solving the remaining integrals gives the final solut ,()(4),,(011212mmppmmtJJFetxrT )(),(sin)1mpJxr (344) 45 PAGE 46 CHAPTER 4 RESULTS AND DISCUSSION This chapter summarizes the computational results obtained by solving the physical models of the reformer numerically and partially validating the results with analytical solution. The steady state solution of the numerical model equation is solved using Matlab. The transient solution is obtained by using SIMULINK, which uses the explicit method to find solution in the next time instant. A major purpose of this study was to be able to develop a model that would correctly predict the reformer behavior in terms of temperature at any particular location. This temperature would then act as an input to a controller that would determine the flow of premix fuel into the system. The validity of the numerical model is established as will be shown in the subsequent sections. For the sake of the study the value of the reformer constants were taken same as the values of a methanol reformer located at UCDavis. The following table shows the values that were estimated to match the UCDavis reformer parameters. Table 4.1 Constants used for finding solution to the reformer model. (Source: Table 4.1, Daniel Betts, 2005) Name Value Used How it was determined Catalyst Bed Heat Capacity, Ce 900 J/kgK Inferred Catalyst Bed Thermal Conductivity, ke 5.00 W/mK Estimated from Data Catalyst Bed Density, e 1983 kg/m3 Measured The catalyst heat capacity was assumed to be equal to an average of the heat capacities of copper and zinc. The void factor used was 50% The catalyst density was measured via a water displacement method for a single pellet 46 PAGE 47 Figure 41 shows the numerical solution of (37) to find the temperature profile for a simple case when the reformer boundaries are kept at constant temperature. The upper boundary or the cylindrical wall of the reformer is kept at 580 K whereas the left side or the face of the reformer is kept at 520 K. Length of reformer [m] 5 10 15 20 25 5 1 0 1 5 2 0 2 5 530 535 540 545 550 555 560 565 570 Radius of reformer [m] Figure 41 Temperature profile of reformer after 230 mins 0.5 Radius of reformer, r [m] 0 Length of reformer, x [m] 0.5 The temperature at different location inside the reformer is shown with color variation. The initial temperature of the reformer was 520 K. As the time increases the temperature inside the reformer gradually increases as can be seen from the plot. The profile shown above is only for one half cross section of the reformer. Now if this transient formulation is run for longer period of time, it gives a steady state solution. Figure 42 shows such time histories for different location along the centerline of the reformer. Figure 42 A shows the temperature plot with respect to time at the inlet of the reformer. It shows constant temperature at the beginning and then it rises and attains a steady state temperature of 567 K whereas for figure 42 B and C the temperature rises from the beginning 47 PAGE 48 (A) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2x 104 523 523.5 524 524.5 525 525.5 526 526.5 Time [s]Temperature [K] (B) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2x 104 520 525 530 535 540 545 550 555 560 Time [s]Temperature [K] (C) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2x 104 520 525 530 535 540 545 550 555 560 Time [s]Temperature [K] Figure 42 Temperature plots at midradius of reformer after 230 minutes at (A) Inlet (B) Center (C) Exit. 48 PAGE 49 Plotting equation (43) in Matlab, we get a temperature profile at r = 0.75, and at different location of distance (x) axially along the reformer. Figure 43 Temperature profile according to analytical solution Tem p erature T [ K ] Time, t [s] Since for the analytical solution we used the Dirichlets boundary conditions as 0K, a numerical simulation with similar boundary conditions was carried out by shifting the wall and inlet gas temperature conditions to 0. The initial reformer temperature is kept at 520 K. As seen from the plot, as the axial length increases the temperature approaches towards the steady state of 0 K at much slower rate. Figure 44 shows the temperature distribution plotted after 230 seconds with zero boundary conditions. Figure 45 shows the temperature profile along the center line of the reformer for the analytical solution. The temperature profile for the inlet shows constant temperature at the beginning before gradually falling and reaching a steady state temperature of 0 K. 49 PAGE 50 R] eformer length [m Reformer Radius [m] 5 10 15 20 25 5 10 15 20 25 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Radius of reformer, r [m] 0.5 Figure 44 Reformer temperature with zero boundary conditions (Numerical solution) 0 Length of reformer, x [m] 0.5 50 PAGE 51 A) 0 5000 10000 15000 0 100 200 300 400 500 600 Time [s]Temperature [K] (B) 0 5000 10000 15000 0 100 200 300 400 500 600 Time [s]Temperature [K] (C) 0 5000 10000 15000 0 100 200 300 400 500 600 Time [s]Temperature [K] Figure 45 Temperature profile of midradius according to Analytical Solution at (A) inlet (B) center (C) exit Figure 46 shows the plot of temperature at r = 0.75 and different locations along the axial center of the reformer with time. This plot represents the analytical solution obtained directly from equation 321 whereas figure 311 shows the same temperature plot at similar location along the axial direction for finite difference solution. The plots match closely with each other when plotted together. Here they are plotted differently for brevity. 51 PAGE 52 Figure 46 Transient analytical solution for axial temperature gradient 0 2 4 6 8 10 12 14x 104 0 100 200 300 400 500 600 t [s]T [k] x=0.1 x=0.25 x=0.50 x=0.75 x=0.90 x increasing Figure 47 Transient numerical solution for axial temperature gradient. Figure 48 shows the temperature plot at different axial location along the reformer when solved analytically. As the axial distance increases, the temperature flattens out. 0 2 4 6 8 10 12 600 14 x 10 4 0 100 200 300 400 500 x = 0.1 x = 0.25 x = 0.5 x = 0.75 Temperature, T [K] Time, t [s] 52 PAGE 53 Figure 48 Analytical heat transfer solution for temperature as a function of axial position and time at radial distance = 0.75R Figure 49 Numerical solution at different location along xaxis at r = 0.75 m. The temperature profiles are shown at different time instances Temperature, T [K] Length of reformer, x [m] Temperature, T [K] Length of reformer, x [m] 53 PAGE 54 Figure 410 shows the results of the nondimensional analysis using the analytical method. The xaxis represents the time 2 L tt the yaxis represents the scaled reformer temperature. Figure 410 Non dimensional temperature profile using analytical method Temperature, T [K] Length of reformer, x [m] So far the results were obtained without taking into consideration the heat generation term and convection term. However, it is important to include the heat of reaction in the current analysis as the reforming reaction is endothermic and apart from the heat transfer between the fluid and reformer wall, energy is consumed to maintain the reaction. The rate at which methanol is consumed is calculated using the Nakagaki rate of reaction as is discussed in detail in Chapter 2. The following results show the comparison between the reformer bed temperatures by addition of the heat generation term. .gEd 54 PAGE 55 A) (B) Figure 411 Comparison of temperature profile of midradius according to Numerical Solution with and without heat generation at (A) exit (B) center (C) inlet 55 PAGE 56 (C) Figure 411. Continued Figure 412 shows the temperature profile of the reformer with the heat generation term. The upper boundary or the cylindrical wall of the reformer is kept at 580 K whereas the left side or the face of the reformer is kept at 520 K. As shown in Figure 411 and Figure 412, it can be seen that due to the endothermic reaction the reformer reaches a steady state at a lower temperature. Thus it can be concluded that the reformer exhibits temperature fluctuations with the change in methanol flow rate. Since the rate of reaction is directly proportional to the number of moles of methanol consumed, any increase of flow of premix fuel will cause a substantial drop in reformer catalyst bed temperature. Thus there is an immediate requirement to ramp up the heat supply into the reformer in terms of burner hot gases. This condition will be especially noticeable in a steep increase in the load demand at the fuel cell stack. 56 PAGE 57 Figure 412 Temperature profile of reformer after 230 mins (with heat generation) Radius of reformer, r [m] 0.5 0 Length of reformer, x [m] 0.5 Betts indicated that the convection term was small compared to the conduction in the heat transfer equation. The authors conclusion was backed by obtaining experimental data from UCDavis reformer which showed that the convection was below 0.4% of the value of the conduction term. These results are expected as the axial temperature gradient is much lesser than the radial temperature gradient and most of the flow occurs in the axial direction [4]. The current analysis includes the convection term and as can be seen from Figure 412, the comparison shows negligible contribution to temperature gradient by the convective heat transfer. 57 PAGE 58 A) (B) Figure 413 Comparison of temperature profile of midradius according to numerical solution with and without convection (A) exit (B) center (C) inlet 58 PAGE 59 (C) Figure 413. Continued 59 PAGE 60 CHAPTER 5 CONCLUSIONS The methodology and techniques used for dynamic modeling of a reformer are presented. The model is further simulated numerically to predict the transient and spatial temperature distribution of the reformer and the composition of the reformate. Steady state steam reformer models have been developed in the past, but most of them are applicable to particular reformer geometry and parameters. The reformer model developed in this research can be used to study the transient behavior of a reformer of any geometry and size. The purpose of this research was to develop a model that could be used in conjunction with a control model to ultimately design a faster responding fuel processing system. The most desired characteristic of a reformer is a quick transient response to load fluctuations. It should be able to ramp up or down the reformate flow according to load demands while maintaining an economic use of primary fuel. The numerical model built in this study helps in understanding the reformer parameters that play an important role in deciding the reaction rate and subsequently the species concentration of reformate. This model in conjunction with a controller model should help build a controller design that will lead to a quick and efficient dynamic response of a reformer. A sensitivity of the various design parameters and their reasonable representation obtained on the basis of the simulation will help in selecting the desired level of control capabilities of various controllers and in evaluating alternative designs of the reformer. From this study, the following conclusions may be drawn: a. At typical boundary conditions and initial conditions of a steam reformer, the model predicts that steady state is reached after approximately 65 minutes. b. Possible location of sensors preferably at the midradius of the reformer. At boundaries, temperatures dont predict accurate reactions rate fluctuations. Aim 60 PAGE 61 is to have lesser number of sensors and less data analysis for the control. Preferably have three sensors at x=0.1L, 0.5L and 0.9L and r=0.75R. Controller input to be average of the three sensors. c. The temperature near the inlet reaches steady state faster than near the exit. However, the steady state temperature reached near the inlet is lower (closer to that of reformate gas temperature boundary condition) than that near the exit. Reactions occurring at the exit are lower as most of the reforming has already taken place near the inlet and heat required for the endothermic reaction is lower at the exit. Thus for higher efficiency and for optimal utilization of premix fuel, it is preferred that the aspect ratio of the reformer closer to 1. d. The effect of preheating the fuel and reformer From the results, it can be concluded that for obtaining faster startup times, it is advisable to preheat both the premix fuel and reformer catalyst bed. e. The convective heat transfer is negligible compared to conductive heat transfer and can be eliminated in further related studies. f. Nonlinear behavior associated with packed bed steam reformers may be noticed during ramp up conditions of the fuel cell stack load conditions due to the endothermic reaction. Since the reformer heat transfer rates are considerably slower than the rest of the system, an immediate increase in premix fuel rate may not produce desired amount of hydrogen in the reformate gas due to a drop in the reformer bed temperature. To avoid this condition of hydrogen starving at the fuel cell stack, designing a feed forward controller is suggested, which will increase the fuel flow into the burner as soon as ramp up condition occurs at the fuel cell stack. Controlling the burner fuel flow with the reformer catalyst bed temperature as the input to the controller will result in sluggish response and incomplete conversion of premix fuel. Further work 1. Input disturbance into the system by making the reformer wall temperature as a Neumann B.C instead of Dirichlets B.C. 2. Gather substantial experimental data to completely verify the numerical model built here. 3. Design a modelbased controller The input variables to the controller can be the reformer catalyst bed temperature, reformate gas temperature and reformate 61 PAGE 62 composition. The manipulated control variables can be the methanol premix flow and the air flow into the burner. 62 PAGE 63 APPENDIX A MATLAB CODES Reformer Input File % reformer_inputfile.m % Mohua Nath 04/23/2006 % Revision 2 % This program specifies all the reformer parameters % This program is called by the mfile temp_cal1 %========================================================== % Catalyst bed Gas mixture effective parameters %========================================================== ro_eff = 1983 ; % effective mass density of catalyst gas mixture in [kg/m3] C_eff = 900 ; % effective specific heat of catalyst gas mixture in [J/kgK] K_eff = 5.00 ; % effective conductivity of catalyst gas mixture [W/mK] %========================================================== % Reformer Dimensions (Cylindrical) %========================================================== Radius = 0.5 ; % Radius of reformer [m] Length = 0.5 ; % Length of reformer [m] %========================================================== %Discretization Parameters %========================================================== M = 30 ; N = 30 ; n = (N2)*(M2) ; del_r = Radius / N ; % Size of finite element along Radius del_x = Length / M ; % Size of finite element along Length %========================================================== % Coefficients of PDE and ODE %========================================================== alpha = K_eff / ( ro_eff C_eff ) ; % coefficient to have Tdot on LHS for i = 1 : M2 p(i) = alpha / ( del_r i* del_r ) ; % Coefficient of del r term of ODE end q = alpha / ( del_r del_r ) ; % Coefficients of del r2 term of ODE 63 PAGE 64 r = alpha / ( del_x del_x ) ; % Coefficients of del x2 term of ODE %=========================================================== % Nakagaki reaction rate constants %=========================================================== k0=1.35e6 ; % reaction rate constant [mol/(g_cat.s.atm)] P=1 ; l=0.13 ; E=1e5 ; % [J/mol] nn=1.3 ; R=8.3143 ; % Universal gas constant %=========================================================== % Ohl's reaction rate constants %=========================================================== pm=10 ; % partial pressure of methanol pw=20 ; % partial pressure of water cm=30 ; % partial pressure of methanol %mc=40 ; % partial pressure of methanol %=========================================================== % Methanol properties %=========================================================== Er=201e3 ; % heat of reaction per mole of methanol consumed %=========================================================== for i = 1 : M2 dia1(i) = (p(i) +2*q+2*r) ; dia2(i) = (p(i) +q) ; end dia11 = 2*(del_r*p(1)+q)2*r ; dia22 = del_r*p(1)+q ; %=========================================================% 64 PAGE 65 Ordinary Differential Equation (Including Heat Generation Term) % matrix_reformer3.m % Mohua Nath 04/23/2006 % Revision 2 % This function executes the Sfunction that solves the N set of ODE's % Invoked by the SIMULINK program matrix_sim_reformer2.sim %=========================================================== % : t= time ; x= vector of stateS ; u= vector of I/Ps ,... % flag= from simulink 