Consequence Analysis of Aqueous Ammonia Spills Using An Improved Liquid Pool Evaporation Model
Consequence Analysis of Aqueous Ammonia Spills Using An Improved Liquid Pool Evaporation Model
Consequence Analysis of Aqueous Ammonia Spills Using An Improved Liquid Pool Evaporation Model
A Thesis
by
VIJAY RAGHUNATHAN
MASTER OF SCIENCE
December 2004
A Thesis
by
VIJAY RAGHUNATHAN
MASTER OF SCIENCE
December 2004
ABSTRACT
Source term modeling is the key feature in predicting the consequences of releases from
hazardous fluids. Aqueous ammonia serves the purpose of a reducing medium and is
replacing anhydrous ammonia in most of the Selective catalytic reduction (SCR) units.
This newly developed model can estimate the vaporization rate and net mass evaporating
into the air from a multicomponent non- ideal chemical spill. The work has been divided
into two parts. In the first step a generic, dynamic source term model was developed that
can handle multicomponent non-ideal mixtures. The applicability of this improved pool
model for aqueous ammonia spills was then checked to aid in the offsite consequence
analysis of aqueous ammonia spills.
The behavior of the chemical released depends on its various inherent properties,
ambient conditions and the spill scenario. The different heat transfer mechanisms
associated with the pool will strongly depend on the temperature of the liquid pool
system at different times. The model accounts for all the temperature gradients within
the contained pool and hence helps us establish better estimation techniques for source
terms of chemical mixtures. This research work will help obtain more accurate and
reliable liquid evaporation rates that become the critical input for dispersion modeling
studies.
iv
DEDICATION
To Amma, Daddy, Ashwin, Rajam Patti & Sethu Patti who have extended great love and
support throughout my life
v
ACKNOWLEDGEMENTS
I would like to express my sincere gratitude to my academic advisor, Dr. Sam Mannan,
for his boundless motivation throughout my master’s program. He has definitely been a
source of inspiration and a role model in my life. I am also grateful to Dr. Michael
Lindell and Dr. Mahmoud El-Halwagi for their cooperation and support. I also wish to
thank Dr. William. J. Rogers for the insightful discussions I have had with him related to
my research. I would also like to express my special thanks to Dr. Harry West for his
timely help and excellent guidance. I also take this opportunity to thank Dr. Peck and Dr.
Ganesan Gopalakrishnan for helping me get financial assistance. I am also grateful to all
my wonderful teammates and the staff of the Mary Kay O’Connor Process Safety Center
for all their help. I would like to extend my special thanks to Srini and Adhi for their
immense help. I also owe a special thanks to Towanna and Ninette for helping me with
all the paperwork throughout these two years. I also take this opportunity to thank all my
nearest and dearest friends around the globe for all their support. Finally I once again
thank my beloved parents and my brother who gave me an opportunity to realize my
dreams.
vi
TABLE OF CONTENTS
Page
DEDICATION ..................................................................................................................iv
ACKNOWLEDGEMENTS ...............................................................................................v
CHAPTER
I INTRODUCTION ................................................................................................1
Page
CHAPTER
IV MODEL DETAILS...........................................................................................27
REFERENCES.................................................................................................................80
APPENDIX A ..................................................................................................................83
APPENDIX B ..................................................................................................................89
APPENDIX C ..................................................................................................................97
VITA ..............................................................................................................................100
viii
LIST OF FIGURES
Page
Fig. 1. Simple cycle system................................................................................................1
Fig. 2. Combined cycle system ..........................................................................................1
Fig. 3. Low dust system .....................................................................................................2
Fig. 4. High dust system.....................................................................................................2
Fig. 5. Catalyst modules .....................................................................................................3
Fig. 6. Reaction mechanism ..............................................................................................3
Fig. 7. Ammonia injection process.....................................................................................4
Fig. 8. Flowchart for identifying the covered processes in the facility ..............................8
Fig. 9. Steps for consequence analysis ............................................................................12
Fig. 10. Generic source term modeling procedure ..........................................................20
Fig. 11. Liquid pool mechanisms .....................................................................................31
Fig. 12. Comparison of non-ideality with Raoult’s law case ...........................................34
Fig. 13. Mass flux of ammonia in air ...............................................................................48
Fig. 14. Cumulative mass of ammonia in air ...................................................................49
Fig. 15. Cumulative mass of water vapor in air ...............................................................49
Fig. 16. Surface temperature ............................................................................................50
Fig. 17. Bulk temperature.................................................................................................50
Fig. 18. Composition of ammonia in pool .......................................................................51
Fig. 19. Partial pressure variation.....................................................................................53
Fig. 20. Mass flux of ammonia in air ...............................................................................54
Fig. 21. Effect of wind speed ...........................................................................................55
Fig. 22. Effect of composition on partial pressure ...........................................................56
Fig. 23. Effect of composition on surface temperature ....................................................57
Fig. 24. Experiment 1-Mass of ammonia in air................................................................59
ix
Page
Fig. 25. Experiment 2-Mass of ammonia in air................................................................61
Fig. 26. Experiment1-Temperature gradients ..................................................................62
Fig. 27. Experiment 1-Composition of ammonia.............................................................62
Fig. 28. Experiment 2-Temperature gradients .................................................................63
Fig. 29. Experiment1-Composition of ammonia..............................................................63
Fig. 30. Mass flux of ethanol............................................................................................65
Fig 31. Mass flux of dichloromethane ............................................................................67
Fig. 32. Matlab code.........................................................................................................73
Fig. 33. Density of vapor cloud........................................................................................87
Fig. 34. Meshed pool surface ...........................................................................................92
Fig. 35. Top surface contour for a 10 cm pool depth .......................................................93
Fig. 36. Ground surface contour for a 10 cm pool depth .................................................94
Fig. 37. Top surface contour for a 5 cm pool depth .........................................................95
Fig. 38. Ground surface contour for a 5 cm pool depth ...................................................96
Fig. 39. Quantity of HCN released into air ......................................................................98
Fig. 40. Quantity of Acrolein released into air.................................................................98
x
LIST OF TABLES
Page
Table 1. Assumptions for modeling scenarios ................................................................16
Table 2. Base spill conditions .........................................................................................47
Table 3. Cumulative mass of ammonia in air..................................................................52
Table 4. Spill conditions for Mikesell et al. data ............................................................58
Table 5. Predicted values for experiment 1.....................................................................59
Table 6. Predicted values for experiment 2.....................................................................60
Table 7. Test conditions for Norman Dowell data ..........................................................64
Table 8. Predicted values for ethanol ..............................................................................65
Table 9. Predicted values for propylene..........................................................................66
Table 10.Meteorological conditions around the U.S........................................................68
Table 11.Test conditions for density validation ...............................................................86
Table 12.Cloud density at different time instants ............................................................88
1
CHAPTER I
INTRODUCTION
1.1 Introduction
Nitrogen oxide (NOx) is one of the most common air pollutants and its removal
technology has gained great importance from an environmental point of view. Selective
Catalytic reduction (SCR) is the most commonly used technology for control of effective
NOx emissions from utility boilers and combustion turbines nowadays (Pritchard et al.
1995). Its applications also include reduction of NOx emissions from diesel engines,
process gas streams like nitric acid plants. The flue gas emitted from the refinery units
need to be treated with a suitable catalyst under the presence of a reduction medium and
ammonia serves that purpose in this process.
SCR
Stack
SCR
This technology has proven to be the most economically viable solution, with cost
estimates in the range of $20-$30 per kilowatt for natural gas units. In the above shown
boiler SCR arrangements (Figures1-4) “ NH 3 ” is the location of the ammonia injection
grid, ESP is the electrostatic precipitator, FGD is the flue gas desulphurization unit, SH
is the superheater, HP, IP, LP evap are the high, intermediate and low-pressure
evaporator. The SCR equipment includes a reactor chamber, catalyst modules,
ammonia storage system, ammonia vaporization and injection system, and monitoring
equipment and sensors. The catalyst modules used in the SCR process is shown below
(Figure 5).
3
catalyst
NO + NO2 + 2 NH 3 ⎯⎯⎯⎯
⎯→ 2 N 2 + 3H 2 O
catalyst
NO + 4 NH 3 + O2 ⎯⎯⎯⎯
⎯→ 4 N 2 + 6 H 2 O
catalyst
6 NO2 + 8 NH 3 ⎯⎯⎯⎯
⎯→ 7 N 2 + 12 H 2 O
NH3
Catalyst
NOx NOx N2
NOx NH3
H2 O
NOx NOx
N2 Gas
Gas NH3
The actual chemistry behind the reduction process is shown in Figure 6. The catalytic
reaction can take place over a wide temperature range (300 deg F-1100 deg F) with
typical applications between 500◦ F- 800◦ F. The ammonia also combines with the SO3
gas to form ammonium sulfate salts and this could cause corrosion or plugging at times.
⎯→( NH 4 ) 2 SO4
catalyst
2 NH3 + SO3 + H 2 O ⎯⎯⎯⎯
Some of the advantages of using a SCR process in NOx removal process are
• Simple process
• Low running cost
5
• Easy operation
• No secondary environmental pollution
• High reliability
• High NOx removal efficiency
• Small pressure drop
Threshold Limit Value -Time Weighted Average, (TLV- TWA): Time weighted average
for a normal 8-hour workday or 40-hour workweek, to which nearly all workers can be
exposed, day after day without adverse effects.
Threshold Limit Value –Short Term Exposure Limit, (TLV- STEL): The maximum
concentration to which workers can be exposed for a period of up to 15 minutes
continuously without suffering (1) intolerable irritation (2) chronic or irreversible tissue
change (3) narcosis of sufficient degree to increase accident proneness, impair self-
rescue, or materially reduce worker efficiency, provided that no more than 4 excursions
6
per day are permitted, with at least 60 minutes between exposure periods, and provided
that the daily TLV-TWA is not exceeded.
Immediate Damage to Life and Health, IDLH: (NIOSH 2001) Occupational Safety and
Health (NIOSH) in mid 1970’s.This refers to a concentration, formally specified by a
regulatory value, and defined as the maximum exposure concentration of a given
chemical in the workplace from which one could escape within 30 minutes without any
escape-impairing symptoms or any irreversible health effects. The IDLH for ammonia is
300 ppm.
CHAPTER II
2.1 Introduction
EPA’s risk management program requirements may be found in Part 68 of Volume 40 of
the code of federal regulations and it works with the industry, local, state and federal
government agencies to assist sources in complying with the regulatory requirements.
According to this rule if the facility handle, manufacture or store any toxic or flammable
substances as listed in the appendix of 48 CFR Part 68,they are required to develop and
implement a risk management program that includes a five-year accident history, an
offsite consequence analysis, an accident prevention program and an emergency
response program.
Applicability of the Risk Management Rule: Part 68 of the CFR in general applies to the
owner of a stationary source (facility) that has more than a threshold quantity of a
regulated substance in a process. The decision approach for identifying the covered
processes is presented in the page below (figure 8).
8
Is your facility a 1
stationary source?
Yes
Any regulated 2 N
substances in your
facility?
Yes
Is the Regulated
substance above the
threshold quantity?
Yes
No
The facility is subject to the
rule
Fig. 8. Flowchart for identifying the covered processes in the facility, (EPA 1998)
9
Stationary Source (EPA 1998): CAA section 112(r) defines “stationary sources” as
“Any buildings, structures, equipment, installations, or substance emitting stationary
activities
• That belongs to same industrial group.
• That is located on one or more contiguous properties.
• That is under control of the same person.
• From which any accidental release may occur.
Regulated Substances: EPA has included a list of 77 toxic chemicals and 63 flammables
that can cause serious health effects or have the potential to form vapor clouds and
explode if released. The rule also includes those flammable mixtures that meet the
criteria for the National Fire Protection Association. The rule applies only to those
facilities that have one of those listed substances over 1 percent concentration or above
their threshold quantities.
• Records showing the estimated quantity of worst-case release, release rate and
duration of release.
This being the first or basic level imposes limited hazard assessment requirements.
Program 2: This is the default level because any covered process that is not eligible for
Program1 or Program 3 is by default subject to Program 2 requirements.
Program 3: A process that is not eligible for Program 1 and meets one of the two criteria
specified below can be subjected to Program 3 requirements.
The eligibility criteria for Program 3 are as follows:
• Process does not meet the eligibility requirements for Program 1
• Process is subject to OSHA PSM or process is one of the nine SIC codes.
If the release of the regulated substance is below its threshold quantity it is not
mandatory to report its release in the five-year accident history.
Information for the Five-Year Accident History: For every reported release the following
data needs to be included.
• Date and Time of the accident/incident
• Release duration
• Quantity Released
• Nature of release (gas, liquid, aerosol, fire or explosion)
• Release source
• Weather Conditions
• Onsite impacts
• Offsite impacts
• Initiating event
• Changes introduced as a result of accident
Select Scenario
Toxics
Worst case and alternative release Flammables
Worst Case Release
Flammables
Estimate Rate
of Release
The first step of data collection depends on the process under consideration. Each
process has its own characteristics that need to be identified to conduct the hazard
assessment. Some of the process conditions are as follows:
13
Nature of the Chemical: The chemical needs to be classified as toxic, flammable or both
based on its flammable properties and toxic exposure levels.
Temperature of the Chemical: The temperature of the chemical can be either at the
inventory or during the process. Both these temperatures need to be noted as some
reactive chemicals might behave completely different under high temperatures.
Moreover at the modeling stage temperature becomes a very critical factor to estimate
other physical properties that can determine the release rate.
Quantity of the Stored Chemical: The quantity of chemical handled again is a very
important input for modeling the system. The quantity of substance can vary according
to the scenario under consideration.
Active and Passive Mitigation Systems Present at the Process Site: Active mitigation is
defined as equipment, devices or technologies that function with human, mechanical or
energy input. Sprinkler systems, water curtains, valves and scrubbers are considered
active mitigation systems. Passive mitigation systems are those that do not require any
human, mechanical or energy input like dikes, building exposures and containment
walls. The scenario selection procedure consists of two elements
vapor cloud, heat from a fire or blast waves from an explosion will travel before
dissipating to the point that serious injuries from short term exposures will no longer
occur. The number of worst-case scenarios that need to be analyzed depends on the
regulatory specifications on the toxicity and flammability of the substance.
Program 2 – Eligibility: If the distance to the endpoint in the worst-case analysis is equal
to or greater than the distance to any public receptor, the process can be included in
Program 2 or Program 3.And for all these processes one worst case analysis for the
regulated toxic substance and one for the regulated flammable substance needs to be
carried out.
Another critical assumption for alternative release scenarios is the presence of active
mitigation systems such as interlocks, shutdown systems, firewater and deluge systems.
It is required to analyze at least one alternative release scenario for each listed toxic
15
Estimation of Release Rate and Distance to Toxic End Point: This step is the most
critical step in estimating the source term and predicting downwind dispersion distance.
These models for a worst-case scenario and alternative case scenario have different
assumptions and these are summarized in Table1 shown below. The appropriate source
term and dispersion models that satisfy the U.S EPA criteria can be used to predict the
intensity of impact of these accidental releases.
Offsite Consequence Analysis for Aqueous Ammonia: The facility using SCR techniques
as discussed in the previous chapter store aqueous ammonia in storage tanks outside at
atmospheric pressure. The ammonia used could be anhydrous or aqueous in nature.
However anhydrous ammonia is being replaced by aqueous ammonia due to its more
hazardous toxic effects compared to the latter and the principal difference between these
chemicals in the context of atmospheric dispersion modeling is that the former chemical
evaporates as a pure vapor, whereas the latter escaped as an aerosol. Aqueous Ammonia
is stored in vertical or cylindrical storage tanks and it is usually placed inside a diked
containment area. The ammonia solution will be a water solution and stored close to
ambient temperature and the worst case release scenario involves the rupture of the
ammonia tank, instantaneous spilling of their contents into the containment dike and the
evaporation of ammonia from this surface. In this case the offsite consequence analysis
needs to be carried out on based on the guidelines discussed above. For comparing the
potential impacts associated with an accidental release of ammonia some of the “bench
mark” exposure levels need to be evaluated.
In the case of aqueous ammonia (concentrations less than 20%) the risks of reaching
these exposure values are much less when compared to anhydrous ammonia due to lower
partial pressures in solution. According to the U.S Environmental Protection Agency’s
(EPA) 112(r) clean air act, all regulated sources are required to conduct an offsite
consequence analysis for a worst-case release scenario involving regulated substances.
The analysis of potential consequences of an accidental release of ammonia needs to
meet the EPA‘s guidance on analyzing accidental chemical releases. The worst-case
assessment will also require an additional dispersion modeling for which the most
important input is the vaporization rate of the chemical. The following chapters will
discuss the modeling technique based on the EPA’s regulatory compliance procedures.
19
CHAPTER III
3.1 Introduction
Source term modeling is the most critical part of any consequence modeling procedure.
Most of the accidents involve release of hazardous chemical from the containment. The
imitating event for these accidents could be a pipeline or vessel rupture, hole in a tank,
pipe or a runaway reaction. The source term model basically gives us information
regarding the rate of discharge, total quantity discharged and the state of discharge. The
units used for source emission term will generally depend on the scenario or the type of
the release. Typically, emission rate units are expressed in terms of mass per unit time
per unit area.
Fig. 10. Generic source term modeling procedure (Hanna and Drivas 1987)
3.2 Background
The concept of detailed source term modeling started as early as the 1970’s.Several
steady state, semi-empirical evaporation models are available which were programmed
with minimum parameters. The models first developed only for pure component liquids
21
and were further developed over the years to suit chemical mixtures. Although many
concepts were derived, only a few of them were actually made into user-friendly
programs or simple equations that could be easily applied to practical situations. In this
chapter a few important pure component and multi-component models will be reviewed
to get a brief overview on the existing work in this field.
− 2/3
⎡ ⎤
⎢ ⎥
0.8 .0.9 − 0.8 ⎢ . + GMV ) ⎥
(31 0.33 2
⎢ T ⎜⎝ + ⎟ ⎥
⎣ 29 GMV ⎠ ⎦
The major limitations of this model are that it does not take into account the effect of
evaporative cooling, radiation and further assumes pool temperature is the same as air
temperature.
Q = K u0.8A (2)
However ALOHA is suited only for pure chemicals, a few select solutions and the
property information in its chemical library is not valid for multicomponent
mixtures. The model also displays the results only for 1 hour after the release
takes place and does not predict vaporization rates at different time intervals.
Moreover when an incorrect property value is used in ALOHA, the model's
release rate and dispersion estimates will not be valid.
assumes that the liquid pool spreads which is not necessary for instantaneous releases as
the liquid quickly occupies the entire diked area. Further the various assumptions made
in this model have not been well justified because of lack of experimental data for each
of the situations considered.
All the above-mentioned models can calculate the source strengths only for pure
chemicals or liquids. There are a few other multi-component models that were developed
after this and some of the important ones are as follows.
CHEMMAP package is basically built of a chemical database, chemical fates model that
estimates the distribution of chemical (as mass and concentrations), a stochastic model
25
that evaluates the probability of impact from a chemical discharge. This software runs
with Applied Science Associates ‘Geographic Information System (GIS). This model is
more suitable for spill modeling to aid the evaluation of consequences of chemical spills
on water.
Some of these assumptions actually point out some drawbacks in the model. One of the
major drawbacks of this model is that it is not capable of treating non-ideal mixtures that
behave very differently when compared to the ideal solutions. However the model again
does not take into account the temperature variation factor. Moreover the model is not
very user friendly. This model currently needs to be run only in DOS mode that makes it
difficult for the end user to interpret.
This model lays great emphasis on studying the effect of simultaneous spread and
vaporization from spills of hazardous liquids water. However this model mostly
concentrates on cryogenic spills like that of LNG. In this model it is assumed that the
26
initial gravitational potential energy is converted into kinetic energy. Fay also assumes
that an ice layer is formed beneath the LNG and this plays a key role in the heat transfer
mechanism. This model does not seem to be very applicable for instantaneous
multicomponent mixture spills on land.
m = 175.7vt 3 (3)
where ‘m’ is the mass of chemical that evaporates and ‘v’ is the mass flow rate and‘t’ is
the duration of spill.
All these models have some important drawbacks and more importantly they are not user
friendly. The new model will have a more comprehensive approach and at the same will
make some valid assumptions to make it generic in nature.
27
CHAPTER IV
4 MODEL DETAILS
4.1 Introduction
The liquid spill needs to be modeled with all the basic assumptions mentioned in the
Risk management plan for industrial facilities. In the case of a chemical spill, the
characteristics of the chemical have to be studied at all stages to develop the model.
• Storage at Low Pressure: The chemicals stored in the tanks are maintained at
very low pressure almost close to atmospheric pressure. Liquids stored under
high pressure above their normal boiling temperature present problems of
flashing (Crowl and Louvar 2001). For all other cases the spill is entirely in
liquid form and there is a very negligible amount in the vapor phase.
• Dynamics of the Release: The release mechanisms may form an important part
of source term modeling if chemical spill occurs due to a leak. However for
instantaneous release, assuming worst-case scenario the dynamics will not
affect the source term calculations.
• Level Liquid Spill Surface: The spill surface is flat without any dents or
irregularities that could affect the surface area term used in calculating the
mass flux.
• Ground Does Not Absorb Spills: Penetration of spilled material into the soil is a
complex phenomenon to be incorporated into the model. Gravitational forces
move the fluid into the soil at the same time capillary forces also influence the
downward movement of the fluid. In this case the concrete is assumed to be
impermeable to spilled chemical.
• Composition of Air Does Not Change Significantly: The composition of air
should actually change as a new chemical is added to it. The properties that
would be of interest will actually include cloud concentration, density and heat
capacity. However the volume of air just above the pool when compared to the
pool volume is in the range of 10-40 times more. Hence these property changes
will not be very significant.
30
NH3 + H2 O → NH4 OH
NH4 OH → NH 4 + + OH −
Moreover ammonia is not found to be very reactive under extreme temperature
conditions.
• Ground Temperatures Remain Constant: In most process industries, nowadays
the diked area is made up of concrete pads. Further in these calculations the
ground temperature is kept constant. This could be a valid assumption for spills
taking place very close to the ambient temperature. In other cases the ground
might actually cool down a bit and the heat flux to the pool will keep changing
with time.
Wind
Surface layer
Bulk
Spilled chemical 2.1 Heat from
Spilled chemical
Conduction from
Liquid Phase Resistance: The liquid phase mass transfer resistance affects the rate of
transfer of chemical from the bulk of the liquid to the interface and this phenomenon
takes into account the eddy and molecular diffusivities in the liquid (Mackay and
Matsugu 1973). For multi-component solutions there are two limiting conditions that
may exist, one in which there is an infinite diffusion rate, the vapor leaving the surface
will be in equilibrium with the bulk liquid. The other limiting condition is where there is
no diffusion and the concentration of the mixture remains constant with respect to time
and depth. In actual cases there exists an intermediate situation that drives the
evaporation.
32
Vapor Phase Resistance: The second phenomenon is the vapor phase resistance that
controls the rate of emission of chemical into the air. In this model vapor phase
resistance is given more importance and an empirical correlation will be used for the
mass transfer coefficient. This coefficient is a function of the molecular diffusion
characteristics and expressed in terms of the Schmidt number. Moreover the theoretical
modeling work by Dr Frie et al (Frie et al. 1992) shows that there is very little effect on
the calculated vaporization rates from liquid phase resistance. Hence the assumption that
vapor phase resistance is the driving force for mass transfer mechanism can be well
justified.
For a multi component mixture each component moves depending on their partial
pressures. Further the partial pressure changes with the temperature of the pool and
hence becomes a dynamic property. The bulk transfer of fluid from the pool to the air
flowing above it is given by
dm i
= K mi * A * MW * ( Psi − Pai ) / R * T pool (4)
dt
R = Gas Constant
mi =Mass flow rate of the components.
33
Here Kmi (mass transfer value) is an empirical value and can be calculated using the
Schmidt number as shown below.
d o =Pool diameter
u w =Wind velocity
The average Schmidt number Sci is a function of the diffusivity of the component in air
Sci = Va/ (da*Di)
The diffusivity Di of each component in turn can be calculated by the Wilkey Lee
method (Reid and Sherwood 1966).
Va =Viscosity of air
da = density of air
The mass transfer equation 4 shown above represents mass transfer limited liquid pool
evaporation. For liquid pools containing mixtures of liquids of varying volatility, the
light components will evaporate first.
Partial Pressure Estimation: The partial pressure of each component in solution will
depends on its interaction with other components in solution. For ideal systems it is
usually done using the Antoine equation where the molecules virtually do not interact
with other particles. In the case of non-ideal multicomponent mixtures like aqueous
ammonia the ideal system assumption will actually introduce a great deal of error in the
estimation of source term. The following graph will show how Raoult’s approximation
can overestimate the partial pressure of ammonia.
34
6
Non ideal
Ideal condition
0
0 20 40 60 80 100 120 140
Temperature (deg F)
The figure above (Figure 12) shows that the ideal solution approximation predicts values
with almost 200% error. This error will be carried on to the source term estimation if
Raoult’s law is applied for aqueous ammonia. Hence partial pressure data from standard
books or an empirical relation needs to be used.
However to calculate each heat transfer property we would need to calculate some
physical properties. All these basic or relevant heat transfer mechanisms shown in figure
11 can be explained in detail as follows.
35
Ambient Heat Convection: The rate of convective heat transfer from the pool to the air
above it is given by the following expression
h = Heat transfer coefficient, calculated from heat and mass transfer analogy
A= Pool Area
Tambient =Ambient temperature
The data on heat transfer coefficient ‘h’ from atmosphere to water is very limited. An
empirical relation for the heat flux due to convection derived from experimental
correlations developed for LSM-90 model (Cavanaugh et al. 1994) has been utilized in
this new model.
h = K mi * da * Cp * (Scav/Pr)0.67
da = Density of air
Cp = Specific heat of air
Pr = Prandtl number
Scav = Average Schmidt number
Where ‘ K mi ’ is the mass transfer coefficient, it is an empirical value that solves a steady
state atmospheric diffusion equation with power-law vertical velocity and eddy
diffusivity profiles.
Heat from Dike: The heat added from the dike can be included depending on the area
and the nature of the dike (insulated or non-insulated). It might not be very significant in
many cases. The energy transfer phenomenon across the solid wall boundary is usually
36
conduction. The equation shown below is basically derived from Fourier’s law of heat
conduction.
N
Qdike = ∑
[2 * k * (T
dike − T pool ) Ai ] (7)
i =0 ( ∏ *α ) * ( t − t i* )
Heat from the Ground: The heat added from the ground can be added depending on the
surface area and material of the ground. In this case energy is transferred across the
ground by means of conduction. The equation shown below is also derived from
Fourier’s law of heat conduction. The rate at which heat is transferred is based on major
assumption that the ground temperature does not vary throughout the vaporization
process.
Q ground = ∑
N [2 * k * (T ground − T pool ) Ai ] (8)
i =0 ( ∏ *α ) * ( t − t i*
The remaining terms in this expression have the same definition as mentioned in the
“Heat of dike” expression mentioned above.
37
Heat from Evaporation: As the pool loses chemical to the air above it, there is also some
substantial amount of energy that is carried away by this vaporizing chemical. This
directly causes the pool temperature to fall and this effect is known as evaporative
cooling.
This heat of vaporization value can be calculated using the Pitzer –Accentric
factor correlation. (Reid and Sherwood 1966)
Heat of the Material Added: This term is added to incorporate the energy transfer
associated with the material that spills into pool. This energy term is more
applicable for a system where there is continuous liquid spill over a substantial
period of time.
where,
mpool-input = mass entering the pool at any time.
Tstorage = Temperature of stored chemical
38
Radiative Heat Transfer: The heat transferred from the atmosphere, as a long wave
radiation and short wave radiation need to be incorporated into the energy balance
equation. However a constant value 130 Btu/hr-ft2 can be used (Studer et al. 1988).
Moreover it has been found that this term is very small in energy balance calculations
when compared to other heat flux terms.
Sensible Heat: The sensible heat added to the pool represents the accumulation term in a
general heat balance. This term balances the remaining heat flux terms in the net energy
balance equation. For highly volatile compounds, due to evaporative cooling effect this
term will be negative indicating the drop in pool temperature.
The pool evaporation model that will be developed in the next chapter will incorporate
these basic heat and mass transfer principles and will be improved step by step.
39
5 CHAPTER V
MODEL DEVELOPMENT
5.1 Introduction
At this stage, with all the basic mass and heat transfer mechanisms associated with the
liquid pool being identified, a systematic approach has been adopted to first develop a
basic pool evaporation model and then improvise it.
Qair + Qdike + Qevap + Qground + Qmass − add + Qrad + Qsensible + Qsun + Qtank = 0
Any small change in this term will affect the energy transferred across the air-pool
interface and this in turn will affect the mass of chemical vaporized from the pool.
Hence the model should track the pool temperature at every time instant.
40
All the three mentioned transfer mechanisms play a very important role in determining
the temperature of the pool at different time instants once the evaporation begins. In this
bulk evaporation method the liquid pool is treated as a semi-infinite slab and evaporation
is assumed to take place from the pool as a whole. There is no mechanism that accounts
for the mixing inside the pool. The heat transfer is based on average pool properties with
no temperature gradients across the pool.
The above expression is derived from the basic heat transfer equation written for the
bulk evaporation method. The new term ‘Qbulk’ represents the energy transferred across
the bulk – pool boundary. As small amounts of the components vaporize from the top of
the pool, they are immediately replaced by the contents of pool below it .The process
involves some amount of energy transfer from the bulk to the surface. The surface of the
pool acts like a blanket and is at different temperature from the bulk of the pool and this
bulk heat transfer term is also critical in determining the actual energy transfer that
occurs. The expression for bulk heat transfer is represented as
‘hbulk’ is the heat transfer coefficient and can be estimated using the following relation
derived by Kawamura and Mackay.(Kawamura and Mackay 1987)
hbulk = k liq / (h * ϕ )
The liquid resistance factor takes into account the relative role of conduction heat
transfer vs. eddy transfer. If conduction dominates ϕ tends towards 0, otherwise ϕ
tends towards 1 and conduction becomes important. For most non-cryogenic materials it
has been observed conduction starts dominating after the first few minutes. (Raj 1980)
In this model, solving the equation for mass transfer, surface and the bulk heat transfer e
simultaneously will solve the energy balance equation. All three mechanisms that take
place are coupled and need to be solved iteratively to estimate the mass vaporization rate
and the temperature gradients in the pool. The most important common term in all three
equations will be the temperature of the pool and this varies with time. A small change
in the bulk temperature will affect ‘Qbulk’ and will in turn affect the mass of the contents
43
vaporizing from the top. These equations being time dependent can be solved for any
amount of time provided the whole of the spill does not evaporate before the iteration
time.
The equation 14 has been derived from the basic heat transfer equations for stationary
sources (2000). However the pool surface has been defined as an infinitesimal layer or as
a fluid layer of negligible thickness and this assumption will help us eliminate the
∂Ts
−k term and the above written equation reduces to:
∂z
∂Ts
ρc p. h. = Qsensible + Qbulk + Qair + Qsol − mH v (15)
∂t
44
All the energy terms other than the flux term are the constant heat source terms. This
differential equation can now be solved with a pseudo equilibrium assumption. This
assumption will hold good for very small time steps.
k =Thermal conductivity
h = Initial depth of pool.
The above written equation for the bulk of the pool can be solved with suitable boundary
and initial conditions.
Initial Condition
At t=0, Tbulk (z, 0)= Tsurface
Boundary condition
∂Tbulk
At z=0, − K = Qground
∂z
∂Tbulk
At z= h, − K = Qbulk
∂z
45
These partial differential equations will be solved with respect to time and pool depth
and the bulk temperature will be represented as a function of both. Now the plane just
below the surface will be considered to average out the bulk temperature instead of the
whole pool. This will further enhance the flux calculations across the surface –bulk
boundary and will help predict better source terms.
6
46
CHAPTER VI
6.1 Introduction
The improved model will now be demonstrated and validated to predict source term
values under different set of conditions. The whole validation process will be divided
into four stages.
The time dependent unsteady state approach will help analyze the results in a more
detailed fashion. The dynamic code was run under the following conditions and the
following outputs were generated for the pool.
25
Mass flux of ammonia( g/sq.cm sec)
20
15
10
0
0 0.5 1 1.5 2 2.5 3 3.5
Tim e (hours)
The vaporization rate of ammonia is clearly diminishing with time (Figure 13) and this
behavior is consistent with the fact that vapor pressure decreases as evaporation
proceeds. The plot of cumulative mass of both chemicals ammonia and water (Figures
14 and 15) are shown below to determine the chemical composition in air.
49
7.00
6.00
5.00
Mass of ammonia in air(Klbs)
4.00
3.00
2.00
1.00
0.00
0 1 2 3 4 5 6 7
Time(hours)
120.0
Mass of water vapor in air (lbs)
100.0
80.0
60.0
40.0
20.0
0.0
0 0.5 1 1.5 2 2.5 3 3.5
Time(hrs)
299
298
297
Surface Temperature(deg K)
296
295
294
293
292
291
290
0 0.5 1 1.5 2 2.5 3 3.5
Time(hours)
It can be seen from the above plots (Figures 16 and 17) that surface temperature and the
bulk temperature drop with time and this can be directly attributed to the evaporative
cooling effect.
30
25
Composition of ammonia in pool (%)
20
15
10
0
0 1 2 3 4 5 6 7 8
Time(hours)
The graph shown above (Figure 18) depicts the composition of ammonia inside the
liquid pool. The ammonia percentage gradually decreases in the pool. However
throughout the run very little water actually evaporated to the air. This could be
attributed to humidity of the air. For really moist conditions the vapor pressure of water
vapor in the air is much greater than the vapor pressure in the pool. However the
52
moisture will condense and drop back into the pool only if the air temperature equals the
dew point temperature at that pressure.
The results for vapor pressure of ammonia and mass vaporization flux of ammonia with
time are as follows. The first graph (Figure 19) shows the trend of partial pressure of
ammonia in solution with time. The initial partial pressure of ammonia in this case can
be compared to the second case with a lower spill temperature. For a spill at 300 K, the
initial partial pressure of ammonia is very high and this in turn causes more ammonia to
escape from the solution.
53
0.45
0.40 Spill
Temperature:290 K
0.35
Partial Pressure (atm)
Spill
0.30 Temperature:300 K
0.25
0.20
0.15
0.10
0.05
0.00
0 0.5 1 1.5 2 2.5 3 3.5
Time (hours)
0.250
Spill Temperature:290 K
0.200
Mass Flux *E+03 (g/sq.cm sec)
Spill Temperature:300 K
0.150
0.100
0.050
0.000
0 0.5 1 1.5 2 2.5 3 3.5
Time (hours)
However as time progresses, as the ammonia content comes down drastically in the pool
in the higher temperature case, the partial pressure drops down drastically (Figure 20).
For the lower temperature the mass flux is steady and hence the drop in the partial
pressure is also steady.
7.00
6.00
5.00
Ammonia in air (Klbs)
Wind:2 m/s
4.00
Wind:4 m/s
3.00 Wind:5 m/s
2.00
1.00
0.00
0 0.5 1 1.5 2 2.5 3 3.5
Time (hours)
In this model, Sutton’s (Sutton 1953) power law will be applied to the velocity profile.
According to this model, velocity is a function of ground roughness and ambient
temperature profile .In the graph (Figure 21) shown above, we can see a clear
comparison of ammonia entrained at 3 different wind speeds. As the wind speed
increases the mass of ammonia entrained also increases. This is consistent with the
power law equation 4 for mass transfer.
In fact the wind speed becomes more critical for actual dispersion estimates or toxic
endpoint distance. When the toxic plume is being modeled it is observed that as the wind
increases the plume becomes longer and narrower. The substance travels faster in air but
is also diluted by a larger quantity of air (Crowl and Louvar 2001).
56
0.25
Ammonia:20%
0.20
Ammonia:10%
Partial Pressure (atm)
0.15
0.10
0.05
0.00
0 0.5 1 1.5 2 2.5 3 3.5
Time (hours)
The partial pressure graph or Figure 22 shows how partial pressure changes with time
for a mixture of 20% ammonia and 10% ammonia. The comparison shows that the
partial pressure value for the more concentrated mixture is always higher than that of
dilute one. This result is different from the one discussed above for effect of initial spill
temperature. This indicates that for aqueous ammonia both energy and mass transfer
effects are equally important in the initial phase, however with time mass transfer
becomes the driving force for evaporation.
299.0
298.0
297.0
Temperature (deg K)
Ammonia:20%
296.0
Ammonia:10%
295.0
294.0
293.0
0 0.5 1 1.5 2 2.5 3 3.5
Time (hours)
The second graph or Figure 23 depicts the effect of composition of temperature of the
pool. The evaporative cooling effect is more prominent in the case of a concentrated
mixture. This result can be justified by the fact that aqueous ammonia has high enthalpy
of vaporization and with increase in concentration these value further increases. Hence
the strong influence of non-ideality is again shown at this stage.
At half hour intervals the solution temperature, concentration, air velocity across the
front and back of the pan and the average humidity were recorded (Mikessell et al.
1991). The liquid pool evaporation models were run under the same conditions and
results (Tables 5 and 6) from both experiments were compared.
59
6.4.1 Experiment-1
800
700
600
Mass of ammonia in air (lbs)
500
400
Experimental
Bulk evaporation
300
Surface evaporation
New Model
200
100
0
0 1 2 3 4 5 6
Time (hours)
6.4.2 Experiment-2
Table 6. Predicted values for experiment 2
700
600
500
Mass of ammonia in air (lbs)
400
Experimental
New Model
300 Surface evaporation
Bulk evaporation
200
100
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time (hours)
From the above graphs (Figures 24 and 25) and tabulated values it is evident that the
improved liquid pool model predicts the closest values to experimental values and the
percentage deviation or error is clearly decreasing. In this case the error is defined as
follows:
The bulk and surface evaporation models provide over conservative results that might
drastically affect the dispersion calculations.
The improved model can now be picked and the validation can be further extended to
compare the surface temperature and the composition of the solution at different time
instants for both experimental trials.
62
305.0
Experimental
295.0
Temperature (deg K)
290.0
285.0
280.0
275.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
Time (hours)
25
Experimental
20
New Model
Weight Percentage (%)
15
10
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time (hours)
288
286 Experimental
New Model
284
Temperature (deg K)
282
280
278
276
274
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time (hours)
25
Experimental
New Model
20
Weight of ammonia (%)
15
10
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time (hours)
The graphs (Figures 26-29) depict the general behavior of the model. The composition
of ammonia in the pool and surface temperature were tracked throughout the simulation
to make sure that the system is consistent with the model theory.
Norman and Dowell (Woodward et al. 2002) conducted indoor and outdoor tests (Table
7). In the indoor test pans were suspended from a tripod and the pure chemical was
allowed to evaporate from the pan with an initial depth of 10 cm. The pure chemical,
ethanol has been used for this comparison. The chemical was tested for 90 minutes and
the source term was recorded. In this particular case, the improved model was run with a
single component in the input to compare it with available data.
0.7
0.6
0.5
Mass Flux (g/sq.cm sec)
0.4
0.3
Experimental
0.2 New Model
ALOHA
0.1
0
0 10 20 30 40 50 60 70
Time (minutes)
On comparing the experimental values (Table 8 and Figure 30) with the improved model
and ALOHA, it is very clear that ALOHA over predicts the values with an average
deviation of 42% .The new model predicts values much closer to the experimental data.
The experimental values are compared with the improved model and ALOHA (Table 9
and Figure 31), and again it is very clear that ALOHA over predicts the values with an
average deviation of 39%. The model predicts values much closer to the experimental
data.
At this stage the model is predicting the mass flux values with good accuracy. However
the validation of aqueous ammonia still shows that the average error can be further
reduced to predict values closer to the actual ones. The model’s prediction under more
realistic conditions depends on the atmospheric conditions during the spill and few
important points in this context also need to be addressed. Moreover some of the
limitations of the model will also be listed in the following paragraphs.
model predictions and the error associated with the dispersion model predictions might
just not be from the source term model alone. A few key meteorological parameters that
need to be analyzed are ambient temperature, wind speed, and relative humidity and
each of these parameters will affect the source term and the dispersion model used. The
applicability of our model under extreme conditions will be discussed with reference to
the extreme weather conditions in different parts of the U.S.
Relative 92 82 55 54
humidity (%)
Table 10 presented above gives us the extreme weather conditions picked from the
National weather service data.
Ambient Temperature: The model takes into account the ambient temperature conditions
in the heat and mass transfer calculations. For aqueous ammonia, it is usually stored at
ambient temperature and for all the above-mentioned temperatures; it will not have any
significant flashing effects. Hence the model should be able to estimate the chemical
release rate with good accuracy.
69
Relative Humidity: In this case for extreme conditions like Houston with more than 90%
humidity, the vapor pressure of water vapor can be greater than the partial pressure of
water at many stages during the evaporation and this could cause the condensation of
some water molecules back into the pool. This could affect the performance of the model
to some extent, as the model does not account for mass transfer into the pool. For places
like Los Angeles, with very mild climatic variations throughout the year, the humidity
never exceeds 60% and the model will work perfectly fine under these conditions.
Wind Speed: The model uses the Sutton’s power law equation. But at very high speeds it
has been proved with experiments that an empirical correlation with a logarithmic
profile for wind velocity is more suitable. The high wind speeds like that in Midland
could result in greater turbulence and more mixing in the air. The vapor cloud in this
case cannot actually be treated like a cylindrical vapor cloud and the. Moreover at high
wind speeds the plume becomes narrower and longer.
7 CHAPTER VII
MODELING SOFTWARE
This chapter is intended to help the user to get an overview of the code written to
develop this model. The non- linear and other differential equations handled in the
model need to be solved systematically and efficiently. The Matlab software has a lot of
built in functions that can be used to solve these mathematical equations.
MATLAB is a high-level technical computing language that provides tools for algorithm
development, data visualization, data analysis, and numerical computation. MATLAB
can be used to solve technical computing problems faster than with traditional
programming languages, such as C, C++, and FORTRAN. Some of the key features of
this software that helped develop this dynamic model are
The ‘eqns’ are symbolic expressions or strings specifying equations. The ‘vars’ are
symbolic variables or strings specifying the unknown variables. SOLVE seeks zeros
of the expressions or solutions of the equations. If no analytical solution is found
and the number of equations equals the number of dependent variables, a numeric
solution is attempted.
• INLINE - Object
The command ‘INLINE (EXPR)’ constructs an inline function object from the
MATLAB expression contained in the string EXPR. The input arguments are
automatically determined by searching EXPR for variable names. This function was
used to feed the input to dynamic partial pressure equation.
‘PLOT (X, Y)’ plots vector Y versus vector X. If X or Y is a matrix, then the vector is
plotted versus the rows or columns of the matrix, whichever line up. If X is a scalar and
Y is a vector, length (Y) disconnected points are plotted. The plots to show the time
dependency of vapor pressure, mass flux are done using this function.
• BVP4C
This function is used to solve boundary value problems for ODEs by collocation.
SOL = BVP4C (ODEFUN, BCFUN, SOLINIT) integrates a system of ordinary
differential equations of the form y' = f (x, y) on the interval [a, b], subject to general
72
7.2 Algorithm
A source code was written in MATLAB 6.5 to calculate the vaporization rate at
different time instants after a chemical spill and the cumulative mass for each
component. The following algorithm will give a brief overview of the code.
Start
Step-1:
• The input and other program variables are declared and initialized.
• The variables used in the loops are initialized to zero at the beginning.
• The actual input will include initial concentration, temperature, area, pool
mass and some physical properties.
• All the other physical properties are also loaded into the Matlab file.
Step-2:
• The iteration (time loop) for calculating all those properties that changes with
time starts here.
• Some of the average properties are defined, including vapor pressure.
Step-3:
• Iteration for tracking all the components in solution begins here
Step-4:
• The mass and heat transfer equations are defined here and an initial solution
is determined.
Step-5:
73
• The temperature and the mass flux terms are interdependent and need to be
solved in a non-linear fashion.
• The routines newbvp, twobc and twode are separate functions that are called
in this main program.
Step-6:
The feedback is then given to start of the time loop to recalculate the properties.
Step-7:
• In this step all output plots are generated using the 2-d or 3-d plot function to
display the time varying mass flux and other key results.
End
The input to the code can be entered manually or using an interactive software like
Microsoft Excel. The main code that was written in MATLAB is presented below
(Figure 32).
global check;
global condensationtemp; %% used in routine “cloud” as condensation temperature
74
pn=pnh3(Ts1,wnh3)/760;
pntemp(j)=double(pn);
aa(j)=double(ph);
bb(j)=double(pn);
rv=double(ph)
x=double(ph-.0167)% flag to check humidity difference.
Da(1)=0.2161;
Da(2)=0.2164;
viscosityair=.158;
Sc(1)= viscosityair/Da(1);
Sc(2)=viscosityair/Da(2);
MW(1)=18;
MW(2)=17;
R=82.05;
u= 5*3600;
Length=((areanew*4/3.14)^.5)/100;
km(1)=(2.7*10^-2*(u^.78)*(Length^-.11)*(Sc(1)^-.67))/36;
km(2)=(2.7*10^-2*(u^.78)*(Length^-.11)*(Sc(2)^-.67))/36;
if x>0
flag=flag+1;
mm=double(km(1)*MW(1)*(ph-0.0167)/R/Ts1+km(2)*MW(2)*(pn-0)/R/Ts1);
mm= double(1.516*10^(-1)*(pn-0)/Ts1)+double(1.645*10^(-1)*(ph-.0167)/Ts1)
mass(j)=double(mm);
else
mm= double((km(2)*MW(2)*(pn-0)/R/Ts1));
mass(j)=double(mm);
end
syms Ts
S=(2.079*10^-6*(302-Ts)+3.17*10^-2*(Tb-Ts)+(9.81*10^-3)-deltasol*mm);
k=solve(S,Ts);
Ts1=k;
m1=double(km(2)*MW(2)*(pn-0)/R/Ts1);
m2=1.58*10^-3*(ph-.0312)/Ts1;
end %% This end corresponds to the 'h' loop.
global Tsurface
Tsurface=double(Ts1) %%Tsurface is just copied to a new variable
Tsurface
ts1(j)=Tsurface;
weightnh3(j)=wnh3;
newbvp % New bvp function call for differential equation solution
matrixrv(j)=vijay;
for loop=1:5
newarray(TIME,loop)=array(TIME,loop);
end
tb1(j)=Tb;
END
8 CHAPTER VIII
The Selective Catalytic Reduction (SCR) technology is now being applied to a wide
variety of emission systems and has in turn the increased demand of aqueous ammonia.
All the facilities using aqueous ammonia need to submit their updated Risk Management
Plan to the EPA. However most of them are still using the anhydrous ammonia source
term models to model the release scenario for aqueous ammonia. Appendix 1 discusses
the difference in behavior of the dense gas clouds for both these compounds and justifies
the use of different source term models.
This new improved model is capable of handling multi-component non- ideal mixtures.
The relevant mass and heat transfer mechanisms have been incorporated into the model
with particular emphasis on the temperature gradients within the pool. The transient
model was built in several stages and improved at every stage to predict source term
values with more accuracy. The model has also been validated extensively with
experimental data and is found to perform well even with pure chemicals .The
comparison of the model with publicly available models like ALOHA further proves its
reliability.
With improving dispersion modeling techniques, the need for more realistic source term
modeling is growing. The performance of the source term models clearly affects the
dispersion model in predicting downwind concentrations. Hence this robust, user-
friendly model will help process engineers predict the consequences from chemical
releases with more confidence.
There is a reasonably good potential to improve this model further with lesser
assumptions by using a slightly different approach. In this case the liquid pool can be
broken down into small meshes and the principle of Rayleigh Benard convection be used
79
to model the system. This work will help us predict the source term from all points in the
pool with more accuracy and feed this as an input to Computational fluid dynamic
(CFD) models. A more detailed discussion of this idea has been done in Appendix B and
requires further work to develop the idea into a fully working dynamic model.
80
REFERENCES
AWS. (1978). "Calculation of toxic corridors." AWS Pamphlet, Air Weather Service,
105, 57.
Byron, R. B., Stewart, W. E., and Lightfoot, E. N. (2000). Transport Phenomena, John
Wiley and Sons, New York.
Cormack, D. E., Leal, L. G., and Imberger, J. (1974). "Natural convection in a shallow
cavity with differentially heated end walls .1. Asymptotic theory." Journal of
Fluid Mechanics, 65(8), 209-229.
Crowl, D. A., and Louvar, J. F. (2001). Chemical Process Safety, Prentice Hall PTR,
Upper Saddle River, New Jersey.
EPA. (1998). "General guidance for risk management programs (40 CFR Part 68)." U.S.
Environmental Protection Agency, Washington, District of Columbia.
Frie, J. L., Page, G. A., and Buckland, A. C. (1992). "Improved evaporation model for
spills of multicomponent non-ideal liquids." Heat Transfer, AIChE Symposium
Series, 88, 344-349.
GEC. (1999). "NOx removal system by SCR process." Global Environment Center
Foundation, Available at http://nett21.gec.jp/CTT_DATA/AIR/AIR_3/html/Air-
115.html
Hanna, S. R., and Drivas, P. J. (1987). Guidelines for Use of Vapor Cloud Dispersion
Models, Center for Chemical Process Safety, New York.
Ille, G., and Springer, C. (1978). "The evaporation and dispersion of hydrazine
propellants from ground spills." CEEDO 712-78-30, Civil and Environmental
Engineering Development Office, Asheville, North Carolina.
81
Kawamura, P. I., and Mackay, D. (1987). "The evaporation of volatile liquids." Journal
of Hazardous Materials, 15(3), 343-364.
Klosse, K., and Ullersma, P. (1973). "Convection in a chemical vapor transport process."
Journal of Crystal Growth, 18(2), 167-174.
Kunkel, B. A. (1983). "A comparison of evaporative source strength models for toxic
chemical spills." Atmospheric Sciences Division, Air Force Geophysics
Laboratory, Hanscom, Massachusetts.
Mackay, D., and Matsugu, R. S. (1973). "Evaporation rates of liquid hydrocarbon spills
on land and water." Canadian Journal of Chemical Engineering, 51(4), 434-439.
Mary, E., Jones, R., and Overstreet, R. (1993). "Modeling hydrochloric acid evaporation
in ALOHA." National Oceanic and Atmospheric Administration, Seattle.
Mikessell, J. L., Buckland, A. C., Diaz, V., and Kives, J. J. "Evaporation of contained
spills of multicomponent non ideal solutions." American Cyanamid, Inc., Wayne,
New Jersey.
Pritchard, S. G., DiFrancesco, C. E., and Von Alten, T. R. (1995). "SCR catalyst
performance under severe operation conditions." CORMTECH, Inc., Durham,
North Carolina.
Raj, P. K. (1981). "Models for cryogenic liquid spill behavior on land and water."
Journal of Hazardous Materials, 5(1-2), 111-130.
Reid, R., and Sherwood, T. (1966). Properties of Gases and Liquids, McGraw-Hill Inc.,
New York.
82
Shaw, P., and Briscoe, F. (1978). "Vaporization of spills of hazardous liquids on land
and water." SRD R 100, United Kingdom Atomic Energy Authority Safety and
Reliability Directorate, Abingdon, Oxon.
Studer, D. W., Cooper, B. A., and Doelp, L. C. (1988). "Vaporization and dispersion
modeling of contained refrigerated liquid spills." Plants/Operations Progress
(Process Safety Progress), 7(2), 127-135.
Woodward, J. L., Krishna , K., and Bagais, A. (2002). "Validation and extension of pool
evaporation and flashing blowdown models." International Conference and
Workshop on Risk and Reliability, Jacksonville, Florida, 1-15.
83
APPENDIX A
A.1 Introduction
The modeling of ammonia release in anhydrous form and aqueous form needs to be
treated separately. However in many cases the models applied to anhydrous ammonia is
also used to model aqueous ammonia releases. The difference in modeling techniques
can be discussed with respect to the vapor cloud behavior. The vapor cloud resulting
from ammonia release will differ in its density depending on the amount and phase of
ammonia actually being transferred into air. The denser vapor cloud will slump more
easily to the ground compared to the buoyant cloud and this has a direct effect on the
toxic endpoint calculation. In most cases the dense cloud will become slowly buoyant as
it mixes with the huge volume of air over time. The density and temperature of the cloud
in the first few minutes after the release will help us check the nature of cloud and the
applicability of the model.
Though it’s a mass balance equation, moles will be used throughout the expression.
A.1.2.1 Mass Balance
[N vap
k + Nk
liq
]initial
[
= Nk
vap
+ Nk
liq
] final
(17)
where ,
vap
Nk = Moles of component ‘k’ in the vapor phase.
liq
Nk = Moles of component ‘k’ in the liquid phase
In this case the initial condition is when no chemical has evaporated into the cloud and
the final condition is fixed based on the time at which measurements are required.
In the first stage the cloud is assumed to be made up of fully vapor. This assumption will
liq
make the ‘ N k ’ term zero at the initial and final stages.
In this part of the problem, an ammonia- water vapor cloud above the liquid pool is
considered. Further the cloud is assumed to be homogeneous. An enthalpy balance for
each of the constituents entering the cloud is performed assuming pseudo equilibrium for
a given small amount of time ∆t. This equation will also take into account the initial
water vapor present in the cloud, hence the final temperature of the cloud at every time
instant is noted.
∑ ( N kvap hkvap + N kliq hkliq ) final + ∑ ( N kliq ∆ Hmix ) final = ∑ ( N kvap hkvap + N kliq hkliq ) initial + ∑ ( N kliq ∆ Hmix ) initial
(18)
85
where,
vap
hk = Vapor Phase enthalpy
liq
hk = Liquid Phase enthalpy
∆ Hmix = Enthalpy of mixing of ammonia and water per unit quantity of final solution.
The mass transfer and the heat transfer equations can be solved simultaneously for the
temperature of the cloud and using this temperature the density of the mixture can be
found out under ideal conditions.
If the temperature exceeds the dew point measurement at any time instant, then the
respective chemical will condense out of the mixture and there is a great possibility of
forming a dense gas mixture. However if there is no condensation, the full vapor cloud
mixture will be less dense and has the potential to lift off quickly.
A subroutine was written in Matlab to calculate the final density of the cloud and this
was called in the main code to plot a graph of density at different times. The following
section will show the results at two different compositions of aqueous ammonia.
the volume or the cloud initial mass and its surface area at the ground, where ‘ρ’ is the
initial mass density of the cloud. In our case we can consider the volume of cloud just
above the diked containment area. Now the results will show the effect of both the
ammonia concentration and cloud volume.
Test Conditions
0.95
Ammonia weight:28%
0.9
0.8
0.75
0.7
0.65
0.6
0.55
0.5
0 0.5 1 1.5 2 2.5 3 3.5
Tim e(hours)
The spill parameters used in this run are presented in Table 11 .The density of the cloud
in most cases was found not to exceed the density of the moist air. As the ammonia
solution gets weaker the density of vapor cloud goes down further justifying the fact that
lesser ammonia goes up into the air. This is clearly shown in the above figure. (Figure
33)
Further the cloud volumes were changed and the density of the cloud changed as
follows. The cloud volume is estimated based on a ratio between the pool volume and
the air above it and this ratio is called the cloud volume factor .So a cloud volume factor
of ‘5’ means that the air above the pool is five times the pool volume. As this factor
increases the dilution of the chemical in air should also increase.
88
It can be clearly seen from the above tabular column (Table 12) that the cloud becomes
less dense due to rapid mixing with the air above it. On the whole for most aqueous
ammonia spills the density of the vapor cloud stays below that of air density and this
depicts the need for different modeling technique for aqueous ammonia.
89
APPENDIX B
The transient model described above accounted for the temperature variation along the
depth but treats the pool of liquid as a stagnant layer of fluid .In this case conduction is
assumed to be the controlling factor. This assumption holds good for most of the liquids
modeled under EPA’s worst-case release scenario. However for cryogenic liquids, it is
believed that energy transfer across the pool occurs more by convection. The cryogenic
liquid problem goes beyond the scope of this work, however a new approach based on
Raleigh Benard convection theory will be demonstrated here to show its effect on our
case.
According to this theory, the fluid is enclosed between two boundaries maintained at
different temperatures. Actually when the components of the pool evaporate the same
amount is replaced to the surface as a result of natural convection. Natural convection
flow patterns in rectangular enclosures with low height to length ratios (aspect ratio)
have been experimentally determined. Cormack et al have studied the phenomena of
buoyancy driven convection in shallow lakes as a shallow cavity problem (Cormack et
al. 1974). In most of these shallow cavity problems formation of unsteady cells have
been observed which they attribute to turbulence (Klosse and Ullersma 1973) also did
some work with a gas in two-dimensional low aspect ratio rectangular enclosures with
two vertical walls at different uniform temperature and a corresponding linear
temperature gradient along the horizontal walls. In both these it has generally been
observed that in low aspect ratio enclosures the flow over most of the length is parallel
counter flow with fluid in the upper half of the cavity flowing from warm end to cold
end.
90
For our case, the chemical spilled is enclosed in rectangular or cylindrical enclosures and
the aspect ratio is usually in the range of 10-3 or less than that. In these cases steady
convection begins in the form of two-dimensional rolls and as time progresses the
vertical boundaries will have a stabilizing effect due to additional viscous shear that
dampens thermal instabilities.
In this Raleigh Benard mechanism, the free convection may be actually opposed by
viscous forces and may therefore occur only if the temperature difference for vertical
boundaries is greater than some critical value. The parameter used to check the criticality
is called the Raleigh number expressed as follows (Byron et al. 2000):
where,
Gr= Grashof Number
Pr = Prandtl number
Tg =Temperature of the ground.
For Raleigh numbers below the critical value of 1101, the fluid is stationary and
conduction is the mode of heat transfer through the static fluid, otherwise convection
cells exists inside the pool. In this new approach the three dimensional Navier Stokes
equation needs to be solved to obtain the mass flux from the pool.
The equations of mass, momentum and heat transfer can be derived from the Navier-
Stokes equation and these can be represented as follows:
91
→
∇ .u = 0
∂u u x .∂u x u y .∂u xy uz .∂uz
+ + + = − ∇ . P + Pr. Ra. Tpool + Pr. ∇ 2 u
∂t ∂x ∂y ∂z
∂Tpool u x . ∂Tpool u y . ∂Tpool uz . ∂Tpool
+ + + = ∇ 2 Tpool
∂t ∂x ∂y ∂z
The problem can be solved with many modified assumptions compared to the improved
surface evaporation model. The pool can be divided into fine meshes using GAMBIT as
shown below. The structured grid (Figure 34) can then be used in fluent after setting the
necessary boundary conditions.
92
The ground temperature is assumed to be varying and the dike temperature is assumed to
be equal to the ambient temperature. The ground surface is set at a higher temperature
compared to the ground temperature. And an initial heat flux was also specified for the
groundside. The whole system was solved for energy and momentum balance with an
initial kick off velocity of 0.001 m/sec. The program was run for 250 seconds to check if
there is any significant temperature gradient along the horizontal section of the pool. The
results are presented as follows for two different runs.
93
B.1.4 Run 1
The following picture (Figure 35) shows the temperature variation across the top surface
and clearly shows that there is some temperature gradient initially under the set
conditions.
This picture again (Figure 36) shows that there is not any significant temperature
gradient across the surface in the bottom section of the pool.
B.1.5 Run 2
In this run the conditions were modified slightly by increasing the heat flux and
decreasing the thickness of the liquid pool. The result for the bottom surface clearly
shows that the Benard cells are starting to appear. This indicates that there could be
some small changes in temperature at the bottom with time.
95
Hence for the first case (Figure 37) there is no vigorous mixing and hence the
temperature difference across the pool can be assumed to be constant. But the second
case (Figure 38) needs to be tested for more time instants .A separate routine needs to be
written and attached to FLUENT with a specified profile for each of these properties and
this will make the system more dynamic. This new system can then be tested for other
chemicals mostly cryogenics and the 3-d model can be used to estimate the source at
every point inside the pool. This in turn will be a very realistic input for the
Computational Fluid Dynamic models.
97
APPENDIX C
The different water solution streams solvent water, stripper overheads, deep well flow and
rich water comprise of many hazardous chemicals in extremely low concentrations. The
spills of such dilute contaminated solutions might result in a liquid pool with some
volatile/toxic components. The objective of this study was to identify the best practices
used in the industry to deal with low concentration chemical spills and check the
applicability of these models to such chemical spills.
The various models that were studied as a part of the literature survey are the following.
• ESL evaporative Model
• Army Model
• Shell Model
• AWS Model
• CHEMMAP
• LSM-90/LPOOL
The LPOOL model theory was found to be the most suitable one to the problem and hence
the model was reprogrammed in Matlab with some simplified assumptions. The Matlab
program was then linked to an excel file which is very simple to handle. The whole package
is very interactive and the program is also very robust, as it has been tested with a variety of
streams.
C.1 Nomographs
A set of nomographs were also developed which could be used by a process engineer to
read off the amount chemical that could have escaped from a given quantity of spill. A few
sample nomographs for some of the process streams are given as follows.
98
In the above graphs (Figures 39 and 40) the Y-axis represents the quantity of spill in gallons
and the X-axis represents the amount of chemical released into air. These graphs can be
used to check if the released chemical has exceeded the reportable quantity. These
nomographs could be used directly for a known spill quantity or if the conditions vary the
model can be rerun and a new set of graphs can be prepared.
100
VITA
Vijay Raghunathan was born to Mr.Raghunathan and Mrs. Rama Raghunathan on March
3, 1981. Vijay completed his schooling at P.S.Senior Secondary School, Chennai in the
year 1998. He pursued his career in chemical engineering at University of Madras and
earned his B.Tech degree in May 2002. Later Vijay entered Texas A&M University in
the year 2002 to obtain his master’s degree in chemical engineering. Vijay
Raghunathan’s permanent address is: