## THERMAL NOx FORMATION IN A TACONITE PELLET ROTARY KILN

## THERMAL NOx FORMATION IN A TACONITE PELLET ROTARY KILN

by RichardA. Davis DepartmentofChemicalEngineering University ofMinnesotaDuluth Duluth, Minnesota 55812-2496

Preparedfor Ron Visness MinnesotaDepartmentofNatural Resources Division ofMinerals Box45, LafayetteRoad St. Paul, Minnesota 55155-4045 March 22, 1995

## ABSTRACT Thereaction mechanisms, thermodynamics, and kinetics ofoxides ofnitrogen 

### (NOx) formation during combustion are summarized. The formation of 

### NOx in the kiln of a taconite pellet induration furnace was modeled as a series of ideal reactors. The combustion zone was treated as a perfectly stirred reactor. The**post-co~bustion** zone was modeled as a plug flow reactor. The temperature in the combustion zone was assumed to bethe adiabatic flame temperature. Complete combustion ofthe fuel was also assumed in order to decouple the 

### NO and combustion reactions. The maximum temperature was calculated to reach 1900 K while the average temperature was found to be 1700 K. The mean residence time in the kiln is estimated to be 3 seconds. The results ofthe model indicate that the primarypathway for the generation of

### NOx is through the mechanism of non-equilibrium thermalNOformation. Themodel also predicts a largeincreaseintheNO formation as the process airflow rate is decreased due to the increased residence time and higher flame temperature. A small decrease in 

### NO formation is predicted for increases in the process air due to the cooling effect on the flame temperature despite the increased availability ofoxygen.



### i

### CONTENTS ABS1R.ACT

................................................................................. i 

### INTR.ODUCTION
..........................................................................1 KILN DESCRIPTION
.....................................................................1 NOxFORMATIONINCOMBUSTIONPROCESSES
...............................1 TIIERMODYNAMICS OFTHERMAL NOxFORMATION
.........................2 MECHANISM AND KINETICS OFTHERMALNOx FORMATION
.............2 IDEALREACTORMODEL
...............................................................4 MODEL RESULTS
.........................................................................? CONCLUSIONS
...........................................................................9 NOMENCLATURE
........................................................................11 REFERENCES
..............................................................................12

List of Tables Table I. 

### USX Minntac plant, Step 3, line 7 kiln
........................................

### 14 Table II. Predicted equilibrium concentrations
..........................................14 Table lll. Predicted residence times to achieve 90% of.
...............................15 Table

### N. Heats offormation (J/gmole) atT = 25 OC
..................................15
List of Figures Figure 1. Taconite pelletinduration process............................................16 Figure 2. Ideal reactors in series
..........................................................16 Figure 3. Qualitative gas and pellettemperature profiles in 

### thekiln
..................17 Figure 4. Dimensionless NO composition and temperature profiles
.................17 Figure 5. The percent change in 

### NO formation vs air

### flow rate......................18 

### APPENDIX Program KILN
..............................................................................19
ii

### 1 INTRODUCTION Oxides of nitrogen (NOx) are common gas phase pollutants formed in high temperature combustion reactions common to industrial devices such furnaces and boilers. The predominant species are nitric oxide (NO) which reacts further with oxygen to form nitrogen dioxide (N02). These compounds contribute to the formation of undesirable ozone, photochemical smog and acid rain (Wark and Warner, 1981). NOx emission rates vary between combustion devices and have been found to depend on a number ofdesign, operating and fuel conditions. Most importantly, NOx formation is a function of the temperature of combustion, process gas residence times above the formation threshold temperature, and turbulence during combustion (Wood, 1994). The focus of this work is NOx formation in the rotary kiln of a taconite pellet induration furnace at the USX Minntac plant, Mountain Iron, Minnesota (Potts, 1991). The essential features of the kiln are presented along with a concise description ofthe thermodynamics and kinetics of NOx formation. Theories of transport, chemical equilibrium, and kinetics ofNOx formation are used in a simple model to investigate the effectsoftemperature andresidence time on NOx formation inthekiln. Empiricalvalues of terminal N02 concentrations, process gas compositions, flow rates, and temperatures are incorporated into the model in ordertoestablishthe principal source ofNOx inthe kiln. 

### KILN DESCRIPTION The major steps of the taconite pellet induration process are depicted in Figure 1. Pellets are flrst dried, sentto a preheater, reach theirflnal--induration temperaturein a rotary kiln and finally cooled withprocess air (Potts, 1991; Josephson, 1992). The pellets move through thekiln countercurrentto the flow ofprocess air and combustion gases. A single burner, fueled primarily with natural gas (CI4), is located atthe pellet discharge point in

. the kiln. Approximately 8% ofthe stoichiometric airrequired for complete combustion is

premixed with the burnerfuel. Hotprocess airfrom the pellet cooler also enters the kiln under turbulent conditions in the vicinity ofthe burner. Pellets are heated in the kiln by directcontactheattransferfrom thekiln walls, thermal radiation from the combustion gases (C02, H20) and forced convection from the process gases. Physical and operating characteristicsofa kiln are summarizedinTable I. Theseparameters areconsidered typical for this studyofNOxformation in a rotary kiln. 

### NOx FORMATION IN COMBUSTION PROCESSES There are three generally accepted pathways to NO formation in combustion processes: thermal, prompt and fuel NO. Thermal and Prompt NO result from high temperatures near the combustion flame forming atomic oxygen and free radicals offuel

### 2 fragments, respectively. These highly reactive species overcome the relatively strong bond in a N2 molecule to form NO. Prompt NO is only formed in the burner flame and is usually considered negligible under conditions ofexcess oxygen where the combustion reaction goes to nearcompletion (Chigier, 1981; Monnotet al, 1985; Wood, 1994). Fuel NO results from oxidation offuel bound nitrogen sources. At cooler temperatures, NO reacts further with oxygento form**N~.** 

### Consideringthekiln burnerprimarily uses natural gas underconditions ofexcess air, fuel andpromptNO formation are assumed negligible. Therefore, this work centers on thermal NOformation as the sole pathway to NOx from the kiln. 

### THERMODYNAMICS OF THERMAL NOx FORMATION The overall reversible reaction ofN2 and 02to form thermal NOis (1) 

### The equilibrium constant for this reaction was correlated from data compiled by Hanson and Salimian (1984).


### K _

### C~

### 0 _ 

### 22 (-21805) No-
- 

### exp 

## eN Co T 2 z (2)


### Table IIcontains calculated equilibrium concentrations ofth.ermal NO in airat a total pressure of one atmosphere for a range of**temperatur~s** 

### typically found throughout the induration process. The results.ofthese calculationsindicate thatequilibrium thermal NO concentrations only begin to become importantfor temperatures above =1200 K (1700 OF). 

### Atprocess gas temperatures exiting the kiln, the thermodynamics ofthermal NO formation predict NO concentrations 

## =800 ppm. Although this is on the order of empirical NO concentrations in gases exiting thekiln, temperatures ofthekiln burnerflame are expected to reach 2000 K (Barnard, 1985) which would result in theoretical equilibrium NO concentrations >8000 ppm. Without considering the kinetics ofthermal NO formation, equilibrium thermodynamics alonecannotexplain the observed NOx formation in a kiln. 

