EPA-650/2-74-011
February 1974
Environmental Protection Technology Series
RADIATION
MODELING
*X-Xv
m
1
'!^
-------
EPA-650/274-011
THERMAL RADIATION
MODELING
FOR POLLUTION
PREDICTIONS
by
G. R. Whitacre, R. A. McCann, and A. A. Putnam
Batelle Memorial Institute
Columbus Laboratories
505 King Avenue
Columbus, Ohio 43201
Grant No. R-800842
ROAP No. 21ADG-16
Program Element No. 1AB014
EPA Project Officer: David W. Pershing
Control Systems Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D. C. 20460
February 1974
-------
This report has been reviewed by the Environmental Protection Agency and
approved for publication. Approval does not signify that the contents
necessarily reflect the views and policies of the Agency, nor does
mention of trade names or commercial products constitute endorsement
or recommendation for use.
11
-------
ABSTRACT
An adequate treatment of thermal radiation is essential to a
mathematical model of the combustion process. Accurate temperatures are
especially important for pollutant predictions since the chemical kinetics
are extremely temperature dependent.
Test cases were run which combined the Hottel Zone Method with a
completely analytical flow and mixing solution. However, the complexities
that arise when this approach is extended to real systems seem to eliminate
it as working tool. The results were used to check the accuracy of several
four-flux radiation approximations which are simplier than the zone method.
Although temperature predictions were reasonable the four-flux approximation
failed in predicting the wall radiation fluxes at many locations.
An alternate approximate method which offers great promise is
presented and tested. It employs an expansion of the radiation transport
equation by spherical harmonics. When the first three terms of the associated
Legendre polynominal series are retained, a single second-order differential
equation results. This equation applies over the entire field in contrast
to the four-flux approaches which use separate equations for each direction.
This report was submitted in fulfillment of Grant Number R-800842
by Battelle under the sponsorship of the Environmental Protection Agency.
Work was completed as of September 8, 1973.
-------
TABLE OF CONTENTS
Page
CONCLUSIONS 1.
RECOMMENDATIONS 2
INTRODUCTION 3
DERIVATION OF APPROXIMATE EQUATION OF RADIATIVE TRANSFER BY
HARMONIC EXPANSION METHOD 6
The Equation of Radiative Transfer . 7
Approximation of the Equation of Radiative Transfer 9
Flux 12
Conservation of Energy 13
CALCULATION APPROACH 13
Harmonic-Expansion Approximation • 17
Spalding 4-Flux Approximation 18
Solution Technique 20
Boundary Conditions • • 21
LOW-BTU-FUEL TEST CASE • • 23
POSSIBLE VALIDATION APPROACHES 27
COMBINED HOTTEL-SPALDING MODEL 28
Problems with Combined Hottel-Spalding Approach 33
ETHYLENE TEST CASE 35
REFERENCES 49
NOMENCLATURE 51
APPENDIX
LISTING OF PROGRAM ENCLOS 54
-------
LIST OF FIGURES
Page
FIGURE 1. COORDINATE SYSTEM 8
FIGURE 2. COMPARISON OF RADIATION CALCULATION METHODS 25
FIGURE 3. COMPARISON OF CALCULATION METHODS WITH SWIRL VELOCITY . . 26
FIGURE 4. NETWORK USED FOR COMBINED HOTTEL-SPALDING TEST CASE ... 31
FIGURE 5. COMPARISON OF RADIATION HEAT FLUXES FOR VARIOUS METHODS . 36
FIGURE 6. COMPARISON OF RADIATION FLUXES ALONG CYLINDRICAL WALL . . 38
FIGURE 7. TEMPERATURE DISTRIBUTIONS AT 10-1/2 INCHES 39
FIGURE 8. PREDICTED TEMPERATURE DISTRIBUTION AT 4-1/2 INCHES. ... 46
FIGURE 9. PREDICTED TEMPERATURE DISTRIBUTIONS ALONG CENTERLINE. . . 47
FIGURE 10. PREDICTED TEMPERATURE DISTRIBUTION AT 5/8-INCH RADIUS . . 48
-------
LIST OF TABLES
Page
TABLE 1. PREDICTED TEMPERATURES FOR ETHYLENE TEST CASE USING
HOTTEL ZONE METHOD 41
TABLE 2. PREDICTED TEMPERATURES FOR ETHYLENE TEST CASE USING
BATTELLE HARMONIC-EXPANSION METHOD 42
TABLE 3. PREDICTED TEMPERATURES FOR ETHYLENE TEST CASE USING
USING SPALDING 4-FLUX METHOD WITH a = a. ...... 43
TABLE 4. PREDICTED TEMPERATURES FOR ETHYLENE TEST CASE USING
SPALDING 4-FLUX METHOD WITH a = 2a 44
-------
CONCLUSIONS
The harmonic-expansion radiation approximation method developed
in this study appears to offer great promise as a means of realistically in-
cluding radiation transport in combustion-chamber calculations. This is
especially true when radiation is directly combined with flow and mixing
solutions, because it adds only one equation which applies over the entire
field. The single equation form is also advantageous when considering non-
gray effects, either by frequency banding or sum-of-gray gases. Inc ompari-
son with the standard test-case run using the combined Hottel radiation
Spalding flow code, the harmonic expansion method predicted the radiation
flux to the wall quite well at all points. The temperature predictions were
somewhat low, however, with the greatest difference at the centerline.
The 4-flux radiation approximation did quite well in predicting
the temperature distribution, especially when the flux attenuation coefficient
was increased to twice the absorption coefficient for this test case. How-
ever, it is not clear what method could be used in general to pick this
factor for other cases. (The harmonic-expansion method employs only the
actual absorption coefficient). In addition, the 4-flux approximation did
quite badly in predicting wall radiation fluxes at several points. It is
possible, however, that some other formulation of the discrete-flux equa-
tions could improve this efficiency.
It was demonstrated that it is possible to combine the Hottel zone
radiation analysis with the Spalding flow and mixing code and that results
can be quite useful as a standard for comparative purposes. However, the
difficulties and complexities that arise when this approach is extended to
real systems make it questionable whether this is a practical approach at
the present time. In addition, with some improvements, either the harmonic
expansion or a 4-flux method may have sufficient accuracy to remove the need
for the more complex combined solution.
The harmonic expansion could be used in its present form for cur-
rent design studies involving radiation. It can be used separately from
the flow and mixing code for many applications, and it is not limited to
combustion calculations. Its simplicity, compared to the complexity of
other methods, has much to offer in design studies.
-------
RECOMMENDATIONS
Several remaining tasks should be completed before a final conclusion
can be reached on the choice of a radiation-calculation method. First, the
harmonic expansion approximation should be extended to include frequency banding
to determine whether it is practical and how much effect on results is actually
obtained. An initial band model has already been formulated. Second, frequency
banding could be compared to the sum-of-gray gases approximation used with the
4-flux or Hottel approach. Third, additional test cases with other boundary
conditions should be considered to determine if the initial conclusions
reached in this study hold in general. In this comparison the coupled 4-
flux method developed by Lowes, et al should also be included, if possible.
Fourth, comparisons with actual experimental data should be made provided
reliable data of sufficient detail become available. Fifth, consideration
should be given to the extension of the radiation calculation methods to
three dimensions. Since the harmonic-expansion approximation can still reduce
to a single equation, the advantage over the current 6-flux methods or the
Hottel zone method would be even greater. Also, if a higher level of accuracy
is desired, the harmonic expansion approximation can be improved using higher
order terms.
After a final choice and implementation of the radiation calculation
method is made, the best available kinetic model should then also be added
to form a combined turbulent combustion model for pollutant prediction.
-------
INTRODUCTION
The overall objective of this Grant Program is to develop a mathe-
matical model of the combustion process within enclosures that comprises (a)
the turbulent flow and mixing aspects, (b) the heat-transfer aspects, (c)
the radiation-emission aspects, and (d) the chemical kinetic aspects, all
in sufficient detail to permit the prediction and minimization of pollu—
tant outputs.
The combustion process is a complex interaction of a number of
physical phenomena. Thus, a complete mathematical model would require a
combination of detailed analytical models of the various phenomena. These
include a flow and mixing model, a heat transfer and heat balance model,
a radiation emission model, and a chemical kinetic model.
At the present time, computer models have been developed in all
these areas but have not been combined into one unified model. Therefore,
the effect of phenomena that are not simulated in a particular model must
be supplied as input or implicitly by various assumptions. For example,
existing heat transfer and heat balance models require the specification
of a heat release distribution which is the result of the flow and mixing
process coupled with chemical kinetics. The actual analytical prediction
of flow and mixing patterns is a very recent advance and is still in the
early stages of development. This capability seems to be the key to bringing
together the previously separated models into a single, complete analytical
model of combustion which includes the prediction of pollutant levels.
-------
The flow and mixing in a combustion chamber is usually highly
turbulent and usually involves zones of internal recirculation. Therefore,
neither laminar or boundary layer flow models are adequate to determine
the two-dimensional flow and mixing patterns which are of interest. In
addition, the effects of a swirl velocity is probably not required for
chambers which are at least approximately cylindrical. A turbulent recircu-
lating flow and mixing model which is two dimensional with a swirl velocity
has been developed at Imperial College in England. This model is in use
at many locations worldwide including Battelle. The effect of turbulence
is handled by an effective turbulent viscosity which is presently obtained
from a mixing length that is a function of the air/fuel inlet kinetic energy.
The work covered in this report concentrated on the radiation
emission and heat transfer using the Imperial College model for flow and
mixing. The addition of chemical kinetic aspects was left for future work.
An adequate treatment of radiation is essential to a mathematical model of
the combustion process both to determine heat flux to the walls and/or work
and to predict the actual temperature distribution. The temperature dis-
tribution is especially important since the chemical kinetics of pollutants
are extremely temperature dependent.
A technique called the Hottel Zone Method has for some time
been available to handle the radiation exchange problem. It is essentially
the complete solution for a gray gas system and approximations have been
developed to handle nongray effects. However, although it has been used
for predictions of heat transfer where the heat release and flow pattern
are either measured or assumed ' ' ' ' , it has been questioned whether
it is practical or even possible to combine it with a completely analytical
flow and mixing solution such as that developed at Imperial College by Spalding
et al. This is mainly due to the differences in the method of solution re-
quired. The flow equations are differential and are usually solved by a
fine-network finite-difference approach. On the other hand, the Hottel
radiation solution is integral and requires the matrix inversions of radiation
interchange factors. In most practical problems, a much coarser network
-------
will usually suffice for the Hottel radiation solution alone, and, when a fine
network and/or nongray effects are included,the input data requirements be-
come quite cumbersome. As part of this study, a test case was run with a
combined Spalding flow and Hottel radiation solution. This demonstrated
that it is possible to combine the methods, although it may not be practical
or even necessary in most cases. This combined result also provided a stan-
dard by which various approximate treatments could be compared.
Various approximations to the description of radiant energy trans-
fer have been evolved beginning with the thick and thin approximations.
These early approximations appear adequate only for special cases. However,
flux methods have been employed recently and appear quite promising. The
flux methods employ diffusion-type differential equations for radiation-
flux parameters and are therefore quite compatible with the flow solution,
although they do add additional equations. Simple flux methods are well
(8)
summarized and compared by Siddall in a recent paper. A four-flux model
with two additional second-order differential equations for axisymmetric
coordinates has been developed at Imperial College and incorporated into
(9)
the Spalding flow and mixing code. This model was programmed as part of
this study and added to the turbulent,recirculating-flow computer program
for combustion problems described in Reference 7. Results from this model
are compared with the results from the combined Spalding flow and Hottel
radiation solution. Quite recently an alternative four-flux model has been
(10)
developed by Lowes, et al. It utilizes four first-order differential
equations in which the axial and radial fluxes are coupled. This coupling
is not present in the Imperial College four-flux model. Time did not permit
this model to be programmed and compared with the combined program results
in this present study.
An alternative approximate method for radiative transfer was formu-
lated, programmed, and tested as part of the current program. It employs
an expansion of the radiative transport equation by spherical harmonics
into a series of associated Legendre polynomials. When only the first three
terms of the series are retained, a single second-order differential equa-
tion results which encompasses both axial and radial effects. The results
from this single equation are quite encouraging. In addition, extension to
three dimensions is straightforward and, if necessary, additional terms can
-------
be retained for higher order approximations. The method is particularly
adaptable to frequency banding* since only one additional equation is re-
quired per band. The form of the equation also fits the general solution
in the present flow and mixing codes so that direct coupling to the flow
solution is simple and straightforward.
This report contains the complete derivation of the approximate
equation of radiative transfer by the harmonic expansion method. Next, the
calculation approach by which this method was added to the existing flow
and mixing code is described. In addition, several flux approximations
were also added to the existing code and all the approximate methods were
tested and compared using a low-Btu-fuel test case. This case did not pro-
vide a real validation, however, so that after examining possible validation
approaches it was decided to compare the approximate techniques against the
combined Hottel-Spalding model. The procedure used in combining the models
is presented and the problems associated with this approach are discussed.
The results of the validation test cases employing the combustion of ethy-
lene are compared and discussed in detail. Tentative conclusions are pre-
sented and future work needed to complete this phase of the complete model
development is outlined.
DERIVATION OF APPROXIMATE EQUATION OF
RADIATIVE TRANSFER BY HARMONIC EXPANSION METHOD
The equation of radiative transfer has been solved exactly for
certain cases but, because of the difficulties involved, numerous approxi-
mate formulations have been advanced. It seems safe to say that an approxi-
mate formulation is required for a practical study of furnace heat transfer.
Hopefully, such a formulation will offer a reasonable balance between
accuracy and computational efficiency.
The radiative-transfer situation involves nongray gases and par-
ticulates which can absorb, emit, and scatter radiation. Although scattering
* Frequency banding divides the radiation spectrum into a discrete number
of bands with different radiation properties. The total radiation is ob-
. tained from the sum of each band solution.
-------
is of negligible importance in this investigation, isotropic scattering is
retained in the derivation for more completeness.
It is the intent of the following sections to illustrate the
systematic derivation of a set of approximate equations from the equation
of radiation transfer. The degree of approximation can be reduced by
expanding the set of approximate equations but, naturally, with an in-
crease in computational effort. The equations are developed in cylindrical
coordinates and specialized for circular symmetry.
The Equation of Radiative Transfer
The equation of radiative transfer for a medium which emits, ab-
sorbs, and isotropically scatters is given by
Q'VI = - Hi + ~ I" Idw + (1 - W)H I, . (1)
4TT J . D
The spectral intensity of black-body radiation, I, , is
2
Co
and Q is a unit vector through the point to which the intensity, I, is re-
ferred. The extinction coefficient, K, is the sum of the coefficient of
absorption, a, and the coefficient of scattering while the albedo for scat-
tering, w, is given by the coefficient of scattering divided by H. The
assumptions included in the above equations are that the index of refrac-
tion is unity, polarization need not be considered, there is local thertno-
dynamic equilibrium, and that the process is quasi steady.
Figure 1 depicts the system of coordinates.
-------
FIGURE 1. COORDINATE SYSTEM
The components or direction cosines of the unit vector, Q, are
0 - sin 9'cos
- W)HI
b
(4)
In the next section, the intensity is expressed in terms of
position (r,cp, z) and, utilizing spherical harmonics, angular dependence
(9 '»
-------
Approximation of the Equation of Radiative Transfer
Numerous methods are available for approximating the equa-
tion of radiative transfer or, equivalently, evaluating the radiative trans-
fer in enclosures. One approach is the zone method of analysis.'1' The
zone method has the advantages of good generality and potentially high ac-
curacy. But, unfortunately, computational effort becomes excessive when
the radiation field must be determined in concert with flow fields and
chemical reactions. Another approach is the discrete-ordinates method which
specifies the intensity in an arbitrary number of directions. This leads
to a set of differential equations, the total number of which corresponds
to the number of discrete directions chosen. The logic of this method is
common to many of the "flux methods" which have been applied extensively
(12)
to plane geometries. Extension of this method to other geometries is
difficult.
The method of spherical harmonics can be applied systemat-
ically to geometries other than plane and can give results of arbi-
trarily high accuracy. The computational effort, of course, increases
(14)
along with accuracy. Other methods have utility in special cases
or have fairly general applicability but cannot be extended to arbitrarily
high accuracy.
The spherical-harmonics method has been selected for this inves-
tigation primarily because of its generality, rigorousness, and potentially
high accuracy. Practically, higher accuracy at the expense of additional
labor may not be desirable. For this reason, it should be kept in mind
that some of the other methods may in certain cases offer superior results
for low levels of effort
The intensity can be represented by a series of spherical har-
monics
11
I = Y rA™(r,cp,z) cos m cp' + Bm(r,—* u n n _i n '
n=o m=o
-------
10
where the Pm are associated Legendre polynomials. For the case of circular
symmetry, the intensity must be an even function of the angle cp'; hence,
the B must be zero. Equation (5) may be substituted into Equation (4).
Making use of the relations
cos 9' « P° = u. , (6a)
sin 6' - P* , and (6b)
dw = - d^dcp , (6c)
the result of the substitution is
>.m .m
n /OA . nuv
. .
S 2 (-— cos cp' cos mcp' P: Pm + — - sin cp' sin mcp' P, Pm
n=o m=o ^ dr Y Ylnr Y ^ • 1 n
N
-jH cos mcp' P° Pm + HAm cos mcp' Pm] - wHA°
9 Ylnn Yn/ o
z
(1-w) xib « 0 . (7)
The above relation is equivalent to an infinite set of differen-
tial equations involving the Am. Each equation in this set may be obtained
in turn by virtue of the orthogonality of the tesseral harmonics. That is,
r2* r1 ' k '
** o 1 *
Performing the integration for k=0 and &=Q gives
1 1
and for k=0 and j6=0 gives
-------
11
I TT " ""1 ~ " ' (1Q)
5 Oz
and for k=l and &«1 gives
Q
+ ^ ' 5 -§T + 5
For the first approximation, the terms A , A.., and A, will be re-
tained. Combining Equations (9), (10), and (11) provides the first approxi-
mation to the radiation transfer equation.
. (120
Equation (12a) expressed in more compact form as
7'(n 7Ao) " 3<1-w>KAo = " 3(l-w)Hlb » .
(14)
is also applicable to cartesian coordinates . The intensity, to this
level of approximation, is given by
_ i **** M \x*» •
O A OO/ A*X 1 O r J.
A -77 —r— P. (COS 9 ) - — —r— COS CD P, (COS 6 ]
oHdzl H or • Y 1
-------
Flux
The flux passing through a unit area whose outward normal is given
by the unit vector N may be found from
q • P I N-Q du> . (15)
JQ
The expression for the intensity is given by Equation (14) and for the di-
rection cosines by Equations (3a, b, c). Consequently, the flux directed
in the positive r direction is
and in the negative r direction is
9A°
resulting in a net radial flux given by
"met
Similarly, fluxes in the z direction are
v • "Ao
. 3A°
4rr o fi-7»\
net = - 3^ ~ol ' (17C)
-------
13
Conservation jaf Energy
An expression relating the divergence of the radiative flux to
either known or calculable quantities is required as one factor in the over-
all energy balance. Returning to the radiative transfer equation and inte-
grating both sides of Equation (1) as indicated
J (Q-Al)du> = J* [- HI + ~ j* Idu> + (1 - W)H Ib] dcu , (18)
and noting that
V'qnet = V ' J* m dU) = J Q ' VI da> ' (19)
the flux divergence may be written
V-q = 4rr(l - W)HI - (1 - W)H J* I duo . (20)
Ilw L> U •* wCTT
Substitution of Equation (14) for the intensity into the above equation
gives the result
= 4u(l - w)H (Ib - A°) (21)
CALCULATION APPROACH
In all the computations performed in this study, the additional
equations and terms resulting from the various treatments of radiative trans-
port were put into the same general format as the flow and mixing model
equations developed by Spalding and described in Reference 7. In most cases
the entire set of flow, mixing, and radiation equations were solved simul-
taneously. In this way, compatibility of the radiation technique with the
flow solution could be assured and the behavior of the solution of the com-
plete set of equations could be studied. In certain of the test-case com-
parisons, however, it was decided to input a flow and mixing solution and
-------
14
allow only the energy and radiation equations to be solved with the various
methods. In this way, a more direct comparison of the radiation treatments
could be made without the added complexity of coupling effects. The details
of the flow and mixing model are net repeated here, since it is well
documented in Reference 7. It is recognized that certain improvements have
been made to both the model and the method of solution.* However, since
these changes were still not finalized, it was decided to simply use the
original version of the flow and mixing code as the starting point in this
comparative study of radiation-prediction methods, Improvements can be
added later and they should not alter this general conclusion of the radia-
tion study.
In the two-dimensional flow and mixing model used, a simultaneous
set of differential equations for stream function, vorticity, and mixture
fraction are solved iteratively by the method of relaxation. In addition,
an equation for swirl velocity in the azimuthal direction for axisymmetric
flow can be added. The effect of turbulence is treated by a local turbulent
viscosity that is related to local conditions by the semiempirical formula
_ _2/3 Tl/3 2/3 ,. .. 2 . 2.1/3 ,„.
u -- = K D L n (m,- V.. + m V ) . (22)
Keff v ^ fu f\i ox ox J ' '
where
fi -f is the local turbulent viscosity
D is the combustion chamber diameter
L is the combustion chamber length
p is the local density
m is mass flow rate
V is velocity
and the subscripts fu and ox refer to conditions in the fuel and air
inlet, respectively. A value of K equal to 0.012 recommended by Pun and
Spalding^ was used for all calculations in this study. Recent work '
* Some work was attempted in this study with a revised solution technique.
The results were encouraging, however, additional difficulties were en-
countered and the approach was dropped for the current efforts.
-------
15
has indicated that a value of 0,006 may be more representative for a large
furnace. The mixing and combustion model used assumes a physically con-
trolled single-step reaction with equal diffusion coefficients.
In the original combustion code given in Reference 7, the addi-
tional assumptions of no radiation heat transfer, adiabatic wall-boundary
conditions, and negligible kinetic heating were made. Under these condi-
tions the energy-conservation equation and its boundary conditions becomes
exactly similar to the concentration equation, and the enthalpy and tempera-
ture can be obtained directly from the mixture fraction by auxiliary rela-
tions. When considering radiation, however, it becomes necessary to solve
both the mixture fraction equation (which remains the same) and an equation
for the specific enthalpy defined by
h = C T + h. m, , (23)
p fu fu ' N '
where
h is the specific enthalpy with the base state taken as
absolute zero
C is the specific heat at constant pressure, which is
" taken as identical for all species and constant with
temperature, I
mf is the mass fraction of the fuel
h. is heat of reaction.
This definition can be modified to allow C to vary with temperature, but
this was not done in the present study. The temperature must be calculated
from both the specific enthalpy and mixture fraction rather than from only
the mixture fraction as in the original code. To calculate tem-
perature, two different regions must first be defined, divided by the
stoichiometric mixture fraction
where i is the mass of air which combines with a unit mass of fuel. In the
fuel-products region, f < f £ 1, the temperature is given by
S t
-------
16
h-h 1(1 + 1) £.Ij
—i - . (25a)
In the air-products region, 0 £ f £ f
81
~ . (25b)
P
The specific enthalpy or energy- conservation equation, even in-
cluding all the radiation approximations considered in this study, can be
expressed in the same form as the general elliptic differential equation of
Reference 7. This form is
IT
• a? Ib2 ir (c$)1 + d = ° » (26)
where
z is the axial coordinate
r is the radial coordinate
$ is the conserved dependent variable
¥ is the stream function
a, is the function multiplying the convection terms
b-, b , and C are functions connected with the diffusion terms
d is the source term.
For the specific enthalpy equation
«1 - ! • (27a)
bx " r Th , (27b)
b2 = r Th , and (27c)
C - 1 , (27d)
-------
17
where F. is the effective exchange coefficient for heat, which is equal to
the effective thermal conductivity divided by the specific heat. For tur-
bulent flow, it is obtained by dividing the effective turbulent viscosity
by the Prandtl number. The term d in the enthalpy equation contains the
appropriate source terms for the various radiation approximations. This
is discussed in detail later.
Harmonic >Expans ion Approximat ion
In programming the harmonic expansion approximation derived in
the previous section, the scattering coefficient was set equal to zero.
Therefore, the extinction coefficient, K, becomes simply the absorption
coefficient, of, and the albedo for scattering, w, becomes 0. (There is
no special problems in including scattering, but for comparative purposes
it did not appear necessary.) The independent radiation variable was de-
fined as
U » TT A° . (28)
Therefore, the source term in the energy equation from Equation (21) becomes
denergy = 4 (1 ' w) H (E " U) • (29a)
4
where E = cr T when the entire frequency spectrum is'considered in a single
equation. For no scattering
d » 4 a (E - U) . (29b)
energy
The differential equation for U from Equation (12a) becomes
sr $
or for no scattering
-------
18
This fits the general differential Equation (26) with
a = 0 , (31a)
bx - I/a , (31b)
b2 = r/a , (31c)
C - 1 , (31d)
d = 3 & (U - E) , (31e)
plus an additional 1/r multiplier on the r-direction diffusion term. The
radiation variable, U, is a scalar potential since the net flux in the r di-
rection
and the net flux in the z direction
«.•-&£ • <32b>
Spalding 4-Flux Approximation
To compare the harmonic expansion approximation with a dis-
(n\
crete ordinate flux approximation, the 4-flux model developed by Spaldingv '
was also programmed. The derivation of this model is presented in the
Symposium paper as well as in an earlier paper by Lockwood and Spalding
and is not repeated here. The derivation is different in that it con-
tains an additional term in the radial-direction equations that is not
*
included in most other derivations with the form of the inward flux divided
-------
19
by the radius. The addition of this term causes the model to reduce to a
known exact solution in the special case of radiation between coaxial
cylinders without an absorbing gas.
The derivation results in four first-order differential equations
which can be reduced to two second-order differential equations and can then
be fit into the general form of the flow and mixing equations. The two in-
dependent radiation parameters are F , the average of the positive and
Z
negative radiation fluxes in the z direction, and F , the average of the
fluxes in the r direction. The model can handle scattering, although
the present study does not include scattering. There is an uncertainty in
the model since the flux absorption coefficient, a, appearing in the equa-
tions is not necessarily the conventional absorption coefficient, a, defined
as the constant of proportionality in the attenuation law for radiation in-
tensity which appears in the harmonic expansion method. Following the sug-
gestion of Spalding, two values of "a" have been used for all calculations
with the A-flux model in this study. The values used are: a - ot, and
a = 2 a, which represent approximate limits on "a" and also demonstrate the
sensitivity of the predictions to the value of "a".
The final equations for F and F for the no-scattering case
r z
become
i A dF
-££
-------
20
There are similarities with the results from the harmonic expansion method
[Equations (29b) and (30b)]. However, there are significant differences
also since the harmonic expansion method results in a single partial dif-
ferential equation over the entire field while the 4-flux method results in
separate total differential equations in the orthogonal directions only.
The net fluxes in these two directions only can be obtained in the 4-flux
method by
dF
- ' (36a)
and
dF
-3T ' (36b)
Solution Technique
When solving both the harmonic expansion and the 4-flux approxi-
mations with the method of relaxation used in Reference 7, it was found
quite desirable to over-relax the radiation equations for U or F and FZ by
a relaxation factor as large as 1.9. This means that the radiation variable
was changed in each iteration. (The theoretical maximum relaxation factor
conditions at that iteration. (The theoretical maximum relaxation factor
is 2.0.) This speeded up the convergence both when complete flow mixing
and radiation equations were solved simultaneously and when the flow and
mixing solution was fixed and only the energy and radiation equations were
solved. It was also found desirable in many cases to under- relax the energy
equation with a factor as low as 0.5 to prevent divergent conditions that
could drive the solution to negative values of the temperature, particularly
near the fuel inlet.
To solve the additional radiation equations required added time
and, in general, the number of iterations required to reach the same
level of convergence increased by about 50 percent when radiation was
added to the flow solution. (The vorticity was still generally the
last to converge, however.) For the cases studied there was no significant
-------
21
difference in calculation time between the harmonic expansion and the 4- flux
methods.
Boundary Conditions
The walls of the combustion chamber are generally either water-
cooled or refractory surfaces or composites of both. Therefore, a conven-
ient wall boundary condition formulation is one that can handle such a
composite wall, at least in an approximate manner. Such a treatment has
(9)
been used by Gosman and Lockwood with the 4-flux method and was used in
this study for both the 4-flux and harmonic expansion cases. The assump-
tions made are: the water-cooled surfaces have negligible emissive power
and are gray, and the refractory portion of the wall is gray, diffuse, and
reradiates all radiation received. The net radiation heat transfer, q ,
across the wall then becomes
where
ew R V ' (37)
e is the emissivity of the cooled surface
R is the ratio of the water-cooled area to the total
surface area of the boundary node
q . is the radiation flux approaching the wall in the
normal direction,
For the cylindrical outer wall with the harmonic -expansion approximation
from Equation (16a)
U - r— r— i , and (38)
w 3n or '
'w
w
Combining Equations (37), (38), and (39) yields
-------
22
w w
Similarly for the end walls,
2- e T
& A
w w
in
± uw
In the Spalding 4-flux approximation a similar analysis yields
2 - e R dF
w r 'w w
and
2 - e R dF
w z w w
The preceding wall-boundary conditions were programmed for both
methods with e R as an input parameter. However, for all the computations
Vr
performed in this study e R was 1.0, i.e., completely water-cooled black
walls. In addition, the wall -boundary condition used for the energy or en-
thalpy remained zero gradient at the wall, i.e., no convective heat trans-
fer at wall. Although these conditions are idealized, they were applied
consistently for all methods so that comparison of results should be valid.
At the centerline of the axisymmetric coordinate system symmetry
requires that q $,= 0.
rnet
For the harmonic expansion from Equation (32a)
,r ---0 . (44)
net
therefore
-------
23
for the 4- flux approximation
dF
' - 2 F
net
However, since F at r c 0 vanishes, no condition can be placed on dF /dr
at r = 0 since symmetry is assured by the presence of F . One possible so-
lution to this problem might be to consider the centerline nodal points as
special points in the field rather than boundary points and use a special
finite-difference formulation to solve for F at the centerline. However,
this did not fit conveniently into the method of solution of Reference 7
and was dropped. Zero gradients at boundaries are handled in Reference 7
by fitting a quadratic to the independent variable at the first two interior
points plus the boundary point^and this approach was used for the harmonic
expansion U at the centerline. Since the gradient could not be specified
for F ,a quadratic was fit at the first three interior points and then ex-
trapolated to the centerline. In a coarse network this procedure may intro-
duce some inaccuracy at the centerline for the F solution and might increase
differences between the harmonic expansion and 4-flux methods.
LOW-BTU-FUEL TEST CASE
To compare results of the harmonic expansiom method and
the flux methods,the test case for the combustion of a low-Btu fuel given
in Reference 7 was solved with the addition of the radiation approximations.
The same assumptions and input and boundary conditions were used as de-
scribed in Reference 7. A constant absorption coefficient, a, of 2.0 ft~
was used for all the calculations; this is representative of a very luminous
flame. The idealized fuel used in the test case was a gas having the same
molecular weight as air and requiring 2 pounds of air for complete combus-
tion of each pound of fuel. The chemical reaction rate was assumed to be
-------
24
rapid enough for combustion to be controlled by mixing alone. The heat of
combustion of the fuel was taken as 1789 Btu/lb of fuel. The combustor
geometry used had a diameter of 0.5 ft and a length of 10 ft to a semi-
restricted outlet.
Figure 2 shows a comparison of various radiation-calculation
methods for this test case without swirl velocity. The curves show the cal-
culated temperatures across the turbulent diffusion flame 1.6 combustor
diameters downstream from the fuel and air inlets with both inlet velocities
100 ft/sec. Results are shown for the Spalding 4-flux method with the flux
absorption coefficient, a, equal to the conventional absorption coefficient,
(9)
a, as well as twice a, Gosman and Lockwood imply that the true answer
should fall between these limiting values but give no method for actually
picking a single value. For this case, the results using the single-equation
harmonic-expansion method with the conventional absorption coefficient
fall between the two different 4-flux results.
It has been suggested that it may be possible in cylindrical com-
•
bustion chambers to consider only the radial 2-flux equation as an approxi-
mation and to neglect the axial radiation flux. To check this ap-
proximation, the test case was also run using only the radial equation of
the Spalding 4-flux model, and the results are shown. It appears that this
would not be an acceptable approximation in the present cases since the
predicted results differ from the more complete solutions as great as 200 R
can occur.
Figure 3 shows another comparison of radiation-calculation methods
for the same test case, except that the inlet air is assumed to have a tan-
gential velocity component of 500 ft/sec. All other conditions remain the
same as for the previous case. To include the effect of the swirl, an
addition differential equation in swirl velocity is added as discussed in
Reference 7. Then with the 4-flux method a total of seven differential
equations were solved simultaneously. Because the effect of swirl causes
a much smaller flame envelope, calculated temperatures closer to the inlet
(at 0.4 combustor diameters) are compared. The trends are quite similar
to those for the previous case; the harmonic-expansion results generally
fall between the two results from the Spalding 4-flux method. The 2-flux
-------
25
2300
2300
— a-Flux (Radial Only)
Spalding 4-Flux a = a
Battelle Harmonic Expansion
Spalding 4-Flux a = 2a
Without Swirl Velocity
1200
MOO
0.05
0.10 0.15
Radius, ft
0.20
0.25
FIGURE 2. COMPARISON OF RADIATION CALCULATION METHODS
-------
2700
2500
With Swirl Velocity at 0.4 Diameters!
Radiation
lux a=a(Radial Only)
Spalding 4-Flux a=a
Battelle Harmonic Expansion
Spalding 4-Flux a = 2a
900
0.05
0.10 0.15
Radius, ft
0.20
N>
0.25
FIGURE 3. COMPARISON OF CALCULATION METHODS WITH SWIRL VELOCITY
-------
27
approximation again appears inadequate,and, in fact gives results approxi-
mately halfway between the no-radiation results and the full 4-flux solu-
tion. This would indicate that the axial flux is just as important as the
radial flux (for this test case at least).
It is recognized that these limited results at relatively low
temperatures do not establish the accuracy of either approximate method.
They do demonstrate that either the harmonic expansion or the 4-flux method
can be added to the flow and mixing solution, and that the results from the
single equation harmonic expansion method generally fall between the
limiting values of the two-equation 4-flux method.
POSSIBLE VALIDATION APPROACHES
To determine the accuracy of the radiation approximations,
a standard test case or, better still, cases are needed to which they can be
compared. One approach would be to compare directly against experiments;
thus, a survey of available data was made. However, it was found that there
are many uncertainties in the data resulting from the inherent difficulties
in measurements in a high-temperature combustion system. Also, real systems
generally have absorption coefficients varying with location and with fre-
quency (non-grayness) which are difficult to specify accurately. In addi-
tion, if the flow and mixing solution is obtained analytically the differ-
ence in results between the various radiation treatments may be completely
overshadowed by the difference between the predicted and actual flow and
concentration profiles. On the other hand, rarely are these profiles mea-
sured so that they could be supplied as fixed input with only the radiation
treatment changed.
Because of all these difficulties and uncertainties it was de-
cided to compare the radiation approximations by a different approach. At
least for a gray system with constant absorption coefficient, the Hottel
zone approach is recognized as an essentially exact solution to the radia-
tion transport. Therefore, the Hottel zone method was chosen as a standard
by which to compare radiation approximations under the condition of constant
absorption coefficient. Hopefully, a single radiation approach can then be
-------
28
chosen which can be extended to handle real effects including non grayness,
either by frequency banding or a sum of gray gases approximation.
COMBINED HOTTEL-SP ALPING MODEL
The Hottel zone approach to radiation transfer is an integral
approach to the problem rather than one of differential approximation, in
that it «6nsiders all possible interchanges of radiation between gas and surface
nodes. If it is used in conjunction with a flow pattern and a mixing or
heat-release distribution, it can be coupled with a series of heat balances
(energy equations) to predict the temperature field and resulting heat
fluxes. This approach can be useful in certain problems and has been used
(2 3 4) (5)
by several investigatorsv ' ' ' including the present authorsv '. It can
also be used in a more restricted sense by specifying the temperature field
and calculating the radiation fluxes. This latter approach has been used
by Lowes, et al to provide a basis for comparison of several flux
/0\
models. Siddall also compared the zone method to several flux models
and to the exact solution for a simplified one-dimensional problem.
In the present study, however, a new approach was taken in that
the Hottel zone technique was coupled directly to the flow and mixing pro-
gram of Reference 7. This was done for several reasons. First, it would
provide a standard of comparison which assures that only the treatment of
radiation is changed while allowing complete freedom of all the independent
variables. This follows from the fact that the various approximations to
be compared were also coupled to the same flow and mixing program.
Second, the coupled approach would demonstrate that it is possible
to solve the combined system, even though the flow program is differential
and a zone temperature is linked only to its immediate neighbors while the
zone method considers all possible linkages.
Third, the coupled approach, if found to be possible, would give
an indication of whether it would actually be practical for use in a
working code. If it is practical, including economic considerations, to
use the zone approach expanded to include non gray effects there would
actually be no real need for the flux or harmonic expansion approximations.
-------
29
The changes to the flow and mixing code to add the Hottel zone
radiation method are actually quite simple, in fact, less complicated than
the changes to add the approximate methods. The difficulties are mainly
involved with generating the input required for the calculating the Hottel
total-exchange areas and devising a system to convert the nodal points of
the flow code network into the gas and surface zones of the Hottel approach.
In the solution technique no additional differential equation is required,
although both the specific-enthalpy equation and the mixture fracture equa-
tion must be used because they are no longer linearly related as in the
special case of Reference 7. In the Hottel approach the term in the energy
equation [Equation (26-) becomesj
S G. G. E . + Z S. G. E .
, j ig g, J ,k ig s, k
d = 4 or. E . - J — , (47)
energy ig g, ig vig
where
Of is the absorption coefficient of the i gas zone
° which surrounds the nodal point of the flow code
network
E . is the black body emissive power of the i gas zone
g) ^g equal to a TJ4
G. G. is the total exchange area between gas zone j and gas
J zone ig
S. G. is the total exchange area between surface zone k and
° gas zone ig
E , is the black body emissive power of the k surface
S • K
zone
v. is the volume of the i gas zone.
The summation over j includes all the gas zones including the i zone and
the summation over k includes all surface zones. The total exchange areas
used in the Hottel zone method are defined as the factor by which the
difference in emissive powers of any pairs of zones must be multiplied to
obtain the total net radiation flux between the pair of zones. The total
radiation flux between zones includes consideration of the diffuse wall re-
flections from all surface zones and absorption of gas zones in addition
to the direct exchange between the zones. The derivation of the procedure
-------
30
for obtaining the total exchange areas is given in Reference 2 and is not
repeated here. The procedure involves assigning in turn unit radiation
emissive power on one zone and zero on all others and solving simultaneously
the radiation balances on the various surfaces. The resulting solutions for
the total exchange areas, although quite complex, can be expressed rather
easily in Fortran using DO-loops and matrix notation. A listing of the pro-
gram ENCLOS, developed to calculate the total exchange areast is given in the
Appendix* For the gray-test-case run in the present study this program
was run once external to the flow and mixing program and the calculated total
exchange areas were then supplied to the main flow code which contained the
energy equation with the additional source term given by Equation (47).
The input required for the calculation of the total exchange areas
include the absorption coefficients for each gas zone, the surface emissi-
vities of each surface zone,and the direct exchange areas for all pairs of
zones. The direct exchange areas are available in tabular or graphical form
for uniform cubic or cylindrical networks and other simple configurations.
Also, separate codes are available to calculate geometric view factors
for general configurations. The direct interchange areas include the effect
of the attenuation by gas absorption between the zones in addition to the
geometric view factor. For the present test case, the'tables for cylindrical
/ION
enclosures are calculated by Erkku and given in the Appendix to Chapter
7 of Reference 1. The direct exchange is nondimensionalized by the square
of the zone dimension in Erkku's tables so that program ENCLOS was written
to accept these nondimensional factors as input. For the test case, a zone
dimension of 0.25 ft was used together with a constant absorption of 1.0 ft
for a nondimensional attenuation factor of 0.25.
Figure 4 shows the network used for the combined Hottel-Spalding
standard test case run. Achieving compatibility and conversion between the
flow code nodal network and Hottel zones is a major problem in using the
combined code. In order to use available data,the maximum zonal network
for which direct exchange areas were calculated by Erkku was used. This
network contains five radial rings of zones with twelve zones each in the
axial direction for a total of 60 gas zones. This results in 22 Hottel sur-
face nodes and a physical size with a radius of 15 inches and a length of
36 inches since each zone is 3 inches by 3 inches. A nodal flow network
-------
8
Hottel Surface Zones
IO II 12 13
JN-7—
6sH
t
Air
Fuel
54-
•4 3-
32-
2 i-
T • ^ «--4-
14
.JL
15 16 17
T • T •—£—••
• 5 i •10) «I5 «20 »25 | «50 | «35 | «40 | «4S | «50 | »55
• 14 . «I9 | «Z4
*60 <
34 »39 *44 »49
O*3, *8 I *I3 »I8 *23 I »28
—j--t—;—
• 2 *7 I *12 I *17 I *22 I *27
--I--I-—I — 4— ->--
• 33 | *3S | *43 I »48 I «53 | *58 <
|--+-t- + ~t ' —
| «32 I *37 I «42 | «47 I «52 | «57 '
!_
*I6
I
J 1
• 41 | »46 *5I I »56 <
H8
20
u>
Outlet
•22
1234
I ^^ Nodal Network
8 9 IO II 12 13 14
IN
•c.—
FIGURE 4. NETWORK USED FOR COMBINED HOTTEL-SPALDING TEST CASE
-------
32
was then superimposed on this zonal network with 7 grid lines of constant
radius and 14 grid lines of constant axial distance. It was felt that this
network would give a realistic answer to the finite-difference solution of
the flow and mixing equations although perhaps not the exact same patterns
as a finer nodal network would yield.
In order to convert from the flow nodal network numbering system
to the Hottel gas zone number the following transformation is used.
ig = J - 1 + (1-2) (JN-2) , (48)
where
ig is the Hottel gas-zone number
I and J are flow network indices
JN is the total number of constant radius grid lines.
This transformation is general so long as the boundaries form a right cir-
cular cylinder. Similar transformations are also required to convert to
the Hottel surface zones along the boundaries. However, the centerline is
not a boundary in the Hottel zone approach. In fact, the Hottel surface
zone heat balances do not enter into the energy-equation solution explicitly
unless the radiation flux is assumed to be related to the convective boun-
dary condition. If, however, the wall radiation fluxes are desired outputs,
they can be obtained for each surface zone by the following expression
= 5ItiE«.j + j[VSiE..j
J A7Ji si Vi • <49>
where
e. is the emissivity of the i surface zone
A is the surface area and the remaining terms and summations
have been defined with Equation (47).
For the special boundary condition of the test case in which it
is assumed that there is no reradiation from the cold wall, Equations (47)
-------
33
and (49) are simplified. In Equation (47) the summation with the S, G. E ,
E K Ig S ,K
terms is dropped and in Equation (49) both the summation with S, S. E . and
tC 3- S j K.
the term 6. E . are dropped.
x s j i
Problems with Combined Hottel-Spalding Approach
Although it was possible to run the combined Hottel-Spalding Model
as a standard test case there are a number of problems in its use as an
actual working code. The network used for the test case shown in Figure 3
was chosen mainly to fit the zone approach. Check runs were made with two
and three times the number of constant radius grid lines and the flow solu-
tion did change, indicating a finer network should be used for the finite
difference to obtain a converged solution. (Because the same network was
used for all the comparisons in the test case, results should be valid for
comparative purposes, however).
The increase in the number of zones required with a finer network
would be a major problem. The additional direct interchange factors that
would be required would further complicate the input and, in fact, make
automatic input generation essential. In addition, the zone method must
invert a square matrix whose size is determined by the total number of gas
zones. For the test cases with 60 gas zones, a single inversion and related
calculations to compute a set of total exchange areas required 115 seconds
on a CDC Cyber 73 computer. Using these areas as input, the combined Hottel-
Spalding run required an additional 42 seconds. In comparison, the complete
flow mixing and harmonic expansion radiation solution required a total of
42 seconds. For finer networks, the time required to calculate the total
exchange areas would increase at a rate almost proportional to the square
of the number of gas zones. The time for other calculations would also in-
crease but only approximately proportional to the number of zones.
For the test case, the absorption coefficient was assumed
constant. However, for most practical cases the absorption coefficient
will vary with temperature (non grayness due to frequence dependence), and
also will vary with local composition. When coupled with a flow and mixing
solution, temperature and composition not only vary with position but can
-------
34
also change from one iteration to the next at the same point. Although in
theory the zone method can handle variations in absorption coefficient, in
practice it becomes extremely difficult. In the first place, the direct
exchange areas themselves are dependent on the absorption coefficients or
attenuation along the path between the pair of zones. Hottel has devel-
oped an approximate method whereby an average absorption coefficient is
used to obtain a single set of direct exchange areas and a correction is ap-
plied later during the calculation of the total exchange areas. Such an
approximation is probably required if the zone method is to be considered
at all practical. Even with a single set of direct exchange areas, in
theory, the total exchange areas should be recalculated each time the ab-
sorption coefficients change. Therefore, if the absorption coefficient is
expressed directly as a function of temperature and/or composition the
ENGLOS program must be run as a subroutine directly with the main flow pro-
gram and the complete matrix inverted each time absorption coefficients
change. Preliminary runs were made on a simplified system in which a was
expressed directly as a function of temperature. It was found that ENCLOS
did not need to be rerun for each iteration but could be called only when
temperatures had changed by some specified amount. This technique reduced
the running time by approximately one-half compared to an inversion with
each iteration, but for practical cases the running time would still appear
to be excessive, at least for a fine network.
A possible alternative has been suggested by Hottel through
the use of the sum-of-gray-gases approach. In this technique the gas ab-
sorption and emissivity are described empirically as a weighted sum of gray
gases. The weighting factors are functions of temperature to handle non-
grayness and could also be function of composition. The matrix inversion
need be performed only once for each gray gas (usually three) plus one clear
gas, and the variation in at is handled by the change in the weighting factors
in a relatively straightforward manner. However, the difficulty comes in
choosing the gray gases and the weighting-factor functions. Although it
has been done by Hottel for several fixed-composition gases and extended by
(3)
Johnson for a particular luminous flame, it does not appear as yet really
a practical approach for an actual working code.
-------
35
Because of all the difficulties associated with the combined-code
approach for an actual design code, in this study the combined code was used
only to supply standard test-case results to compare approximate techniques.
The difficulties were avoided by using constant or, a fairly coarse network,
and simplified boundary conditions. If an approximate treatment appears
adequate for practical use there is really no need to attempt to overcome
all the combined code difficulties.
ETHYLENE TEST CASE
In order to compare the various approximate radiation calculation
methods with the combined Hottel-Spalding treatment, a test case with high
radiation was desired to accentuate the differences. However, it should
still be physically realistic rather than a completely hypothetical case.
Therefore, the combustion of ethylene (C^) was chosen as the standard test
case. The high adiabatic flame temperature combined with a luminous type
flame maximize radiation and tend to fit the constant=absorption-coefficient
assumption. A value of 1.0 ft was used for the absorption coefficient,
throughout the combustion chamber. The chamber was of 15 inches radius by
36 inches long with the uniform network shown in Figure 4. The basic gas
zone dimension is 0.25 ft thus giving a constant nondimensional attenuation
factor of 0.25 which is used to obtain the direct exchange areas.
The heat of reaction used was 21,636 Btu/lb of fuel with i, the
mass of air which combines with a unit mass of fuel, taken as 14.28. This
gives a stoichiometric mixture fraction, f defined by Equation (24) of
s t
0.0633. The test case was run with an overall air-to-fuel ratio of 16 to 1
which is, therefore,slightly on the lean side. Another advantage in choosing
ethylene with a molecular weight of 28 is that the simplifying assumption
can be made that the air and fuel have equal densities. A constant mean
specific heat of 0.33 Btu/lb F was also used.
As previously discussed, all the radiation methods tested were
combined with the same flow and mixing code and utilized the same general
solution technique. The same input conditions were used for all methods.
Figure 5 shows the comparison of the calculated radiation heat
flux to the wall for the various calculation methods tested. The predicted
-------
80
70
CJ
8J60
CM
CO
150
x40
ZJ
o
X 30
c
o
o
o20*
10
0 3 6 9 12
Inlet End Wall
Radius Inches
. I I I ,1
IETHYLENE TEST CASE|
Legend
Hottel Zone Method
Bottelle Harmonic Expansion
Spalding 4-Flux a=a
Spalding 4-Flux a =2a
i rn r
15
0
12 18 24
Distance Along Wall, inches
15 12 9 63
30 36 Out|et End Wa,|
Radius Inches
0
FIGURE 5. COMPARISON OF RADIATION HEAT FLUXES FOR VARIOUS METHODS
-------
37
fluxes are shown for each surface zone along the inlet and outlet end walls
as well as the zones along the cylindrical wall. The harmonic-expansion
method agrees quite well with the results from the Hottel zone method with
differences usually less than 10 percent. On the other hand, the results
from the Spalding 4-flux method differ markedly from the zone result. In-
creasing the flux-attenuation coefficient to twice the absorption coefficient
causes some improvement but there are still large disagreements. These ap-
pear to be caused by the two-directions-only nature of the 4-flux treatment
which does not include coupling between the r and z directions. The very
high results in the second zone from the center on both end walls reflect
the high-temperature region occurring along the entire chamber at that radius.
Both the zone method and the harmonic-expansion method allow radiation from
the region to reach the walls in all directions. For the same reason the
4-flux treatment underpredicts the flux in the inlet corner.
Figure 6 shows a more detailed comparison of predicted radiation
fluxes along the cylindrical wall only. In addition to the methods shown
on the previous figure, 2-flux method results are also shown. In the 2-flux
method used the axial equation was dropped and the same equation used for
the radial direction as was used in the Spalding 4-flux. For this particu-
lar case the 2-flux methods actually are in better agreement than the 4-flux
method with the zone method along the cylindrical wall. However, the 2-flux
method gives no prediction at all on the end-wall radiation fluxes. The
t
harmonic expansion prediction is better than any of the flux methods.
Figure 7 shows the predicted temperature distributions for the
various methods at a constant distance of 10-1/2 inches. Although the 4-flux
methods did not do well in predicting the radiation flux to the wall, for
this case the 4-flux methods do reasonably well in predicting the tempera-
ture distribution. The prediction from the zone method generally falls be-
tween the values predicted using the two different values of absorption
coefficient with the 4-flux method. The value of a = 2ot appears to be a
somewhat better choice for this case. As the 2-flux method does not appear
adequate for temperature prediction, particularly near the peak, it was
therefore dropped from further consideration. No curve can be shown for
the zero-radiation case because no converged solution could be obtained.
-------
20
18
16
CM
CD
s
3
o
0)
X
O
10
8
O
or
6
4-
2-
Legend
Hottel Zone Method
Battelle Harmonic Expansion
Spalding 4-Flux a = a
Spalding 4-Flux a = 2a
2-Flux a = 2a(Radial Only)
2-Flux o = a (Radial Only)
6 12 18 24 30
Distance Along Cylindrical Wall,-inches
36
FIGURE 6. COMPARISON OF RADIATION FLUXES ALONG CYLINDRICAL WALL
-------
39
3800
3600
3400
3200
3000
o:
u
O
i_
Q)
2800
2600
2400
2200
2000
1800
IETHYLENE TEST CASE I
Legend
Hottel Zone Method
Battelle Harmonic Expansion
Spalding 4-Flux a = a
Spalding 4-Flux a=2a
2-Flux a=a (Radial Only)
2-Flux a = 2a(Radial Only)
1600
0.25
0.50 0.75
Radius, ft
1.00
1.25
FIGURE 7. TEMPERATURE DISTRIBUTIONS AT 10-1/2 INCHES
-------
40
For this particular case, the presence of any of the radiation methods
actually stabilized the complete solution.
The temperature prediction from the harmonic-expansion method is
somewhat lower than the zone method results near the centerline and at peak
conditions. The agreement near the outer wall is excellent, however, The
consistency of the difference indicates that a higher order approximation
or possibly even some simple correction might greatly improve the tempera-
ture prediction of the harmonic expansion. Some initial work along these
lines was performed. However, even in the present single-equation simple
form the temperature prediction of the harmonic expansion is probably ade-
quate for most purposes.
All the previous test case results were run using the combined
flow and mixing code with the various radiation methods added. This did
confirm that all the treatments were compatible in a complete solution. In
addition, it meant that the flow and mixing solution was free to change as
the radiation treatment affected the temperature. While this comparison is
quite useful it does tend to exaggerate the effect of differences in radia-
tion treatment, particularly in the fuel-rich region. As a more direct com-
parison of radiation treatments alone, the test case was rerun with the
flow and mixing solution fixed with a pattern obtained from the complete
solution with the harmonic expansion. The same code and solution techniques
were employed but only the energy or specific enthalpy equation and the ap-
propriate radiation equation(s) were solved using the flow and mixing solu-
tions as input.
The predicted radiation heat fluxes to the wall with the fixed
flow solution were almost identical to the results with the variable-flow
solution. Therefore, Figures 4 and 5 and the comments concerning them apply
for this case also. The spread in the temperature distribution between the
various methods was reduced somewhat as compared to the variable flow
results.
Tables 1 through 4 show the predicted temperatures for nodal
points shown in Figure 4 for the ethylene test case for the Hottel zone
method, the Battelle harmonic expansion method, and the Spalding 4-flux
method with two values of a ( a = ot and a = 2o?), respectively. The tempera-
tures are given in R. All results are given for the fixed flow and mixing
-------
T«E
t - A(I,J,1E*1>
NUMBER OF 1IEKATIUN =
7
6
b
4
3
2
1
1
1
1.0.
5.
5«
5.
5.
1.264E* 3
1.245E+C3
9.67bt»r2
2.iiit*.:3
1.720t*03
1.85iE*03
l.BVOE+03
1.953E.+03
"
3.5V2E+03
3.860E+03
2.352E*03
2.73oE*03
2.967|*C3
3J473C.* 03
3.112t*63
2.037E*03
2.393E»fl3
3.5Q1E+03
3.4b7E*03
3.i34E*fl3
3.099E*o3
2.036E*03
2.076E*i)3
2.411E4.03
3.391E*03
3.451E+03
3.139E+03
3.100E+03
2.Q77E+03 2.11*E*03
3.282E+03
3.133E*03 3.125E+Q3
10
7
6
5
4
3
2
1
2.173t*';3 2.186t*03 2.1B8E+C3
2.177L*U3 2.197E»i3 2.2v6t+OJ 2.207E*03
2.42.E*i)3 S.S^St*:.^ 2«3feot*03 «i.356E*03
2.913£*^3 2.Bc,4t*u3 2.790E*03
3.47Rt*-';3 3.46ot + 6-i 3.453E+C3
3.y29E*(!3
11
12
13
14
3.'849t*02
1.000E*02
RtYNOLOS NO BASED ON MEAN MASS FLO* RAjE AND MAX OlA = 115*8
TABLE 1. PREDICTED TEMPERATURES FOR ETHYLENE TEST CASE USING HOTTEL ZONE METHOD
-------
^JMBER OF ITERATION * *a
-
7
6
3
4
3
2
1
1.219E+G3
1.196£*03
1.0lGt»03
S.379E»C2
5.379E»02
5'.383£*02
5.383E>02
T
1.274E*03
1.2'i5E»03
1.1C1E*03
9.699E*o2
2.11&E*C3
1.3695*03
9.3o*E*02
2
1.7UE*o3
1.725E*03
1.823E*n3
2.l87E*o3
*.369E«n3
2.109E»03
1.827£»n3
3
1.853£»03 1.915E*03. l.a5lE«63 2.007E*03 2.051II»03 2.093c*03 2«131E*D3
1.896E*03 1.96aE*03. 2.TD5£*03 2.90iK03 2,12?E*03 2,164E*B3
2.2D5£»03 2.32*E*03- 2.^$lt*o3 2.3a6E»Q3 2.405£|»03 2,£»03 2*95o£»03
*• 5 S< 7 89 10 .
7
6
5
4
3
2
1
£,1£4E*03
2.192£«03
2.416£*C3
3.015E«03
3.389E«03
2.989t*C3
2.^40£*03
11
2.1695*03
2.212E*03
2.3P1E»03
2.881E»03
3.4o5t*03
2.969E*03
2.91)E«03
12
2.2o3E*n3 2.235?t»03
2.221E*fi3
2.362E*n3
2.779E*n3
3.39*E*n3
2.92SE*n3
2.870E*n3
13
2.222£»03
2.359E»03
2.7b&£*03
3.3?2E»03
2.92*1*03
2.8&8E4Q3
4>v
ro
!*•
v?»
VEXIT
1.00C£*02 . VS= 1
a; 3.753E»02'
.OOOE»02
FT/SEC» REYNOLDS ^0' BASED ON MEAN MASS F.P* ^ATtl AVO MAX DIA * 115,9
TABLE 2.
PREDICTED
TEMPERATURES FOR ETHYLENE TEST CASE USING BATTELLE HARMONIC -EXPANSION METHOD
-------
NUMBER OF ITERATION
THE DISTRIBUTION Of TEMPERATURE • A
T
6
5
4
3
5.379E*;?
T.719E*p3 1.862E*03 T.915E*03
1.733E*63 1.902E*03 1,961E»03
T.844E»03 2.222E*03 ?.336E*03
?.ai5E*o3 3.266E*03 3.609E*03
4.4?3E*03 3.894E«03 3.609£*03
1.93lE*03 2.473E«03 !?.793E*03
1.955E*03
2.66lE*03
2.366E+Q3
3.523E*03
2.94BE*03
I.99SE*S3
?.038E*n3
?.385E»n3
2.633E*03
2.073E+03
3,399E*C3
3.0*2?:*03
•2t935E*03
2.106E+03
2.405E*63
2,Jo2E*03
2.136£*03
3.5t)6E«03
3.499E*03
3.647E*03
2.988E*f>3
10
7
6
5
4
3
2
1
2.131E*-3
2.3AOE«(>3
2.16RE'63
2»187E»03
2.8n6E*03
3«5oTE*03
2.993E*03
2.929E»03
2.170E*03
2.335E*13
2.792F*03
3.499E*03
2.9295*03
12
13
14
1.55?E*?2 vSs l.556E»02 FT/SECt
= 3.»31E*02
REYNOLDS NO BASER ON MEAN MASS FLOW PATE AND MAX DIfl = 115.8
TABIE 3. PREDICTED TEMPERATURES FOR ETHY1ENE TEST CASE USING SPALDING 4-FLUX METHOD WITH a = Of
-------
OF JFESATIOW =
THE DISTRIBUTION OF.
53
- /UI»J,IE*I>
7
6
5
4
3
2
1
1.207E*C3
1.184E*C3
i.u3
5'.379t*C2
5'.379E*02
5'.388E*02
5'.38SE*tZ
1.
1.2fc3E*C3
1.245E+C3
l.099E*03
9. 7S"E*02
2.1?2£»u3
l.b16E*03
1.095E+03
2
1.7l3E*n3
1.72&E*n3
I.b33£*fi3
2.l8SE+ii3
4-.279E»?3
2.551E»n3
2.335E*03
3
1.857E»03
1.B?6E*03
2.208E»03
3. 216^*03
3*756^*03
2lb37£i»03
2.792E*03
V-
1.91*E*03.
1.96QE*03
2.325E+03.
3.53lE*03
3.53*E*03
3.053£*C3.
2.992E*03.
5
l.Q59E*63
2.?o4E+o3
2.T&3E*63
3.=,21E*03
3.47&E»63
3.flBE*03
3.r7»E*63
5'
2.0o3E*03
2.0*5E»63
2.3B8E»63
3.439E»03
3.4&7E+03
3.146E»03
3.111E»63
7
2»0*5-|»03
2,o85£i*o3
2,406-|»03
3.342£l*03
3,460^1^03
3, l52r|>03
3.113^03
8
2.08*£*03
2.121E«03
2»4l3£i»03
3.245E*03
3.457E+03
3,144£*03
3,10»£-»03
9
2.113E*03
2.132£*D3
2.422£«t>3
3.1V7E03
3.450E*03
3.128£«03
3.096£*33
10
7
6
5
4
3
2
1
2.167E*o3 2.175E*ft3 2.17fa£*03
2.177t*c3 2.192E*[j3
2.t59SE*l;3
3.445£*U3 3.*36E»C3
3.1i"ie.*t3 S.'-SSE+uS
3.059t»t.3 3.C'".8E*03
2.35U«03 2.3V7£*03
2.7?1E*03 2.77S£*03
3.400E*n3 3»395£*03
2.997E+S3 i
2.9'+2£-»03
11
12
13
tf?= 1.000E*C2 VSs
VEXIT = 3t3l5t»i/a
1.0001*02 FT/SEC*
NO. 3ASED OM MEAN MASS f .0* 3ATEI ASO MAX 3d A =: 115.8.
TABIE 4. PREDICTED TEMPERATURES FOR ETHYIENE TEST CASE USING SPALDING 4-FLUX METHOD WITH a = 2ot
-------
45
case. In order to examine the differences it is helpful to study several
plots at constant r or z.
Figure 8 shows the predicted temperature distribution at 4-1/2
inches from the inlet. At this location the harmonic expansion matches the
zone result quite well, especially at the peak. The two Spalding results
bracket the zone result.
Figure 9 shows the predicted temperature distributions on
the centerline of the combustion chamber. At this location, the harmonic-
expansion result is consistently low, whereas the 4-flux method with a = 2ot
agrees best with the zone result.
Figure 10 shows the predicted temperature distributions at a con-
stant radius of 5/8 inch. Here the harmonic expansion result is in reason-
able agreement although somewhat low. This result is quite typical of most
of the combustion chamber.
In summary, for this standard test case for the combustion of
ethylene, the temperatures predicted by the Hottel zone method nearly always
falls between 4-flux,a = a,and a « 2a results but generally closer to the
a » 2a case. The harmonic-expansion results are consistently low with the
greatest difference on the centerline, improving to excellent agreement at
the outer wall. On the other hand the radiation heat fluxes to the walls
predicted by the harmonic expansion method are everywhere in good agreement
with zone-method wall-heat fluxes. However, the 4-flux methods do not pre-
dict the wall-heat fluxes very well for this case, especially at several
locations on the end walls. Therefore, the harmonic-expansion method ap-
pears to be a useful approximation, although some further checking and pos-
sible refinements to the method to improve temperature prediction should
_be investigated.
-------
46
4600
4200
3800
3400
0)
o 3000
&
2600
2200
1800
1400
0
ETHYLENE TEST CASE]
Legend
Hottel Zone Method
1 Battelle Harmonic Expansion
— - ——— Spalding 4-Flux a = a
Spalding 4-Flux a = 2a
0.2 0.4 0.6 0.8
Radius, ft
i.O
1.2
1.4
FIGURE 8. PREDICTED TEMPERATURE DISTRIBUTION AT 4-1/2 INCHES
-------
47
3300
3000
IETHYLENE TEST CASE
Hottel Zone Method
Battelle Harmonic Expansion
Spalding 4-Flux a = a
Spalding 4-Flux a = 2a
1800
1600
12 18 24
Distance From Inlet, inches
FIGURE 9. PREDICTED TEMPERATURE DISTRIBUTIONS ALONG CENTERLINE
-------
48
3800
3600
3400
3200
o:
- 3000
o
-------
49
REFERENCES
1. Hottel, H. C. and Sarofim, A. F., Radiative Transfer. McGraw-Hill,
New York, 1967.
2. Hottel, H. C. and Sarofim, A. F., "The Effect of Gas Flow Patterns on
Radiation Transfer in Cylindrical Furnaces", International Journal
of Heat and Mass Transfer, Vol. 8, 1965, pp. 1153-1169.
3. Johnson, T. R., "Application of the Zone Method of Analysis to the
Calculation of Heat Transfer from Luminous Flames", Ph.D. Thesis,
Sheffield University, England, 1971.
4. Johnson, T. R. and Bee'r, J. M., "The Zone Method of Analysis of Radiant
Heat Transfer: A Model for Luminous Radiation", Paper No. 4
Presented at the 4th Flames and Industry Symposium, London,
September, 1972.
5. Whitacre, G. R., "An Analytical Approach to Thermal Design of a Low-
Emission Combustor for Rankine-Cycle Automobile Engines", Report to
National Air Pollution Control Administration, November 13, 1970,
Contract EHS-70-117.
6. Fitzgerald, F. and Sheridan, A. T., "Prediction of Temperature and
Heat Transfer in Gas-Fired Pusher-Reheating Furnaces", Paper No. 12
Presented at the 4th Flames and Industry Symposium, London, Septem-
ber, 1972.
7. Gosman, A. D., Pun, W. M., Runchal, A. K., Spalding, D. B., and Wolfshtein
M., Heat and Mass Transfer in Recirculating Flows, Academic Press,
London, 1969.
8. Siddall, R. G., "Flux Methods for the Analysis of Radiant Heat Transfer",
Paper No. 16 Presented at the 4th Flames and Industry Symposium,
London, September, 1972.
9. Gosman, A. D., and Lockwood, F., "Incorporation of a Flux Model for
Radiation into a Finite Difference Procedure for Furnace Calculations",
Fourteenth Symposium (International) on Combustion (1972) Penn State
University, Published by The Combustion Institute, 1973, pp.661-671.
10. Lowes, T. M., Bartelds, H., Heap, M. P., Michelefelder, S., and Pai, B. R.,
"Prediction of Radiant Heat Flux Distribution", International Flame
Research Foundation Document No. G02/a/26, Paper presented at the 1973
International Seminar Heat Transfer from Flames, Trogir, Yugoslavis,
August, 1973.
11. Chandrasekhar, S., Radiative Transfer. Dover Publications, Inc.
New York, 1960.
-------
50
REFERENCES
(Continued)
12. Chu, C. M., and Churchill, S. W., "Numerical Solution of Problems in
Multiple Scattering of Electromagnetic Radiation", J. Physical
Chemistry, Vol. 59, 1955, pp 855-863.
13. Davison, B., Transport Theory of Neutrons. Oxford University Press,
1952.
14. Giovanelli, R. G., "Radiative Transfer Non-Uniform Media", Australian
Journal of Physics, VoL 12, 1959, pp 164-170.
15. Pun, W. M., and Spalding, D. B., "A Procedure for Predicting the
Velocity and Temperature Distribution in a Confined, Steady, Trubu-
lent, Gaseous, Diffusion Flame", Imperial College, Mech. Eng. Dept,
Rept. SF/TN/11, 1967.
16. Pai, B. R., and Lowes, T. M., "The Prediction of Flow, Mixing and Heat
Transfer in the IJmuiden Furnace", International Flame Research
Foundation Document No. G02/a/22, October, 1972.
17. Lockwood, F. C., and Spalding D. B., "Prediction of a Turbulent
Reacting Duct Flow with Significant Radiation", Thermodynamic
Sesson, Proceedings Colloques d'Evian de la Societe Francaise, 1971.
18. Erkku, H., Sc. D. Thesis in Chemical Engineering, M.I.T., Cambridge,
Mass., 1959.
-------
51
NOMENCLATURE
A coefficients in spherical harmonic series
A, surface area
a flux absorption coefficient
a.. function multiplying the convection term in Equation (26)
B coefficients in spherical harmonic series
b. functions connected with diffusion terms in Equation (26)
b functions connected with diffusion terms in Equation (26)
C functions connected with diffusion terms in Equation (26)
C specific heat at constant pressure
D diameter of combustion chamber
d source term in Equation (26)
E black-body emissive power
F average radiation flux
f mixture fraction
G G total exchange area between gas zones
h specific enthalpy or heat of reaction
I intensity of radiation
i mass of air which combines with unit mass of fuel
K constant in Equation (22)
L length of combustion chamber
m mass fraction
m mass flow rate
P associated Legendre polynomial
q heat flux
R ratio of water-cooled area to the total surface area
r radial coordinate
-------
52
NOMENCLATURE
(Continued)
S G total exchange area between surface zone and gas zone
S S total exchange area between surface zones
T temperature
U radiation potential
V velocity
v volume
w albedo for scattering
z axial coordinate
Greek Symbols
ot absorption coefficient
F transfer or exchange coefficient
V vector gradient operator
e emissivity
9 angle in coordinate system (see Figure 1)
H extinction coefficient
(j, viscosity
V frequency
p density
<* Stefan-Boltzmann constant
$ conserved dependent variable in Equation (26)
cp angle in coordinate system (see Figure 1)
cp angle in coordinate system (see Figure 1)
•if stream function
n unit vector
U) solid angle, steradians
-------
53
NOMENCLATURE
(Continued)
Subscripts
b black body
eff effective
fu fuel
g gas zone
h heat or energy
I flow network index
1 index on surface zone
ig index on gas zone
J flow network index
j running indices over all gas zones
k running indices over all surface zones
m positive integer
n positive integer
o initial or base
ox oxidizer
r r-direction
s surface zone
st stoichiometric conditions
w wall
z z-direction
tp cp-direction
cp' cp'-direction
-------
54
APPENDIX
LISTING OF PROGRAM ENCLOS
-------
COC-64.00 JTJN.V.»^ 0*P351 OPT=J
t INPUT»OUTP_UT,PUNCH,TAPE60 =
SSOH?<22»22) » A (£2.32) ,SSO (22,22> »UG (22t22) »E (22.22) MM
1* Sfi')K--13(2?»60) ,SGO(22«(>0) . SGOH2 (22» 60 ) HN
_5 _ ,?. niAK^-Ht (f>0 .60 ) , fiBO (frfl «<*O > _ ; _______ '. ___ M N
3. EMS (?,»>. AEt2?).AtuB2*<22>»LTEMP(22>.LPR<22>,LPC(22> M
... ______ 4, BM60)-«Y*K.lb.(>.).*S.i60J ________________________________________ N
NN=i>? *tPb=0.0
... ...... .100 HLAO lQ 20MeS»SX*K,G=»I2»2I9//« CASE NO.«F10.3,5X»K%(J)
3 -AdSOt-'PTlON .F_ACTOR&.FOH.-EACH_GAS_ZONE?_.9X«B =»F9.5//1 IVXlOFil.b) )
12? FO^MAT(//10X*SU«FACE-rO-SUHFACE DIRECT EXCHANGE FACTORS - SS(I.J)/
m^« /> ___^__-
(10 1?3 1 = 1,M SIF(IH.EO.O) GO TO 123
25 .... . KEAO lOdt tSS08i.L!?*IJ_tK=l.»ll.
P*!NT 12, (SSO& SIF.(IGG.LO.O) GO TO 128
_ ..... .Jfc>JQ_127..J=l»N ___ 4IF.MHJ.tU.OJ_ GO. T0..127 _________ ______
FORMAT (//10X«Gab-TO-GAS ni«ECT EXCHANGE FACTOHS - Gfri I* J) /K28** />
___ _35
127 P^INT !£:, «i(jK2K4(K» J> »K = 1» J)
12>- *£Af> lWHt(fc.MSU)..»I = l.t»M ________ ..... »PRINT 1E9 » ( EhS ( I ) » J
129 FOWM&T(//1'»)'»EM1SS1V1TIES FOH EACH SURFACE ZONE*// ?20.1 = 1.,M ___ J ___________ $AOB2=O.OL .. .. .. __________ . ...
40 00 . - -
l) $DM( I, I )=D»« ( I» I >-AOB2R
220 *EtI)=EMSiK tK)<>BK (J)
0.0 $S(J)=0.0 $00 240 K=J»N
$00 243 K=1»M_
50 ?* * V4KOn2=V4Kfie2 *GC«0(K»J) iOO 245 1 = 1»M
-------
J?JiOGHAH__£JiCl_aS CHC_A*fla-JETN_V^.O»Pibl_Q?I=l-..06/20/73 _20^S8.55. PAGE i-
24S
250 V4K {j)=V4Kfi42*t<2
._ DO .310 1*1.1'* , 5DO_J1Q_J=1»_M ._._..
310 MItJ)=DMUtJ) »I=0 *J=0
55 rail r>rTi(fl.M. NNtFPs.t TFMP.TFHR'.D. NPTV.PiVtLPft.LPo
TF(IE^.\-E.O) (rO TO 500
00 330 1 = 1,M *[)G 330 J=1«I. _....._. SK1 = 0
TO 325 K = 1.M *L1 = 0 JIF(K.NE.I) GO TO 321 $K1=1 SCO TO 325
321 10 3?* L=lt(Ji tIF(L.N£.J)_GO_TD 323 SL1=1 ... iGO TO 32*
60 323 A(K-<1,L-L1>=DM(KtL>
CONTINUE iMl=M-l
CALL OFT1 (fl»Ml,NN»t"^S»LTEKP»IEK«»OET»NPIV,PIV»LPH»LPC)
IF(IEKK.NE.O) bO TO bOO $•£ ( I » J) =- (-1 . ) «• ( I » J) *OET/D
65 SS -~ ( I.J)-:-A£Of-.2K _ *IF(IP.EU.1>PUNCH 333» fSSO (I»J)»J=ltM)
___ 340 PrflNJT j"l»(SSO (IbU^FACE-TO-GAS TOTAL EXCHANGE AREAS S&(IfJ)
lX*aM) SUMueTI.ONS. .OVEH J» / VX»ANO SUMMATIONS OVEH !»/)
no "40 1 = 1. M SSJ = 0. 100 435 J=1»N 5>S ( J) =S ( J) +SGO ( 1 1 J)
^35. SJ =SJ *ifiO (I«J) ilF (IP.tG.l)tJiJNCH 333»(SGO (ItJ)»J=l»N)
440 PMINT 3*1«(SGU (I»J) »J=1»N) »SJ SPRINT 481
«5 ____ L^rtTNT 34.lt (StJ)»J=l>H) _______ . __________ _ ________
r>o 470 1 = 1. N S.S(I)=0. tOQ 470 J=ltl $S6«=0.0
K = l «-4
..'G(ItK) SGG=GGO ( I » J)
90 470 filiOUtl) =G(,0{1»J) $PKINT 471
_______ _*_H._^l»«ttT(//wx»GAS-IO-GAS TOTAL EXCHANGE AKEAS ____ GG.(I»J) ____________ «/ .
1 «awo SUMMftflONS OVEW I» /)
no 4^U 1 = 1. N .. SDO 475 J=ltN
S(J)=S(J)* P(50(I»J) SIF(IP.fc-0.1)PUNCH 333» (6GO ( I » J) t J=l .N)
°^]MT 3*lt (bGO(ItJ) tJ=l»N) S^HINT 481
481
_
49] F
-------
_____ EtoCLOS __ CUC-AASLd-ETN- V4.0*.fc3;>I- OPT=1 ___ U6/20X73_2 O.Sri. 55. _ RAG&
(//9X*CAL.CUIATED ..A£(I)*/) ....... ... ... ._. . . .
T41, ( AE ( I) »I = 1»M)
ft(j TO 100 ___________ _____
SOO fKINT SOI , T .J»PlV»MPIV,LPR»LPC
105 ___ soi FORMAT (»i PKOGKAM STDHPFI-) TN
.END.
-------
58
TECHNICAL REPORT DATA
^ (Please read Instructions on the reverse before completing)
'•REPORT NO. 2.
ETPA-650/2-74-011
"•TITLE AND SUBTITLE
Thermal Radiation Modeling for Pollution
Predictions
'• AUTHOR(S)
G. R. Whitacre , R. A. McCann, and A. A. Putnam
8- PERFORMING ORGANIZATION NAME AND ADDRESS
^attelle Memorial Institute, Columbus Laboratories
505 King Avenue
Columbus, Ohio 43201
12. SPONSORING AGENCY NAME AND ADDRESS
EPA, Office of Research and Development
NERC-RTP, Control Systems Laboratory
Research Triangle Park, N. C. 27711
3. RECIPIENT'S ACCESSIOItNO.
5. REPORT DATE
February 1974
6. PERFORMING ORGANIZATION CODE
8. PERFORMING ORGANIZATION REPORT NO.
10. PROGRAM ELEMENT NO.
1AB014, ROAP 21ADG-16
11. CONTRACT/GRANT NO.
Grant No. R- 800842
13. TYPE OF REPORT AND PERIOD COVERED
Final
14. SPONSORING AGENCY CODE
15. SUPPLEMENTARY NOTES
IB. ABSTRACT ^g report gives results of a study of the mathematical modeling of the
combustion process. It indicates that adequate treatment of thermal radiation is
essential and that accurate temperatures are especially important for pollutant
predictions since the chemical kinetics are extremely temperature dependent. It
describes test cases that were run which combined the Hottel Zone Method with a
completely analytical flow and mixing solution: the complexities that arise when this
approach is extended to real systems seem to eliminate it as a working tool. The
results were used to check the accuracy of several four-flux radiation approximations
Which are simpler than the Zone Method: although temperature predictions were
reasonable, the four-flux approximation failed in predicting the wall radiation fluxes
at many locations. The report also describes tests of an alternate approximate
method which offers great promise: an expansion of the radiation transport equation
by spherical harmonics. When the first three terms of the associated Legendre
polynomial series are retained, a single second-order differential equation results.
This equation applies over the entire field in contrast to the four-flux approaches
which use separate equations for each direction.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b. IDENTIFIERS/OPEN ENDED TERMS
c. COSATI Field/Group
Air Pollution Combustion
Thermal Radiation Chambers
Mathematical Models Nitrogen Oxides
Predictions Combustion Control
Heat Transfer
Aerodynamics
20M
13B
'8. DISTRIBUTION STATEMENT
19. SECURITY CLASS (ThisReport)
Unlimited
21. NO. OF PAGES
63
20. SECURITY CLASS (Thispage)
22. PRICE
EPA Form 2220-1 (9-73)
-------
INSTRUCTIONS
1. REPORT NUMBER
Insert the EPA report number as it appears on the cover of the publication.
2. LEAVE BLANK
3. RECIPIENTS ACCESSION NUMBER
Reserved for use by each report recipient.
Title should indicate clearly and briefly the subject coverage of the report, and be displayed prominently. Set subtitle, if used, in smaller
type or otherwise subordinate it to main title. When a report is prepared in more than one volume, repeat the primary title, add volume
number and include subtitle for the specific title.
Each report shall carry a date indicating at least month and year. Indicate the basis on which it was selected (e.g., date of issue, date of
approvcl, date of preparation, etc.).
6. PERFORMING ORGANIZATION CODE
Leave blank.
Give name(s) in conventional order (John R. Doe, J. Robert Doe, etc.). List author's affiliation if it differs from the performing organi-
zation.
8. PERFORMING ORGANIZATION REPORT NUMBER
Insert if performing organization wishes to assign this number.
9. PERFORMING ORGANIZATION NAME AND ADDRESS .
Give name, street, city, state, and ZIP code. List no more than two levels of an organizational nirearchy.
10. PROGRAM ELEMENT NUMBER ,..,.,_,.
Use the program element number under which the report was prepared. Subordinate numbers may be included in parentheses.
11. CONTRACT/GRANT NUMBER
Insert contract or grant number under which report was prepared.
12. SPONSORING AGENCY NAME AND ADDRESS
Include ZIP code.
13. TYPE OF REPORT AND PERIOD COVERED
Indicate interim final, etc., and if applicable, dates covered.
14. SPONSORING AGENCY CODE
Leave blank.
15. SUPPLEMENTARY NOTES . ... _ , .
Enter information not included elsewhere but useful, such as: Prepared in cooperation with, Translation of, Presented at conference of,
To be published in, Supersedes, Supplements, etc.
Include a brief (200 words or less) factual summary of the most significant information contained in the report. If the report contains a
significant bibliography'or literature survey, mention it here.
17. KEY WORDS AND DOCUMENT ANALYSIS .,„.,.
(a) DESCRIPTORS - Select from the Thesaurus of Engineering and Scientific Terms the proper authorized terms that identify the major
concept of the research and are sufficiently specific and precise to be used as index entries for cataloging.
(bj IDENTIFIERS AND OPEN-ENDED TERMS - Use identifiers for project names, code names, equipment designators, etc. Use open-
ended terms written in descriptor form for those subjects for which no descriptor exists.
(c) COS ATI FIELD GROUP - Field and group assignments are to be taken from the 1965 COS ATI Subject Category List. Since the ma-
jority of documents are multidisciplinary in nature, the Primary Field/Group assignment(s) will be specific discipline, area of human
endeavor, or type of physical object. The application(s) will be cross-referenced with secondary Field/Group assignments that will follow
the primary posting(s).
18. DISTRIBUTION STATEMENT . ,
Denote reusability to the public or limitation for reasons other than security for example "Release Unlimited." Cite any availability to
the public, with address and price. /
19.&20. SECURITY CLASSIFICATION
DO NOT submit classified reports to the National Technical Information service.
21. NUMBER OF PAGES
Insert the total number of pages, including this one and unnumbered pages, but extiuae aisiriDuuon usi, 11 any.
22 PRICE
Insert the price set by the National Technical Information Service or the Government Printing Office, if known.
EPA Form 2220-1 (9-73) (Reverie)
------- |