Nothing Special   »   [go: up one dir, main page]

Consequence Analysis of Aqueous Ammonia Spills Using An Improved Liquid Pool Evaporation Model

Download as pdf or txt
Download as pdf or txt
You are on page 1of 110

CONSEQUENCE ANALYSIS OF AQUEOUS AMMONIA SPILLS USING AN

IMPROVED LIQUID POOL EVAPORATION MODEL

A Thesis

by

VIJAY RAGHUNATHAN

Submitted to the Office of Graduate Studies of


Texas A&M University
in partial fulfillment of the requirements for the degree of

MASTER OF SCIENCE

December 2004

Major Subject: Chemical Engineering


CONSEQUENCE ANALYSIS OF AQUEOUS AMMONIA SPILLS USING AN

IMPROVED LIQUID POOL EVAPORATION MODEL

A Thesis

by

VIJAY RAGHUNATHAN

Submitted to Texas A&M University


in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Approved as to style and content by:

M. Sam Mannan Mahmoud El-Halwagi


(Chair of Committee) (Member)

Kenneth Hall Michael K. Lindell


(Head of Department) (Member)

December 2004

Major Subject: Chemical Engineering


iii

ABSTRACT

Consequence Analysis of Aqueous Ammonia Spills Using an Improved Liquid Pool


Evaporation Model. (December 2004)
Vijay Raghunathan, B.Tech, University of Madras
Chair of Advisory Committee: Dr. Sam Mannan

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

ABSTRACT ..................................................................................................................... iii

DEDICATION ..................................................................................................................iv

ACKNOWLEDGEMENTS ...............................................................................................v

TABLE OF CONTENTS ..................................................................................................vi

LIST OF FIGURES........................................................................................................ viii

LIST OF TABLES .............................................................................................................x

CHAPTER

I INTRODUCTION ................................................................................................1

1.1 Introduction ...............................................................................................1


1.2 Chemistry Involved in NOx Reduction.....................................................3
1.3 Ammonia Injection System.......................................................................4
1.4 Regulatory Implications of Ammonia Used in SCR Process....................5

II OFFSITE CONSEQUENCE ANALYSIS ...........................................................7

2.1 Introduction ...............................................................................................7


2.2 Risk Management Plan for Industrial Facilities ........................................7
2.3 Five-Year Accident History ....................................................................10
2.4 Offsite Consequence Analysis Procedure ...............................................11
2.5 Emergency Response Program................................................................15
2.6 EPA’s Effort on OCA for Anhydrous Ammonia/Aqueous Ammonia....16

III SOURCE TERM MODELING.........................................................................19

3.1 Introduction .............................................................................................19


3.2 Background .............................................................................................20
3.3 Pure Component Models.........................................................................21
3.4 Multi-component models ........................................................................24
3.5 Other Concepts........................................................................................25
vii

Page
CHAPTER

IV MODEL DETAILS...........................................................................................27

4.1 Introduction .............................................................................................27


4.2 Liquid Spill Characteristics .....................................................................27
4.3 Characteristics of Chemical Emitted.......................................................28
4.4 Modeling Assumptions ...........................................................................29
4.5 Theoretical Basis for Models ..................................................................30

V MODEL DEVELOPMENT ...............................................................................39

5.1 Introduction .............................................................................................39


5.2 Bulk Evaporation.....................................................................................39
5.3 Surface Evaporation ................................................................................40
5.4 Improved Surface Evaporation Model ....................................................43

VI MODEL VALIDATION AND TESTING .......................................................45

6.1 Introduction .............................................................................................46


6.2 Base Case ................................................................................................46
6.3 Sensitivity Analysis.................................................................................52
6.4 Validation Against Experimental Data ...................................................58
6.5 Validation for Other Chemicals ..............................................................64
6.6 Applicability of the Model Under Different Meteorological Conditions67
6.7 Limitations of the Model.........................................................................69

VII MODELING SOFTWARE .............................................................................70

7.1 Programming and Key Functions Used ..................................................70


7.2 Algorithm ................................................................................................72

VIII CONCLUSION AND FUTURE WORK........................................................78

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.

Simple Cycle system

SCR

Gas Turbine catalyst module Stack

Fig. 1. Simple cycle system (Pritchard et al. 1995)

Combined Cycle system

SH and HP SCR IP LP Evap and FWH


Evap

Gas Turbine catalyst module Stack

Fig. 2. Combined cycle system (Pritchard et al. 1995)

This thesis follows the style of Journal of Environmental Engineering.


2

Low Dust System


Heater NH3
SCR

Boiler Air Heater ESP FGD

Stack

Fig. 3. Low dust system (Pritchard et al. 1995)

High Dust System

SCR

Boiler Air Heater ESP FGD Stack

Fig. 4. High dust system (Pritchard et al. 1995)

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

Fig. 5. Catalyst modules (Pritchard et al. 1995)

1.2 Chemistry Involved in NOx Reduction


The chemical reactions that take occur in the presence of the SCR catalyst, NOx
reduction is as follows: (GEC, 1999)

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

Fig. 6. Reaction mechanism (GEC, 1999)


4

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 ⎯⎯⎯⎯

⎯→( NH4 )SO4


catalyst
NH3 + SO3 + H 2 O ⎯⎯⎯⎯

1.3 Ammonia Injection System


The aqueous ammonia stored in tanks is vaporized and mixed with air before it is
brought into contact with the flue gas. The untreated gas and the oxygen rich ammonia
mixture react with the gas in presence of a catalyst in the SCR chamber. The whole
process can be represented as follows (Figure 7).

Fig. 7. Ammonia injection process (GEC, 1999)

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

1.4 Regulatory Implications of Ammonia Used in SCR Process


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. For comparing the potential impacts associated with an accidental
release of ammonia some of the “bench mark” exposure levels need to be evaluated.

Threshold Limit Values: The American Conference of Governmental Industrial


Hygienists (ACGIH) has established threshold doses called threshold limit values
(TLV’s) for a large number of chemical agents. The TLV usually refers to airborne
concentrations that correspond to conditions under which no adverse effects are
normally expected during a worker’s lifetime. The exposure occurs only during normal
working hours, eight hours per day and five days a week. The following are the
definitions for all the TLV’s used to commonly evaluate toxicity levels (Crowl and
Louvar 2001).

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.