### MECHANISM AND KINETICS OF THERMAL NOx FORMATION The generally accepted mechanism for the reaction of02and N2 to form thermal NO is due to Zeldovich (Bowman, 1973; Zeldovich et al, 1947). At relatively high temperatures common in combustion processes, oxygen dissociates via three body collisions to form atomic 0 (Westenberg, 1971; Warnatz, 1984):

### 3
02

### +M 

### <=> 20

### +M (3) where M is a third body in the airsuch as argon. The bond energy ofN2 is considerably larger than that of02 such that the thennal dissociation ofN2 to form free N atoms is not important (Gupta, 1985). The highly reactive free oxygen atoms breakthe N2 bond to fonn NO and free N, which in

### tum reacts with molecularoxygen. (4) (5)
Thereaction rate constants have been correlated from several experimentsby an Arrhenius equation for temperature dependence. Hanson and Salimian (1984) have surveyed the literature and recommend the following correlations:

( -425) 3.8

### x 107exp 

# T m3 gmol·s k3 

## =1.8

## x 104Texp(-4680) m3 T _gmol·s k4 

## =3.8

### x 103Texp(-20820) 

### m3 T gmol·s 

### (6)
"(7)
(8)
(9)
The large activation energy and temperature dependence in eq 6 indicates that the forward reaction ofeq4 is rate controlling and thatNOwill not approach equilibrium in a combustion process. The overall reaction also depends on the availability of 

### 0 atoms. Hence, the NO formation reaction is expected to become significant only at high temperatures. (Palmer, 1974). Upon cooling, the consumption of NO by reverting to N2 and 02 involves overcoming the large activation energy ofthe reverse reaction ofeq5. As a consequence, essentially all NO forms N02 at lower temperatures according to the following stoichiometric equation (Palmer, 1974;Wark and Warner, 1981): (10)

### 4 Assuming elementary reaction kinetics and application of the pseudo-steady state hypothesis to the intermediate free N atom in eq 

### 5 results in the following reversible rate law for thermal NO formation from 02and N2. gmol 

### m3

### ·s (11)
where 

### Ci is the concentration ofspecies i. 

### Forideal gas behavior, the concentration ofeach species can be 

### represented in terms ofthe partial pressure in the gas phase. C.= 

### YiP

I 

### RT g 

### gmol 

# 7 (12)


### where 

### Yi is the mole fraction ofspecies i, 

### P is the total pressure ofthe gas mixture, 

### Rg is the gas constant and Tisthe absolute temperature. Application of this rate law requires an assumption that the oxygen dissociation reaction approaches equilibrium much faster than the NO formation reaction (Wark and Warner, 1981; Clarke and Williams, 1991). The composition ofatomic oxygen is found from the equilibrium relationship for eq 

### 3 (Warnatz, 1984): 

# C ( 27123) 

### (gmmo 3

# l) 112 

### K

0 

### = 

# S 2 

### =34.64T112 

### exp ---- 

### Coz T (13)


### where three body reactions of02are considered. Eq 11 was used to calculatethermal NO production in air at atmospheric pressure, for several different residence times at various temperatures. The results are listed in Table 

### ill. At temperatures below 1800 K, these calculations reveal that significant thermal NO production requires unreasonably large residence times (Milani and Nelli, 1992; Wood, 1994). The rate law in eq 11 is also used in a model ofthe kiln in order to investigate temperature and residence time effects on thermal NOx formation. IDEAL REACTOR MODEL The rotary kiln is modeled as a series ofideal reactors operating under steady-state conditions (Figure 2). 

### This modeling technique has been employed for several combustion processes to qualitatively predictperformance while avoiding the complexity and difficulty ofthe rigorous field equations for momentum, heat and mass transfer (Essenhigh, 1974; Khalil, 1982; Gupta, 1985). One significant limitation of this type of model is that a thorough investigation ofmixing effects is notpermitted. Thekiln is divided into two ideal

### 5 regions: the burner zone (BZ), comprising the region ofthe combustion flame, and the post-burner zone (PBZ) for the remaining volume ofthe kiln in which the process gas is assumed to achieve fully developed flow. Sudden expansion of the fuel and process air, accompanied by entrainment and recirculation ofthe process gases in turbulentflow around the burnerflame, is anticipated upon introduction of the fuel mixture to the kiln (Gupta, 1985; Carvalho 

### et al, 1990). Completecombustion ofthe fuel is also assumed. Therefore, the BZis assumed to behave as a perfectly stirred reactor (PSR). Under ideal conditions ofperfectly mixed flow, the concentration ofNO in the BZ canbe found from a mass balance assumingno NOx entering the kiln in the process air.



### - rN09Vt YNo- F 

### (14)
where YNO 

### is the exiting mole fraction ofNO from the BZ, 

### Fis the total molarflow rate of all gases through the kiln, Vt is the total volume ofthe kiln, 9 is the fraction ofthe kiln volume represented by the PSR, and r 

### is the rate ofproduction ofNO, from eq 11, in terms ofthe gas phase composition in the BZ. Eq 14 is non-linear in 

### y due to the reaction rate law and must be 

### solved by an iterative trial and errormethod. Thetemperature in the BZis assumed to be 

### the adiabatic flame temperature (AFf) 

### in the presence ofexcess air. Since the actual flame temperature is difficult to calculate in turbulent flames, the AFf is often used under these circumstances with little loss of accuracy (Barnard, 1985; Kunz et al, 1991). The 

### AFT is calculated from an enthalpy balance around the gases entering and exiting the burner zone.

n 

### L,(FioHio- FiHi) 

## =0 (15) i=l 

### where 

### Fi and 

### Hi are the molar flow rate and enthalpy of species i 

### respectively. The subscript 

### o is for entering conditions. The enthalpy ofeach species can be described in terms ofits heat offormation and molarheatcapacity.
T 

### Hi=Hir 

### +

## fcpidT (16) Tr 

### where 

### Cpi is the temperature dependent, constantpressure heatcapacity and

### Hir is the heat offormation atthe reference temperature, 

### T,. Values for 

### cp and

### Htused in the calculations are listed in Table IV. The upper limit ofthe integral is the inlet temperature for the feed6 species and is the AFT, or Tfi for the exiting species. Thisresults in a non-linearequation in Tt. Eqs 14 and 15 are coupled through the composition of 02 and N2 and must be solved simultaneouslyfor Ttandy. The combustion reaction in the presence ofexcess airis relativelyfast compared with NO formation reactions and is assumed to go irreversibly to completion (Monnot eta/, 1985). The overall combustion reactionis (17) The air supply to the kiln, as listed in Table I, is more than 400% in excess of the air required for completecombustion. Also note thatthe overall NOformation reaction in 

### eq 1, and the combustion reaction in 

### eq 17 have equal moles ofreactants and products such thatthe total molarflow rate ofgases in thekiln remains constant. The NO formation reaction is considered to be much slower than the combustion reaction (Chigier, 1981). Combustion reaction times are on the order of 0.0001-0.001 seconds (Essenhigh, 1974), several orders ofmagnitude smaller than the NO formation times calculated for Table Ill. Consequently, NO formation is treated as uncoupled from the combustion reactionin this model. The PBZis modeled as a plug flow reactor (PFR) assuming the process gases reach fully developed turbulent flow in the remaining portion ofthekiln.
. The species balance for NOin this ideal differential reactorresults in a first orderordinary differentialequation. (18) subjectto the initial condition 

