18-20 June, 2003, Budapest, Hungary OSSKI Center (Törley Palace) | 
with Exhibition and Pre-Session on Thermal Energy in Hungarian
"THERMO-BRIDGE"
between East and West for technology transfer and
information exchange
| Scientific Society of
Measurement, Automation and Informatics (MATE)
Branch of Thermal Engineering and Thermogrammetry (TE and TGM)
Méréstechnikai, Automatizálási és
Informatikai Tudományos Egyesület
Hőtechnikai és Termogrammetriai (HT és TGM) Szakosztály
MATE Secretariat: H-1372 Budapest, POB. 451.,Hungary
House of Technology, Budapest
V.,Kossuth Lajos tér 6-8.III.318.
Phone: +361-332-9571, Fax:+361-353-1406
E-mail: mate@mtesz.hu benko@energia.bme.hu |
| | | |
|
49/5 |
Development of the Monte-Carlo method to predict radiative heat transfer within the boilers and furnaces
Prof., Dc.Sc.(Eng.) B. SOROKA, Dipl.Eng. V. ZGURSKYY, Dipl.Eng. K. PYANYKH
Gas Institute, National Academy of Science of Ukraine |
S7IR07 |
|
| |
SummaryComparison of various methods of
radiative heat transfer calculations with the
Monte Carlo procedure (MCP) is carried out in
frame of some principal geometries of the
emitting medium. The detailed conditions of MCP
performance to conform the results with those in
accordance with the zone and flux methods are
discussed for infinite flat layer, boiler
furnace, industrial and process furnaces.
Adequacy of the model of the weighted sum of
grey gases (WSGG) is tested to compute the gas
medium radiation in combination with proposed
MCP version.
Comparison of various methods of radiative heat transfer calculations with the Monte Carlo procedure (MCP) is carried out in frame of some principal geometries of the emitting medium. The detailed conditions of MCP performance to conform the results with those in accordance with the zone and flux methods are discussed for infinite flat layer, boiler furnace, industrial and process furnaces. Adequacy of the model of the weighted sum of grey gases (WSGG) is tested to compute the gas medium radiation in combination with proposed MCP version.
The method of statistical tests or Monte-Carlo procedure has taken a wide spread by furnace equipment calculations and by high-temperature firing units (boilers and furnaces) improvement in frame of solution of following tasks [1,2,3]:
calculation of local and average angular radiation factors in diathermic medium;
definition of generalized and resolving angular factors of radiation, including those accounting anisotropic scattering being used at zone calculation methods of radiative heat exchange;
Calculation of mentioned angular factors corresponds to determination of geometrical factors and the surfaces of a mutual radiative exchange: direct and total – between the computation zones of both surface and volume types within the framework of Hottel’s zone method [4].
In generalizing P. Lybaert’s report [5] Monte-Carlo method is mentioned among other methods of definition of direct surfaces of a mutual exchange, including analytical (being offered for the simple geometries of the system) and alternative procedures – the discrete transfer methods (of the fixed directions of the beams being suitable practically for any system). Unlike the Monte-Carlo method, the simplified procedures cannot be used for calculation of the generalized angular factors and of emissivity factor of radiating volumes.
|
2. Comparison of calculations for a flat radiating layer |
Our earlier calculations have been submitted and generalized in the monograph book [6] basing on some procedures: of multiple (repeated) reflections, of integro-exponential E-functions, the approximate method of effective fluxes – are compared below with the new computations providing an examination of adequacy of Monte-Carlo procedure (MCP).
In our presentations [7,8] the Monte-Carlo procedure was developed for the case of general calculation of the combustion chambers and the boiler furnaces filled in of the emitting and absorbing medium. The last was considered in frame of approximation of non-grey radiation being considered as the weighted sum of grey gases (WSGG). This approach is used within present paper.
Selectivity and non-isothermal state of the gas volume of non-grey gases were taken into account by use of H. Hottel’s model [7,8] and by representation of the integral emissivity factor
 and absorption abilities Ag of real gas (combustion products) by means of weighted sum of radiation of n grey gases each of them being characterized by absorption factor kn not depending on temperature. Contribution of the n-th gas radiation depends on gas Tg or emitter Ts temperature, the corresponding fraction an being submitted to a condition of normalization:
|
 ;
(1) |
|
 ;
(2) |
|
 , (3) |
|
where
 ;
 |
Various approximations for non-grey radiation are collected in the literature where the number n makes value between 3 up to 6 components [5] and even up to 19 [9]. The an and kn values for the combustion products of natural gas ( kPa,
 kPa) were taken accordingly Taylor and Foster data [10] by number n = 4 in our calculations have been performed earlier and within the framework of the present presentation.