Occupational Safety and Health Administration -Permissible Exposure Level, (OSHA-


PEL): (OSHA, 1992) OSHA sets permissible exposure limits (PELs) to protect workers
against the health effects of exposure to hazardous substances. PELs are regulatory
limits on the amount or concentration of a substance in the air. They may also contain a
skin designation. PELs are enforceable. OSHA PELs are based on an 8-hour time
weighted average (TWA) exposure. This value for ammonia has been set to 50 ppm.

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.

Emergency Response Planning Guidelines, ERPG: The ERPG-2 represents the


concentration below which it is believed nearly all individuals could be exposed for up
to one hour without irreversible or serious health effects (Crowl and Louvar 2001). The
American Industrial Hygiene Association issues these ERPG values and for ammonia
the standard ERPG-2 value is 200 ppm.
7

CHAPTER II

2 OFFSITE CONSEQUENCE ANALYSIS

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.

2.2 Risk Management Plan for Industrial Facilities


Environmental Protection Agency’s (EPA) risk management program rule requires all
the industrial facilities to perform a hazard assessment to provide information to the
government and public about the potential consequences of a chemical release. Section
112(r) of the Clean Air Act (CAA) section directed the EPA to issue regulations for
facilities with large quantities of very hazardous chemicals to prepare and implement
programs to prevent the accidental release of these chemicals and mitigate the
consequences of any releases that can occur (EPA 2004). In short, the objectives of this
program are to prevent any serious damage to human health or environment and reduce
the impact of such accidents that could otherwise have a severe consequence.

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

Define the processes

The facility is not


covered by the rule

Is the Regulated
substance above the
threshold quantity?

Yes

No
The facility is subject to the
rule

Assign program levels to


covered process

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.

Process: Process again is defined as “any activity involving a regulated substance,


including any use, storage, manufacturing, handling, or onsite movement such
substances (EPA 2004).The flowchart needs to be followed to check if the facility is
covered by the rule and at this stage program levels can be assigned to the processes.
The program levels are defined as follows:
Program1: Processes that do not have public receptors within the distance to an endpoint
from a worst-case release and no accidents with specific offsite consequences within the
last five years are covered under Program1.The following points summarize the
eligibility criteria for Program 1 coverage.
• A brief description of the worst-case release scenario, which must specify the

vessel and substance selected as worst-case and modeling assumptions.

• Methodology used to determine the distance to endpoints.

• Data used to determine that no public receptor would be affected.


10

• 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.

2.3 Five-Year Accident History


The five-year accident history comprises of a detailed study of the effects of any
accidental releases of one or more regulated substances from a covered process in those
five years prior to the submission of a Risk Management Plan. The information on
accidents should be collected for all processes that take place at the facility irrespective
of whether its covered in Program1 or not.

2.3.1 Accident Reporting Criteria


There are some key points that need to be considered before actually including the
accident/incident in the five-year accident history list. (EPA 1998)
• The release must be from a covered process and involve a regulated substance
held above its threshold quantity in the process.
• The release must have causes at least one of the following
o Onsite deaths, injuries or significant property damage.
o Offsite deaths, injuries, property damage or environmental damage.
11

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

2.4 Offsite Consequence Analysis Procedure


The basic steps involved in an offsite consequence analysis based on EPA guidelines are
shown in a flowchart (Figure 9).
12

Gather Basic Data

Select Scenario

Toxics
Worst case and alternative release Flammables
Worst Case Release
Flammables

Estimate Rate
of Release

Define Distance to Endpoint of Concern

Fig. 9. Steps for consequence analysis (EPA 1998)

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.

Pressure of the Toxic/Flammable Chemical Used: Liquid pressure becomes very


important as most liquids when stored under pressure above their boiling point
temperature present substantial problems of flashing. The liquid leaking from a pipe or a
tank can partially flash into vapor and can also be explosive at times.

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

A Worst-Case Release Scenario: It is defined as the release of largest quantity of a


regulated substance from a single vessel or process line failure that results in the greatest
distance to the toxic endpoint. The distance to the toxic endpoint is the distance a toxic
14

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 1 – Eligibility: The worst-case release scenario analysis is applicable to all


covered process. For a process to be eligible for Program 1 the worst-case analysis
should indicate that there are no public receptors within the distance to an endpoint. If
there is more than one process that may qualify under this category a separate analysis
needs to be performed to verify its eligibility.

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.

Alternative Release Scenario: An alternative scenario must be picked with reasonable


assumptions and the credibility of the scenario must be justified. It is not necessary to
demonstrate greater likelihood of occurrence or carry out any analysis of probability of
occurrence while choosing the alternative scenario.

Some examples of alternative release scenario will be as follows (EPA 1998)


• Transfer hose releases due to splits or sudden uncoupling.
• Process vessel or pump releases due to cracks, seal failure or plug failure.

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

substance in Program 2 and Program 3, however no analysis is required for regulated


substances in Program 1.

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.

2.5 Emergency Response Program


For a facility having at least one Program 2 or Program 3 process, then it is
recommended to implement an emergency response consisting of an emergency
response plan, emergency response equipment procedures, employee training and
procedures to ensure that the program is up- to-date. If a chemical facility is planning to
respond to release of substances with their own employees, the response program must
consist of the following elements (EPA 1998).
• Standard procedures to inform the public and other agencies about the release
• Documentation of proper first aid and emergency medical treatment necessary to
treat human exposures, and
16

Table 1. Assumptions for modeling scenarios (EPA 1998)

Parameter Worst Case Alternative


Endpoints Toxic substances endpoints Toxic substances endpoints
• Specified in appendix A • Specified in appendix A of
of 40 CFR Part 68. 40 CFR Part 68.
Flammable substances Flammable substances endpoints
endpoints • Overpressure of 1 psi for
• Overpressure of 1 psi vapor cloud explosions.
for vapor cloud • Radiant heat of 5 KW/sq.m
explosions. • Lower Flammability Limit
for vapor cloud fires.

Wind speed A wind speed of 1.5 m/sec and F Site-specific meteorological


