Mathematical modeling of toxic pollutants dispersion from large tank res

and assessment of acute effects for re ghters

N.C. Markatos
*
, C. Christolis, C. Argyropoulos
Computational Fluid Dynamics Unit, School of Chemical Engineering, National Technical University of Athens, 9, Iroon Polytechniou Str., Zografou Campus, GR - 157 80 Athens, Greece
a r t i c l e i n f o
Article history:
Received 5 December 2008
Aceepted by 23 March 2009
Available online 8 May 2009
Keywords:
Tank re
Pollutant dispersion
Smoke plume
Turbulence
CFD
Risk zones
a b s t r a c t
The paper presents a Computational Fluid Dynamics (CFD) approach to major hazard studies, by apply-
ing a nite-domain technique to predict the dispersion of combustion products (CO, SO
2
, smoke, Volatile
Organic Compounds, VOC, Polycyclic Aromatic Hydrocarbons, PAH, etc.) from fuel-tank res, as well as
the toxic plume rise.
Furthermore, a methodology is presented for the identication of risk zones for the rst respondents
(re ghters). Numerical simulations were performed with the use of the SIMPLEST algorithm, a sec-
ond-order accurate MUSCL scheme with deferred correction for the convective terms in the momentum
equations, CUPID for the convective terms in the scalar conservation equations, and a modied, RNG k-e,
two-equation turbulence model. Radiation was computed by a discrete transfer model. Various scenarios
were examined and the results are presented and discussed.
2009 Elsevier Ltd. All rights reserved.
1. Introduction
Many analysts, such as Beck [1], claim that we live in a society
with many risks, which spread in all facets of our life. In the last
40 years, intense industrialization in the developed countries and
the constant economic growth have brought new technologies,
new products, new opportunities, better conditions of life, fast
transportation and modern communications. However, there are
also some indirect disadvantages from these improvements,
including industrial accidents. Every year hundreds of work acci-
dents, some with human losses and high cost for the economy
and the environment, take place in the world.
Petrochemical industries normally use large storage tanks,
which contain large amounts of ammable and hazardous chemi-
cals. Hence, the occurrence of a tank accident is possible and it usu-
ally leads to re and explosions. Experience has shown that the
continuous dispersion of toxic pollutants from large tank res,
such as smoke, SO
2
and CO is responsible for potential environ-
mental and health problems.
Industrial accidents such as that of February 24, 1986, in Thes-
saloniki, Greece, and the relatively recent large industrial accident
of December 11, 2005 at Oil Storage Depots (Bunceeld, Hertford-
shire, England) demonstrated unanticipated consequences. The
latter one produced a massive re which engulfed over 20 large
fuel-storage tanks. The continuous generation of smoke from these
res presents a potential environmental and health problem that is
difcult to assess. In order to try and manage it, it is important to
be able to estimate the concentration of the re-plume pollutants
over a wide range of conditions.
Hence, the need for new techniques and methodologies, in
order to predict the consequences and to reduce the frequency of
industrial accidents, is imperative. Reneries and chemical plants
comprise the major part of storage-tank accidents, because of the
large amounts of chemicals and fuels which are stored, mostly in
atmospheric external oating-roof and cone-tap tanks.
According to the research of Chang and Lin [2], among the 242
tank accidents that have been reported in the last 40 years, the re
incident is the most frequent type of accident, caused by lightning
and maintenance errors. Thus, the study of tank res is of great
importance, since, unfortunately, there is always a chance for a
new accident. Under these circumstances, the re plume produced
from the burning fuel generates large amounts of toxic pollutants
which are transported by the wind over a large area. The transpor-
tation of pollutants, such as sulfur dioxide (SO
2
), carbon monoxide
(CO), hydrogen sulde (H
2
S), as well as soot and particulates, have
adverse effects on human health and the ecosystems.
Intensive research has been performed by many scientists and
engineers for the investigation and explanation of the buoyant
re-plume dispersion. Ghoniem et al. [3] and Ghoniem and Zhang
[4] have devised Lagrangian numerical techniques for the solution
of governing equations, based on the extension of the vortex meth-
od to variable-density ows. Of a different numerical approach is
the work of McGrattan et al. [5], with the use of Large Eddy Simu-
lation (LES), for the simulation of turbulent smoke plume from
large oil res. Finally, Miles et al. [6], Christolis et al. [7], Sideris
0017-9310/$ - see front matter 2009 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijheatmasstransfer.2009.03.039
* Corresponding author.
E-mail address: [email protected] (N.C. Markatos).
International Journal of Heat and Mass Transfer 52 (2009) 40214030
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer
j our nal homepage: www. el sevi er . com/ l ocat e/ i j hmt
et al. [8] and Kefalas et al. [9] determined the lift off condition of
buoyant re plumes released from warehouses and buildings,
using CFD techniques and the dimensionless buoyant-ux number.
The purpose of the present effort is to estimate the dispersion of
combustion products and the consequences to the environment
from large hydrocarbon-tank res, as well as the height of the toxic
plume (plume rise), and it is a continuation of the work by
Argyropoulos et al. [10]. Furthermore, the risk zones for the rst
respondents, the re ghters, are identied, by comparing the
ground-level concentration of the pollutants with the safety limits.
The mathematical model used in the present work consists of the
full RANS (Reynolds-Averaged NavierStokes) differential equa-
tions that describe turbulent ow and heat and mass transport
[11]. A MUSCL scheme [12] with deferred correction for the con-
vective terms in the momentum equations, and the CUPID scheme
[13] for the convective terms in the scalar conservation equations,
were employed to discretize the equations and the solution was
obtained using the iterative SIMPLEST algorithm of Spalding
[14,15]. The RNG k e model [16] was employed with the buoy-
ancy modications of Markatos and Pericleous [15], for the model-
ing of turbulence. Radiation is accounted for by a discrete transfer
model [17,18].
Parametric analysis is performed for eight scenarios, all for an
adiabatic atmosphere: for two different fuels with heat-release
rates of 1 and 1.3 MW/m
2
, and for four values of wind velocity, 7,
8, 9 and 11 m/s.
2. The physical problem considered
For the numerical simulation of the toxic contaminants disper-
sion and of the plume rise, an external oating-roof tank has been
selected, with dimensions of 85 m diameter and of 20 m height.
The tank is surrounded by bunds of 4 m height and 0.5 m width.
The geometry of the tank with the bund is illustrated in Fig. 1.
This type of tank is characterized by the ability of the roof to rise
and fall on the stored-fuel surface, according to the latter level,
with a maximum storage capacity of 113,432 m
3
. Numerical simu-
lations were performed with the use of CFD techniques for a phys-
ical domain of 30000 m length (Z-axis), 2200 m width (X-axis) and
2500 m height (Y-axis). These dimensions were chosen following a
parametric study, in order to evaluate the extent of the domain
Fig. 1. The tank geometry considered.
Nomenclature
A
f
surface area of pool (m
2
)
C
e1
; C
e2
; C
e3
; C
l
; g
o
; b
1
empirical constants inthe turbulencemodel
g gravitational acceleration (m/s
2
)
G
B
buoyancy source/sink term
h enthalpy (J/kg)
I radiative intensity
IDLH Immediate Dangerous to Life and Health
k turbulence kinetic energy (m
2
/s
2
)
LC
1
Lethal Concentration
_ m
00
mass-loss rate (kg/s)
P pressure (N/m
2
)
P
k
turbulence kinetic-energy generation term
_ q volumetric heat-release rate (J/(m
3
s))
S mean rate of strain
S
ij
deformation tensor
S
h
enthalpy source term (J/(kgs))
SP source term of pollutants
T temperature (K)
T
ref
reference temperature (K)
U
i
mean velocity component in the i-direction (m/s)
U
j
mean velocity component in the j-direction (m/s)
U average wind speed (m/s)
X, Y, Z Cartesian co-ordinates
z
o
surface roughness (m)
Greek symbols
b thermal expansion coefcient (1/K)
d height of atmospheric boundary layer
e turbulence dissipation rate (m
2
/s
3
)
n ratio of turbulence time scale to mean strain-rate scale
u time-average of instantaneous ow variable
u
P
value of ow variable at the central grid node, P
u
I
value of ow variable at the neighbouring grid-node, I
m kinematic laminar viscosity (kg/(m
2
s))
m
t
kinematic turbulence viscosity (kg/(m
2
s))
q air density (kg/m
3
)
r
h
turbulence Prandtl number for enthalpy
r
k
turbulence Prandtl number for kinetic energy of turbu-
lence
r
e
turbulence Prandtl number for dissipation rate of turbu-
lence
r
T
turbulence Prandtl number for temperature
W Van Leer ux limiter
Subscripts
H higher-order approximation
i, j i, j directions
L lower-order approximation
P central grid node
ref reference
t turbulent
W the node West of the central grid-node P
WW the node west of the west-node of P
4022 N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030
after which no major change of the variables was computed close
to the tank. Z-axis is the parallel to the wind direction axis, X-axis
is the vertical to the wind direction on the ground-level axis, while
Y-axis is the vertical axis. The tank is situated 543 m from the ori-
gin of the Z-axis and 1100 m from the X-axis (symmetry).
3. The mathematical model
3.1. Governing equations
The mathematical model is based on the set of elliptic, partial-
differential equations (PDE) expressing the conservation of mass,
momentum, chemical species, enthalpy and two turbulence vari-
ables in steady, three dimensional ow. Assumptions adopted for
the ow equations are that the Coriolis force is negligible and that
the Boussinesq approximation is valid. The steady-state, 3D equa-
tions with buoyancy forces can be written in Cartesian-tensor
notation as follows:
Continuity equation:
@U
i
@x
i
0 1
Momentum equations:
U
j
@U
i
@x
j

