Estimating Rates of Spreading and Evaporation of Volatile Liquids

Follow these guidelines to calculate how fast a pool of spilled liquid will spread across a surface, evaporate and potentially form a flammable mixture with the air.

Despite operators’ meticulous efforts to avoid spills during the handling of volatile liquids, accidents can, and do happen. In such cases, the ability to predict the rate at which the liquid spreads and how fast it evaporates prove valuable. The former would be instrumental in planning and designing containment. The latter would be useful in finding the vapor concentration of the substance in ambient air, which would help one to determine the electric al area classification, fraction of lower explosive limit achieved, and address other similar safety-related issues. In this article, methodologies will be presented for calculating spill spreading and evaporation rates. Examples featuring these methods are used to find the liquid mass remaining at any given moment and the time required to evaporate the entire spill.

First, one must assume that: the volume of spilled liquid is known (e.g., de rived from batch darn as the larges t vol u me used in a process) or can be derived (e.g., using the flowrate from a leak point and approximate du ration of the leak); the liquid is well characterized i n terms of density (p) , surface tension (cr), viscosity (µ) and vapor pressure (P ,.,) ; the li quid is at ambient temperature and barometric pressure , which are known; and the spill progresses as a liquid spreading  across a smooth, level sur


How far will a spill spread?)

A liquid that is spilled on a flat surface will progress in three regimes ( I ):

  • the gravity-inertia regime, in which gravity tends to spread the fluid and is opposed by the inertia of the fluid
  • the gravity-viscous regime, in which gravity tends to spread the fluid and is opposed by the viscosity of the fluid
  • the viscous-surface tension regime, in which the liquid viscosity is opposed by the surface tension of the fluid

The most volatile liquids have viscosities less than or close to that of water (~1 cP), and will enter the viscous-surface tension regime in a few seconds. Eq. 1 yields the time required to enter the viscous-surface tension regime: 

tvs = 0.023462(g   VpJ.J/ O)        ( I )

where g is the gravitational constant (ft/s2), V is the spill volume (ft3), and p (lb/ft3), u (cP) and O (dyne/cm) are evaluated at the ambient temperature (TA) or the temperate of the air above the liquid pool. With tvs one can calculate the radius of the spill at zero time (aO) (i.e., at the onset of the evaporation): 

o0 = I .4 1 3 1 42(crV1,./ µ) t1•

The shape of the volume of spilled liquid should be modeled in such a way as to enable the calculation of the area exposed to the atmosphere. In the real world, the spill assumes the shape of a spherical cap. If one determines the proportions of the spherical cap {and those of the corresponding hypothetical sphere), one can find the exposed surface area of the spill. The volume of a spherical cap is calculated as (2 ):

v,,.,,= (1th2t3)(3 r.,. ,, -  h)

where, h is the depth of liquid at the center of the spill, a is the radius of the spill, and rsph is the radius of the hypothetical sphere of 1ich the cap is part. The initial radius of the spill (i.e., the radius measured immediately following the brief interval tvs after the liquid is first spilled) is noted as "ao· It is cal­culated with Eq. 2. where the volume of the spherical cap is V0, and the time with respect to the evaporation process is zero. Note that during the spreading phase, no evaporation takes place.

Correspondingly, the height of the spherical cap at the center, h0, the maximum depth of the spill, is calculated as:

This cubic equation can be solved analytically, or, more conveniently, by the use of a spreadsheet solver function. If one assumes that the central half angle (13) of the cap (Figure I ) remain s constant during the course of the spill, the n:

where r0 is the initial radius of the spherical cap, and:

Then, substituting Eq. 7 into q. 5:

A collection of terms leads to the expression:

Eq. 9 is a cubic equation and may be solved for tanB using a spreadsheet. Reader seeking a rigorous solution should consult Ref . 3 and use the key words, "cubic equation."

With tanB (and therefore 13), /, and" in hand, one can calculate the surface area of the cap using (3): 

Weisstein defines rsph as (2):

Where the a is the compliment of B:

Other important relationships include: 

These equations will come in handy when seeking a solving equation for h. It is now possible to write the unsteady state mass balance on the spill, assuming that the evaporative mass flux (E) - or the evaporative rate, normalized over the area exposed to the environment - remains constant:

Here, Acap  and Vcap  are the surface area (ft2) and volume (ft3) of the spherical cap at any time t (mi n) ; W is liquid mass (lb); p is liquid density (lb/ft3); £ is mass flux (lb/ft2-min) ; and rph and h arc measured in ft. It i assumed that B is constant throughout the course of the spill. Eq. 4  may be rearranged to solve for rsph:
 