### YNo 

# =YNO,PSR at 

# v=

# evt (19) The temperature profile in the PFRis calculatedfrom a differential energy balance. dT 

# -[q +(fllir,NO)rNO] dV

## = FL,y;cpi (20)

with the initial condition (21) where the heat transferperunitvolume from the gases is

### 7 q=UaAT (22) 

### Ua is the overall heat transfer coefficient times the effective heat transfer area per unit volume, ATis the temperatures difference between the hotprocess gasses and the cooler surroundings ofthekiln, &/ r is the heat ofreaction permole ofNOformed (calculated as thedifferencein heats offormation between reaction products and reactants), and 

### c p is the heat capacity of the process gas components. Eqs 18 and 20 are coupled through the compositions in the reaction rate law and temperature. They were solved simultaneously for theexitmole fraction ofNOfrom thekilnbyintegrating numerically overtheremaining volume ofthekiln, 

### 8Vt

### < V

### < Vt. 

### MODEL RESULTS The BZ model equations were solved numerically using the iterative secant method and the PBZ model equations were solved numerically by a 4-5th order Runge-Kutta method with variable step sizesforerrorcontrol. A FORTRANprogram was developed to solve for the NO composition and 

### AFfin the BZjust before integrating the PFR design equations. Thecomputerprogram code, KILN, is found in the Appendix. The

### AFfalong with the two adjustable parameters, 

### 8and qwere found byforcing the modelpredictions to matchexperimentalresults for the NOcomposition and temperature oftheexitgasesfrom thekiln. Theprogram performs these step automatically bysequentially upgrading guesses for 

### 9 and q using a modified form ofNewton's method for two unknowns. As discussed earlier, the pellets move counter current to the process gasses through the kiln. The temperature profile ofthe pellets and gasses may behave in a qualitativefashion similarto the plot in Figure 3 such thatATand q are considered constantin the model (although this isnot required). For the operating conditions of Table I, the flame temperature in the BZ was calculated to be 1900 K. At atmospheric pressure, Barnard (1985) reports a value for the stoichiometric flame temperature of2222 K for Cf4in air. Monnotetal (1985) report the stoichiometric Cf4flame temperature of2600 K with airpreheated to 1500 K. Theflame temperature in the kiln is expected to be less than the stoichiometric flame temperatures cited here due to the additional cooling effect oflarge quantities ofexcess air. Thefraction ofthe totalkiln volume attributed to the BZis estimatedfrom the numerical solution to be 31% while the heat transferrate perunitvolume from the process gas was calculated to be 65,000WJm3. A plotofthe model predictions forthe dimensionless NOconcentration and process gas temperature with dimensionlessposition in the PBZis shown inFigure4. The

### 8


### dimensionless variables for NO concentration, gas temperature, and volume are defmed as follows: y*

# = YNo
- YNo,Bz (23) YNO,BZ
- YNO,PBZ (24)


### V*= 

# v-evt (1- 8)V, (25)


### The subscripts BZ and PBZ refer to the terminal conditions ofeach zone. Recall that the temperature ofthe gas leaving the BZ is the AFT. Inspection of the NO concentration profile in Figure 4 indicates that the majority of NO is expected to form in the BZ and quickly reach the terminal concentration early in the PBZ where the gas temperature drops below the threshold temperature for NO formation. The mean residence time ofthe process gases in the kiln can be calculated from the reactor volume divided by the volumetric flow rate. The mean residence time in the BZ, assuming ideal gas behavioris: 8VtP 1: - _..;.___

BZ- FR T g f (26)


### The meanresidence time in the PBZis determined by integrating the temperature dependent volumetric flow rate overthe volume ofthe PBZ.



### _ [(1-

# e)vkrp 1:PBZ- V.t (27) FR8 jTdV 8V.t 

### Since the temperature profile in the PBZ, as shown in Figure 4, is nearly linear, the integrated average temperature in the PBZ is approximately the algebraic average of the adiabatic flame and PBZ exit temperatures, which yields the result in eq 28 for the mean residence time. (28)