We’ll call the n-th grey gas contribution in case of non-grey gas as n-th “quasi-level” by analogy with the structure of the selectively emitting gas. Using Monte-Carlo procedure accordingly [11] the different random numbers must correspond to each of "quasi-levels" of non-grey radiation. This approach simulates subjecting of the selective radiation characteristics to stochastic computation procedures.
By our present performance the approach is changed as contrasted to the paper [11]: the complete cycles of statistical trials are fulfilled for each of n grey gases in the same manner. The similar procedure was recommended in the manual [12] for separate spectral intervals of selectively radiating gas.
The difference of the method, originated by us, from those in the book [12] is, that we have fulfilled the straight calculation of resulting heat fluxes by radiative transfer between all volume and surface zones while in [12] the spectral resolvent angular coefficients (SRAC) were computed basing upon optical and geometrical performance between volume and surface zones inherent to the enclosure. The direct interchange areas represent analogue of SRAC by use of H. Hottel's zone method grounded upon WSGG model of radiating gas. The computation area represents two infinite parallel plates: the radiator (emitter) – an adiabatic brickwork of unknown temperature
 and the receiver of given temperature
 - enclosing the flat gas layer of
 thickness, the last being filled in of the completed combustion products of natural gas. Calculation of the resulting specific heat flux
 and refractory temperature
makes the problem to be solved.
Table 1. Temperature pattern of the
gas layer sections across the gas flat
volume within the furnace by various heat
transfer conditions: uniform distribution,
direct, indirect, direct + indirect
heating (in accordance with classification
[6]).
|
MCP computations were carried out by
substitution of infinite flat layer for finite
sizes one surrounded by mirror reflected side
surfaces. The results of comparison of the
resulting heat fluxes qA for
various temperature profiles across the layer
are demonstrated by Fig.1 and are fairly
coincided between two computation procedures
under consideration.
|
 |
|
Fig.1. Dependence of
resulting heat flux density
absorbed by receiving surface of flat
layer of natural gas combustion products
in dependence on reflectivity of
refractory surface ρR. |
The number of the curves
corresponds to the variant of temperature
pattern across the flat layer of natural gas’
combustion products (Table 1)
a – our calculations by Monte-Carlo
procedure. Relative energy remainder ξ = 0.002;
b – our data being computed by multiple
reflections procedure and adopted from [6].
| 3. Comparison of calculations for
the boilers’ and furnaces’ combustion
chambers
|
Reliability of performance of
the Monte-Carlo procedure with reference to heat
exchange calculations within the boilers and
furnaces is of great practical interest. The
following typical designs of the furnaces have
been chosen to testify MCP
- the parallelepipeds have been extended in
height, schematising the furnaces of the steam
generating units and of the tubular process
furnaces (of oil refining and of chemical
branches of industry),
- the parallelepipeds have been extended in
horizontal planes, schematising the reheating
furnaces.
2.1. The following
geometrical ratios are accepted for the units
calculated by H. Kremer:
 ,
 ,
the range of the size b variation – up to 15
m, the step of change – 5 m – for the
boiler or the process furnace (Fig.2);
 , ,
range of b variation – up to 10 m,
step of the size change – 2 m – for the
industrial reheating furnace.
|
 |
|
Fig. 2. Resulting radiative heat flux
density q absorbed by the receiving
surface within the boiler or process
furnace in dependence on the combustion
chamber width b.
|
-
Solid lines – our calculations
by Monte-Carlo procedure, dotted lines – by
H. Kremer [13, Fig.15]
Number of computation zones Ncz
=1.
Zone method computations parameter – furnace
gas temperature Tg, ˚C, is
indicated near the curves.
O2 content in oxidant [O2],%:
a– 30%; b– 50%
The following initial
conditions are assumed by calculations: natural
gas (methane) is considered as a fuel, air or
air enriched with oxygen (О2
fraction in an oxidiser )
- as an oxidizer, supplied for combustion by
oxidant excess factor till 1.1 (the residual
О2 contents in dry combustion
products of 2 %). A range of a variation
[O2] in an oxidizer
= 0.21;
0.3; 0.5; 1.0, that corresponds to the sum of
partial pressures of СО2 and of
Н2О in dry combustion products
 ,
resulted in Table 2.
Table 2. Partial pressure of the
emitting components in combustion products
of natural gas
|
 |