Taking the cubic root of both sides yields the following:

Combining Eqs. 22 and 27, and performing extensive rearrangement and substitution leads to:

The constants b and c are expressed in units of ft/lb1/3 when English units are used. Subsequently, an expression for Acap in terms of mass of the liquid may be derived: 

If one assumes that B and E are the constant, the cubic root of the mass in the spherical cap decreases linearly over time.


Determining the Evaporate Flux

There are three methods for estimating £ . Two apply to a spill exposed to a moving air stream. The third method, Stiver-MacKay, can be extended to cover the case of a spill exposed to still air.

  1. U.S. Air Force method. This empirical method is based on the evaporation of hydrazine at ambient temperatures (4). The evaporative flux for other liquids is estimated using the following equation, which is normalized for the effects of the temperate and the properties of liquid other than hydrazine:

In the equation above, M is molecular weight (lb/lbmol), PvpS and Pvp,H are vapor pressures of the spilled substance and hydrazine respectively (torr) and TF is a temperature correction factor defined conditionally as follows:

In the original work, Pvp,S, and Pvp,H are expressed in torr, but when using Eqs. 41 and 42, any consistent set of units is applicable, since the vapor pressure contribution is dimensionless. The original work also evaluates the vapor pressures at TA . Although TA  is not equal to Tp,, it is reasonable to assume they are equivalent, barring special situations (e.g., a cold liquid spilled in a warm environment).

  1. U.S. EPA method. Below is a slightly modified form of the empirical equation developed by the U.S. Environmental Protection Agency (EPA) to define a evaporative flux (5, 6):

Where u is the air velocity (ft/min) and Pvp,S is expressed in the units of torr, since the vapor-pressure contribution term is not dimensionless.

  1. Method of Stiver-MacKay. This method employs a mass transfer coefficient explicitly. As such, it lends itself to situations other than that of a liquid pool exposed to a moving air stream (7, 8, 9, 10):

In this case, k is the mass transfer coefficient measures in ft/min or ft/s, and R’ is the ideal gas constant measures in ft3torr/lbmooR. One can define k using the following empirical relationships. 


Beyond Constant Flux – Forced Convection

The derivations presented thus far are predicated on the assumption that the evaporative flux is independent of the geometry (and thus the characteristic dimension) of the spill. However, the mass transfer coefficient - and therefore flux – is usually a function of some characteristic length of the geometry in question. 

E will vary with the changing geometry of the spill because, in the real world, there is usually movement of air above the spilled liquid. This creates a pressure differential, causing evaporative mass transfer to occur by force d convection. To account for the effects of forced convection, a mass transfer coefficient that depends upon a characteristic dimension of the spill is introduced into the evaporative flux equation.

It is assumed here that the term "flow," except for the transient case of the spreading of the spilled liquid. refers to the air above the spill. Typically, the radius of curvature of the spill is sufficiently large such that the flow of air behaves like air flowing past a flat plate. This flow can be turbulent or laminar

Bennett and Myers state that for flow past a flat plate, the laminar-to-turbulent transition occurs at about Re 3 x 105 (11), where Re is the Reynolds number calculated for a plate of length L as:

L is a characteristic length of the geometry in question, and u, p and µ are the velocity, density and viscosity of the moving fluid, respectively. L takes the form of 2a, where a is the radius of the spherical cap. The velocity is assumed to have been measured (e.g., by a local anemometer) or determined otherwise (e.g., by a local weather report). The p and µ of air are obtained from table s of physical properties in a standard reference (e.g. , Perry' s Chemi cal Engineers' Handbook). 

Bennet and Myers show a dimensionless expression for the mass transfer coefficient for laminar now using the Sherwood (Sh), Schmidt (Sc) numbers:

where µ is liquid viscosity (lb/ft-min) and DAB is the diffusivity of substance A diffusing through substance B (ft2/min) and may be found by the methods described by Reid, et al.(12). Equating Eqs. 47 and 50, and solving for k.

The characteristic dimension a can be expressed in terms of the mass of the spherical cap using a combination of Eqs. 16 and 28: 

The result of substituting Eq. 57 in the Stiver-MacKay relationship for flux yields:

Use of Eq. 59 in the unsteady-state mass balance, followed by integration, results in this empirical equation for forced convection evaporative mass flux in the laminar flow regime: 

Thus, for laminar flow, when accounting for a change in the mas s transfer coefficient (and therefore £) as a function of the changing dimensions of the liquid pool, the square root of W decreases linearly with time.