### 9
The total mean residence time ofthe process gasses is the sum ofthe mean residence times inboth zones. V,P[ 

### 8 

### 2(1- 8) ] t'

# =t' BZ

### + t' PBZ

# =-- -- + ( ) FR8 

# r, r, + TPBz 

### (29)

Theresults from the model predictions in eq 

### 29 give an overall mean residence time of3.1 seconds and an average temperature of 

### 1710 K. Confirmation of turbulent gas flow through the kiln requires the dimensionless Reynolds number, shown here as a function ofthe residence time and average temperature ofthe process gases: 4V 

### !R.f,

## =__ k_

### >4

### x 103 

### for turbulentflow tcDVt' (30)


### where D 

### is the inside diameter ofthe rotary kiln and 

### vis the kinematic viscosity of the process gas at the average temperature. Approximating the process gas viscosity to bethat of air at the average temperature, v

### = 3 x I0-4 m2fs (lncropera and DeWitt, 1990), the average Reynolds number equals 2.5 x 105, indicating fully developed turbulent flow of the process gas in the kiln. It

### is instructiveto use the modelresults to investigateprocesschanges which result in a reduction of the NO formation in the kiln. Reductions in temperature, oxygen concentration and residence times in the kiln are expected to dramatically reduce NO formation. Most remedies for cutting 

### NOx formation work by lowering the peak combustion temperature (Wood, 1994). Increasing the flow rate ofprocess air from the cooler should have the simultaneous effect oflowering the burner flame temperature and decreasing the process gas residence time while increasing the available oxygen. Keeping 

### q and 

# econstant, the calculated percent change in the final 

### NO concentration, AFT

### and residence time are plotted in Figure 5 as a function ofthepercentchangein theprocess air flow rate. The most dramatic effect ofchanging the air flow rate is manifest in the sharp rise in the NO formation at lower air flow rates. However, the assumption ofconstant q 

### implies that the pellets would notreach the required temperature ofinduration in the kiln when the 

### air flow rate is increased. 

### CONCLUSIONS According to the literature on 

### NOx formation, fuel and Prompt

### NO formation should not significantly contribute to the effluent 

### NOx concentrations under conditions of combustion ofnatural gas in excess air. A simple model ofthe kiln based on ideal reactors

### 10


### in series was developed to estimate the temperatures and residence time of the process gases in the kiln. The model predictions indicate that the majority ofthe NO is formed in the burnerzone·around the combustion flame by the thermal NO mechanism. Although the exitNO concentrations appearto represent equilibrium values atthe exittemperatures, this is clearly not the case in light ofthe residence time calculations. The residence times to reach equilibrium values are several orders of magnitude larger than the actual mean residence time ofthe process gases. The model predicts a maximum temperature of 1900 K and an average temperature of 1700 K. The estimated process gas residence time in the rotary kiln is 3 seconds. Based on approximate reaction times for thermal NOformation, the concentration ofNO in 

### the kiln 

### effluent stream can only 

### be explainedby the thermal NO formation rate at temperatures above 1800 K. The results of the model give an upper bound on the temperature and the percent volume due to the burner zone. Although it was not specifically investigated, these results can also be considered a lower bound on the effect ofmixing inasmuch as any deviation from perfect mixing will increase the volume and residence times in the burner zone while decreasing the cooling effect ofthe process air.

### NOMENCLATURE



### a 

## = heattransfer area perunit volume, 

### m2 

### Im3 

### BZ 

## = burner zone (PSR) 

### cP 

## = heatcapacity, J

### Igmol· K cj 

## = concentration ofspeciesj, 

### gmol

### Im3 D 

## = kiln diameter, 

### m 

### F 

## = molarflow rate ofgas, 

### gmol

### Is k 

## = reaction rate constant, 

### m3 

### Igmol ·

### s H 

## = enthalpy, J

### Igmol KNO 

## = NOformation equilibrium.constant 

### Ko 

## = 0 formation equilibrium constant, 

### (gmol

### I

### m 3 ) 112 p 

## = total pressure, 

### atm 

### PBZ 

## = postburner zone 

### PFR 

## = plug flow reactor 

### PSR 

## = perfectly stirred reactor 

### q 

## = heat transfer rate perunit volume ofgas, 

### WI m3 

### rNO 

## = rate offormation ofNO, 

### gmol

### Im3 •

### s**~** 

## = Reynolds number 

### Rg 

## = ideal gas constant, 8.205 X 10-5 

### m3 

### ·atmlgmol·K T 

## = temperature, 

### K Tt 

## = adiabatic flame temperature, 

### K 

# u = overall heattransfer coefficient, 

### WIm2 •

### K 

## v 

## = volume variable, 

### m3 

# vt = total kiln volume, 

### m3 yj 

## = mole fraction ofspeciesj 

### Greek Letters 9 

## = fraction oftotal volume attributed to PSR



### r v 

## = = = pi, 3.14159
... residence time, 

### s kinematic viscosity, 

### m2 

### Is 11

### 12 

### REFERENCES Barnard, J.A.(1985) Flame and Combustion 2ndEd., Chapman andHall, London, 52-3. Bowman, C.T. (1973) "Kinetics ofNitric Oxide Formation in Combustion Processes," 14th Symposium on Combustion, TheCombustion Institute, 729. Carvalho, M.G., Semiao, V., Lockwood, F.C. and C. Papadopoulos (1990) "Predictions of Nitric Oxide Emissions from an Industrial Glass-Melting Furnace," J. 

### Inst. Energy, 63, 39-47 Chigier, N.A. (1981) Energy, Combustion and Environment, McGraw-Hill, New York, 328-334. Clarke, A.G. and A. Williams (1991) "The Formation and Control ofNOx Emissions," Chemistry &

### Industry, Dec. 16, 917-920. Essenhigh, R.H. (1974) "An Introduction to Stirred ReactorTheory Applied to Design of Combustion Chambers," Combustion Technology: Some Modern Developments, Palmer, H.B. and J.M. Beer, Eds, Academic Press, New York, 376-385. Gupta, A.K.(1985) FlowfieldModelling andDiagnostics, Abacus Press
.. Kent, Ohio, 129- 336. Hanson, R.K. and S. Salimian (1984) "Survey of Rate Constants in H/N/0 

### Systems," Combustion Chemistry, Gardiner, W.C. Ed., 361. Himmelblau, D.H. (1981) Principles and Calculations in Chemical Engineering, 4th Ed. Prentice-Hall, Inglewood Cliffs, NJ. Incropera, F.D. and D.P. DeWitt (1990) Fundamentals ofHeat and Mass Transfer, Wiley, New York ,p.Al5. Josephson, G.R. (1992) Air Pollution Engineering Manual, Buonicore, A.J. and W.T. Davis, Ed.s, Van Nostrand Reinhold, New York, pp. 679-681. Khalil, E.E. (1982) Modelling ofFurnaces and Combustors, Abacus Press, Kent, Ohio, 250-252. Kunz, R.G., Smith, D.D., Patel, N.M., Thompson, G.P., and G.S. Patrick (1992) "Control NOx from Furnaces," Hydroc. Proc., 71, 57-62
. Milani, A. and D. Nelli (1992) "Low-NOxCombustion Techniques Applied to Steelworks Plants Firing Gas and Oil," J. 

### Inst. Energy, 65, 35-40. Monnot, G., Feugier, A., Le Bouc, F., Mauss, F., Meyer, C., Perthuis, E. and C. De Soete (1985) Principles ofTurbulent Fired Heat, Gulf Publishing Co., Houston, 135-224. Palmer, H.B (1974) "Equilibrium and Chemical Kinetics in Flames," Combustion Technology: Some Modern Developments, Palmer, H.B. and J.M. Beer, Eds, Academic Press, New York, 19-20.

### 13 Potts, R.D. (1991) "Induration of Taconite Pellets Using a Grate-Kiln**System~"** 

### AIME/University ofMinnesota Mining Symposium, Duluth, MN, Jan. 8-9. Siddiqi, A.A. and J.W. Tenini (1981) 

### "NOx Controls in Review," Hydroc. Proc., 60, 

### 115-124. Warnatz, J (1984) "Rate Coefficients in the C/H/0 System," Combustion Chemistry, Gardiner, W.C. Ed., 336. Wark, K.; Warner, C.F. (1981) AirPollution: Its Origin and Control, 2nd Ed., Harper and Row, New York, 376-377. Westenberg, A.A. (1971) Combustion Science and Technology, Vol. 4, 59. Wood, S.C. (1994) 

### "Select the Right NOx Control Technology," Chern. Eng. Prog., 90, 

### 32-38. Zeldovich, Y.B., Sadovnikov, P.Y., and D.A. Frank-Kamenetskii (1947) Oxidation of Nitrogen in Combustion (Translations by M. Shelet) Academy ofSciences ofUSSR (Institute ofChemical Physics) Moscow.

### Table

### I. USX Minntac plant: Step 3, line7 kiln operatingConditions (Potts, 1991) Inside Volume 1100m3 

### Inside Diameter 5.94 m 

### Length 39.6 m Mass flow rate ofprocess inlet air 

### 2400 gmolls Mass flow rate ofairto burner 47 gmol/s Mass flow rate ofCI4toburner 61 gmolls Temperature ofinlet air

### (fromcooler) 1370 K (2100 OF) 

### Temperature ofexiting air 

### 1370 K (2100 OF) 

### Average N0:2 concentration in exiting gas =500 ppm

Table

## n. Predicted equilibrium concentrations ofNO in air atP == 

### 1atm. T(K) T(OF) 

### NO (ppm)

600 620 0.02 800 980 2.3 1000 1340 36 1200 1700 220 1400 2060 810 1600 2420 2110 1800 2780 4600 2000 3140 8300 2200 3500 14000 2400 3860 21000 14

### 15

Table

### Ill. 

### Predicted residence times to achieve 90% 

### of equilibrium NO concentrations in air 

### atP

# =1atm. T(K) TCF) t 

### (sec)



### 1000 1340 4 X 

### 1012 1500 2240 6 X 

### 105 1600 2420 85000 1700 2600 14000 1800 2780 3300 1900 2960 840 2000 3140 250 2100 3320 79 2200 3500 25 2300 3680 12 2400 3860 4.7 2500 4040 

### 2.3

Table IV. 

### Heats offormation (J/gmole) at T

### =25 tJC 

### and heatcapacity data (Himmelblau, 1981): c, 

### =a+bT

## +cTz 

## +dT 3 ( 

### 11• ) 

### gmo· C 

### Species 

## H,,, X 

### 10=4 

### a b X 

### 102 cX106 dx 109 (>2 

### 0 29.1 1.16 w6,08 1.31 N2 

### 0 29.0 0.22 

### 5.72 w2,87 

## NO 9.0370 29.5 0.8188 w2,925 0.3652 

### CH4 -7.4840 34.31 5.47 3.66 w 11 

### C02 -39.351 36.11 4.23 -28.9 7.46 

### H20 -24.183 33.46 0.688 7.6 -3.59Pellets In Dryer



### Gas Out Preheater Kiln

AirOut 

### Figure 1. 

### Taconite pellet induration process.



## PFR PSR Pellets Figure 

### 2. Ideal reactors in series. 16

Pellets Air

Out In



### Air

## In

### •**~** 

### • > 

### TT0

(1-0V) 0V 

# y._j Figure 3. Qualitative gas and pellettemperature profiles in the kiln.

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0 0 0.2 0.4 0.6 0.8 V* 

### Figure 4. Dimensionless NO composition and temperature profiles in the PBZ. 17

### 18



### 450 o/oA NO 0 

### /o AAFT

- - ·- - o/o A res. time Cl) 