H. Kremer has executed
alternative furnace calculations by means of
zone method approach with acceptance of the
optical characteristics of the solid surfaces
and of gaseous combustion products, being
corresponded to grey radiation. Combustion
products temperature Tg = 800…1800 °
C is settled as uniform value by furnace
space. The step of variation of the last by
computations was stated
of 200 ° C.
It was supposed, that in case of the boiler or
the process furnace all walls are arranged with
the tube heat receivers of the set temperature
 °
C.
In case of the industrial furnace
heat receiver surface of the assigned
temperature
 °
C occupies the area
in the bottom base of the furnace presented as a
parallelepiped. The step of 400 ° C of
variation of receiver temperature was assumed by
computations. Other five side surfaces and the
roof of the furnace are covered of refractory
bricks. Specific heat loss through the walls (density
of heat flux) was assumed
of 5 kW/m2. |
The MCP was realized by
account of "history" of 104 beams
being emitted by each of surface zones and of 106
beams let off by volume zones on each of n
quasi-levels of radiation (model of the weighted
sum of grey gases) by fulfilment of the boiler (or
process furnace) calculations in case of the
furnace space consideration similar to 1-zone
system. Number of the beams have been taken into
account makes: 103 – for each of
surface zones, 105 – for each of
volume zones - in case of division of the boiler
(furnace) space into 54 volume zones.
Number of the beams taken
into account by emission on each of n
quasi-levels for any of volume and each of
surface zones made 80× 103 while
using of MCP by calculation of the industrial
reheating furnace.
Comparison of the
calculations results obtained by means of
Monte-Carlo method with those computed by zone
method for two objects – the boiler (the process
furnace) (Fig.2) and the industrial reheating
furnace – shows an increase in divergences of
the resulting heat fluxes density q as
the characteristic size of the unit or О2
fraction in an oxidizer is rising. It is
possible to assume, that the certain explanation
of such character of a divergence of the results
gives an account of the limited range of optical
thickness of the gas volume, at which
and
factors were chosen by generalization of the
data on emissivity factor
 .
Really, area of
presentation (and 
- the same) falls at parameter
values according to [10]. Meanwhile the
calculations were carried out for the furnaces
of the characteristic sizes
 .
By considered conditions of
combustion of natural gas with nitrogen-oxygen
oxidizers, including an air ([O2]=20.95%)
till pure [O2] (100%) at l
up to 1.1,
under consideration varies within the range of
values 0.065…3.255 atm·m. Thus the
parameter
in our computations exceeds the top limit of the
parameter values being met to Taylor and
Foster’s recommendations [10] on use of the
and
values providing the real
and
values. As a result the WSGG model doesn’t fit
to maximum sizes of the furnace under
consideration.
2.2.
E. Scholand has generalized the results of
calculation of the tube cracking furnace (Fig.3)
in the paper [14], executed by three methods: a
– accordingly to the flux model, b – by Hottel’s
zone method and, c – by means of his own
Scholand’s simplified one-dimensional zone model.
The last assumes division of the furnace under
consideration by separate volume zones. These
volume cells are submitted to computation by the
perfect stirring reactor (PSR) model. All of the
zones are calculated in combination by solution
of the matrix formed by the system of the linear
equations concerning effective fluxes from the
surfaces enclosing the volume elements.
Reciprocal heat and mass transfer between the
vertical zones hasn’t taken into account.
The initial data taken for
calculation: fuel – the gas mixture of 70 %
СН4, of 25 % Н2,
of 5 % N2, being burnt
with an oxidizer – an air of excess air factor .
The tubes forming the reception surface are of
an emissivity factor
and of temperature 1000 K. All other surfaces
are considered as adiabatic ones. The
temperature profile within the vertical planes
of the furnace gas is assigned of six isothermic
sections of the same height, while the
horizontal cross-sections are completely met to
the isothermal pattern of the combustions
products. By the zone method calculations the
furnace space was divided by height resulting in
formation of six horizontal zones, of the size
 .
Another type of the zones formation supposes an
additional division for calculations of each of
the mentioned horizontal zones by three ones
each of the last being of
in the sizes. Half of furnace space divided by
an axial vertical plane
was examined due to symmetry of the furnace
space (Fig.3).
The data submitted in [21]
testify to rather good conformity of the
calculation results basing on zone method and on
simplified one-dimensional method. Unlike
mentioned methods the calculations by flux
method have caused the strong deviation from
their results.
Processing of the results
submitted in [14] provides at the same time more
than three times difference in the values
obtained in separate sections of the furnace by
calculations by flux method from those taken by
zone method. The most reliable from the compared
approaches is seemed apparently the H. Hottel’s
zone method.
Comparison of the mentioned
results of calculation obtained in the paper
[14] with our calculations, which have been
carried out by Monte-Carlo method, is submitted
in Fig.4. Division of the furnace space into the
zones of two sizes mentioned above, as well as
problem decision in two statements has been used:
for grey and non-grey gas radiation. In the
first case the single value of medium absorption
factor
has been accepted for all volume of combustion
products. In the second case the model of gas
radiation as weighted sum of grey gases was used
with acceptance of empirical factors
and
- accordingly data by P.B. Taylor, P.J. Foster
[10].
|
 |