ln many standard texts (II), the analogy between heat and mass transfer developed origin ally by Chilton and Colburn I used to derive mass transfer relationships for known geometries and flow conditions based on dimensionless numbers for heat transfer.

As a corollary, if one has a relation for heat transfer for a given geometry (e.g., flow past a flat plate), then by analogy, one has a relation for mass transfer for that same geometry. For turbulent now, this analogy between heat and mass transfer is used Lo find £. For heat transfer involving turbulent flow past a flat plate, Bennett and Myers give this correlation for the Nusselt number(13):

Nu = 0.0365Re4/5Prn1/2                                                               (62)

The exponent on the Prandtl number (Pr) is open to some debate. Based on other work (14), an exponent of 1/2 is used here. Thus, the mass transfer analog is assumed to be:

Sh = 0.0365Re4/5Sc1/2                                                                  ( 63)

Eq. 63 is equated with Eq. 50 to solve for k. As for laminar flow, subsequent expressions fork are derived and substituted into the solving equation for E in the Stiver-MacKay relationship. The integration of the unsteady-state mass balance yields the following equation for turbulent convection:

In the equations defining q and x, constants b and c are calculated as per Eqs. 29 and 32, respectively.

Extending the Method to Free convection

Next, the evaporative flux is examined under the conditions of free convection. The air above the pool i assumed to be completely still, and the driving force for mass transfer is the difference in concentration of the volatile compound between the liquid pool and the air above the pool.

To adapt the method explained for forced convection to the case of free convection, one needs a free-convection mass-transfer coefficient, which is again derived using the Chilton and Colburn analogy between heat and mass-transfer, as applied to free convection past a flat plate of liquid.

For free convection in the laminar flow regime (i.e., when 105 <; GrL Pr < 2 x 107, where Gr is the dimensionless Grashof number for laminar flow), the heal transfer coefficient may be expressed as:

and, for free convection in the turbulent flow regime : (i.e., when 2 x 107< GrLPr < 3 x 1010)                                       

Using the Chilton and Colburn analogy, the mass transfer coefficient for free convection in the laminar flow regime: (i.e., when 105< GRABSC<107)

Sh = 0.54(GrABwSc)1/4                                                                    (69)

and, for the free convection in the turbulent flow regime, (i.e., when 2 x 107 < GrABSc)1/3                                                                                                                  (70)

In the equations above, Gr and Sc are defined as (15):

For the purposes of this article, L = 2a. In addition, l\S refers to the difference in concentration of the evaporating species between the boundary layer of liquid and the bulk fluid above it. Usually, the concentration of the evaporating species in the bulk fluid is zero or effectively zero.

Following a procedure similar to that used previously, one obtains for laminar flow:

Follow the procedure used for laminar flow to assess turbulent free convection:

Where:

Rarely is the outdoor atmosphere completely still for any appreciable period of time. Therefore, for spills that occur under the condition of light winds to calm air, it is suggested that the estimated time for evaporation is calculated based on the average of the forced convection and free convection cases, since the actual situation lies somewhere between these two extremes. Furthermore, the upper limit on the product of the Gr and Sc numbers may limit the applicability of this analysis to small spills.

Example Problems

Physical and transport properties, where required, are cal­culated from empirical correlations given by Yaws (16).

Example I. Assume 50 gal of methanol spills onto a level surface outdoors. A local thermometer read TA= 59 ° F, and a local anemometer gives an average wind speed of u = 5 mi/h. Estimate the greatest depth of the spill (h) and the time it will take the spill to evaporate (t2-t1).

Summarize the known conditions and the physical properties of methanol: Pvp = 69 .0 58 mmHg, M = 32.044 lb/lbmol, W = 332.24 lb, R = 555 mmHg-ft3/lbmoloR, µ= 0.61 9 cP, cr = 24.869 dyne/cm, p = 49.707 lb/ft3, and DAB = lb/ft.

As a first pass, assume that the evaporative flux is independent of the dimensions of the spill ( i.e., E remains constant during the evaporation process). A preview of the cal­culation s reveals that the EPA method yields the shortest evaporation time, while the Stiver-MacKay method 

yields the longest evaporation time. Therefore, for a conservative estimate, the Stiver-MacKay method will be used.

Calculate the initial spreading time, tvs using Eq. I:

tvs = 0.023462[(32.I 74 ft/s2)(50   gal x 7.48 gal/fl3) (49.707 l b/ft3)(0.6l9  lb/ft-s)/24.869 dyne/ cm] = 6.24 s.

Calculate the pool radius at tvs using Eq. 23:

a0 = 1 .413 1 42 [ (24.869 dyne/cm)(6.684 ft3)(6.243 s)/

(0.619cP)]1/4 = 9.04 ft.

In this calculation, the unit conversion factors for µ and cr have been worked into the coefficient. The liquid pool is assumed to take the form of a spherical cap, due to the effects of surface tension. Given the volume and the radius at time zero, solve Eq. 5 for the maximum depth of the pool at its center:

1 2.77 ft3; therefore, h = 0.0S2 ft

B is found by rearranging g Eq. 9 and using a spreadsheet solver function:

Tan3 P + 3tanB = (1/ao3)(6Vo). Thus , tanB = 0.005755 rad, and B = 0.00576.

Per Eqs. 43 and 45:

E = 0.1 758(5 mi/h)(69.058 mm Hg)(32.044 lb/lbmol)/(555 m mHg-ft3 lb/mo1°R)(59+ 453.49°R)) = 6.76 x 10-3l lb/ft2min.

Use this result in Eq. 37 to find the evaporation time, t2 - t1. Solve for t2 with t1 = 0 and W2 = 0. This leads to:

t2 = 3W1 1/3 /2bc£ where:

b = (6W/((1 + 3cot2B))) 1/3 = 7.51 x 101/3 ft/lb1/3

and

c = 3(1 + 3cot2B)2/4p)1/3 + (6/(p(1 + 3cot2P)))1/3

= 113.42 ft/lb1/3.

Thus:

t2 = 3(332.24lb)1/3/(2 x 3.14 x (7.51 x 10-3 ft/lb1/3) x

(113.42 ft/lb1/3 )(6.76 x 10-3 lb/ft2min))= 574.1 7 min.

Example 2: Repeat Example 1, but this time, assume that the evaporative flux is a function of the pool radius (a0) under conditions of forced convection. Since the flux varies throughout the evaporation process, one needs an integrated mass balance that accounts for the effect of the pool's shrinkage on the flux. The Stiver-MacKay method is the only one that includes an explicit term for k, and will be used to perform the calculations. All of the physical properties and constants (e.g., b and c) are consistent with those cited in Example 1.

First, determine whether convection is turbulent or laminar using Eq. 46:

Re = (5 mi/h)(5,280 ft/mi)(0.076 lb/ft3)(2X 9.04 ft)/ ((0.018 cP)(2.419 lb/ft-h)/cP)= 8.33 x  105.

Since Re is greater than 3 x 105, now is turbulent and Eq. 64 should be used. This equation requires the determination of several constants.

DAB is determined using physical properly estimation methods described in Ref. 1 2, 11-4.4

and Table 11 - 1 10 be 0.160 cm2/s = 0.010 ft2/min . The central half- angle is calculated as P = 0 .0 0576 rad. Also

per Eq. 66:

x = 0.0365((0.010 ft2/min ) ln )((5 mi/h x 88

fi/mi n/(mi/h))413(0.619 cP(0 .04032 lb/fl

mi n/cP))"10(2(7.5 I x JQ- l ft/lb1il )co 1(0.00576 rad))-113

11.62 fl-lb111S/ min.

And, per Eq. 65:

q = 2rcbcxP"l' _sMI R'TA = 2rc(7.5 l x lQ- l ftllb til) x (113 .4 2 ft/lb1il)(I 1.62 ft-lb1l15/min )(69.058 mmHg) x (32.0422 lb/lbmol)/((555 mmHg-ft3/lbmol0   R)(59 +

453.49°R)) = 0.484 lb215/min.

Assume W2 = 0 and 11 = 0 , a od solve for r2

using Eq. 64, which is rearranged as:

12 = 5W .2'512q = 5 x (332 .24 l b)215/(2 x 0.484) =

52.68 min . 

As may be expected, the predicted time required to evaporate the entire spill decreases significantly when one ac­ count s for a change in the evaporative flux with the decreasing size of the pool.

Example 3. Consider a smaller spill ( V0 = 5 gal) of methanol. Once again, assume that the evaporative flux varies during the evaporation process. Assume that 11 = 0 ft/s and thus, only free convection takes place. Also, as­

sum e that the air above the spill contains a negligible concentration of vapor. Calculate the amount of time ii will take to evaporate the entire pill.

The Stiver-MacKay method will be use d because it includes an explicit term for k. This case exhibits turbulent flow free convec1ion, since ScGr = - 3.9 x 109. Thus, Eq. 76 is used with W2 = 0  and t1  = 0  to calculate /2:

where:

11  = 2rcbcmPvpSMIR'T ,. = 0.0887 (lb/min )11l