stability class unless the local conditions will be used.
conditions are more severe.
Ambient For toxic substances, highest Average temperature and humidity at
temperature and daily maximum temperature the chosen meteorological conditions
during the past 3 years and will be used.
Humidity average humidity for the site.
Height of release Ground level release for toxic Depends on the release scenario.
substances
Topography Urban or rural topography will be Urban or rural topography will be
used used.
Dense or neutrally Tables or models used for Tables or models used for dispersion
buoyant gases dispersion of regulated toxic gas of regulated toxic gas must account
must account for gas density. for gas density.
Temperature of For liquids, the highest daily Substances may be considered to be
substance maximum temperature or process released at the process or ambient
temperature whichever is higher. temperature.

• Procedures and measures for emergency response


• Training for all employees in relevant procedures.
• Procedures for using, inspecting, testing and maintaining emergency response
equipment.

2.6 EPA’s Effort on OCA for Anhydrous Ammonia/Aqueous Ammonia


EPA gives some guidance to perform the OCA for regulated substances like anhydrous
and aqueous ammonia.
17

Anhydrous Ammonia: In general anhydrous ammonia is stored as a liquid under pressure


to be used in refrigeration systems. If the temperature and pressure are sufficiently high,
and if there is a sudden release of ammonia, the release mixture will be two phases, a
mixture of vapor and fine liquid droplets. The EPA thus conducted a comparative study
for a worst-case release of anhydrous ammonia at rural site. The toxic endpoint for
ammonia as specified in any RMP 200 ppm was taken as reference.

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.

Exposure Values for Ammonia: (NIOSH 2001)


• Immediate Damage to Life and Health, IDLH: 300 ppm
• Threshold Limit Value-Time Weighted Average, TLV- TWA: 25 ppm
• Threshold Limit Value-Short Term Exposure limit TLV STEL: 35 ppm
Emergency Response Planning Guidelines, ERPG-2: 50 ppm
• Occupational Safety and Health Administration -Permissible Exposure Level,
OSHA- PEL: TWA 50 ppm
18

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 SOURCE TERM MODELING

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.

The accidental releases of chemicals can further be classified as gas or liquid,


instantaneous or continuous, from storage tanks or pipelines, refrigerated or pressurized,
on land or water, confined or unconfined (Hanna and Drivas 1987). The four basic steps
involved in determining the source emission rate are as follows:
• Determining the time dependence of release scenario
• Identifying the most applicable source model
• Gathering specific input data and physical properties necessary for modeling.
• Calculating the source emission rate.
The following flowchart in figure 10 will depict the important steps mentioned above
and any source term modeling procedure will have to follow these guidelines for a
systematic approach.
20

Determine type of emission scenario

Instantaneous Continuous Continuous


spill time varying steady state

Selecting most appropriate source


term model

Gas jet Liquid jet Pool Cryogenic Two- phase


evaporation evaporation flow

Gathering input for the chosen model

Physical & Ambient conditions Source geometry Surface


chemical characteristics
properties

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.

3.3 Pure Component Models


3.3.1 Ille and Springer Model (Ille and Springer 1978)
This model is considered to be the most sophisticated of all the models developed during
the 1980’s.In this model the evaporation rate is computed by the method of Mackay and
Matsugu and is function of the concentration driving force as determined from vapor
pressure. The liquid pool temperature is also assumed to be equal to the ambient
temperature initially and is determined from a steady state energy balance equation at
any other time. The model also assumes evaporation of a pure liquid and ideal gas
behavior of the vapor film. These assumptions give a highly overestimated vaporization
rate for the spilled liquids. The only special feature of this model was that the
temperature inside the pool was tracked throughout the evaporation time and mass
transfer coefficients were calculated based on that.

3.3.2 Army Model (Kunkel 1983)


This model was adapted from the Chemical Engineer’s Handbook .The methodology
used here computes source strength as a function of the Reynolds number, the Schmidt
number, velocity of air and the vapor pressure of the pure chemical. The model gives a
very approximate equation for determining the source term for liquids with Reynolds
numbers less than 20,000.
22