### 300 E += (/)**~** 

# t 150 <( 

### 0 

## z <:]**~** 

# - -· - 

### 0

- 

## - - - -



### -150 -50 -25 

### 0 25 

### 50 o/o A Air Flow Figure 5. The percent change in 

### NO formation 

### vs air flow rate.

### 19



### APPENDIX Program KILN ************************************************************************ * * By: Richard Davis Language: FORTRAN 77 DATE: 6/5/94 * * This is program KILN.F by R. A. Davis. This program solves for the * exiting temperature and NO concentration in the process gas exiting * the kiln. The kiln is modeled as a PSR and PFR in series. The PSR * is used to model the •burner zone•. The temperature is assumed to be * the adiabatic flame temperature of the burner fuel. The combustion is * assumed to go to completion (no prompt NO reaction)
. The fuel * combustion reaction is fast relative to the NO formation reaction. * Thus, it is assumed that the NO reaction is decoupled from the * combustion reaction (occurs after the combustion)
. The PFR is used to *model the "post burner zone". The process gas and pellets pass * counter-currently, hence, the rate of heat loss from the process * gas is assumed to be constant (q [=] W/m3). The two unknowns q and * V(PSR) are determined by matching the model calculations with the * process conditions for the exiting temperature and NO concentration * using Newton's method with first order derivative approximations. * The model assumes ideal gas behavior due to the low P and high T in * the process. * The PFR design equations are solved by a 4-5th order variable step * size Runge-Kutta numerical method for error control. * ********************************** Nomenclature ************************ * * Species identification for array values * CH4 1 * C02 = 2 * H20 

### = 3 * 02 4 * N2 

## = 5 * NO 

## = 6 * a = flag for adiabatic vs non adiabatic PSR * b() = boundary conditions; 1 = NOx concentration. 2 = Temperature * c(,) 

### = heat capacity, kJ/kgmol-K * CO, C02, CN2, CNO = concentrations of 0, 02, N2 and NO in moles/m3 * del = delta g factor * delHl, delH2 = Ho

- Hi for Secant method in the adiabatic flame * calculation * dg() = increment in g for derivative approximations * du() =derivative of residuals with respect to guesses * dwm,dwt = temporary derivative evaluations for RK4 calculations * eps = error limit for variable step size calculation * F() =molar feed rate (inlet conditions), kgmol/s * F2() =molar exit rate from the PSR or feed to the pfr * g() = initial guess for the fitting paramters; 1 = q, 2 = V(PSR)/Vt * g1 = temporary value of g(1) * h = step size in the ·PFR solution * hh = half step for RK4 * h6 

### = h/6 for RK4 * Hf() 

### = heat of formation at reference temperature, kJ/kgmol * k() =reaction rate constants; 1=k(l), 2=k(2), 3=k(-1), 4=k(-2) * m = number of Newton iterations* n = number of derivative iterations * NO = flag to signal reaction in PSR * nstep = integration step counter 

### * P absolute pressure, atm 

### * PSRl, PSR2 = temporary variable in the secant iteration method for 

### * yNO calculation from PSR. * r = rate of reaction, kgmol/m3-s * Rg = ideal gas constant, 8.205e-5 m3-atm/mol-K * sumcl,sumc2 = summation of heat capacity * T = temperature in adiabatic flame iteration calculations, K 

### * TF adiabatic flame temperature, K 

### * Tfl, Tf2 = iteration variables for secant method in adiabatic flame 

### * temperature calculation * Tfc = flame temperature in degree C * Ti = inlet temperature to PSR from coolers, K * Trc = reference temperature in degree C * TK 

## = 273.15 * tol = convergence tolerance * Tr 

## = reference temperature, 298.2 K * u() = residuals of calculated boundary conditions * V = dimensionless volume of PFR * Vt = PFR volume, m3 * w() 

## = dimensionless yNO and T * wt() 

### = temporaty w value for RK4 calculation * y() mole fraction of species i * yNO =mole fraction of NO in the process gas. * 

### 20

********************************** Driver ****************************** ***

***

***
***
***
***
***

*** This is the main driver program for subroutines FLAME, implicit double precision (a-h,o-z) dimension b(2),c(6,4),dg(2),F(6),F2(6),g(2),Hf(6), Ti (6) I u ( 3 I 2) I w (2) I y (6) Initialize data call init(a,b,c,del,eps,F,F2,Ft,g,Hf,hrnin,hl, m,n,NO,P,Rg,Tf,Ti,TK,tol,Tr,Vt,y,yNO) Get adiabatic flame temperature call flame(c,F,F2,Hf,Tfi,Ti,TK,tol,Tr) Newton iterations do 2000 i 

### = l,m PSR and PFR.

Increment in guesses for residual derivative approximations do 1000 j 

### = l,n Increment the guesses for the partial derivative approximations. call pdg(del,dg,g,j) Model "burner zone" as a PSR, assume T = adiabatic flame temperature call PSR(a,c,F,Ft,F2,g,Hf,NO,P, Rg,Tf,Ti,TK,tol,Tr,Vt,y,yNO) Model "post burner zone" as PFR. Init. condtns=exit condns from PSR. call pfr(b,c,eps,Ft,g,hrnin,hl,Hf,j, P,Rg,Tf,TK,Tr,Vt,w,y,yNO)*** Get the residuals in the calculated solution call resid(j,u,w) 1000 ***
*** continue Get the upgraded values for g1 and g2 call upgrd(dg,g,u) Print out the upgraded results for the guessed values. call output(g,i,Tf,Tfi,u,Vt,yNO) rt=Vt*P*(g(2)/Tf+2.0*(1-g(2))/(Tf+b(2)))/Ft/Rg Tave=Vt*P*(g(2)+(1-g(2)))/Ft/Rg/rt write(*,*) write(*,*) •Residence time =•,rt, •sec• write(*,*) •Average T =•,Tave, •K• *** Check for convergence in Newton iterations. call check(tol,u) 2000 continue write(*,*) •No Newton convergence achieved. • stop end 