m = 0. 14 ( D,.8)2f.l(g 1'16/ µ)1il = 7.93 ft/mi n. Thus, 12 = 63 .6 min .

Example 4. Building o n Example 3, in which the evaporative flux varies, calculate lhc mass of liquid remaining, along wi1h the volume and radius of the spill, a evaporation progresses, until all of the liquid is evaporated.

Us e Eq. 76 10 solve for W2 with 11 = 0 and 12 varying from I min to 60 min. T o solve for a, find V us in g W/ p.

T he n, us i ng Eq. 2, solve for a. The results are shown in

t he Table and Figu re 2.

Literature Cited

  1. Pulorli, A. D., •I al., '"Flamma ble and Combustible Liquid Spill/Bum Panems:• U.S. Dept. of Justice; Office of Justice Programs. National lnsiitutc of J ustkc Report 604-00, Washington, DC (2001).

2. Weisslein, E., "Spherical Cap," hnp://math world.wolf111J11. com/ Spheri­ calCap,h lml.

J. Jeffrey, A., lf rmdbook of Matlrematicaf Fon11ulas   tmd hutgm ls, Aca•

demic Press, §0,8, San Diego, CA (1995).

4. U.S. Envlrorunentul P rotection Agency (EPA), Federal Emergency Management Agency (FEMA) and U.S. Dept. of Tnmsporlallon, Handbook of Ch, mica/ Hau1'-d Analysis Procttl11res, Washingt o n, DC (1989 ).

S. U . EPA and U.S. F EMA, '"l'cchnicol Guidance For Hazards Analysis," Equation 7, Section G-2, Appendix G, Washington, DC (Dec. 198 7); www.epa.gov/swercc pp/pubs/1cch.pdf.

  1. U.• EPA, ''Risk Management Progmm Guidance For Offsite Con e­ qucnce Annly is;• publication EPA-550-B-99-00'.I, Section D.2.3, Ap• pendix D, Equation 0-1 , Wn.,hington, DC (Apr. 1999); www.cpa,govlceppo/pubs/oca/oca-all.pdf.
  2. liver, W., and D. MacKay, "A Spill Hozard Ranking System For Chemicals,"   Environ1n e.n1  Canodo   F'i    t  Technical   Spil ls   Seminar, Toronto, Canada (1993).
  3. Clewell, H. J., "A Simple Method For Estim ating the Source Strength of Spill of Toxic Liquids," Energy Systems Laboratory, ESL-TR-83- 03, ovoiloblc at Air Force Weather Technical Library; Ash eville, NC (1983).
  4. Ille, G., and C. Springer, "The Evaporation and Dis pen;ion of Hy• dro1Jne Propellants from Ground Spills," Civil and Environmental En­ gineering Development Office, CEEOO 712-7 -30; available al Air Force WemherTechni cal Library, Asheville, NC (1978).
  5. Kohler, P.1 , al., "Calc ulating Toxic Co rridms", Air Force Weather Ser­ vice. AWS TR-80/003, 1980; nvailable at Air Poree Weather Technical Library, Asheville, NC (1980).
  6. Bennelt, . O., and J. E. Myers, "Momentum. Heut. and MassTrans­

fer," 2nd Ed.. McGmw-Hill. New York (1974).

  1. Reid, Robert C.,et al., 77,o Propertie, of Gases and liq11ids, 4th Ed.. McGn11v-Hill, New York (1987).
  2. Benn ell, 0. , And J. E. Myers, Equation. 23-51, p. 367.
  3. Hanna , O . T., and ,I. E. Myers, Engineering 6Jrpcrimem Stution. Bul­ letin No. 148, Purdue Univ. (1962).
  4.             wm v.svce.oc.in/- msubbu/Lec1u reN01es/M n.ss'lr ansfer/MT-2.htm , Eq.4.30.

16, Yaws, Carl L., Chemical Properti,s Handbook, McGmw- liill, New York (1999).

 

JOHN BARRY is a lead process engineer, currently working under contract for DuPont Corp. at the Chambers Works site in Deepwater, NJ (Phone: (856) 540-2364; Fax: (856) 540-3080; E-mail: john.barry.engr@comcast.net). Prior to consulting with DuPont, Barrv was a contract process engineer at Jacobs Engineering's Premcor refinery (Delaware City, DE) and a lead process engineer for CDI Engineering (Philadelphia, PA). Barry received a BS In chemistry from the Univ. of Delaware and a MS In chemical engineering from the Univ. of Maryland. He Is a licensed professional engineer In the slates of DE and NI, and Is a member of AtChE.