− 2/3
⎡ ⎤
⎢ ⎥
0.8 .0.9 − 0.8 ⎢ . + GMV ) ⎥
(31 0.33 2

Q = 0.3u A T (GMW )(Vp) (1)


⎢ 0.5 ⎛ 1 1 ⎞ ⎥
0.5

⎢ T ⎜⎝ + ⎟ ⎥
⎣ 29 GMV ⎠ ⎦

Q = Evaporation rate (kg/hr)


u = Wind speed (m/sec)
A = Area of spill (sq .m)
GMW = gram molecular weight of liquid
GMV = gram molecular volume of liquid at normal boiling point (cu.cm/g mole)
Vp = Vapor pressure of liquid (mm-Hg)
T = Temperature of the liquid (deg K)

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.

3.3.3 Shell Model (Fleischer 1980)


This is an unsteady state model that deals with the estimation of evaporation rates of
pure chemicals. The model estimates the vaporization rate as a function of time. Even
though it is a dynamic approach there is no heat transfer effect due to evaporative
cooling and radiation. And moreover this model cannot be reduced to a simple
expression like the Army model.

3.3.4 AWS Model (AWS 1978)


The air weather service (AWS) gives an empirical equation based on the laboratory
studies for dinitrogen dioxide.

Q = K u0.8A (2)

Q = Source strength (kg/hr)


23

u = Wind speed (m/sec)


A = Area of spill (sq .m)
K = A constant which is equal to the vapor pressure of the chemical.
The major deficiency of this model is its lack of dependency on temperature.

3.3.5 ALOHA (Mary et al. 1993)


ALOHA (Area Locations of Hazardous Atmospheres) is a computer program jointly
developed by National Oceanic and Atmospheric Administration (NOAA) and U.S
Environmental Protection Agency (EPA) that uses information you provide it, along
with physical property data from its extensive chemical library, to predict how a
hazardous gas cloud might disperse in the atmosphere after an accidental chemical
release. This model uses the Kawamura Mackay relation for the mass transfer
mechanism. The ALOHA package can be used to model a chemical puddle scenario and
in this case the actual chemical, atmospheric conditions including surface roughness and
the site location are taken into consideration to calculate the emission term. This
information is then used to perform the dispersion the dispersion calculations.

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.

3.3.6 Raj Spreading Model (Raj 1981)


This model was developed in the early 80’s to address the cryogenic liquid spills on
water and land. In modeling spills on land the heat transfer rate is obtained from quassi
one-dimensional theory. This model addresses both continuous and instantaneous spills.
In both cases the model is derived by considering the hydrodynamics of the spread, the
varying heat flux to the spreading liquid and the relation between them. The model also
24

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.

3.4 Multi-component Models


3.4.1 Mikesell et al (Mikessell et al. 1991)
This model addresses the multi-component mixture spill and the model is transient in
nature. The heat and mass balance equations are solved simultaneously to estimate the
mass vaporization rate and the temperature. In this model various heat transfer
mechanism are incorporated, however there is no proper justification for some of the
assumptions. The model included condensation effects but it does not seem to show the
behavior of vapor cloud in air to decide if condensation actually occurs. Moreover this
model has not been programmed or made user friendly to enable to suit the end user.

3.4.2 CHEMMAP (McCay et al. 2003)


This model predicts the fate of a wide variety of chemicals products including floating,
sinking soluble and insoluble chemicals and product mixtures. In this model the state
and the solubility of the chemical are the most important factors. CHEMMAP is used for
chemical spills on water surface to simulate slick spreading transport, evaporation and
volatilization. The model actually uses the physical and chemical properties like density
to predict the fate of a chemical spill.

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.

3.4.3 LSM-90 (Cavanaugh et al. 1994)


This is a computer model and helps run simulations for vapor emissions from a
multicomponent spill. This model is actually capable of handling boiling and non-
boiling spills. The major assumptions of this model are that:

• All liquids and vapors are ideal


• The physical properties of the liquid will be averaged
• Multicomponent will always remain well mixed.

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.

3.5 Other Concepts


Some other important concepts, which have been adapted into a few computer models,
are described below.

3.5.1 Fay Model (Shaw and Briscoe 1978)

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.

3.5.2 Lind Model (Shaw and Briscoe 1978)


Lind worked with instantaneous spills of 4000, 10,000 cu.m of LNG and studied the
effects of chemical explosion that occurs after a large spill of LNG from a tanker. This
model used an empirical pool spread rate and mass vaporization rate.

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.

3.5.3 GASP Model


GASP is an acronym for Gas accumulation over spreading pools and was developed by
SRD/HSE. This model incorporates all the heat and mass transfer coefficients. However
the temperature dependence of the source term is not discussed in detail in this model.
Moreover the model has been converted to a computer program but is not very user
friendly or interactive with current day applications.

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.

4.2 Liquid Spill Characteristics


The liquid spilled on land will be emitted into the atmosphere depending on its physical
properties and some extraneous factors. The spilled fluid characteristics depend on the
following external factors:

4.2.1 Storage Tank


The chemical stored in facilities is usually stored in vertical or horizontal tanks. In some
cases it may also be transported in pipelines in between the actual process in which it is
involved. Source models are developed from fundamental equations that depict the
physiochemical process occurring during the release of the material. The phase of the
release also mostly depends on the storage tank conditions and the nature of release from
the tank.

4.2.2 Presence of a Containment Dike


Passive mitigation systems like a containment dike serve a very important purpose while
storing hazardous chemicals. As the liquid spills on the ground, they have a tendency to
spread randomly as the gravitational force tries to balance the inertial force. A
rectangular or circular dike around the storage tank can contain the spreading liquid.
Moreover the release rate of the chemical is also reduced as the surface area for
evaporation is reduced.
28

4.2.3 Characteristics of Spill Surface/Spill Geometry


The spill surface can be made of dry soil, concrete, insulated concrete and depending on
the nature of the surface the chemical pool behavior changes with time. In the case of a
contained liquid spill, the geometry of the bund is also important to model the heat
transfer associated with the pool.

4.3 Characteristics of Chemical Emitted


4.3.1 Chemical Composition
The vaporized chemical can be a pure chemical, binary mixture, multicomponent
mixture or even water solutions with very low concentration of hazardous chemicals.
Most of the physical properties of these chemicals keep changing with their composition.
For non-ideal solutions the physical property like vapor pressure are greatly influenced
by the composition of the mixture.

4.3.2 Ambient Condition


Atmospheric conditions also play a very important role in predicting the vapor cloud
behavior .The ambient temperature is a critical parameter for modeling the mixing of
hazardous chemical with air. The cloud density can vary depending on the ambient
temperature and this in turn will affect the end point concentration.

4.3.3 Spilled Chemical Conditions


Inherent properties like vapor pressure, density, and enthalpy of vaporization very
clearly depend on the conditions of the pool. The temperature especially becomes the
driving factor at one stage and the model needs to keep track of the temperature to
predict accurate mass vaporization rates.
29

4.4 Modeling Assumptions


The chemical spill to be modeled needs some basic assumptions to reduce the
complexity of the system. Most of these assumptions have been made in reference to the
U.S EPA worst- case scenario modeling conditions.

• 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

• No Physical or Chemical Reactions: For materials like ammonia, spills on land


usually do not give rise to any reactions. However if the same spill were to take
place on some water surface, ammonia undergoes a fast irreversible ionization
reaction. (Raj 1980)

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.

4.5 Theoretical Basis for Models


The basic modeling approach involves simultaneous heat and mass balances around the
pool. The chemical spilled on land will usually have a huge volume of air above it and
the ground at a particular temperature below it. Hence the mass emitted to the air above
clearly depends on the mass transfer mechanisms and the heat transfer mechanisms
involving the liquid pool, air and the ground. However in this case the energy and mass
transfer mechanisms are coupled and hence their equations also need to be solved
simultaneously. The figure 11 shown below clearly shows the various heat transfer
mechanisms around the pool that need to be analyzed.
31

Wind

Convection with air


Solar insolation Evaporation

Surface layer

Bulk
Spilled chemical 2.1 Heat from
Spilled chemical

Conduction from

Fig. 11. Liquid pool mechanisms

4.5.1 Mass Transfer Mechanism


Bulk movement of each volatile component will occur as air flows over the pool. The
mass transfer process can be explained in two parts.

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

K mi = Mass transfer Coefficient


A = Diked Containment area.
MW = Molecular weight of the components.
Psi =Partial pressure of the components in solution.

Pai =Vapor pressure of the components in air.

Tpool =Temperature of the pool.

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.

0.0292 * u 0w.76 * d o−.011 * Sci−0.67


K mi =
RTpool

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

Partial Pressure (atm)


5

0
0 20 40 60 80 100 120 140
Temperature (deg F)

Fig. 12. Comparison of non-ideality with Raoult’s law case

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.

4.5.2 Energy Transfer Mechanism


The volatile components of the pool evaporate as heat is added to the pool continuously
from various energy sources. The initial storage conditions of the pool (temperature,
pressure) suggest that the pool does not boil and mass loss takes place only due to
vaporization. The net energy balance around the pool can be represented by the
following equation.

Qair + Qdike + Qevap + Qground + Qrad + Qsensible + Qsun + Qmass-add = 0 (5)

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

Qair = h * A * (Tambient − T pool ) (6)

h = Heat transfer coefficient, calculated from heat and mass transfer analogy
A= Pool Area
Tambient =Ambient temperature

T pool =Pool 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* )