### 21

********************************** pdg *********************************
***

***

1000 ***

*** subroutine pdg(del,dg,g,j) Returns the guesses plus increments in g. implicit double precision (a-h,o-z) dimension dg(2),g(2) Get the delta g values if (j.eq.1) then do 1000 k = 1,2 dg(k) = del*g(k) + 1.0d-12 continue Get g1 + delg1 else if (j.eq.2) then g(1) = g(1) + dg(1) Restore g1 and get g2 + delg2 else if (j.eq.3) then g ( 1 ) = g (1 ) - dg ( 1 ) g(2) 

## = g(2) + dg(2) end if return end ********************************** upgrd *******************************
*** * subroutine upgrd(dg,g,u) Calculates residuals from over the boundary conditions at end of calculations. implicit double precision (a-h,o-z) dimension dg(2),du(2,2),g(2),u(3,2) *** Restore g2 to initial value at start of Newton iteration. g ( 2 ) 

### = g ( 2 ) - dg ( 2 )*** Get the derivatives of residuals with respect to guesses.

*** du(1,1) (u(2,1)

- u(1,1))1dg(1) du(1 12) (u(3 1 1)

- u(1,1))1dg(2) du ( 2 I 1 ) ( u ( 2 I 2 )

- u (11 2 ) ) Idg ( 1 ) du(2 12) (u(3 12)

- u(1 12))1dg(2) Solve Newton equations for upgraded values of g(). g1 (u(1 1 1)

- g(1) *du(1,1) + (g(2) *du(2 12) + g(1) *du(2 11)

-**u(1****2))*du(1~2)1du(212)** - g(2)*du(1 12))1(du(2,1)*du(1,2)1du(2,2) du(1,1)) 

### 22

g (2) 

### = ( g (2) *du (2, 2) + g (1) *du( 2 1 1)

- 

## u(11 2)

- g1*du (2, 1) ) Idu ( 2 1 2) g(1) g1 *** Keep the guesses within the physiscal bounds if (g(1) .lt.O.O) g(1) = 1.0d-12 if (g(2) .lt.O.O) g(2) = 1.0d-12 return end ********************************** check *******************************
*** subroutine check(tol,u) Check for convergence in g() 1 i.e., u < tol implicit double precision (a-h 1o-z) dimension u(3 12) u2 

## = (u(1,1)**2 + u(1,2)**2)**0.5 if (u2.lt.tol) stop return end ********************************** output ****************************** subroutine output(g 1i,Tf1Tfi,u,Vt,yNO) *** Print out the intermediate results implicit double precision (a-h,o-z) dimension g(2),u(3,2) write(*,*) write(*,1000) i 1000 format(' Iteration# ',I2) write(*,2000) Tfi,Tf,yN0*1.0d6 2000 format(' Adiab. Flame T 

## = ',F7.2,' K; PSR T 

## = ',F7.2,' K; y(PSR) = ' 1F12.2,' ppm') write(*,*) 'PFR Residuals:' 1U(1,1) 1 '; ' 1U(1 12) write ( * 1 * ) I g ( 1 ) = I 1 g ( 1 ) 1 1 ; g ( 2 ) = I 1 g ( 2 ) write ( * I 3000) g ( 1 ) , g ( 2 ) *Vt I ( 1 . 0 - g ( 2 ) ) *Vt 3000 format (' q 

## = ',F9.2,' Wlm3; V(PSR) 

## = ',F7.2,' m3; V(PFR) 

## = ' , F7
. 2 , ' m3 ' ) write(*,*) return end ********************************** Flame *******************************
*** subroutine flame(c,F,F2,Hf,Tf,Ti,TK,tol,Tr) This solves for the adiabatic flame Tf by the secant method. implicit double precision (a-h,o-z) dimension c(6,4),F(6),F2(6),Hf(6),Ti(6)*** Initialize variables for secant iterations Tf1 = 2000.0 Tf2 = 0.99*Tf1 delH2 = tol Trc 

### = Tr

- TK *** Secant loop for adiabatic flame temperature do 2000 i 

### = 1,100 delH1 

### = 0.0 Tfc = Tf1

- TK do 1000 j=1,5 Tic 

## = Ti(j)

- TK 

### 23

delH1 = delH1 + (F2 (j)

- F(j)) *Hf (j) + c(j,1)*(F2(j)*(Tf1- Tr)

- F(j)*(Ti(j)

- Tr))+ c(j,2)*(F2(j)*(Tfc**2 Trc**2) F(j)*(Tic**2

- Trc**2))/2.0 +

1000 continue c(j,3)*(F2(j)*(Tfc**3

- Trc**3)

- F(j)*(Tic**3

- Trc**3})/3.0 + c(j,4)*(F2(j)*(Tfc**4

- Trc**4)

- F(j)*(Tic**4

- Trc**4)}/4.0 Tf = Tf1

- delH1*(Tfl

- Tf2)/(de1Hl

- delH2) *** Check for convergence if(dabs((Tf-Tf1)/Tf1) .lt.tol) then return else· *** Upgrade previous solutions Tf2 

### = Tfl Tf1 = Tf delH2 

### = delH1 end if 2000 continue write(*,*) 'No convergence in adiabatic flame temperature secant calculation' stop end ********************************** ·PSR ********************************
***

*** subroutine PSR(a,c,F,Ft,F2,g,Hf,NO,P,Rg,Tf,Ti,TK,tol,Tr,Vt,y,yNO) Calculates T and yNOx from the PSR (Burner Zone) implicit double precision (a-h,o-z) dimension c(6,4),F(6),F2(6),g(2),Hf(6),Ti(6),y(6) Initialize temporary variables for secant method Tfl Tf Tf2 = 0.99*Tf1 Trc 

### = delH2 yNOl yN02 PSR2 Tr

- TK 

### = tol

yNO 0.99*yNOl 0.0 *** Secant iterations do 2000 i 

### = 1,100*** Initialize the NO temperature and mole fraction in the PSR Ti(6) = Tfl y(6) = yN01 24

*** Upgrade the molar flow rates of 02, N2 & NO F2(4) = (y(4)

- 0.5*y(6))*Ft

***

1000

*** F2(5) = (y(S)

- 0.5*y(6) )*Ft F2 (6) y (6) *Ft Calculate the adiabatic flame temperature delH1 = 0.0 Tfc = Tf1

- TK do 1000 j=1,6 Tic = Ti(j)

- TK delH1 = delH1 + (F2(j)

- F(j))*Hf(j) + c(j,1)*(F2(j)*(Tf1- Tr)

- F(j)*(Ti(j)

- Tr))+ c(j,2)*(F2(j)*(Tfc**2 Trc**2) F(j)*(Tic**2

- Trc**2))/2.0 + c(j,3)*(F2(j)*(Tfc**3 Trc**3) F(j)*(Tic**3

- Trc**3))/3.0 + c(j,4)*(F2(j)*(Tfc**4 Trc**4) F(j)*(Tic**4

- Trc**4))/4.0 continue delH1 = a*g(1)*g(2)*Vt + delH1 Tf = Tf1

