THERMAL NOx FORMATION IN A 

TACONITE PELLET ROTARY KILN 
THERMAL NOx FORMATION IN A 

TACONITE PELLET ROTARY KILN 

by 

Richard A. Davis 

Department of Chemical Engineering 

University of Minnesota Duluth 

Duluth, Minnesota 55812-2496 

Prepared for 

Ron Visness 

Minnesota Department of Natural Resources 

Division of Minerals 

Box 45, Lafayette Road 

St. Paul, Minnesota 55155-4045 

March 22, 1995 
ABSTRACT 

The reaction mechanisms, thermodynamics, and kinetics of oxides of nitrogen (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 

be the adiabatic flame temperature. Complete combustion of the 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 of the model 

indicate that the primary pathway for the generation of NOx is through the mechanism of 

non-equilibrium thermal NO formation. The model also predicts a large increase in the NO 

formation as the process air flow 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 of oxygen. 

i 
CONTENTS 

ABS1R.ACT ................................................................................. i 

INTR.ODUCTION .......................................................................... 1 

KILN DESCRIPTION ..................................................................... 1 

NOx FORMATION IN COMBUSTION PROCESSES ............................... 1 

TIIERMODYNAMICS OF THERMAL NOx FORMATION ......................... 2 

MECHANISM AND KINETICS OF THERMAL NOx FORMATION ............. 2 

IDEAL REACTOR MODEL ............................................................... 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 of formation (J/gmole) at T = 25 OC .................................. 15 

List of Figures 

Figure 1. Taconite pellet induration process ............................................ 16 

Figure 2. Ideal reactors in series .......................................................... 16 

Figure 3. Qualitative gas and pellet temperature profiles in the kiln .................. 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 of design, 

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 of the 

thermodynamics and kinetics of NOx formation. Theories of transport, chemical 

equilibrium, and kinetics of NOx formation are used in a simple model to investigate the 

effects of temperature and residence time on NOx formation in the kiln. Empirical values of 

terminal N02 concentrations, process gas compositions, flow rates, and temperatures are 

incorporated into the model in order to establish the principal source of NOx in the kiln. 

KILN DESCRIPTION 

The major steps of the taconite pellet induration process are depicted in Figure 1. 

Pellets are flrst dried, sent to a preheater, reach their flnal--induration temperature in a rotary 

kiln and finally cooled with process air (Potts, 1991; Josephson, 1992). The pellets move 

through the kiln counter current to the flow of process air and combustion gases. A single 

burner, fueled primarily with natural gas (CI4), is located at the pellet discharge point in 

. the kiln. Approximately 8% of the stoichiometric air required for complete combustion is 

premixed with the burner fuel. Hot process air from the pellet cooler also enters the kiln 

under turbulent conditions in the vicinity of the burner. Pellets are heated in the kiln by 

direct contact heat transfer from the kiln walls, thermal radiation from the combustion gases 

(C02, H20) and forced convection from the process gases. Physical and operating 

characteristics of a kiln are summarized in Table I. These parameters are considered typical 

for this study of NOx formation 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 of fuel 
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 of excess oxygen where the combustion 

reaction goes to near completion (Chigier, 1981; Monnot et al, 1985; Wood, 1994). Fuel 

NO results from oxidation of fuel bound nitrogen sources. At cooler temperatures, NO 

reacts further with oxygen to form N~. Considering the kiln burner primarily uses natural 

gas under conditions of excess air, fuel and prompt NO formation are assumed negligible. 

Therefore, this work centers on thermal NO formation as the sole pathway to NOx from the 

kiln. 

THERMODYNAMICS OF THERMAL NOx FORMATION 

The overall reversible reaction of N2 and 02 to form thermal NO is 

(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 II contains calculated equilibrium concentrations of th.ermal NO in air at a total 

pressure of one atmosphere for a range of temperatur~s typically found throughout the 

induration process. The results. of these calculations indicate that equilibrium thermal NO 

concentrations only begin to become important for temperatures above =1200 K (1700 OF). 

At process gas temperatures exiting the kiln, the thermodynamics of thermal NO formation 

predict NO concentrations = 800 ppm. Although this is on the order of empirical NO 

concentrations in gases exiting the kiln, temperatures of the kiln burner flame are expected 

to reach 2000 K (Barnard, 1985) which would result in theoretical equilibrium NO 

concentrations >8000 ppm. Without considering the kinetics of thermal NO formation, 

equilibrium thermodynamics alone cannot explain the observed NOx formation in a kiln. 

MECHANISM AND KINETICS OF THERMAL NOx FORMATION 

The generally accepted mechanism for the reaction of 02 and 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 air such as argon. The bond energy of N2 is considerably 

larger than that of 02 such that the thennal dissociation of N2 to form free N atoms is not 

important (Gupta, 1985). The highly reactive free oxygen atoms break the N2 bond to fonn 

NO and free N, which in tum reacts with molecular oxygen. 

(4) 

(5) 

The reaction rate constants have been correlated from several experiments by an Arrhenius 

equation for temperature dependence. Hanson and Salimian (1984) have surveyed the 

literature and recommend the following correlations: 

(

-425) 3.8 x 107 

exp 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 of eq 4 is rate controlling and that NO will 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 of the reverse reaction of eq 5. 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 02 and N2. 

gmol 

m

3 

·s 

(11) 

where Ci is the concentration of species i. For ideal gas behavior, the concentration of each 

species can be represented in terms of the partial pressure in the gas phase. 

C.= YiP I RT 

g 

gmol 

7 (12) 

where Yi is the mole fraction of species i, P is the total pressure of the gas mixture, Rg is 

the gas constant and Tis the absolute temperature. 

Application of this rate law requires an assumption that the oxygen dissociation 

reaction approaches equilibrium much faster than the NO formation reaction (W ark and 

Warner, 1981; Clarke and Williams, 1991). The composition of atomic oxygen is found 

from the equilibrium relationship for eq 3 (Warnatz, 1984): 

C ( 27123) (gmmo3l)112 

K0 = S2 =34.64T112 exp ---- Coz T 

(13) 

where three body reactions of 02 are considered. 

Eq 11 was used to calculate thermal 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 of the 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 of ideal reactors operating under steady-state 

conditions (Figure 2). This modeling technique has been employed for several combustion 

processes to qualitatively predict performance while avoiding the complexity and difficulty 

of the rigorous field equations for momentum, heat and mass transfer (Essenhigh, 197 4; 

Khalil, 1982; Gupta, 1985). One significant limitation of this type of model is that a 

thorough investigation of mixing effects is not permitted. The kiln is divided into two ideal 
5 

regions: the burner zone (BZ), comprising the region of the combustion flame, and the 

post-burner zone (PBZ) for the remaining volume of the 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 of the process gases in turbulent flow around the burner flame, is anticipated 

upon introduction of the fuel mixture to the kiln (Gupta, 1985; Carvalho et al, 1990). 

Complete combustion of the fuel is also assumed. Therefore, the BZ is assumed to behave 

as a perfectly stirred reactor (PSR). 

Under ideal conditions of perfectly mixed flow, the concentration of NO in the BZ 

can be found from a mass balance assuming no NOx entering the kiln in the process air. 

- rN09Vt 

YNo- F (14) 

where YNO is the exiting mole fraction of NO from the BZ, F is the total molar flow rate of 

all gases through the kiln, Vt is the total volume of the kiln, 9 is the fraction of the kiln 

volume represented by the PSR, and r is the rate of production of NO, from eq 11, in terms 

of the 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 error method. 

The temperature in the BZ is assumed to be the adiabatic flame temperature (AFf) in 

the presence of excess 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 of each species can be described in 

terms of its heat of formation and molar heat capacity. 

T 

Hi= Hir + f cpidT (16) 

Tr 

where Cpi is the temperature dependent, constant pressure heat capacity and Hir is the heat 

of formation at the reference temperature, T,. Values for cp and Htused in the calculations 

are listed in Table IV. The upper limit of the integral is the inlet temperature for the feed 
6 

species and is the AFT, or Tfi for the exiting species. This results in a non-linear equation 

in Tt. Eqs 14 and 15 are coupled through the composition of 02 and N2 and must be 

solved simultaneously for Ttand y. 

The combustion reaction in the presence of excess air is relatively fast compared with 

NO formation reactions and is assumed to go irreversibly to completion (Monnot eta/, 

1985). The overall combustion reaction is 

(17) 

The air supply to the kiln, as listed in Table I, is more than 400% in excess of the air 

required for complete combustion. Also note that the overall NO formation reaction in eq 

1, and the combustion reaction in eq 17 have equal moles of reactants and products such 

that the total molar flow rate of gases in the kiln 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 of magnitude smaller than the NO formation 

times calculated for Table Ill. Consequently, NO formation is treated as uncoupled from 

the combustion reaction in this model. 

The PBZ is modeled as a plug flow reactor (PFR) assuming the process gases reach 

fully developed turbulent flow in the remaining portion of the kiln. . The species balance 

for NO in this ideal differential reactor results in a first order ordinary differential equation. 

(18) 

subject to the initial condition 

YNo = YNO,PSR at v = evt (19) 

The temperature profile in the PFR is calculated from a differential energy balance. 

dT -[q + (fllir,NO )rNO] 

dV = FL,y;cpi (20) 

with the initial condition 

(21) 

where the heat transfer per unit volume from the gases is 
7 

q=UaAT (22) 

Ua is the overall heat transfer coefficient times the effective heat transfer area per unit 

volume, AT is the temperatures difference between the hot process gasses and the cooler 

surroundings of the kiln, &/ r is the heat of reaction per mole of NO formed (calculated as 

the difference in heats of formation 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 the exit mole fraction of NO from the kiln by integrating numerically over the remaining 

volume of the kiln, 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 sizes for error control. A FORTRAN program was developed to 

solve for the NO composition and AFf in the BZ just before integrating the PFR design 

equations. The computer program code, KILN, is found in the Appendix. The AFf along 

with the two adjustable parameters, 8 and q were found by forcing the model predictions to 

match experimental results for the NO composition and temperature of the exit gases from 

the kiln. The program performs these step automatically by sequentially upgrading guesses 

for 9 and q using a modified form of Newton's method for two unknowns. As discussed 

earlier, the pellets move counter current to the process gasses through the kiln. The 

temperature profile of the pellets and gasses may behave in a qualitative fashion similar to 

the plot in Figure 3 such that AT and q are considered constant in the model (although this 

is not 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 of 2222 K for Cf4 in air. Monnot et al (1985) report the 

stoichiometric Cf4 flame temperature of 2600 K with air preheated to 1500 K. The flame 

temperature in the kiln is expected to be less than the stoichiometric flame temperatures 

cited here due to the additional cooling effect of large quantities of excess air. The fraction 

of the total kiln volume attributed to the BZ is estimated from the numerical solution to be 

31% while the heat transfer rate per unit volume from the process gas was calculated to be 

65,000 W Jm3. A plot of the model predictions for the dimensionless NO concentration and 

process gas temperature with dimensionless position in the PBZ is shown in Figure 4. The 
8 

dimensionless variables for NO concentration, gas temperature, and volume are defmed as 

follows: 

y* = YNo - YNo,Bz (23) 

Y NO,BZ - Y NO,PBZ 

(24) 

V*= v-evt 

(1- 8)V, (25) 

The subscripts BZ and PBZ refer to the terminal conditions of each zone. Recall that the 

temperature of the 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 of the 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 behavior is: 

8VtP 1: - _..;.___ 

BZ- FR T 

g f 

(26) 

The mean residence time in the PBZ is determined by integrating the temperature dependent 

volumetric flow rate over the volume of the PBZ. 

_ [(1- e)vkr p 

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 of the process gasses is the sum of the mean residence times 

in both zones. 

V,P [ 8 2(1- 8) ] t' = t' BZ + t' PBZ = -- --+ ( ) FR8 r, r, + TPBz 

(29) 

The results from the model predictions in eq 29 give an overall mean residence time of 3.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 of the residence time and average temperature 

of the process gases: 

4V !R.f, = __ k_ > 4 x 103 for turbulent flow 

tcDVt' 

(30) 

where D is the inside diameter of the rotary kiln and v is the kinematic viscosity of the 

process gas at the average temperature. Approximating the process gas viscosity to be that 

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 instructive to use the model results to investigate process changes 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 of process air from the 

cooler should have the simultaneous effect of lowering the burner flame temperature and 

decreasing the process gas residence time while increasing the available oxygen. Keeping 

q and e constant, the calculated percent change in the final NO concentration, AFT and 

residence time are plotted in Figure 5 as a function of the percent change in the process air 

flow rate. The most dramatic effect of changing the air flow rate is manifest in the sharp 

rise in the NO formation at lower air flow rates. However, the assumption of constant q 

implies that the pellets would not reach the required temperature of induration 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 of natural gas in excess air. A simple model of the 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 of the NO is formed in 

the burner zone ·around the combustion flame by the thermal NO mechanism. Although the 

exit NO concentrations appear to represent equilibrium values at the exit temperatures, this 

is clearly not the case in light of the residence time calculations. The residence times to 

reach equilibrium values are several orders of magnitude larger than the actual mean 

residence time of the 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 NO formation, 

the concentration of NO in the kiln effluent stream can only be explained by 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 of mixing inasmuch as any deviation from perfect mixing will increase the volume 

and residence times in the burner zone while decreasing the cooling effect of the process 

air. 
NOMENCLATURE 

a = heat transfer area per unit volume, m2 I m3 

BZ = burner zone (PSR) 

cP = heat capacity, J I gmol· K 

cj = concentration of species j, gmol I m3 

D = kiln diameter, m 

F = molar flow rate of gas, gmol I s 

k = reaction rate constant, m3 I gmol · s 

H = enthalpy, J I gmol 

KNO = NO formation equilibrium .constant 

Ko = 0 formation equilibrium constant, (gmol I m3

)

112 

p = total pressure, atm 

PBZ = post burner zone 

PFR = plug flow reactor 

PSR = perfectly stirred reactor 

q = heat transfer rate per unit volume of gas, WI m

3 

rNO = rate of formation of NO, gmol I m3 • s 

~ = Reynolds number 

Rg = ideal gas constant, 8. 205 X 10-5 m3 ·atml gmol·K 

T = temperature, K 

Tt = adiabatic flame temperature, K 

u = overall heat transfer coefficient, W I m2 • K 

v = volume variable, m3 

vt = total kiln volume, m3 

yj = mole fraction of species j 

Greek Letters 

9 = fraction of total 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 2nd Ed., Chapman and Hall, London, 52-3. 

Bowman, C.T. (1973) "Kinetics of Nitric Oxide Formation in Combustion Processes," 

14th Symposium on Combustion, The Combustion 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 of NOx Emissions," 

Chemistry & Industry, Dec. 16, 917-920. 

Essenhigh, R.H. (1974) "An Introduction to Stirred Reactor Theory 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) Flowfield Modelling and Diagnostics, 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 of Heat 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 of Furnaces 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-NOx Combustion 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 of Turbulent 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 of Minnesota 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) Air Pollution: 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 of Sciences of USSR 

(Institute of Chemical Physics) Moscow. 
Table I. USX Minntac plant: Step 3, line 7 kiln 

operating Conditions (Potts, 1991) 

Inside Volume 1100 m3 

Inside Diameter 5.94 m 

Length 39.6 m 

Mass flow rate of process inlet air 2400 gmolls 

Mass flow rate of air to burner 47 gmol/s 

Mass flow rate ofCI4 to burner 61 gmolls 

Temperature of inlet air (from cooler) 1370 K (2100 OF) 

Temperature of exiting air 1370 K (2100 OF) 

Average N0:2 concentration in exiting gas =500 ppm 

Table n. Predicted equilibrium concentrations 

of NO in air at P == 1 atm. 

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 at P = 1 atm. 

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 of formation (J/gmole) at T = 25 tJC 

and heat capacity data (Himmel blau, 1981): 

c, = a+ bT + cTz + dT3 

( 

1

1

• ) 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.59 
Pellets 

In Dryer 

Gas Out 

Pre heater Kiln 

Air Out 

Figure 1. Taconite pellet induration process. 

PFR PSR 

Pellets 

Figure 2. Ideal reactors in series. 

16 

Pellets Air 

Out In 

Air In 
• 

~ 

• > 

T 

T 

0 

(1-0V) 0V 

y._j 

Figure 3. Qualitative gas and pellet temperature 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 A AFT 

- - ·- - 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 11) - u(1,1))1dg(2) 

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

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

Solve Newton equations for upgraded values of g(). 

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

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) ) I du ( 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,Tf 1Tfi,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 3 0 0 0 ) 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) 

1 0 0 0 format ( I v I I t 15 1 I W ( 1 ) I 1 t 3 0 I I W ( 2 ) I ) 

write(*,2000) 

2 0 0 0 f 0 rma t ( 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,2 
wsav(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 if 

h = 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 
*** 

*** 

*** 

2000 

Total 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 

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