K = Thermal conductivity of ground


α = Thermal diffusivity
t = Time in seconds
t i* = Time section ‘i’ was first wetted in seconds

Ai = Area of the section of dike in contact with the liquid.

Tdike=Temperature of the dike

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*

Tground =Temperature of the ground.

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.

Qevap = ∑ (m j * ∆H vapj ) (9)

mj = Rate of evaporation of component j.


∆Hvapj = Heat of vaporization of component j.

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.

Qmass− add = m pool − input (Tstorage − Tpool ) (10)

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.

Qsensible = mpool cp (Tpool–Tpool-prev) (11)

Tpool-prev = Temperature of the pool at the previous time instance (t-1)

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.

5.2 Bulk Evaporation


As the chemical spills on to the land, it tries to equilibrate within a few minutes and after
which the various energy transfer mechanisms can be applied. At this stage most of the
critical heat transfer methods involve the temperature of the pool. The liquid pool can be
treated as a stagnant layer of fluid that interacts with the atmosphere and the ground. The
basic energy balance expression from equation 5 will be used again to discuss this
section

Qair + Qdike + Qevap + Qground + Qmass − add + Qrad + Qsensible + Qsun + Qtank = 0

5.2.1 Convection in Air


The pool temperature directly affects this term as shown in the following proportionality
expression
Qair ∝ (Tambient − T pool )

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

5.2.2 Ground Energy Transfer Term


Similar to the ambient convection term the interaction with the ground also has direct
relation with pool temperature. However, in this expression the ground temperature
remains constant and the heat flux will again depend only on the pool.
Q ground ∝ (Tground − T pool )

5.2.3 Evaporation Energy Term


This mechanism is the most important since it accounts for most of the heat loss from
the pool. The evaporative energy causes the pool to cool down continuously and this
effect becomes very prominent for non-ideal liquids where the heat transfer properties
are a strong function of temperature.

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.

5.3 Surface Evaporation


In this approach the whole pool will be divided into two layers, one being the surface
and the other known as the bulk. The surface of the pool is an infinitesimal layer of pool
from which evaporation actually takes place. The layer below this surface is of
considerable thickness and actually interacts with the ground and also supplies heat to
the surface above.
41

5.3.1 Pool Surface Conditions


The surface on top acts like a cap for the pool bulk and is at a slightly lower/higher
compared to the layer beneath it. The heat transfer equation for the surface can be
written as shown below.

Qsun + Qair + Qbulk + Qsensible = Qevap (12)

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

Qbulk = hbulk (Tbulk –Tsurface)

‘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 * ϕ )

Tbulk =Bulk pool temperature


Tsurface= Temperature of the surface.
kliq = Average thermal conductivity of the liquid.
h = Depth of pool at any instant of time.
ϕ = Liquid resistance factor.
The liquid resistance factor can be defined by the following expression.
42

ϕ = 1 + exp( − 0.6(Tboiling ,av − 34316


. ))) −1

Where Tboiling, av=Average boiling point of the mixture.

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)

5.3.2 Pool Bulk Conditions


The bulk of the pool, just below the surface layer will also have some energy
mechanisms associated with it. When the pool contents evaporate, they are replenished
to the top from this layer of the pool. However in this surface evaporation model, the
∂Tbulk
temperature in the bulk of the pool is averaged out and hence =0. The heat
∂z
balance in the bulk of the pool is given as follows in equation 13.

Qdike + Qground = Qbulk (13)

Qbulk = Bulk heat transfer due to conduction through liquid.

Tbulk =Temperature of the Bulk.

z = Depth of the pool.

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.

5.4 Improved Surface Evaporation Model


The slightly modified “surface evaporation approach” described above clearly shows the
temperature sensitivity of the model. Now the bulk of the pool beneath the surface can
be further divided into a number of planes of very small thickness. This approach will
assume that there is a temperature gradient along the depth of the pool. This allows us to
estimate the exact temperature of the bulk layer just below the surface layer and the heat
transfer across bulk- surface boundary can be further enhanced improved by using a
better temperature average.

5.4.1 Surface Conditions


In this approach, initially the energy balance for the surface can be represented with a
more basic equation where temperature of surface changes with time and depth.

∂Tsurface ∂Ts urface


ρc p. h. = −k + Qsensible + Qbulk + Qair + Qsol − mH v (14)
∂t ∂z

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.

5.4.2 Bulk Conditions


For the bulk case, the temperature can be tracked by accounting for the different
mechanisms as mentioned in the surface evaporation method. The major assumption in
this model is that conduction dominates after a certain amount of time. However this
model is different from the surface evaporation in that it accounts for bulk temperature
variation in the vertical plane.
∂Tbulk ∂ 2Tbulk
ρc p.h. = −k (16)
∂t ∂z 2

ρ = Bulk density of the pool


c p = Specific heat capacity of liquid

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

MODEL TESTING AND VALIDATION

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.

• Aqueous ammonia base case demonstration.


• Sensitivity analysis will be conducted to test the performance of model under
different conditions.
• Validation against experimental data will be performed to check the confidence
of the predicted values.
• The model will also be tested for other chemicals and their experimental results

6.2 Base Case


Aqueous ammonia spill on to a diked containment area will be modeled with a particular
set of parameters and the various results obtained from this simulation will be presented
to discuss the applicability of this model to our case. The spill temperature will be
assumed to be equal to the ambient temperature. The Table 2 presented in the next page
will summarize all the spill parameters.
47

Table 2. Base spill conditions

Aqueous Ground Ambient Air Relative Duration


Ammonia temperature temperature velocity Humidity of spill
weight (deg K) (deg K) (m/sec) (%) (hours)
(%)
28 297 302 5 50 3

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.