|
Fig.
4. Profile of resulting heat flux
density qt by furnace
height h / Comparison of original
predicted data with former calculations. |
Number of computation zones/
zone sizes (m3):
a – 6 / 1×1×3; b– 18 / 1×1×1
1 – Monte-Carlo procedure, grey gas radiation,
K=0.2 m-1;
2 – Monte-Carlo procedure, non-grey gas
radiation (WSGG)
pCO2 + pH2O
= 0.22 atm, pCO2 : pH2O = 2.36
: 1;
3 – Zone method by H. Hottel, data summarized by
E. Scholand and taken from [14].
The dynamics of power
decrease by movement of 240000 beams by each
zone at each of n radiation quasi-levels
(WSGG model) was investigated by performance of
Monte-Carlo procedure.
For recalculation of the computed values of heat
fluxes directed to conditional reception surface
in respect to true tube surface
the next equations were used:
|
 ,
(5) |
where
 ,
- heat flux density at the surface
of the tubes and within conditional plane
-
correspondingly.
It is possible to draw a
conclusion from Fig.4 consideration, on good
concurrence of the results of our calculations
by Monte-Carlo procedure both for grey and
non-grey gas emissivity models to those obtained
by zone method, the results of our calculations
for case of non-grey gas being in better
response to the data resulted in the E. Scholand’s
paper for “grey” gas model of absorption
coefficient
.
The maximum divergence of the results by
Monte-Carlo model and by zone method use makes
14 % (min – 0 %).
Let's note in conclusion, that additional
division of the isothermal zones does not result
in change of the computation values of the heat
fluxes (compare the data on qt
in Fig.4,a,b), that is connected apparently the
advantages of developed generator of the random
numbers by stochastic procedures.
1. It’s stated that
Monte-Carlo method represents an operative
computation means to predict both integral (and
averaged within the furnace) and local thermal
condition of the boiler furnace as well as of
the separate facilities of the last with account
of furnace design features and specifies of
operational mode: by substitution the fuel and/or
oxidant composition, of by input of
supplementary components – recirculating
combustion products, water vapour etc.
2. The Monte-Carlo
method and procedures are developed for the
furnace chamber, being divided into the
computation zones. Selectivity of combustion
products radiation within the furnace space is
accounted by means of use of a weighted sum of
emissivities of grey gases as a gas emissivity
model, the weighting coefficients being
dependent on temperature of emitting zone.
It was succeeded with this to
investigate contribution into the resulting heat
flux of each of heat transfer components: of the
initial radiation of each zone (gas volume,
surface) and of numerical reflections of the
rays from the surfaces enclosing the furnace
chamber or space.
3. Validation of the
Monte-Carlo heat exchange computation method had
been performed
- by use of the models of infinite flat
layer – grounded on comparison with the
results obtained by the method of multiple
reflections. The method of presentation of the
finite sizes flat layer surrounded by the side
walls of mirror reflection to simulate of
infinite flat layer by Monte-Carlo procedure
calculation is developed.
- by computation of the boiler furnace and
the oven chamber facilities – grounded on
comparison with the zone method calculations.
The results of calculations
of the boiler furnaces and of the industrial
furnaces obtained by Monte-Carlo method were
compared with the computations presented in some
references, by the results of various existing
methods: flux, zone 3D (three-dimensional) and
simplified 1D (one-dimensional).
H. Hottel’s zone method is
the most perfect procedure of the existing
methods of calculation of industrial furnaces of
arbitrary geometry with account of gas radiation
selectivity.
Monte-Carlo approach with
respective procedures of its realization being
proposed in frame of the paper provides
qualitative and quantitive similarity to the
results of 3D zone method computations but has
some advantages from the standpoint of the
method simplicity and of opportunity to
determine the detailed heat pattern within the
boiler furnace.
[1.] Zhuravlev Yu.A.
Radiative heat exchange within firing
engineering plants. – Krasnoyarsk: Krasnoyarsk
Univ. Publ. House. – 1983. – 256 pp. (In Russian).
[2.] Arutyunov V.A., Bukhmirov V.V.,
Krupennikov S.A. Mathematical modeling of
thermal operation of industrial furnaces // Sc.
editor V.A. Arutyunov. – Moscow: “Metalurgiya”
Publ. House. – 1983. – 256 pp. (In Russian).
[3.] Lisienko V.G., Volkov V.V., Malikov
Ju.K. Fuel use improvement and heat exchange
control within metallurgical furnaces. – Moscow:
“Metalurgiya” Publ. House. – 1988. – 231pp. (In
Russian).
[4.] Boineau P., Copin C., Aquile F. Heat
transfer modelling using advanced zone model
based on a CFD code // 6th European
Conference on Industrial Furnaces and Boilers:
Preprints. – Estoril – Portugal. 02-05 April
2002. vol.III: Modelling of Furnaces and
Combustion Systems. – 10 pp.
[5.] Lybaert P. From CFD to zonal
furnace models: an overview of head treatment
furnace modeling // 6th European
Conference on Industrial Furnaces and Boilers:
Preprints – Estoril – Portugal. 02-05 April
2002. – 14 pp.
[6. Soroka B.S. Intensification of thermal
processes within fuel furnaces. – Kiev: “Naukova
dumka” Publ. House. – 1993. – 416 pp. (In
Russian).
[7.] B. Soroka, К. Pyanykh, V. Zgursky, M.
Khinkis, H. Abbasi, J. Rabovitser. Mathematical
modeling of low-emission combustion processes
basing upon Monte-Carlo procedures // Preprints
of 5th European conference on
industrial furnaces and boilers (11-14 April
2000). – Espinho-Porto (Portugal) – Vol. II. –
p. 12.
[8.] Soroka B. et al. Heat - and Mass
Transfer Processes in Low-Emission by NOx
Boiler Furnaces under Condition of Two-Stage
Gasburning and Reaction Products Recirculation.
// Proceedings of Heat and Mass Transfer IV
Minsk Intern. Forum , 22 – 26 May, 2000, - V.10.
–P.436-445. (in Russian).
[9.] Strohle J., Coelho P.J., Schnell U.,
Hein K. A non-grey radiation model for the
simulation of coal-fired furnaces // 6th
European Conference on Industrial Furnaces and
Boilers: Preprints. – Estoril – Portugal. 02-05
April 2002. vol.III: Modelling of Furnaces and
Combustion Systems. – 10 pp.
[10.] Taylor P.B., Foster P.J. The total
emissivities of luminous and non-luminous flames//
Int. Journ. Heat & Mass transfer. – 1974, 1974,
vol.17. - №12 - P.1591-1605.
[11.] Steward F.R., Cannon P. The
calculation of radiative heat flux in a
cylindrical furnace using the Monte-Carlo method
/ Int. Journal of Heat & Mass Transfer. - 1971.
- Vol.14, No.2. - P.245-261.
[12.] Lisienko V.G., Volkov V.V., Goncharov
A.L. Mathematical modeling of heat exchange
within furnaces and assemblies. Kiev: “Naukova
dumka” Publ. House. – 1984. – 230 pp. (In
Russian).
[13.] Kremer H. Moglichkeiten und Grenzen
der Strahlungeswarmeubertragung in industriellen
Gasfeurungen // GasWarme Int. – 1994. – J.43,
H.10. – S.482-495.
[14.] Scholand E. Modern procedures for the
calculation of radiant heat transfer in
direct-fired tubed furnaces // Int. Chemical
Engineering, 1983. - Vol. 23, № 4. - P.600-610.
| | |
Contact details:
Prof., Dc.Sc.(Eng.) B. SOROKA
Gas Institute, National Academy of Science of
Ukraine
Kiev,
39 Degtyarivska St., Ukraine, 03113,
E-mail: soroka@elan-ua.net |
| |
|
| 18-20 June, 2003, Budapest,Hungary |
OSSKI
Center (Törley Palace)
"Fodor
József" National Center of Public Heath
"Frédéric Joliot-Curie"
National Research Institute for
Radiobiology and Radiohygiene. (OKK-OSSKI)
www.osski.hu |
|
|
|
profile across the gas layer is submitted in the
Table 1. The gas layer is divided upon N = 20
discrete sections of the same thickness.
Numbering of the sections (i = 1) is
beginning from the layer, being the nearest to
the refractory wall. The layer contacting with
the receiver wall is numbering as i = N = 20).
Various temperature profiles were compared from
the efficiency standpoint (estimation of the
maximal
(4)
and convective
components in frame of additivity ability:
,
and
.
The surfaces A and R, enclosing
the gas layer, were assumed as grey bodies by
performance of calculation. Other conditions and
limitations are similar to those assumed in book
[6].