1
q
@p
@x
i

@
@x
j
m m
t

@U
i
@x
j

@U
j
@x
i
_ _ _ _
g
i
T T
ref
T
ref
_ _
2
Energy equation:
U
j
@h
@x
j

@
@x
j
m
P
r

m
t
r
h
_ _
@h
@x
j
_ _
3
Chemical species concentration equation:
U
j
@C
@x
j

@
@x
j
m
m
t
r
C
_ _
@C
@x
j
_ _
4
where q is the uid density, U
j
the velocity vector, p the pressure, g
i
the gravitation acceleration, T
ref
the reference temperature, equal to
298 K, T the temperature, v the kinematic viscosity, v
t
is the Bous-
sinesq eddy viscosity, h is the enthalpy, P
r
the Prandtl number, r
h
the turbulence Prandtl number for h, C the concentration of chem-
ical species and r
C
is the turbulence Schmidt number for C. The
subscripts i, j(=13) denote the three space coordinates.
Since for the cases under consideration buoyancy plays an
important role in promoting turbulent mixing, the RNG k e mod-
el (Yakhot and Orszag, [16]) was employed with the modications
of Markatos and Pericleous [15], in order to account for buoyancy
forces in the production/destruction of the kinetic energy of turbu-
lence. The selection of RNG k e model was preferred to the stan-
dard one as, according to the work of Kim and Patel [19] but also
according to our experience, it was shown to give the best predic-
tion of pollutant dispersion under neutrally stratied atmosphere.
Hence, the turbulence kinetic energy k and its dissipation rate e,
obey the following transport equations:
Equation for kinetic energy of turbulence (k):
U
j
@k
@x
j