• Mass flux of ammonia vapor


• Cumulative mass of ammonia and water in air
• Temperature of the surface
• Bulk temperature gradients along the depth of the pool
• Composition of the pool
48

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)

Fig. 13. Mass flux of ammonia in air

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)

Fig. 14. Cumulative mass of ammonia in air

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)

Fig. 15. Cumulative mass of water vapor in air


50

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)

Fig. 16. Surface temperature

Fig. 17. Bulk temperature


51

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)

Fig. 18. Composition of ammonia in pool

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.

6.3 Sensitivity Analysis


The model will now be tested under different conditions with all the sensitive parameters
and suitable discussion will be presented at every stage to justify the working of the
model.

6.3.1 Initial Spill Temperature


The temperature of the chemical is probably the most important parameter that could
affect the nature of spill and the behavior of the spilled liquid once it begins to
evaporate. The model was run at two different initial spill temperatures and the
following results were observed. Table 3 shown below summarizes the spill parameters.

Table 3. Cumulative mass of ammonia in air


Weight % of Ground Spill Air Duration of
aqueous temperature temperature velocity spill (hours)
ammonia (deg K) (deg K) (m/sec)
28 297 302 5 3

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)

Fig. 19. Partial pressure variation


54

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)

Fig. 20. Mass flux of ammonia in air

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.

6.3.2 Wind Speed


Wind speed will play a very important role in the mass transfer equation. The velocity
profiles of the lower atmosphere have been studied in detail and it is found that for oil
spills, velocity of air above it is very critical in controlling the eddy diffusivity.
55

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)

Fig. 21. Effect of wind speed

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

6.3.3 Composition of Chemical


The weight percentage of the toxic chemical determines most of the physical properties
and the behavior of the mixture in case of a spill. For aqueous ammonia, if ammonia is
present in very small quantities the non-ideal equation used for the partial pressure
determination will not be applicable. Moreover the mass transfer correlation developed
by Mackay also needs to be modified for very dilute solutions. In the case of
concentrated aqueous ammonia solutions the mixture exhibits strong non-ideal
properties. The model was run for different composition of aqueous ammonia and the
mass flux and temperatures gradient trends were observed.

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)

Fig. 22. Effect of composition on partial pressure


57

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)

Fig. 23. Effect of composition on surface temperature


58

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.

6.4 Validation Against Experimental Data


Mikesell et al have carried two experiments for aqueous ammonia on a small scale and
the source term has been determined under a specified set of conditions. The general
spill conditions for these two experiments are given below in a tabulated form (Table 4).
In both these experiments 2240 grams of aqueous ammonia solution was spilled into a
pan and allowed to equilibrate for about 2 minutes.

Table 4. Spill conditions for Mikesell et al. data

Experiment Average Wind Average Initial spill Weight Area of


Ambient speed humidity temperature percentage spill
temperature (m/s) (%) (deg C) ammonia (sq cm)
(deg C) (%)
1 24.5 1.59 59 24.5 28.8 522
2 28.9 2.2 68 28.9 28.8 522

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

Table 5. Predicted values for experiment 1


Time Experimental Predicted values for
(hrs) data for Ammonia (g)
Ammonia (g) Bulk Error Surface Error New Error
± 15% model (%) model (%) model (%)

0.5 190 345 -81 270 -42 228 -20


1.0 240 498 -107 391 -63 339 -41
1.5 300 572 -91 462 -54 405 -35
2.0 350 616 -76 510 -46 449 -28
2.5 400 647 -62 546 -36 482 -20
3.0 420 669 -59 575 -37 507 -20
3.5 450 682 -52 599 -33 529 -17
4.0 490 683 -39 618 -26 547 -11

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)

Fig. 24. Experiment 1-Mass of ammonia in air


60

6.4.2 Experiment-2
Table 6. Predicted values for experiment 2

Time Experimental Predicted values for


(hrs) data for Ammonia (g)
Ammonia (g) Bulk Error Surface Error New Error
±15% model (%) model (%) model (%)

0.5 240 377 -57 325 -35 267 -11


1 300 503 -68 454 -51 378 -26
1.5 380 559 -47 518 -36 439 -16
2 420 593 -41 559 -33 480 -14
2.5 470 617 -31 590 -26 510 -9
3 500 634 -27 610 -22 534 -7
3.5 520 648 -25 628 -21 555 -7
4 550 654 -19 645 -17 571 -4
61

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)

Fig. 25. Experiment 2-Mass of ammonia in air

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:

Error = (Experimental Value – Predicted Value)/Experimental Value

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

300.0 New Model

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)

Fig. 26. Experiment 1-Temperature Gradients

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)

Fig. 27. Experiment 1-Composition of ammonia


63

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)

Fig. 28. Experiment 2-Temperature Gradients

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)

Fig. 29. Experiment 2-Composition of ammonia


64

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.

6.5 Validation for Other Chemicals


The model prediction also needs to be validated for chemicals other than aqueous
ammonia to check the generic nature of the model. At this stage two completely different
sets of experimental data will be used for this test.

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.

6.5.1 Test Conditions

Table 7. Test conditions for Norman Dowell data

Ambient Initial spill Air Relative Initial Area of


temperature temperature velocity humidity pool pan
(deg C) (deg C) (m/sec) (%) depth (sq.cm)
(cm)
20 20 0.1-1 35 10 113
65

6.5.2 Ethanol Data Comparison

Table 8. Predicted values for ethanol

Time Experimental data for Predicted values for


(minutes) Ethanol (g/sq cm sec) Ethanol (g/sq cm.sec) *1E-04
*1E-04 Model Error ALOHA Error
±10% (%) (%)

10 0.48 0.51 -6.7 0.63 -30


20 0.42 0.51 -22 0.58 -38
30 0.40 0.5 -28 0.58 -44
40 0.38 0.5 -35 0.56 -47
50 0.36 0.5 -42 0.54 -50
60 0.36 0.49 -39 0.52 -44

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)

Fig. 30. Mass flux of ethanol


66

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.

6.5.3 Propylene Data Comparison

Table 9. Predicted values for propylene

Time Experimental data Predicted values for