- delH1*(Tf1

- Tf2)/(delH1

- delH2) Calculate the PSR yNO if (NO.eq.1) then call rxn(P,r,Rg,Tf,y) PSR1 = yN01*Ft

- r*g(2)*Vt yNO 

## = yN01

- PSR1*(yN01

- yN02)/(PSR1 end if - PSR2)
*** Check for convergence if (NO.eq.1) then check= (((yNO- yN01)/yN01)**2 + ((TfTf1)/Tf1)**2)**0.5 else if (NO.eq.O) then check= dabs((Tf- Tf1)/Tf1) end if if (check.lt.tol) then y(6) = yNO return else *** Upgrade the temporary variables.

end if 2000 continue Tf2 = Tf1 Tf1 = Tf delH2 = delH1 yN02 yN01 yN01 yNO PSR2 PSR1 write(*,*) 'No convergence in yNO and Tf calculation.' stop end

### 25

********************************** rxn ********************************* subroutine rxn(P,r,Rg,T,y) *** Returns the reaction rate at the temperature supplied implicit double precision (a-h,k,o-z) dimension k(4),y(6) *** Set the rate constants k(1) (1.8d8)*dexp(-38370.0/T) k(2) = (1.8d4)*T*dexp(-4680.0/T) k(3) (3.8d7)*dexp(-425.0/T) k(4) = (3.8d3)*T*dexp(-20820.0/T) !m3/mol-s !m3/mol-s !m3/mol-s !m3/mol-s *** Set the concentrations CNO = y(6)*P/Rg/T !mol/m3 C02 = y(4)*P/Rg/T- 0.5*CNO !mol/m3 CN2 = y(5)*P/Rg/T- 0.5*CNO !mol/m3 CO= 36.64*({T*C02)**0.5)*dexp(-27123.0/T) !mol/m3 *** Calculate the reaction rate of NO production r 2.0*CO*(k(1)*k(2)*C02*CN2-k(3)*k(4)*(CN0**2.0))/{k(2)*C02 + k (3) *CNO) !mol/m3-s return end ********************************** resid *******************************
***

1000 subroutine resid(j,u,w) Get the residuals in the calculated boundary conditions. implicit double precision (a-h,o-z) dimension u(3,2),w(2) do 1000 L = 1,2 u(j,L) = 1.0- w(L) continue return end ********************************** pfr *********************************
***

***

1000

*** subroutine pfr(b,c,eps,Ft,g,hmin,h1,Hf,j,P,Rg,Tf,TK,Tr,Vt,w,y,yNO) This subroutine calculates the exiting yNO and final T in the kiln. implicit double precision (a-h,o-z) dimension b(2),c(6,4),d(2),g(2),Hf(6),w(2),wscal(2),y(6) Set initial conditions for dimensionless V, yNO and T v 

### = 0.0 h = h1 do 1000 i = 1,2 w(i) = 0.0 continue nstep = 0 call pfrout(j,nstep,V,w) This loop is adapted from "Numerical Recipes," pp. 559 do 3000 nstep = 1,10000 call deriv(b,c,d,Ft,g,Hf,P,Rg,Tf,TK,Tr,Vt,w,y,yNO)2000 do 2000 i 

### = 1,2 wsca1(i) 

### = dabs(w(i}) + dabs(h*d(i)) + 1.0d-30 continue if ((V + h- 1.0)*(V +h) .gt.O.O) h 

## = 1.0- V call rkqc(b,c,d,eps,Ft,g,hmin,h,hdid,hnext,Hf, P,Rg,Tf,TK,Tr,V,Vt,w,wscal,y,yNO) call pfrout(j,nstep,V,w) if ((V- 1.0) .ge.O.O) then return end if if (dabs(hnext) .lt.hmin) then
else end if write(*,*) •stepsize smaller than minimum. • h hmin h hnext 3000 continue write(*,*} 'Too many steps. • return end 26

***************************** pfrout *********************************** subroutine pfrout(j,nstep,V,w) *** Print out intermediate results implicit double precision (a-h,o-z) dimension w(2) if (j.eq.1) then if (nstep.eq.O) then write(*,1000) 1000 format ( I v I I t 15 1 I W ( 1 ) I 1 t 30 I I W (2 ) I ) write(*,2000) 2000 f 0 rmat ( I
--------- - I 1 t 15 1 I
---- -- --- - I 1 t30 1 I---------- I) write(*,3000) V,w(1) ,w(2) 3000 format(t1,F10.6,t15,F10.6,t30,F10.6) else end if end if return end write(*,3000) V,w(1), w(2)

***************************** rkqc *************************************

*** subroutine rkqc(b,c,d,eps,Ft,g,hmin,htry,hdid,hnext,Hf, P,Rg,Tf,TK,Tr,V,Vt,w,wscal,y,yNO) This subroutine calculates 5th order RK Step for error control implicit double precision (a-h,o-z) dimension b(2),c(6,4),d(2),dsav(2),g(2),Hf(6), w(2),wsav(2),wscal(2),wtemp(2),y(6) pgrow 

### = -0.2 pshrnk 

### = -0.25 fcor 

### = 1.0/15.0 safety 0.9 errcon 

### = 6.0d-4 Vsav = V do 1000 i 

### = 1,2wsav(i) dsav(i) 1000 continue w( i) d (i) h 

## = htry 1 hh 

## = 0.5*h call rk4(b,c,dsav,Ft,g,hh,Hf,P,Rg,Tf,TK,Tr,Vt,wsav,wtemp,y,yNO) V 

### = Vsav + hh call deriv(b,c,d,Ft,g,Hf,P,Rg,Tf,TK,Tr,Vt,wtemp,y,yNO) call rk4(b,c,d,Ft,g,hh,Hf,P,Rg,Tf,TK,Tr,Vt,wtemp,w,y,yNO) V 

### = Vsav + h if (V.eq.Vsav) then write(*,*) 'Stepsize not significant in rkqc.' h hmin V = Vsav + h end if call rk4(b,c,dsav,Ft,g,h,Hf,P,Rg,Tf,TK,Tr,Vt,wsav,wtemp,y,yNO) errmax 

## = 0.0 do 2000 i 

## = 1,2 wtemp(i) 

## = w(i)

- wtemp(i) errmax 

### = dmax1(errmax,dabs(wtemp(i}/wscal(i))) 2000 continue errmax 

## = errmax/eps if (errmax.gt.1.0) then
else

end ifh 

### = safety*h*(errmax**pshrnk) go to 1 hdid 

### = h if (errmax.gt.errcon) then hnext safety*h*(errmax**pgrow) else hnext 4.0*h end if do 3000 i 

## = 1,2 w(i) 

### = w(i) + wtemp(i)*fcor 3000 continue return end 27

********************************** rk4 *********************************
***

1000

2000 subroutine rk4(b,c,d,Ft,g,h,Hf,P,Rg,Tf,TK,Tr,Vt,w,wo,y,yNO) This subrouine calculates one 4th order Runge-Kutta step implicit double precision (a-h,o-z) dimension b(2),c(6,4),d(2),dwm(2),dwt(2},g(2),Hf(6), w ( 2 ) ,wo ( 2 ) ,wt ( 2 ) , y ( 6 ) hh 