@
@x
j
m
m
t
r
k
_ _
@k
@x
j
_ _
P
k
G
b
e 5
Equation for dissipation rate of turbulence (e):
U
j
@e
@x
j

@
@x
j
m
m
t
r
e
_ _
@e
@x
j
_ _
c
e1
P
k
c
e2
e c
e3
G
b

e
k

c
l
g
3
1 g=g
0

1 b
1
g
3
e
2
k
6
where r
k
and r
e
are Prandtl numbers for k and e, respectively, and
C
e1
, C
e2
, C
e3
, C
l
, n
0
and b
1
are the model constants. Herein, P
k
sig-
nies the production rate of turbulence kinetic energy, G
b
is the
term of production of turbulence energy due to the action of buoy-
ancy forces, g is the ratio of the turbulence time scale to the mean
strain-rate scale, and v
t
is the Boussinesq eddy viscosity. P
k
, G
b
, g
and v
t
are given by the following equations, and the values of the
model constants are presented in Table 1.
P
k
m
t
@U
i
@x
j

@U
j
@x
i
_ _
@U
i
@x
j
7
G
b
g
1
T
ref
m
t
r
T
@T
@y
8
m
t
c
l
k
2
e
9
g
Sk
e
10
where in Eq. (10) S denotes the mean strain-rate of the ow, dened
as:
S
2
2S
ij
S
ij
11
where the deformation tensor, S
ij
is expressed as:
S
ij

1
2
@U
i
@x
j

@U
j
@x
i
_ _
12
3.2. Boundary conditions
Boundary conditions are specied as follows. A logarithmic pro-
le was used at the inlet boundary for the wind velocity within the
atmospheric boundary layer, and then it was kept constant above
that height. The logarithmic prole is dened by the following
equation (Tennekes [20]):
U
u

k
a
ln
z
z
0
_ _
13
Here U is the average wind speed at height z, k
a
is the von Karman
constant (=0.41), u
*
is the friction velocity and z
0
the roughness
height of the surface. Friction velocity u
*
is estimated by inserting
a known reference wind speed into Eq. (13). The reference wind
speed at 10 m above the surface is applied, being a standard height
for wind measurements. The surface roughness, z
0
, is taken as 0.3 m
for urban environment; and at the top of the boundary the normal
change of all variables is set to zero. Finally, outside the ow eld
the external pressure is considered uniform.
At the inlet boundary the prole of turbulence kinetic energy (k)
is calculated by the equation of Huser et al. [21] and the dissipation
rate of turbulence (e) is estimated by Eqs. (16) and (17), below
Turbulence kinetic energy (k):
k
u

C
l
_ 1
z
d
_ _
for z 0:9d 14
k 0:1
u