(minutes) for Propylene (g/sq cm.sec) *1E-04
Propylene (g/sq cm Model Deviation ALOHA Deviation
sec) *1E-04 (%) (%)
±10%
10 2.6 2.72 -4.6 4 -54
20 2.6 2.72 -4.6 3.7 -42
30 2.6 2.72 -4.6 3.7 -40
40 2.6 2.72 -4.6 3.5 -35
50 2.6 2.72 -4.6 3.5 -33
60 2.6 2.72 -4.6 3.5 -33
67

4.1 New Model

Mass flux (g /sq.cm sec)


3.9 Experimental
ALOHA
3.7
3.5
3.3
3.1
2.9
2.7
2.5
0 20 40 60 80
Time(hours)
Fig 31. Mass flux of dichloromethane

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.

6.6 Applicability of the Model Under Different Meteorological Conditions


The source term model developed in this work predicts emission rates much closer to
experimental data and hence reduced the deviation from actual values. These source
term values are in turn fed into a dispersion model to estimate the distance to the toxic
endpoint. The meteorological conditions might sometimes have a great effect on the
68

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.

Table 10. Meteorological conditions around the U.S.

Parameters Houston, Midland, Jersey City, Los Angeles,


Texas Michigan New Jersey California

Ambient 93.5 16.1 27.5 81.1


temperature
(deg F)

Relative 92 82 55 54
humidity (%)

Wind speed 6.2 11.1 11.4 6.6


(mph)

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.

6.7 Limitations of the Model


• The effect of spreading is neglected and can be important for some highly
viscous fluids.
• The ground surface is assumed to be dry and presence of water would cause
ammonia to react with water.
• The ground surface is assumed to be at a constant temperature with heat flux
depending only on the pool temperature.
• The mixing of the pool contents is not considered explicitly and this effect could
be important for some cryogenic fluids.
70

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

• Interactive tools for complex design problems.


• Mathematical functions for linear algebra, partial differential equation,
statistics, optimization, and numerical integration
• Graphics functions for representing data in multi-dimensions.
• Functions for integrating MATLAB based algorithms with external
applications, such as FORTRAN, Microsoft word, and Microsoft Excel.

7.1 Programming and Key Functions Used


The Matlab code was written to estimate the vaporization rate of the chemical and also
to keep track of the various important parameters involve din modeling. The temperature
inside the pool was solved using some PDE solvers. Some of the important functions
(MathWorks 2004) used in this code are:
71

• SOLVE - Symbolic solution of algebraic equations

A typical example of a SOLVE function is shown below

SOLVE ('eqn1',’eqn2’...'eqnN','var1, var2...varN')

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 - Linear plot

‘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

two-point boundary conditions of the form bc (y (a), y (b)) = 0. ODEFUN is a function


of two arguments: a scalar X and a vector Y. ODEFUN (X, Y) must return a column
vector representing f (x, y). In this ODE solver BCFUN represents the boundary
condition for the differential equation. SOLINIT specifies the initial condition for the
same equations.

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).

GLOBALIZATION AND INTIALIZATION BEGINS

global tt; %% Time variable used for pseudo equilibrium calculations.


global array; %% variable used for passing values from newbvp to main program ( Tbulk
values)
array=zeros(100,20); %% The array initial values are set to zero
global TIME;
global x ; %% Variable is used for moisture check.
global m1; %% mass of vaporization of ammonia
global j; %% Main looping variable
global m2; %5 mass of vaporization of water
xx=linspace(1,5,5); %% Depth is assumed to be 5 cm.
timetime=linspace(1,16000,20);
global Tb; %% Bulk temperature variable
global vijay; %% surface temperature variable.

Fig. 32. Matlab Code

global check;
global condensationtemp; %% used in routine “cloud” as condensation temperature
74

global pnh3final; %% Partial pressure of ammonia


global ph20final; %% Partial pressure of water
global kk; %%
global wnh3; %% weight fraction of ammonia
global check2;
global densitykg; % % Density of cloud in kg/m3
global yvnh3;%% mole fraction of ammonia in vapor
global yvh2o;
global areanew;
index=1;
Mpool=10^7; %% Intial spill mass
mixturedensity=.896;
userarea=(Mpool/mixturedensity)/4.2 ;
areanew= userarea;
matrixrv=300; %%
Ts1=302;% %Surface temperature, can be user input.
Tb=302; %% Bulk temperature, can be user input.
Tg=297; % %Ground temperature, can be user input.
wnh3= 28.80;%% Weight fraction of chemical
deltasol=302;
h2evap=0;
nh3evap=0;
p=1;
tt=600;
flag=0;

GLOBALIZATION AND INTAILIZATION ENDS

MAIN BODY OF CODE BEGINS

for j=1:18 %% Main Time Loop begins here


TIME=j;
Tbav=373.2-2.663*wnh3+.0133*wnh3^2;
Cpav=4.2-.006*wnh3+.000191*(wnh3)^2;
massammoniacumu(1)=0;
pnh3=inline('exp(17.17-4294/Ts+(.137*wnh3)-
(.00347*wnh3^2)+(23.01*wnh3/Ts))');%% vapor pressure of ammonia
ph20=inline('exp(18.37-4552/Ts+(.1477*wnh3)-(.00081*wnh3^2)-(43.53*wnh3/Ts))');
%%vapor pressure of water
for h= 1:10 %%Loop to check for convergence.
ph=ph20(Ts1,wnh3)/760;

Fig. 32. Continued


75

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;

Fig. 32. Continued


76

Mntemp=m1*areanew*tt; %% Mass of ammonia that is evaporated