### = h*0.5 h6 = h/6.0 do 1000 i 

### = 1,2 wt(i) 

### = w(i) + hh*d(i) continue call deriv(b,c,dwt,Ft,g,Hf,P,Rg,Tf,TK,Tr,Vt,wt,y,yNO) do 2000 i 

### = 1,2 wt(i) = w(i) 

### + hh*dwt(i) continue call deriv(b,c,dwm,Ft,g,Hf,P,Rg,Tf,TK,Tr,Vt,wt,y,yNO) do 3000 i = 1,2 wt(i) 

### = w(i) + h*dwm(i)dw.m(i) = dwt(i) + dwm(i) 3000 continue call deriv(b,c,dwt,Ft,g,Hf,P,Rg,Tf,TK,Tr,Vt,wt,y,yNO) do 4000 i=1,2 wo(i) = w(i) + h6*(d(i) + dwt(i) + 2.0*dwm(i)) 4000 continue return end 28*********~**************************** deriv ******************************
***

***

***

1000 subroutine deriv(b,c,d,Ft,g,Hf,P,Rg,Tf,TK,Tr,Vt,w,y,yNO) This subroutine calculates values of the derivatives for T and yNO. implicit double precision (a-h,o-z) dimension b(2),c(6,4),d(2),g(2),Hf(6),w(2),y(6) Calculate yNO and T from dimensionless variables y(6) = (b(1)

- yNO)*w(1) + yNO T = Tf- (Tf- b(2))*w(2) Tc = T - TK Trc = Tr

- TK Get dy*/dV* call rxn(P,r,Rg,T,y) d(1) = Vt*(1.0- g(2))*r/(b(1)

- yNO)/Ft sumc1 = 0.0 do 2000 i = 1,6 sumc2 = 0.0 do 1000 j 

## = 1,4 sumc2 sumc2 + c(i,j)*(Tc**(j-1)) continue sumc1 = sumc1 + y(i)*sumc2 2000 continue *** Get dT*/dV* delct (c(6,1)- 0.5*(c(4,1) + c(5,1)))*(T- Tr) + (c(6,2)

- 0.5*(c(4,2) + c(5,2))) *(Tc**2

- Trc**2}/2.0 + (c(6,3)

- 0.5*(c(4,3) + c(5,3)))*(Tc**3- Trc**3)/3.0 + (c(6,4)

- 0.5*(c(4,4) + c(5,4)))*(Tc**4- Trc**4)/4.0 d(2) (g(1) + (Hf(6) + delct)*r)*Vt*(1.0- g(2))/(Tfreturn end b(2))/Ft/sumc1

********************************** init ********************************

***
*** subroutine init(a,b,c,del,eps,F,F2,Ft,g,Hf,hrnin,h1, m,n,NO,P,Rg,Tf,Ti,TK,tol,Tr,Vt,y,yNO) implicit double precision (a-h,o-z) dimension b(2),c(6,4),F(6),F2(6),g(2),Hf(6),Ti(6),y(6) Adiabatic PSR (a = 0.0 for adiabatic; a = 1.0 for non-adiabatic) a 

### = 0.0 Boundary conditions b(1) SOO.Od-6 b(2) = 1367.0***

***
***
*** Heat Capacity & Heat of Formation Data from Himelblau, c(1,1) 34.31 !Jimol-C c(1,2) 5.47d-2 c(1,3) = 3.66d-6 c(1,4) -1.1d-8 c (2, 1) = 36.11 !Jimol-C c(2,2) 4.23d-2 c(2,3) -2.89d-5 c(2,4) 7.46d-9 c(3,1) = 33.46 !Jimol-C c (3' 2) 6.88d-3 c(3,3) 7.6d-6 c(3,4) -3.59d-9 c(4,1) = 29.1 !Jimol-C

. c(4,2) = 1.16d-2

c(4,3) -6.08d-6 c(4,4) = 1.31d-9 c (5, 1) 29.0 !Jimol-C c(5,2) 2.2d-3 c(5,3) = 5.72d-6 c(5,4) = -2.87d-9 c (6, 1) = 29.5 !Jimol-C c(6,2) 0.8188d-2 c(6,3) = -0.2925d-5 c(6,4) 0.3652d-9 Increment in the guesses del 

## = 1.0d-4 Error control parameter eps 1.0d-6 Initial molar feed rates to the "burner zone" write(*,*) '%of the actual air feed rate:' read(*,*) pair pair=pairl100 F(1) 61.0 F(2) 0. 0 F(3) 

### = 0.0 F(4) 

## = pair*503.0+0.16*F(1) F(5) pair*1890.0+0.602*F(1) F(6) 0.0 !molls !molls !molls !molls !molls !molls 

### 29

App. E & F

*** Calculate molar feed rates to the "post burner zone" as exit values *** from the "burner zone" F2(1) 

## = 0.0 F2 (2) F.( 2) + F (1) F2 ( 3 ) F ( 3 ) + 2 . 0*F (1) F2 (4) F (4) - 2 . 0*F (1) F2(5) F(5) F2 (6) 

## = 0. 0 !molls !molls !molls !molls !molls !molls *** Calculate the total molar feed rate to the PFR Ft 

### = 0.0 do 1000 i=1,5 Ft = Ft + F2(i) 1000 continue***

***

***
***
***
*** Set the initial guesses g (1) :: 6.5d4 g(2) :: 0.31 Heats for formation Hf(l) -7.484d4 Hf(2) :: -3.935ld5 Hf(3) :: -2.41826d5 Hf(4) 0.0 Hf(S) :: 0.0 Hf(6) :: 90.374d3 Initial step size h1 == 1.0d-3 !W/m3 !1/100
!J/mol !J/mol !J/mol !J/mol %

!J/mol
!J/mol

Minimum step size in RK calculations hrnin 1.Od-4 Maximum number of Newton iterations m == 100 Number of derivative approximation iterations n 

## = 3 *** NO 

## = 1 for reaction, NO 

## = 0 for no reaction in PSR NO 1 *** Set m n 1 for sensitivity analysis m 1 n m *** Pressure in kiln
***
***
***

***
***
*** P 

### = 1.0 !atm Ideal gas constant Rg 

## = 8.205d-5 !m3-atm/gmol-K Initial guess for flame temperature. Tf 

## = 1900.0 Initial temperatures Ti(l) = 1367.0 Ti(2) = 1367.0 Ti(3) = 1367.0 Ti(4) 1367.0 Ti(S) 

### = 1367.0 Ti(6) 1367.0 Degrees Kelvin TK:: 273.15 of all

Newton convergence tolerance tal 

## = 1.0d-8 Reference temperature for Hf Tr 

### = 298.2 species 

### 30***
***

***
2000Total volume of the kiln Vt = 1 . 1d3 !m3 Initial guess for yNO if (NO.eq.O) then yNO 

## = 0.0 else yNO end if 380d-6 

### 31

Calculate mole fractions of species entering the
..post burner zone.. do 2000 j 

### = 1,6 y (j) = F2 (j) /Ft continue return end ************************************************************************