C
l
_ for z > 0:9d 15
Dissipation rate of turbulence (e):
e
C
3=4
l
k
1:5
k
a
z
for z 0:22d 16
e
C
3=4
l
k
1:5
k
a
d
for z > 0:22d 17
Table 1
Values of RNG k e model constants [16].
C
e1
C
e2
C
e3
C
l
r
k
r
e
n
o
b
1
1.42 1.68 tanh/uv 0.0845 0.719 0.719 4.38 0.012
N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030 4023
where d is the height of the atmospheric boundary layer. In this
work the atmospheric stability has been selected to be of class D,
which represents neutral conditions. Another important parameter
for the dispersion of pollutants is the height of the atmospheric
boundary layer, d, the value of which has been taken as equal to
800 m.
At solid boundaries of the physical domain, the no-slip condi-
tion is imposed and zero-ux conditions were applied at the sym-
metry plane for all variables. Moreover the momentum ux to the
walls is considered to obey the wall-function relationships of
Launder and Spalding [22].
3.3. The radiation model
The basis of all methods for the solution of radiation problems is
the radiative transfer equation (RTE) [17,18]:
s rIr; s jrIr; s Qr; s 18
which describes the radiative intensity eld, I, within the domain,
as a function of location vector (r) and direction vector (s); Q repre-
sents the total attenuation of the radiative intensity due to the gas
emission and to the in-scattered energy from other directions to the
direction of propagation, and j is the total extinction coefcient.
The discrete transfer model of Lockwood and Shah [18] discret-
izes the RTE along rays. The path along a ray is discretized by using
the sections formed from breaking the path at zone boundaries.
Assuming that the physical properties remain constant inside a
zone, Eq. (18) can be integrated from zone entry to zone exit
(Fig. 2) to yield:
I
n1
I
n
e
sn
L
n
Q
n
1 e
sn
s
n
_ _
; s
n
jL
n
19
where L
n
is the path length in the n-th zone, I
n
and I
n+1
are the inten-
sities at zone entry and zone exit, respectively, and:
Q
n
r; s
k
a
p
rT
4
n
k
s
J
n
r 20
where k
a
, k
s
, are the absorption and scattering coefcients for a gray
medium, J
n
r
1
4p
_
I
n
r; sdX is the mean intensity of the in-scat-
tered radiation, r is the Stefan-Boltzmann constant and dX is the
element of solid angle containing s.
The rays are chosen by xing nodes to all the physical surfaces,
dividing up the interior hemisphere into elements of equal solid
angle and projecting one ray into each solid angle (Fig. 3).
For gray surfaces, integration of Eq. (18) yields the required
boundary conditions:
Ir; s
e
W
p
rT
4
r
1 e
W
p
Rr 21
where e
W
is the wall emissivity, and Rr
_
sn<0
s nIr; sdX is the
radiation ux on a surface, and n is the inward pointing unit vector
normal to the surface at r.
The spectral nature of radiation is generally considered to be
important in combustion processes. The error introduced by the
assumption of a gray gas cannot be a priori estimated; however,
evaluation, of the nal results (see Section 5.4 below) indicates
that in the present case it may be only a few percent.
3.4. Heat and mass sources of re and toxic contaminants
3.4.1. Estimation of heat-release rate and mass burning rate of the
fuels
The computational work required for the present application is
very demanding in terms of CPU time. Therefore, the use of de-
tailed chemistry modeling was avoided. Instead, the heat and mass
source of the re on the surface of the external oating-roof tank is
simulated as a volumetric source, which releases heat at a constant
rate. Hence, the re on the top of the external oating-roof tank is
characterized as a large pool re. Thus, the estimation of mass loss
rate _ m
00
and the total heat-release rate _ q (HRR) are important input
parameters for the model. The calculation of the above quantities is
determined by the following expressions (Babrauskas [23]):
_ m
00
_ m
00
1
1 e
jbD
22
_
q _ m
00
DH
c;eff
A
f
23
where _ m
00
1
is the innite-diameter pool mass-loss rate, k is the
absorption extinction coefcient of the ame, b is a mean beam
length corrector, D is the pool diameter, DH
c,eff
is the heat of com-
bustion and A
f
is the surface area of the pool. In this study the type,
HRR and mass burning rate of fuels and the appropriate values for
the _ m
00
1
, k, b and DH
c,eff
have been adopted from the work of Babr-
auskas [23]. The mass burning rate for both liquid fuels used has
been calculated for steady-state conditions, even though the mass
burning rates vary with time.
The estimated value of HRR for all accident scenarios presented
here (Table 2) has been calculated from the total HRR, which is
estimated from Eq. (23), minus the computed thermal radiation.
3.4.2. Emission factors for the toxic contaminants of the fuels
Combustion products of various fuels are characterized by the
presence of toxic pollutants such as smoke, sulfur dioxide (SO
2
),
fire
zone
I
n
I
n+1
L
zone
boundary
physical
surface
node
Fig. 2. Sub-division of the solution domain into zones.
Fig. 3. Projection of a ray at a node into a number of angular divisions.
4024 N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030
carbon dioxide (CO
2
), carbon monoxide (CO), etc. The prediction of
their ground-level concentration is of great importance for this
study, in order to assess its impact on human health.
In the present study, two fuels have been examined, crude oil
and diesel oil, and the ground-level concentration of smoke, CO,
SO
2
, Polycyclic Aromatic Hydrocarbons (PAH) and Volatile Organic
Compounds (VOC) have been estimated.
Unfortunately, experimental data for large-scale pool res with
diameter 85 m are not available. Many attempts from various
researchers (Notarianni et al. [24]; Koseki et al. [25]; Evans et al.
[26]; Walton et al. [27]) have been performed in the hope of study-
ing the phenomenon but without prospective results, because of
the complexities of the physical problem. The re experiments
have been conducted for various pool re sizes of 5, 10, 15, 17,
20 and 50 m diameter.
Because of the absence of experimental data for the emission
factors of pollutants, available data from the literature have been
used. More specically, laboratory data have been adopted in
meso-scale and offshore burning experiments (Notarianni et al.
[24]; Koseki et al. [25]; Evans et al. [26]; Walton et al. [27]; Booher
and Janke, [28]; Lemieux et al. [29]). Then, with the appropriate
data for emission factors, it has been easy to calculate the source
term of the pollutants for the model. Thus, the source term of
the pollutants is estimated by the following equation:
SP m
00
EF

A
f
24
where SP is the source term of the pollutant, _ m
00
is the mass burning
rate of the fuel, EF is the emission factor of the pollutant for the spe-
cic fuel and A
f
is the area of pool re.
3.4.2.1. Crude oil. Crude oil contains both hydrocarbons and differ-
ent compounds containing nitrogen, oxygen, hydrogen, sulfur and
metals. The liquid quantity used in the numerical simulation was
set equal to the capacity of the tank, and the density was 880 kg/
m
3
(Babrauskas [23]).
According to the evidence of the works of Evans et al. [26] and
Walton et al. [27] the value of smoke yield for crude oil appears to
be from 10% to 15%. Herein, the value of smoke yield is taken equal
to 12.5% which is the average of the above measurements. The
source term for smoke is calculated with the help of Eq. (24) and
is equal to 24.1 kg/s.
The emission factor for CO is taken equal to 0.09 kg CO/kg fuel
according to the evidence of Evans et al. [26] and Evans [30]. For
SO
2
, the value of emission factor depends on the sulfur content
of the fuel and it is taken equal to 0.04 kg SO
2
/kg fuel. PAHs and
VOCs emissions factors are adopted from the research of Booher
and Janke [28] and their values are 0.0004 kg PAHs/kg fuel and
0.0005 kg VOCs/kg fuel, respectively.
3.4.2.2. Diesel oil. Diesel oil is composed of about 75% saturated
hydrocarbons and 25% aromatic hydrocarbons. Similar to the
assumptions for the crude oil, the smoke yield for diesel oil varies
from 15% up to 20% (Walton et al.). Herein, the value of smoke
yield is taken equal to 17.5%, which is the average of the above val-
ues from the literature.
In the present work the emission factor for estimating the con-
centration of sulfur dioxide (SO
2
) for diesel is based on the sulfur
content of the fuel, as for the crude oil, and the value is 0.02 kg
SO
2
/kg fuel. Finally the emission factors for CO, PAHs and VOCs
are adopted from the research of Booher and Janke [28] and their
values are taken equal to 0.03 kg CO/kg fuel, 0.0009 kg PAHs/kg
fuel and 0.0004 kg VOCs/kg fuel, respectively.
4. Method of solution
4.1. Solution procedure and grid sensitivity analysis
The set of the model partial-differential equations, along with
the appropriate boundary conditions and the auxiliary relations
have been solved by means of the Finite Volume Method (Patankar
and Spalding [31]).
The model is implemented in the general computer program
PHOENICS (Spalding [32]).
The treatment of the convective terms in conservation equa-
tions determines the accuracy of the solutions. These terms are dif-
cult to handle numerically because the more accurate schemes
tend to be less robust and/or slower to converge. Consider the con-
trol volume shown in Fig. 4. The scheme used here is illustrated
taking an example of the west face (w), for variable u.
Monotone Upstream-Centered Schemes for Conservation Laws
(MUSCL) [12] are modications to the higher-order upwind
scheme with ux limiters to ensure boundedness and monotonic-
ity of the solution. In this work a second-order MUSCL scheme is
used as follows:
U
w
1
1
2
W
_ _
U
W