nh3evap=nh3evap+Mntemp;
m2=double(km(1)*MW(1)*(ph-0.0167)/R/Ts1);% Cumulative mass
massammonia1(j)=m1;
masswater2(j)=m2;
massammoniacumu(j+1)=m1*areanew*tt+massammoniacumu(j);
if x>0
m2=double(km(1)*MW(1)*(ph-0.0167)/R/Ts1) %% RV needs to check these values and
replace with generic terms
Mhtemp=m2*areanew*tt;
h2evap=h2evap+Mhtemp;
xn=(Mntemp/17)/((Mntemp/17)+(Mhtemp/18));
xh=1-xn;
deltanh3=5562*((1-Tsurface/385)/.377)^.38;
deltah20=10598*((1-Tsurface/647)/.423)^.38;
deltasol=double(deltanh3*xn/17+deltah20*xh/17);
Mpoolnew=Mpool-(Mntemp+Mhtemp);
wnh3=100*((wnh3/100*Mpool)-Mntemp)/Mpoolnew
Mpool=Mpoolnew;
massammonia1(j)=m1;
masswater2(j)=m2;
time(j)=tt*j/60;
if wnh3<0
break;
end
else
time(j)=tt*j/60;
deltanh3=327.22*((1-Tsurface/385)/.43)^.38;
deltasol=double(deltanh3);
Mpoolnew=Mpool-(Mntemp);
wnh3=100*((wnh3/100*Mpool)-Mntemp)/Mpoolnew
Mpool=Mpoolnew;
end
weightpercent(j)=wnh3;
cloudhouston
temparray(index)=densitykg;
tempcloud(index)=kk;
tempyvnh3(index)=yvnh3;
tempyvh20(index)=yvh2o;
tempcheck(index)=check2;
tempcondensationtemp(index)=condensationtemp;
index=index+1;
end

Fig. 32. Continued


77

MAIN BODY OF CODE ENDS

PLOTS & FIGURES


for nb= 1:j
masscumuammonia(nb)=massammoniacumu(nb+1);
end
dlmwrite('alpha.xls',temparray,'\t',0,0)
plot(time,massammonia1,'-.b*','Linewidth',2)
ylabel(' mass flux ()')
xlabel('Time in Minutes')
plot(time,tempcloud)
ylabel(' condensation temp (K)')
xlabel('Time in Minutes')
figure;
plot(time,temparray)
ylabel(' density of air (g/cm3)')
xlabel('Time in Minutes')
SURF(xx,time,newarray)

END

Fig. 32. Continued


78

8 CHAPTER VIII

CONCLUSION AND FUTURE WORK

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.

Cavanaugh, T. A., Siegell, J. H., and Steinberg, K. W. (1994). "Simulation of vapor


emissions from liquid spills." Journal of Hazardous Materials, 38(1), 41-63.

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.

EPA. (2004). "Changes to chemical accident prevention rule in 2004." Chemical


Emergency Preparedness and Prevention Office, Chemical Emergency
Preparedness and Prevention Office, Washington, District of Columbia.

Fleischer, M. T. (1980). "SPILLS-an evaporation/air dispersion model for chemical


spills on land." Shell Development Company, Houston, Texas.

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.

MathWorks. (2004). "Matlab." The Mathworks, Inc., Available at


http://www.mathworks.com/

McCay, D. P. F., Isaji, T. (2003). "Evaluation of the consequences of chemical spills


using modeling: chemicals used in deepwater oil and gas operations."
Environmental Modeling & Software, 19(7-8), 629-644.

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.

OSHA. (1992) "Permissible exposure limits." Occupational Safety and Health


Association, Washington, District of Columbia.

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. (1980). "Emission from liquid pools." Arthur D Little, Boston,


Massachusetts.

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.

Sutton, O. G. (1953). Micrometeorology, McGraw-Hill, New York.

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

ESTIMATION OF CLOUD DENSITY

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.

A.1.1 Properties of Ammonia- Water Vapor Mixture in Air


The ammonia and water vapor evaporate from the pool with time and this will affect the
physical properties of the vapor cloud. The major assumption in this discussion is that air
above the pool is maintained at 1 atmosphere pressure and that the initial and final
clouds are at thermal equilibrium. The thermodynamics of the ammonia –water vapor
mixture in air needs to be studied in detail to estimate its density with time. In a practical
scenario, the air above the pool will be humid and this factor needs to be accounted
while performing these calculations. The physical state of the final cloud will be
systematically determined by solving the mass and energy balance equations of the
cloud mixture.

A.1.2 Stepwise Estimation of Final Cloud Conditions


The first step in these calculations will be the mass balance for the initial and final cloud.
84

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.

A.1.2.2 Energy Balance

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.

A.2 Calculation of Cloud Density


The two important factors that will decide the density of the mixture will be the
composition of aqueous ammonia and the cloud volume used in calculations. In this
problem the cloud formation step is used to define its initial shape. At release time, the
cloud is assumed to take the shape of a vertical gas cylinder characterised by its radius
"R" and its height "H". The main assumption is that its mass density ‘ρ’ is within the
whole cloud. The main input of this step is the cloud volume. It can be given either by
86

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

Table 11. Test conditions for density validation

Ground Ambient Spill Air Relative Density Duration


temperature temperature temperature velocity Humidity of moist of spill
(deg K) (deg K) (kg) (m/sec) (%) air (hours)
(g/cu.cm)
297 303 293 2 50 1.09 3
87

0.95

Ammonia weight:28%
0.9

0.85 Ammonia weight:20%


Density of vapor cloud(kg/cu.m)

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)

Fig. 33. Density of vapor cloud

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

Table 12. Cloud density at different time instants

Cloud Volume Factor Cloud density (kg/cu.m)


30 (mins) 60(mins) 90(mins)
5 0.89 0.80 0.77
10 0.80 0.76 0.74
20 0.75 0.74 0.74

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

CONVECTION DRIVEN EVAPORATION METHOD

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.

B.1 Natural Convection Model in Enclosures

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):

Ra= Gr.Pr= ρ 2 . g.(Tg − Tpool ). h 3 . c p . β / µ. k

where,
Gr= Grashof Number
Pr = Prandtl number
Tg =Temperature of the ground.

β = Thermal expansion coefficient


µ =Dynamic viscosity of the solution.

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.

B.1.2 Basic Equations

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

In the above equations,

u = Velocity of the fluid

Tpool= Temperature of the pool.

x, y, z = The respective rectangular coordinates.

ux,, uy,u z =velocity components in the x,y,z direction.

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

Fig. 334. Meshed pool surface

B.1.3 Boundary Condition and Test Conditions

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.

Fig. 35. Top surface contour for a 10 cm pool depth


94

Fig. 346. Ground surface contour for a 10 cm pool depth

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

Fig. 37. Top surface contour for a 5 cm pool depth


96

Fig. 38. Ground surface contour for a 5 cm pool depth

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

CASE STUDY FOR DILUTE PROCESS STREAMS

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

Fig. 39. Quantity of HCN released into air

Fig. 40. Quantity of acrolein released into air


99

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:

A-5, Lakshmi Enclave,


Apt # 16, First Street Abhirampuram,
Chennai, 600018

You might also like