1
2
WU
WW
25
Where W is the ux limiter which is given in terms of the ratio:
r
U
P
U
W
U
W
U
WW
26
Van Leers ux limiter is used [33]:
W r jrj=1 jrj 27
and deferred correction was applied:
U
w
U
L
w
aU
H
w
U
L
w

old
28
where U
L
w
stands for the usual rst-order upwind approximation
and U
H
w
for the above Van Leer approximation. The term in brackets
is evaluated using values from the previous iteration, as indicated
by the superscript old. Although this term is obviously small com-
Fig. 4. The control volume.
Table 2
Accident scenarios with plume rise estimation.
Scenarios Fuel HRR (MW/m
2
) Wind velocity (m/s) Plume rise (m)
(1) Crude oil 1 8 1746
(2) Crude oil 1 11 1335
(3) Diesel 1.3 8 1781
(4) Diesel 1.3 11 1417
N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030 4025
pared to U
L
w
, the factor a that blends the two schemes is taken here
as 0.4, in order to diminish further the explicit treatment, thus
improving convergence.
An attempt was also made to use CUPID [13], that intuitively
appears well suited, but convergence was not obtained when it
was used for the momentum equations; its use was thus restricted
to the equations of the scalar variables.
In addition, the computational domain was divided into regions
for the better simulation of the phenomenon in each space dimen-
sion. The selection of the appropriate grid for the numerical simu-
lation was based on a grid sensitivity analysis among various
different grids, between 2 and 5 million cells. It was found that
the results of grids 66 156 248, 70 167 260 and 80
190 280 were sufciently close to each other. Therefore, the grid
of 66 156 248(2,553,408 cells) was used for the nal runs, as
the results were virtually grid-independent and the computational
cost signicantly lower than for the ner grids.
4.2. Convergence and time requirements
Runs were performed on a Linux PC with CPU 2.8 GHz and main
memory of 1 Gbytes. Convergence was achieved by applying un-
der-relaxation techniques. More specically, the following crite-
rion for all dependent variables was met:
maxj/
n1
/
n
j 10
4
29
between sweeps n and n+1. The false transient type of relaxation
was used for the velocity components, enthalpy and concentration
of chemical species, the value of false time step taken equal to
0.1. For pressure and turbulence the usual linear type of relaxa-
tion was employed with a value of 0.30. About 10,000 sweeps of
the computational domain were needed to ensure full convergence
for the grid of 66 156 248 (2,553,408 cells). The CPU time re-
quired was about 200 h.
5. Results and discussion
The overwhelming amount of results obtained dictates a judi-
cious choice of some of them for presentation. The results are pre-
sented below, according to the authors judgement of what is
practically meaningful, and are discussed under headings of plume
rise, ground-level concentrations and hazard identication.
Figs. 57 present velocity vectors and enthalpy contours around
the tank on re. These gures show clearly the way the ow devel-
ops (magnitude and direction) and the location of the highest
temperatures.
5.1. Plume rise analysis
Plume rise analysis was performed for four different scenarios.
More specically, crude oil and diesel are examined with two values
of wind velocity, 8 and 11 m/s. Table 2 presents the general charac-
teristics of four accident scenarios with the results from the
calculation of plume rise for all cases and Fig. 8 illustrates the center
line curve of the plume downwind from the center of the tank for
each scenario. The plume axis is specied by tracing the location of
the maximum concentration at each vertical z-slab.
According to the results from the four accident scenarios the
highest plume rise occurs for scenario (3), at 1781 m. This is due
to the highest HRR (1.3 MW/m
2
) strength of the re and lowest
wind velocity (8 m/s), compared with the other scenarios. There-
fore, the buoyancy forces are strong enough and with the help of
medium-speed wind the plume is transported high in the atmo-
spheric boundary layer.
On the other hand, the lowest plume trajectory takes place in
scenario (2), at the height of 1335 m. This attitude is explained
by the lowest HRR (1 MW/m
2
) strength of re and the highest va-
lue of wind velocity (11 m/s) among the other scenarios. Thus, the
buoyancy forces are not strong enough to raise the plume and, due
to the high wind velocity (11 m/s), the plume is transported to the
ground-level with anticipated consequences for the environment
and human health.
5.2. Ground-level concentration
The estimation of the anticipated consequences from the toxic
contaminants is usually based on the ground-level concentration
of these toxic gases after dispersing. For all accident scenarios,
ground-level concentration is examined for two zones at a height
of 1 m from the ground. The rst zone lies closer to the tank re
and has a range of approximately 10001500 m from the tank
and the pollutant concentrations there are rather high. The second
zone is further away from the tank and follows the rst zone up to
the end of the domain. The toxic contaminants concentrations are
Fig. 5. Velocity vectors (YZ plane).
4026 N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030
lower, but it is still possible to be signicant depending on the
toxic characteristics of the fuel.
From the parametric analysis it is concluded that the worst-case
scenario, among the four scenarios which were examined, for
smoke ground-level concentration is scenario (4) and the lowest
smoke concentration appears for scenario (1).
The above results are shown in Fig. 9 and it is seen that the
ground-level concentrations of smoke in the rst zone near to
the tank are in the range of 315 mg/m
3
, which are very high.
The values in the second zone are small and range from 0 to
3 mg/m
3
. Substantial differences in ground-level concentrations
of smoke depend also on the fuel because of the different % smoke
yield.
According to Fig. 10, it is observed that CO emissions for sce-
nario (2) are the highest, while the lowest values appear for sce-
nario (3). In the rst zone, near to the tank, the ground-level
concentrations of CO present high values within a range of
Fig. 6. Enthalpy contours (XZ plane).
Fig. 7. Enthalpy contours (YZ plane).
Fig. 8. Plume rise vs. the downwind distance from the tank for each scenario.
Fig. 9. Ground-level concentrations of smoke for all scenarios at 1 m height from
the ground.
N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030 4027
28 mg/m
3
, while the range of values for the second zone, away
from the tank, is from 0 to 2 mg/m
3
. An important factor of these
differences between the ground-level concentration of crude oil
and diesel is the emission factor of CO which is different for the
two fuels.
Moreover, in Fig. 11 for SO
2
emissions, scenario (2) presents the
highest ground-level concentration and scenario (3) the lowest.
Similar to the above situations, the ground-level concentrations
for SO
2
in the rst zone are extremely high, with concentration val-
ues from 1.5 to 3.8 mg/m
3
. In the second zone the values for SO
2
emissions are very low, from 0 to 1.5 mg/m
3
.
From the above Figs. 911, it is observed that for all ground-
level concentrations of pollutants, except for smoke, scenario (2)
gives the highest values of concentration and scenario (3) the
lowest.
Moreover, Figs. 1214 illustrate concentration of smoke, CO and
SO
2
proles taken at 5 km perpendicular to the ow axis.
It is observed that for low heights the concentrations of smoke,
CO and SO
2
are extremely high at a position 5 km downwind of the
ow axis. On the contrary, with the increase of height the concen-
trations of smoke, CO and SO
2
decrease.
In Fig. 12, it is noticed that for heights in the range of 500
750 m above the ground-level, the value of maximum concentra-
tion of smoke is 9 mg/m
3
. Scenario (4) presents the maximum
value of smoke concentration. The lowest value of smoke concen-
tration is observed for scenario (1) and it is equal to 4 mg/m
3
at a
height of 800 m.
In Fig. 13, it is observed that the maximum value for CO concen-
tration is found for scenario (2) and is equal to 4 mg/m
3
. On the
other hand the minimum value appears for scenario (3) and it is
1 mg/m
3
, approximately.
Finally, in Fig. 14, the same behaviour is exhibited for the con-
centration of SO
2
, the maximum value found for scenario (2) and
Fig. 10. Ground-level concentrations of CO for all scenarios at 1 m height from the
ground.
Fig. 11. Ground-level concentrations of SO
2
for all scenarios at 1 m height from the
ground.
Fig. 12. Concentration of smoke vs. the height at a position of 5 km from the ow
axis.
Fig. 13. Concentration of CO vs. the height at a position of 5 km from the ow axis.
Fig. 14. Concentration of SO
2
vs. the height at a position of 5 km from the ow axis.
4028 N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030
the minimum value for scenario (3), with values of 1.8 and 0.7 mg/
m
3
, respectively.
5.3. Hazard identication analysis
Hazard identication analysis is performed with the character-
ization of risk zones (Fig. 15), by comparing the ground-level con-
centration of the pollutants with the safety limits.
The present methodology for the determination of hazard range
can be described by comparing the ground-level concentrations of
the pollutants (smoke, CO, SO
2
) with the Lethal Concentration (LC
1
)
and Immediately Dangerous to Life and Health (IDLH) values of
them [9]. The (LC
1
) determines the range of the high-risk zone
where only re ghters are allowed wearing the appropriate gear,
and (IDLH) species the hazard range for the general population
[9]. The values for the above safety limits are adopted from the Na-
tional Institute for Occupational Safety and Health (NIOSH) web
site (http://www.cdc.gov/niosh/homepage.html). The LC
1
(CO) =
9190.18 mg/m
3
and IDLH (CO) = 1374.23 mg/m
3
, for SO
2
LC
1
(SO
2
) = 2525.97 mg/m
3
and IDLH (SO
2
) = 262 mg/m
3
. Limit values
for smoke are LC
1
(smoke) = 25,000 mg/m
3
and IDLH (smo-
ke) = 2500 mg/m
3
.
From the parametric runs above, the ground-level concentra-
tions of toxic pollutants for the two zones are computed for all sce-
narios. According to those data, there are no death zones due to
the concentrations of smoke, CO and SO
2
. In the rst zone the con-
centrations which are observed are high, especially for smoke, but
do not exceed the safety limits of LC
1
or IDLH.
5.4. Some further results
It is interesting to point out that comparison with earlier work
by the author and his colleagues [10], reveals that despite the pres-
ent more sophisticated treatment of the convective terms (use of
Van Leer and CUPID Schemes rather than the hybrid one of Spal-
ding [34]), and the detailed calculation of the radiation effects,
the results do not differ by more than 78%, with maximum of
10% for the enthalpy. In [10] the estimated value of HRR for all sce-
narios was calculated from the total HRR (Eq. (23)), minus the frac-
tion of thermal radiation, which constitutes a 30% of the total HRR,
according to literature [35]. A rough calculation of present com-
puted radiation contribution reduces this fraction to 20%, which
may be attributed to the very large re diameter that leads to
smoke obscuration [35].
In total 18 scenarios were run for different re strengths and
wind speeds. Some results concerning the effect of the tank diam-
Fig. 15. Conguration of risk zones.
Fig. 16. Ground-level concentration of benzene at 1 m height from the ground for
various scenarios. Tank diameter is 50 m.
Fig. 17. Ground-level concentration of benzene at 1 m height from the ground for
various scenarios. Tank diameter is 70 m.
Fig. 18. Ground-level concentration of benzene at 1 m height from the ground for
various scenarios Tank diameter is 85 m.
N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030 4029
eter on the pollutants ground concentration are shown in Figs. 16
18. It is seen that the larger the diameter the higher the pollutant
concentration and the more abrupt the slope of its distribution
with distance downwind.
6. Conclusions
The work described herein focuses on the computation of toxic
contaminant releases from large hydrocarbon-tank res. An inte-
grated methodology is developed and applied in order to quantify
the consequences. The cases studied were four accident scenarios
with different HRR and wind velocities. The available data were
inadequate for the detailed description of the accident scenarios
and several assumptions had to be made, concerning the input
data, based on the worst-case scenario.
Fires of great diameter in storage tanks of liquid fuels, such as
crude oil, kerosene, etc., are very difcult to be extinguished. This
means that the re may burn for days until all the fuel in the tank
is burnt. In these cases the dispersion of the plume and toxic pol-
lutants depends on the meteorological conditions and the charac-
teristics of the re. The most serious consequences are exhibited
for high wind velocity and weak re, i.e., when the plume is
trapped very close to the ground with high values of pollutant con-
centrations. On the contrary, with low wind velocity and strong
re, the plume of smoke is dispersed at great heights away from
the ground and the pollutant concentrations decrease. It is also
concluded that, for the accident scenarios considered, there are
no death zones due to the concentrations of smoke, CO and
SO
2
and that near and away from the tank, the safety limits are
not exceeded.
The results have been shown to be physically reasonable. The
author believes that the results shown in this paper, and the far
greater number of similar results that have been obtained but
not presented, for many accident scenarios, conrm that the math-
ematical formulation is a satisfactory one, from the point of view of
conformity with physical reality; and, that the procedure followed
is reliable for gaining useful insight into industrial accidents and
therefore alleviating their consequences.
Acknowledgements
The author wishes to thank Prof. D.B. Spalding and CHAM Ltd. of
London, for permitting the use of the software PHOENICS with
which the reported calculations were performed.
References
[1] U. Beck, The Risk Society, Sage, London, United Kingdom, 1992.
[2] L. Chang, C.C. Lin, A study of storage tanks accidents, J. Loss Prev. Process Ind.
19 (2006) 5159.
[3] A.F. Ghoniem, X. Zhang, O. Knio, H.R. Baum, R.G. Rehm, Dispersion and
deposition of smoke plumes generated in massive res, J. Hazard. Mater. 33
(1993) 275293.
[4] A.F. Ghoniem, X. Zhang, A computational model for the rise and dispersion of
wind-blown, buoyancy driven plumes, Part I Neutrally stratied atmosphere,
Atmos. Environ. 27 (1993) 22952311.
[5] K.B. McGrattan, H.R. Baum, R.G. Rehm, Numerical simulation of smoke plumes
from large oil res, Atmos. Environ. 30 (1996) 41254136.
[6] G. Miles, M. Cox, M.N. Christolis, C.A. Christidou, A.G. Boudouvis, N.C.
Markatos, Modelling the environmental consequences of res in warehouses,
in: T. Kashiwagi (Ed.), Proceedings of the Forth International Symposium on
Fire Safety Science, Ottawa, Canada, 1994, pp. 12211232.
[7] M.N. Christolis, C.A. Christidou, A.G. Boudouvis, N.C. Markatos, Modelling
pollutants dispersion around buildings on re, in: H. Power (Ed.), Proceedings
of the International Conference on Air Pollution, Porto Carras, Greece, 1995, pp.
8592.
[8] G.M. Sideris, M.N. Christolis, N.C. Markatos, Numerical simulation of pollutants
dispersion around warehouses on re, in: D.T. Tsahalis (Ed.), Proceedings of the
First International Conference, From Scientic Computing to Computational
Engineering, Athens, Greece, 2004.
[9] D.A. Kefalas, M.N. Christolis, Z. Nivolianitou, N.C. Markatos, Consequence
analysis of an open re incident in a pesticide storage plant, J. Loss Prev.
Process Ind. 19 (2006) 7888.
[10] C.D. Argyropoulos, M.N. Christolis, N.C. Markatos, Z. Nivolianitou, Assessment
of acute effects for re-ghters during a fuel-tank re, in: Proceedings of the
4th International Conference on Prevention of Occupational Accident in a
Changing Work Environment, WOS 2008, 30 September3 October 2008,
Crete, Greece.
[11] N.C. Markatos, Mathematical modelling of single- and two-phase ow
problems in the process industries, Oil & Gas Sci. Technol. Rev. IFP 48
(1993) 631662.
[12] W.K. Anderson, J. Thomas, B. Van Leer, Comparison of nite volume ux vector
splittings for the Euler equations, AIAA J. 24 (1986) 1453.
[13] M.K. Patel, M. Cross, N.C. Markatos, An assessment of ow oriented schemes
for reducing false diffusion, Int. J. Numer. Meth. Eng. 26 (1988) 2279.
[14] D.B. Spalding, Mathematical modelling of uid-mechanics, heat-transfer and
chemical-reaction processes, A Lecture Course, CFDU Report, HTS/80/1,
Imperial College, London, 1980.
[15] N.C. Markatos, K.A. Pericleous, Laminar and turbulent natural convection in an
enclosed cavity, Int. J. Heat Mass Transfer 27 (5) (1984) 755772.
[16] V. Yakhot, S.A. Orszag, Renormalization group analysis of turbulence I. Basic
theory, J. Sci. Comput. 1 (1986) 351.
[17] R. Siegel, J.R. Howell, Thermal radiation heat transfer, 4th ed., Hemisphere,
Washington, 2002.;
E.P. Keramida, H.H. Liakos, M.A. Founti, A.G. Boudouvis, N.C. Markatos,
Radiative heat transfer in natural gas-red furnaces, Int. J. Heat Mass
Transfer 43 (2000) 351371.
[18] F.C. Lockwood, N.G. Shah, A new radiation solution method for incorporation
in general combustion prediction procedures, in: Proceedings of the
Eighteenth Symposium (Inter.) on Combustion, The Combustion Institute,
Pittsburgh, 1981, pp. 14051416.
[19] H.M. Kim, V.C. Patel, Test of turbulence models for wind ow over terrain with
separation and recirculation, Bound. Layer Meteorol. 94 (2000) 521.
[20] H. Tennekes, The logarithmic wind prole, J. Atmos. Sci. 30 (1973) 234
238.
[21] A. Huser, P.J. Nilsen, H. Sketun, Application of k e model to the stable ABL:
pollution in complex terrain, J. Wind Eng. Ind. Aerodyn. 67 and 68 (1997) 425
436.
[22] B.E. Launder, D.B. Spalding, The numerical computation of turbulent ows,
Comput. Methods Appl. Mech. Eng. 3 (1974) 269289.
[23] V. Babrauskas, Estimating large pool re burning rates, Fire Technol. 19 (1983)
251261.
[24] K.A. Notarianni, D.D. Evans, W.D. Walton, D. Madrzykowski, J.R. Lawson,
Smoke production from large oil pool res, Interam 93, Fire Safety, in: C.A.
Franks (Ed.), 6th Int. Fire Conference, Oxford, England; Interscience
Communications Ltd., London,, 1993, pp. 111119.
[25] H. Koseki, Y. Iwata, Y. Natsume, T. Takahashi, T. Hirano, Tomakomai large scale
crude oil re experiments, Fire technol. 36 (2000) 2438.
[26] D.D. Evans, G.W. Mulholland, H.R. Baum, W.D. Walton, K.B. McGrattan, In situ
burning of oil spills, J. Res. Natl. Inst. Stand. Technol. 106 (2001) 231
278.
[27] W.D. Walton, W.H. Twilley, A.D. Putorti, R.R. Hiltabrand, Smoke measurements
using an advanced helicopter transported sampling package with radio
telemetry, in: Proceedings of the 18th Arctic and Marine Oilspill Program
Technical Seminar, Environment Canada, Ottawa, Ontario, 1995, pp. 1053
1074.
[28] L.E. Booher, B. Janke, Air emissions from petroleum hydrocarbon Fires during
controlled burning, Am. Ind. Hyg. Assoc. J. 58 (1997) 359365.
[29] P.M. Lemieux, C.C. Lutes, D.A. Santoianni, Emissions of organic air toxics from
open burning: a comprehensive review, Progress Energy Combust. Sci. 30
(2004) 132.
[30] D.D. Evans, Combustion of oil spills on water (OCS Study MMS 88-0057)
Environment Canada, Reston, VA, 1991, pp. 169179.
[31] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mats and
momentum transfer in three-dimensional parabolic ows, Int. J. Heat Mass
Transfer 15 (1992) 1787.
[32] D.B. Spalding, A general purpose computer program for multi-dimensional
one- and two-phase ow in mathematics and computers in simulation, vol.
XXIII, North Holland, Amsterdam, 1981. pp. 267276.
[33] B. Van Leer, Towards the ultimate conservative difference scheme. A second
order sequel to Godunovs method, J. Comp. Phys. 32 (1979) 101.
[34] D.B. Spalding, Novel nite-difference formulation for differential expressions
involving both rst and second derivatives, Int. J. Numer. Methods Eng. 4 (4)
(1972) 551559.
[35] K.B. McGrattan, H.R. Baum, A. Hamins, Thermal radiation from large pool res,
NISTIR 6546, National Institute of Standards and Technology, U.S. Department
of Commerce, Washington, DC, 2000.
4030 N.C. Markatos et al. / International Journal of Heat and Mass Transfer 52 (2009) 40214030