Two-layer integral model for calculating the evaporation rate from a liquid surface

Journal of Hazardous Materials 59 Ž1998. 167–187

Two-layer integral model for calculating the evaporation rate from a liquid surface J.P. Kunsch Department of Mechanical Engineering, Swiss Federal Institute of Technology, CH-8092 Zurich, Switzerland ¨ Received 26 March 1997; received in revised form 7 October 1997; accepted 12 October 1997

Abstract An integral model is presented for predicting the evaporation rate from a single component liquid layer exposed to wind. This approach simplifies the solution of the convection–diffusion equation by averaging over the height of the concentration boundary-layer thickness. The model is valid for hydraulically smooth as well as rough surfaces and yields results which are in good agreement with those obtained by mathematically more rigorous methods. The new integral model is conceptually simpler and provides an explicit estimate of the concentration boundary-layer thickness. It allows a compact analytic solution for the evaporation rate with emphasis on the

Abbreviations: A, Area, e.g., area of a pool A p ; b, Width of the pool; c, Concentration; c f , Friction coefficient; c p , Specific heat; cs , Concentration at the saturation point; c l , Specific heat of the pool; D, ŽMolecular. diffusion coefficient; EŽ x ., Mass-flow density flux of the evaporating fluid; E p Ž x ., Mass-flow density flux of the evaporating fluid for high vapour pressure; EhŽ x ., Heat-flow density flux from the wind to X the evaporating fluid; E, Mass-flow rate, averaged over the pool area; E , Dimensionless mass-flow density X X flux of the evaporating fluid E s ErŽ cs u ) .; E , Dimensionless, area averaged mass-flow rate of the X evaporating fluid E s ErŽ cs u ) .; k h , Heat-transfer coefficient; k m , Mass-transfer coefficient; K, Eddy X X viscosity; K , Eddy diffusivity; L, Length, e.g., length of the pool in streamwise direction; L , Dimensionless X length of the pool L s Lr z o ; m l , Mass of the pool; p, Pressure; e.g., vapour pressure p v ; q, ˙ Heat-flow density flux; q˙s , Heat-flow density flux through radiation; re , Latent heat of vaporisation; Re, Reynolds number; Sc, Žmolecular. Schmidt number; St, Stanton number; St m , Stanton number for mass transfer; T, Temperature; Ts , Temperature at the saturation point; T` , Temperature of the ambient; u, Velocity of the windfield; u1 , Velocity at the reference height z1 ; u ) , Friction velocity; x, Coordinate in streamwise X X direction; x , Dimensionless coordinate in streamwise direction x s x r z o ; z, Coordinate normal to the liquid X X surface or pool; z , Dimensionless coordinate normal to the liquid surface or pool z s z r z1 ; z1 , Reference height; z o , Roughness length; z ov , Vapour roughness length; Greek letters: b , Empirical function in the concentration profile Eq. Ž7.; g , g s Ž s r k n b .Ž1q2 n.r n ; G , G s Ž dc r z1 .1q 2 n ; d , Thickness of the velocity boundary layer; dc , Thickness of the concentration boundary layer; k , von Karman ´ ´ constant; m, Dynamic viscosity; n , Kinematic viscosity; r , Density; s , Turbulent Schmidt number; t , Shear stress 0304-3894r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 0 4 - 3 8 9 4 Ž 9 7 . 0 0 1 4 5 - 3

168

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

diffusion process in the laminar sublayer, accounted for by means of empirical wall functions. q 1998 Elsevier Science B.V. Keywords: Evaporation rate; Hazardous liquid pools; Turbulent diffusion layer; Two-layer integral model

1. Introduction Evaporation processes are frequently encountered in our natural environment w1x but they are also of interest in many technical applications, e.g. in the process industry. Analyses of hazards related to loss prevention also require accurate prediction of evaporation rates. The main objective herein will be the derivation of compact analytical expressions for predicting the evaporation rates, in view of the needs of the loss prevention analyst to produce quick but reliable estimates. Past investigations of the evaporation phenomena in the free atmosphere include Webber and Brighton w2x and Boke ¨ and Hartwig w3x. In confined environments Že.g. tunnels., new problems occur, as discussed by Wren et al. w4x. The problem considered herein will be evaporation from wet surfaces, from pools or spills exposed to an airflow or the atmospheric wind. For volatile fluids with high vapour pressure, the evaporation rate will be high and the pool can cool down considerably unless it is deep. In this case, heat transfer from the ground can be important. For pools consisting of more than one component, concentration gradients can be important w5x. Here, only single-component pools will be considered. 1.1. General outline of the present contribution The present study will be restricted to the problem of pool evaporation in the wind field of a neutral atmosphere. The problem will be two-dimensional with an evaporating surface of finite length in the flow direction but of infinite extent transverse to the flow Žsee Fig. 1.. In a general overview ŽSection 2.1., the two representative models of Sutton w6x and Brighton w7x for solving the convection–diffusion equation will be presented. The underlying assumptions and simplifications will be discussed, in particular, the improvements of Brighton’s model with respect to that of Sutton.

Fig. 1. Concentration boundary layer above a liquid pool in the wind.

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

169

In subsequent sections, the integral methods, familiar to engineers interested in mass and heat transfer, following the lead of von Karman ´ ´ or Pohlhausen, are presented. Here, the equations are satisfied only by averaging over the thickness of the concentration boundary layer w8x. In order to demonstrate the main features of the integral method, the integral formulation and the corresponding solution of Sutton’s model will be compared, using the same physical assumptions as in the original model of Sutton. A new method, based on the same physical assumptions as used in the method of Brighton w7x, will be presented in Section 2.3. It is conceptually simpler and more intuitive than the method of Brighton and makes use of the concept of boundary-layer integral methods. As in similar applications in the past, it will be demonstrated that the method offers good accuracy with little computational effort. The simple analytic expressions for the mass-flow rate of the evaporating fluid will be useful for safety studies. The present method allows also an explicit but approximate solution for the thickness of the concentration boundary layer not derived in the mathematically more rigorous approaches of Sutton or Brighton. The power profiles, approximating the logarithmic velocity distribution, will be improved based on the analytic estimate of the concentration boundary-layer thickness.

2. Evaporation from liquid surface layers and pools Fluids, with boiling points well above the ambient temperature, are in a liquid state under standard conditions Že.g. water, diesel fuel, etc.. and their evaporation rates are small due to their small vapour pressures. The evaporation rate and how it changes with time is controlled by a mass-transfer process. 2.1. General oÕerÕiew The solution of Sutton w6x as well as that of Brighton w7x of the convection–diffusion problem will be presented in this section. The improvements of the approach of Brighton with respect to that of Sutton will be outlined. The physical basis of both models will be described as it is used also in the formulation of the integral model proposed herein. The concentration field above an evaporating pool can be described by means of a convection–diffusion equation. When simplified using the boundary-layer approximation, it takes the form w9x: Ec uŽ z .

E s

Ex

Ez

KX Ž z.

Ec Ez

.

Ž 1.

In what follows, a rectangular pool of unit width Ž b s 1. and of length L in the direction of the wind is considered. The gradients of velocity vanish in the flow direction because the wind field is fully developed: EurE x s 0. The solution given by Sutton w6x will be discussed first. The velocity field is approximated by means of a power law: u s u1Ž zrz1 . n, where u1 represents the velocity at a reference height z 1. Sutton presented his solution for a general exponent n, but

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

170

specified n s 1r7 for neutral atmospheric boundary layers. A justification for this value of the exponent is given, e.g. by Schlichting w8x. The exponent n s 1r7 is hereby derived from the resistance formula of Blasius for fully developed turbulent pipe flow. This is the same value used for flat plate flows w8x. Due to the similarity of the turbulent flow in a pipe and over a flat plate, the same power law is used for the velocity profiles in both cases w8x. The mass-flow rate of the evaporating liquid E s yK X EcrE z depends on the eddy diffusivity K X describing the concentration field. The eddy diffusivity K X is proportional to the eddy viscosity K, i.e. K X s Ž1rs . K, with the turbulent Schmidt number denoted s. Sutton used von Karman’s ´ ´ formulation of the eddy viscosity with mixing length l s z k E2EurE 2 . Applying the ‘statistical theory method’, outlined in w6x, Sutton obtained for urE z the eddy viscosity: r28 , K s 0.1386 n 1r4 u13r4 z 6r7 zy3 1

Ž 2.

with n s 1r7 in the power law for the velocity field. The shear stress t is accordingly:

t r

1r4

n

s 0.020

ž /

u12 ,

u1 z 1

when von Karman ´ ´ constant k s 0.40 and a turbulent Schmidt number s s 1 are used w6x. The classical value of the coefficient from PrandtlrBlasius is 0.0225 for hydraulically smooth surfaces, remembered here for comparison w6,8x. The solution of Sutton is based on the method of singularities, and the rectangular pool is modeled by superposition of line sources of infinite length transverse to the wind. For the area–averaged mass–flow rate of the evaporating substance, Sutton w6x obtained: r9 1r9 E s 0.0277cs Rey2 z 1 u1 Ly1r9 , 1

Ž 3.

where cs is the vapour concentration for saturated vapour and Re1 s u1 z 1rn . The frequent use of Sutton’s model is due mainly to the compact analytic expression of E. Some weak points of the model have been pointed out by Brighton w7x. One objection is that the turbulent mass transfer is not linked to the molecular diffusion process, which plays a major role in the laminar sublayer close to the evaporating surface. Moreover, E is not formulated in dimensionless form with respect to the relevant similarity parameters Re and Sc. Finally, parameters, such as the friction velocity u ) and the roughness length z o , that define the inner structure of the boundary layer, are not taken into account. As an improvement, Brighton makes use of the logarithmic velocity distribution, defined by u ) and z o , for the neutral boundary layer: uŽ z . s

u)

k

ln

z

ž / zo

.

Ž 4.

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

171

corresponding to a linear eddy-viscosity profile K s k u ) z Žsee Ref. w10x or w8x.. The logarithmic velocity distribution Eq. Ž4. is not well suited for an analytic solution of Eq. Ž1. and it is approximated by means of a power law: u s u1

n

z

ž / z1

.

Ž 5.

Formally, this is the same as that used by Sutton, but in order to obtain an optimal approximation in the concentration boundary layer, the reference height z 1 should be related now to the thickness of the concentration boundary layer. Therefore, equality of the velocity u and of the gradient d urd z is required at location z 1 for both profiles Eqs. Ž4. and Ž5.. Then u1 and n can be derived in a subsequent step. An estimate of z 1 is given by the theory of line sources w10x. Brighton used a linear eddy-diffusivity profile K X s Krs s k u ) zrs , yielding the logarithmic velocity distribution Eq. Ž4.. By substituting K X and the velocity distribution Eq. Ž5. into the convection–diffusion Eq. Ž1., he obtained the following linear differential equation: u1

z

ž / z1

n

Ec

k u ) z Ec

E s

Ex

s

Ez

Ez

.

Ž 6.

Brighton solved this equation by means of Laplace transforms. He matched the concentration profile, obtained by solving Eq. Ž6., to a wall profile that introduces the effects of molecular diffusivity in a small region close to the liquid surface. By using a linear eddy-diffusivity profile Žsee, e.g. Ref. w10x. E s yK X EcrEz s yŽ Krs .EcrE z s yŽ K m ) zrs .EcrEz can be integrated to obtain such a wall profile: c cs

s1y

E

s

cs u ) k

ln

z

ž / zo

qb .

Ž 7.

The term b can be considered as an integration constant and accounts for molecular diffusion processes in the laminar sublayer. b stems from empirical correlations and has different forms depending on the properties of the surface. This can be either hydrodynamically smooth or rough. An overview of the different forms of b , including a discussion of the physical background, is given by Brutsaert w1x. A ‘heuristic’ representation of the two-layer model for the velocity and the concentration field, with definitions used in the present work, is given in Appendix A. Brighton w7x takes for the hydrodynamically rough surfaces: u) zo b s 7.3Re1r4 Sc 1r2 y 5s where Re o s , o n and for the hydrodynamically smooth surfaces: s s 2 b s Ž 3.85Sc 1r3 y 1.3 . q ln Ž Sc . q ln Ž 0.135 . . k k The profile given by Eq. Ž7. is universal in the vicinity of the wall Žfor Schmidt numbers larger than 0.5.. It includes the molecular diffusion process in the laminar

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

172

sublayer and the turbulent diffusion described by K X . In the following, it will be denoted the ‘wall profile’ in accord with Brighton w7x. The effective roughness for hydrodynamically smooth surfaces is z o s 0.13 nru ) ; whereas only crude estimates are available for hydrodynamically rough surfaces Žcf. Brutsaert w1x.. The integrals occurring in the algorithm of Brighton can be solved by means of developments in series of the appropriate variables in the transformed space. The results of this mathematically rigorous integration of the convection–diffusion equation will be used for comparison with the results of the approximate integral approach, which is developed in the present contribution. 2.2. General features of the integral procedures This section deals with the integral approach of the convection–diffusion equation. Integral methods, in general, are approximate methods where the boundary-layer equations are satisfied only on an average extended over the thickness of the boundary layer. These methods, which can be attributed to von Karman ´ ´ and Pohlhausen, are well documented in Schlichting w8x. Techniques for approximating the profiles in a velocity or concentration boundary layer are addressed and their limitations are discussed in this reference. On integrating Eq. Ž1. over the height dc , one obtains: E

dc

H0

Ex

Ž uc . d z s K

X

Ec

dc

Ez

0

s EŽ x . .

Ž 8a .

Because the diffusive transport is much larger at the liquid surface than at the outer edge dc of the concentration boundary layer, whose thickness increases continuously, one obtains for the RHS of Eq. Ž8a.: E Ž x . s yK X

Ec

K Ec

.

sy Ez

s Ez

zs 0

Ž 8b .

zs0

The concentration vanishes at the outer edge dc of the concentration boundary layer, i.e., cŽ dc . s c` s 0 Žor in case of c` / 0: cX Ž dc . s 0 where cX s c y c` .. Such approximate formulations of the boundary conditions are normally used in the integral approach in boundary-layer theory w8x. With cŽ dc . s 0, the differential operator can be put in front of the integral to obtain: d dx

dc

H0

Ž uc . d z s E Ž x . .

Ž 9.

The integration of Eq. Ž9. over the length L of the pool gives a mass-flow rate per unit width, which is averaged over the length of the pool: Es

1

L

1

d EŽ x . d x s H H L 0 L 0

cŽ L .

Ž uc . d z.

Ž 10 .

Fig. 1 allows a straightforward interpretation of this formula: Eq. Ž10. represents the convective mass-flow rate through the vertical plane at location x s L.

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

173

2.2.1. Integral formulation and the corresponding solution of Sutton’s model The integration of the height–averaged differential equation for the convection-diffusion problem, Eq. Ž9., is the main objective of this section. Agreement between the integral approach and Sutton’s work will be demonstrated. Basic assumptions related to the physics of Sutton’s model, as discussed in Section 2.1, are therefore also used in the present integral formulation. The eddy viscosity Eq. Ž2., used by Sutton, is introduced in Eq. Ž9.. In the integral approach, the profile for the velocity u is approximated using the power law of Eq. Ž5. already used by Sutton. The one-seventh power law for the velocity profile suggests also the form of the concentration profile, which fulfils the following boundary conditions at the surface and at the outer edge, i.e., cŽ z s 0. s cs , resp. cŽ z s dc . s 0: c s cs 1 y

z

1r7

ž /

.

dc

Ž 11 .

The thickness of the concentration boundary layer, that can vary in the streamwise direction, i.e., dc s dc Ž x ., is still unknown. Here, similarity concepts in boundary-layer theory are used by choosing profiles where the ratio Ž zrdc . Žwith dc s dc Ž x .. can be considered as a similarity variable h. Eq. Ž11. is only an approximation of the real concentration profile in the concentration boundary layer up to the upper boundary z s dc . The good agreement of the area–averaged mass–flow rate with Sutton’s result demonstrates the validity of the profile approximation, in particular the choice of the one-seventh power in Eq. Ž11.. For a documentation of the integral approach, including a discussion of different assumptions for the boundary-layer profiles, it is referred to Ref. w8x. The mass-flow density flux EŽ x ., i.e. Eq. Ž8b., can be obtained with the help of the eddy viscosity, Eq. Ž2., and the concentration profile, Eq. Ž11.: r28 y1 r7 dc , E Ž x . s 0.020 sy1 cs n 1r4 u13r4 zy3 Ž 12 . 1 where EŽ x . does not depend on z. With the profile assumptions Eqs. Ž5. and Ž11., the integral on the LHS of Eq. Ž9. can be evaluated. Substituting EŽ x . from Eq. Ž12. in the RHS of Eq. Ž9., a differential equation for the unknown concentration boundary-layer thickness dc Ž x . is obtained. The solution is: 7

7

1

dc Ž x . s 0.318 sy 9 uy 1 36 z1 36 n

7

7

36

x 9.

Ž 13 .

With this result, the mass-flow density flux EŽ x ., Eq. Ž12., is also known. The integral in Eq. Ž10. can be evaluated to yield the area–averaged mass–flow rate: 1 L Es E Ž x . d x s 0.0262 sy8 r9 cs Rey2r9 u1 z 11r9 Ly1r9 . Ž 14 . 1 L 0 The coefficient 0.0262 is quite close to Sutton’s value 0.0277. The numerical value obtained depends on the von Karman ´ ´ constant, the value used here is k s 0.40. For the turbulent Schmidt number, Sutton chose the value s ( 1. In contrast to Sutton’s model, the solution obtained from the integral formulation yields also an estimate of the thickness of the concentration boundary layer dc Ž x . with elementary algebra.

H

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

174

2.3. Present integral model The approximate profiles for velocity and concentration, which are required for the integration of the height–averaged differential equation for the convection–diffusion problem, Eq. Ž9., are identical with those used by Brighton w7x. The profile for velocity is given by Eq. Ž4. and the profile for concentration by Eq. Ž7.. The term b is discussed in Section 2.1. Since the concentration profile includes the molecular diffusion processes in the laminar sublayer as well as the turbulent diffusion described by K X ŽSection 2.1., it represents a reasonable approximation for the whole height of the concentration boundary layer. The boundary conditions are accordingly c s cs at the liquid surface or inner edge and cŽ z s dc . s c` s 0 at the outer edge of the concentration boundary layer. If the evaporating component is present in the outer field, c s c` / 0. In this case, a new variable cX s c y c` is defined. In the ‘exact’ analytic solution of Brighton w7x, the logarithmic velocity profile was approximated by a power law given in Eq. Ž5.. The value of the reference height z 1 should be related to the thickness of the concentration boundary layer. It will be discussed later. The values n and u1 are determined by the requirement that the power law should be a suitable approximation to the logarithmic velocity profile at height z 1: 1 z1 s ln , Ž 15 . n zo

ž /

and: u1

ž /

1 s

kn

u)

.

Ž 16 .

Eq. Ž16. is substituted into Eq. Ž5. to obtain: u

ž / u)

1 s

kn

z

n

ž / z1

.

Ž 17 .

Since Eq. Ž7. is not well suited for an analytic integration over the relevant height of the convection–diffusion equation, an approximation is also required for the concentration field. By combining Eq. Ž4. with Eq. Ž7., one obtains a concentration profile which is consistent with the velocity profile: c

ž / cs

s1y

E cs u )

s

u

ž / u)

qb ,

Ž 18 .

and with the approximation Eq. Ž17. for u: c

ž / cs

X

s1yE s

1

kn

z

ž / z1

n

q b where EX s

E cs u )

.

Ž 19 .

In order to obtain a suitable approximation for the concentration boundary layer, z 1 is set equal to the thickness dc of the concentration boundary layer at the end of the pool. z 1 s dc Ž x s L . .

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

175

ŽThe mean value z 1 s 0.5dc Ž x s L. would be more accurate, but the effect on the final numerical result is small.. The value of n is unknown a priori, but the value n s 1r7, from the power law for the velocity, can be taken as a first approximation, as dc is unknown in the first iteration step. A better value of n can be obtained from Eq. Ž15. in a subsequent step. The boundary condition c s c` s 0 is valid at the outer edge z s dc of the concentration boundary layer. On using this condition in Eq. Ž19., which is solved for EX , one obtains:

s

X

Es

kn

y1

n

dc

ž /

qb

z1

.

Ž 20 .

The boundary condition cŽ z s 0. s cs at the liquid surface or inner edge of the boundary layer is not correctly accounted for. In effect, it is shifted from the liquid surface z s 0 to a finite but small z-value, producing a small error. This error in the concentration profile near the wall has little effect on the final result for two reasons: Ža. the approximation of the logarithmic velocity profile by means of a power law is not very accurate in the vicinity of the wall anyway Žlikewise true for the ‘exact’ solution of Brighton., and Žb. the convective contribution in the neighbourhood of the wall is negligible. Also, the mismatch of the concentration gradient EcrEz at the outer edge of the concentration boundary layer z s dc , which cannot be correctly accounted for, is a weak point here, as in other integral boundary-layer solutions using power-law profiles. The successful comparison of the results obtained, by the present model with those of Brighton, supports the validity of the assumptions made. The convection–diffusion equation in integral form, Eq. Ž9., can be written in the non-dimensional form: z1

ž / zo

d

H d rz dx 0 Ž

c

1.

X

u

c

u)

cs

ž /ž /

d z X s EX ,

Ž 21 .

when the dimensionless variables xX s xrz o and zX s zrz1 are introduced. The coordinate x is nondimensionalized by the roughness length z o , which is assumed constant for the problem considered, in accord with Brighton. The definition of zX is also in accord with Brighton; z 1 corresponds approximately to the thickness of the concentration boundary layer. The ratio Ž z 1rz o ., from Eq. Ž15., can be substituted into Eq. Ž21. to give: e 1r n

d

d rz H dx 0 Ž

c

1.

X

u

c

u)

cs

ž /ž /

d z X s EX .

Ž 22 .

The profiles Eqs. Ž17. and Ž19. are now substituted into the integrand of Eq. Ž22.. Upon evaluation, one obtains: 1

d

C1 d x X

E

X

dc

ž / z1

1q 2 n X

s E with C1 s

k 2 nŽ 1 q n. Ž 1 q 2 n. s e 1r n

.

Ž 23 .

The differential equation Eq. Ž23. and the algebraic equation Eq. Ž20. can be solved

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

176

for the unknown dimensionless mass-flow density flux EX and the boundary-layer thickness dc . The solution procedure, which is based purely on algebraic operations, can be found in Appendix B. The solution for the thickness of the concentration boundary layer is:

dc

ž / z1

1 X

s Ž C1 C 2 x .

1q2 n

with C2 s 1 q

1

n

2 1qn

.

Ž 24 .

When this result is substituted into Eq. Ž20., the mass-flow density flux EX can be obtained: E s bq

kn

y1

n

s

X

X

Ž C1 C 2 x .

1q 2 n

.

Ž 25 .

The concentration boundary-layer thickness, normalized by the known roughness length z o , is:

dc

z1

dc

zo

z1

1

ž / ž /ž / s

zo

se

1r n

X

Ž C1 C 2 x .

1q2 n

,

Ž 26 .

where Ž z 1rz o . is provided by Eq. Ž15. and Ž dcrz1 . by Eq. Ž24.. Now the estimate of the power n can be improved. The first approximation of n will be denoted by n a . The logarithmic profile Eq. Ž4. for the velocity field is approximated by the power law of Eq. Ž5.. In order to describe the convective transport of the evaporated fluid adequately, the power-law approximation should be valid over the height of the concentration boundary layer. This is achieved when z 1 corresponds approximately to the thickness of the concentration boundary layer at the downstream edge of the pool dc Ž x s L.. An improved value of z 1 corresponding to dc Ž x s L., denoted z 1new , can be obtained when xX is set equal to LX in Eq. Ž26., i.e., 1

z 1new

ž / ž s

dc Ž x s L .

zo

zo

/

s e 1r n a Ž C1C2 LX .

1q 2 n a

.

Ž 27 .

When Ž z 1new rz o . is substituted into Eq. Ž15., an improved value of n is obtained: 1 n

s ln

z 1new

ž /

1 s

zo

1 q

na

1 q 2 na

ln Ž C1C2 LX . with n a s 1r7.

Ž 28 .

The values of n at both extremes of the range considered here are nŽ LX s 10 3 . ( 1r4.5 and nŽ LX s 10 7 . ( 1r11.6. The calculations documented in Figs. 2 and 3 are based on the values of the power n given by Eq. Ž28. Žone iteration.. The dimensionless mass-flow rate, averaged over the pool area, is: EX s

1 LX

X

L

H0

EX d x X s

1 LX C1

EX

dc

ž / z1

1q 2 n

, X

X

x sL

Ž 29 .

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

177

Fig. 2. Area–averaged dimensionless mass–flow rate for the evaporation above a hydrodynamically smooth surface ŽParameter Sc s 4; 2; 1; 0.5; from the left to the right. ‘exact solution’ of Brighton w7x: - - -, numerical solution by means of the integral formulation: – P –, approximate analytic solution of the integral formulation: ŽEq. Ž30.., numerical values of Reijnhart et al. w12x: I I, ŽParameter Sc s1.7; 1.4; 1.0; from the left to the right., experimental values of Reijnhart and Rose w11x ŽSc s1.7.: v v.

where the integrand EX is given by the LHS of Eq. Ž23.. The initial condition that the concentration boundary layer starts at the upstream edge of the pool, i.e. dc Ž xX s 0. s 0, has been used here. The dimensionless area–averaged mass–flow rate Eq. Ž29. can be rewritten, when Eqs. Ž24. and Ž25. are substituted into Eq. Ž29.: C2

EX s

s bq

kn

n

Ž C1C2 LX .

.

Ž 30 .

1q 2 n

with the power n given by Eq. Ž28.. The mass-flow rate, averaged over the pool area is now E s EX u ) cs s k m cs where k m is the mass-transfer coefficient. The area–averaged dimensionless mass–flow rate, as obtained by means of the compact expression Eq. Ž30., can be compared to the area–averaged dimensionless

178

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

Fig. 3. Area–averaged dimensionless mass–flow rate for the evaporation above a hydrodynamically rough 2 . surface ŽParameter Re1r o Sc s100; 50; 20; 10; 5; 2; 1; from the left to the right ‘exact solution’ of Brighton w7x: - - -, numerical solution by means of the integral formulation: – P –, approximate analytic solution of the ŽEq. Ž30.., numerical results of Reijnhart et al. w12x: ^ integral formulation: ^, correlations of Brutsaert w1x Ževaporation pans.: ` `, correlations of Brutsaert w1x Žlakes.: I I.

mass–flow rate, determined by means of Brighton’s model, for the hydrodynamically smooth case Žcf. Fig. 2. as well as for the hydrodynamically rough case Žcf. Fig. 3.. The discrepancy is always less than 10% In case of a fixed power n s n a s 1r7, instead of the improved value given by Eq. Ž28., the discrepancy can exceed 20% Žfor Sc s 0.5 Žhydrodynamically smooth surface. and for Re1r2 Žhydrodynamically rough o Sc s 1 surface... A comparison with the most recent experimental and theoretical data for evaporation of single component liquids, as already presented by Brighton, is shown in order to give an idea of the accuracy of the present model. The main features of the different approaches will be resumed here. The experimental data of Reijnhart and Rose w11x ŽFig. 2: v v v, Sc s 1.7. for smooth surfaces were obtained in a wind tunnel. They provided numerical values for u ) and achieved a fully developed logarithmic boundary

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

179

layer, two prerequisites for a reasonable comparison to the present model or that of Brighton. Both requirements are generally not met in most experiments reported in the literature. The numerical solution of Reijnhart et al. w12x ŽFig. 2: I I, Sc s 1.0; 1.4 and 1.7. for solving the convection–diffusion equation, with allowance for molecular effects near the wall, is based on physical arguments similar to those of the present model. The numerical results from Reijnhart and Rose w11x, using the Reijnhart et al. w12x model for rough walls, with no allowance for molecular diffusion at the wall, correspond Ž to the following value of the parameter Re1r2 ^.. The o Sc s 0.34 Fig. 3: ^ correlations for large evaporating surfaces, like lakes, ŽFig. 3: I I. and for small evaporation pans ŽFig. 3: o o., given by Brutsaert w1x, are based on his recommended value z o s 2.28 = 10y4 m. The discrepancies observed have to be considered in the light of the scatter of the original data used by Brutsaert w1x. The existence of a roughness change from the shore to the evaporating surface has also been ignored in the present model. In order to render a comparison possible, Brighton assumed a typical wind speed of 5 mrs at height 10 m, yielding Re o s 2.8 and Re1r2 o Sc s 1.0 under these conditions. When the vapour pressure of the fluid is high, the normal velocity at the pool surface is considerable and the convective contribution to the vertical mass transfer cannot be neglected. With the boundary condition, stipulating that no air penetrates into the pool Žone-sided diffusion., the dimensionless mass-flow rate at high vapour pressure Žcf. Refs. w7,13x. can be obtained: EXp s yEX

p pv

ž

ln 1 y

pv p

/

,

where p is the total pressure and pv is the vapour pressure. The area–averaged dimensionless mass–flow rate can be written accordingly: EXp s yEX

p pv

ž

ln 1 y

pv p

/

.

Ž 31 .

where EX is given by Eq. Ž30.. The area–averaged mass–flow rate E is: Ep s cs u ) EXp s cs k m .

Ž 32 .

2.3.1. Dependence on pool length (comparison with Sutton) The dependence of the area–averaged mass–flow rate E on the pool length L will be analysed next. Eq. Ž30. yields: C2

E s c s u ) EX s c s u )

s bq

kn

n X

Ž C1 C 2 L .

1q 2 n

.

Ž 33 .

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

180

b in the denominator can be neglected if the length of the pool LX s zLo 4 C11C2 Ž bks n . It follows:

1q 2 n n

.

yn X

Ef Ž L .

1q2 n

.

Ž 34 .

With n s 1r7, Sutton’s well known dependence on pool length is recovered as a particular case Žsee Eq. Ž14... It follows that the validity of Eq. Ž34., adopted by Mackay and Matsugu w14x and many other authors as a basis for empirical or numerical models, is restricted. The neglect of b would be equivalent to the neglect of the laminar sublayer. 2.4. Heat transfer A convection–diffusion equation for the temperature field, analogous to Eq. Ž1., can be formulated in case of heat exchange with the ambient w9x. The solution Eq. Ž30. can be adopted accordingly for the area–averaged dimensionless heat-flux: C2

EXh s

bh q

sh kn

n X

Ž C1 C 2 L .

,

Ž 35 .

1q2 n

and the dimensional area–averaged heat-flux: C2

E h s r c p u ) Ž T` y Ts . EXh s

bh q

sh kn

n X

Ž C1 C 2 L .

r c p u ) Ž T` y Ts . ,

Ž 36 .

1q2 n

or, in compact form, in terms of the heat-transfer coefficient k h : E h s k h Ž T` y Ts . .

Ž 37 .

b h corresponds to the b of Eq. Ž7. where the Schmidt number Sc has to be replaced by the Prandtl number Pr. There is still some discussion about the best value of the turbulent Prandtl number s h , with proposed values ranging from 0.74 to 1.0 w1x. A value for the turbulent Prandtl number s h s 0.85, which is consistent with the turbulent Schmidt number used by Brighton, is considered to be reasonable. 2.5. Heat balance of the pool A general formulation of the heat balance of a pool is given by: ml cl

dTs dt

s yEp re q E h q qA ˙ p q q˙s A p ,

Ž 38 .

where m l is the mass of the pool, c l its specific heat, Ep re the loss of heat due to evaporation, Eh the heat addition from the wind or airflow, qA ˙ p the heat addition from the ground and q˙s A p the heat addition by radiation. The system attains an equilibrium

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

181

state when the heat needed for evaporation is balanced by the heat gained from the airflow. In this case, the heat balance of the pool Eq. Ž38. yields, in terms of the area-averaged quantities: 0 s yEp re q E h .

Ž 39 .

The expressions in Eq. Ž32. for Ep and Eq. Ž37. for E h can be substituted into Eq. Ž39. to obtain: 1 kh re k m

Ž T` y Ts . s cs .

Ž 40 .

It has been assumed here that no radiation is present and that the pool temperature is equal to the temperature of the ground. When the vapour pressure is low and when the molecular Schmidt and Prandtl numbers are equal Žincluding the case when the turbulent Schmidt and Prandtl numbers are assumed to be equal, i.e. s s s h ., a well-known formula, e.g., used in psychometry, is reproduced w15x: r cp Ž T` y Ts . s cs . Ž 41 . re When the airflow over the pool is purely laminar, one obtains for nonequal Schmidt and Prandtl numbers:

r c p Sc re

2r3

ž / Pr

Ž T` y Ts . s cs .

Ž 42 .

This expression, often generalized to the turbulent case, does not account for the laminar sublayer, in contrast to Eq. Ž40a.. Both the left and the right hand side of Eq. Ž40. can be represented graphically as functions of the temperature Ts . The temperature at the equilibrium point is found by the intersection of the Žnear linear. contribution of the left hand side with the concentration at the saturation point cs s csŽTs . of the right hand side. It is noted that cs is large for volatile fluids with high vapour pressure, so that the pool can cool down considerably. In this case, heat transfer from the ground becomes important. The temperature, the vapour pressure and the evaporation rate will depend strongly on time when the initial temperature of the pool is near ambient or higher ŽWren et al. w4x.. 2.6. Summary of the procedure When Eq. Ž30. is substituted into Eq. Ž31., which in turn is substituted into Eq. Ž32., the calculation of the area–averaged mass–flow rate can be reduced to the evaluation of the two expressions Eq. Ž43. and Eq. Ž28.: Ep s ycs u )

p pv

ž

ln 1 y

pv p

/

C2

s bq

kn

n

Ž C1C2 LX .

s k m cs

Ž LX s Lrz o . ,

1q2 n

Ž 43 .

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

182

with:

k 2 nŽ 1 q n. Ž 1 q 2 n.

C1 s

se

C2 s 1 q

1r n

1

n

2 1qn

.

The power n s n iq1 is given by Eq. Ž28a., which is written in a form suitable for an iterative procedure: 1

1 s

n iq1

1 q

ni

1 q 2 ni

ln Ž C1C2 LX . ,

Ž 28a.

with the starting value n i s n a s 1r7. Numerical experiments show that one iteration is sufficient. The total mass-flow rate is Ep s AEp with the surface of the pool A s Lb and its length L in the streamwise direction. The calculation of the heat-flux averaged over the surface area of the evaporating pool can be calculated in an analogous manner with Eq. Ž36.: E h s r c p u ) Ž T` y Ts . EXh C2

s

bh q

sh kn

n X

Ž C1 C 2 L .

r c p u ) Ž T` y Ts . s k h Ž T` y Ts . ,

Ž 36a.

1q2 n

with n s n i m given by Eq. Ž28a.. When the system has attained an equilibrium state, the temperature at the saturation point Ts is given by Eq. Ž40.: 1 kh re k m

Ž T` y Ts . s cs Ž Ts . .

Ž 40a.

3. Discussion The present work is focused on an integral model for solving the height-averaged convection–diffusion equation. One main objective is to develop a simple solution of the convection–diffusion problem and a second objective is to show the potential capabilities of the integral approach in general. In this approach, the equations are integrated over the thickness of the concentration boundary layer thereby reducing the number of independent variables by one. Results of the present model are compared to those obtained with Brighton’s model in Figs. 2 and 3. Brighton’s model is considered as a reference here since it has been validated by comparison with carefully compiled experimental data. Fig. 2 is related to the hydrodynamically smooth case with the Schmidt number as parameter and Fig. 3 to the hydrodynamically rough case with Ž Re1r2 o Sc as parameter. The good agreement between the numerical solution Appendix Ž . . A, Eq. B8 , obtained by means of the Runge–Kutta procedure and Brighton’s ‘exact’

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

183

solution confirms the approximations and simplifications made in the present integral model. A feature of the integral model, not encountered in Brighton’s model, is related to the approximation of the concentration profile. The chosen profile, required for the integration of the height-averaged convection–diffusion equation, does not account for all the boundary conditions correctly. Comparison with Brighton’s more rigorous mathematical approach shows, however, only a minor loss in accuracy of the present approximate approach. It can be concluded that the concentration profile, used in the present integral procedure and which accounts for the laminar sublayer by means of the parameter b , is approximated adequately. Brighton’s estimate of the boundary-layer thickness is based on Batchelor’s model with continuous line sources Žcf. Brighton w7x.. The present model, in contrast, allows the determination of the thickness of the concentration boundary layer. A new power n for the velocity profile can be derived, which is consistent with the model, yielding an improved approximation to the logarithmic profile over the thickness of the concentration boundary layer. The approximate analytic solution differs somewhat more from Brighton’s ‘exact’ solution. The accuracy could be improved by extending the analytic solution to higher order, but the additional effort at the expense of simplicity and physical insight is hardly justified in light of the uncertainties in the parameters defining the problem Žroughness, etc... In conclusion, the integral approach, which includes the possibility of estimating the thickness of the concentration boundary layer, provides a more intuitive, in some ways conceptually simpler, view of the problem than a mathematically more rigorous approach. The simplicity of the results obtained by means of the present approach comes at the expense of the numerical accuracy, but the compactness of the corresponding analytical expressions Žsee Section 2.6. meets the requirements of producing quick and reliable estimates for safety studies.

Appendix A A ‘heuristic’ representation of the two-layer model for the velocity and the concentration field, including a laminar sublayer and a turbulent logarithmic boundary layer is given in this appendix. The fully turbulent region, at a certain distance from the surface, can be described by means of the eddy viscosity K s k u ) z, which after integration of t s r KEurEz, yields the well-known logarithmic velocity profile w8x: uŽ z . s

u)

k

ln

z

ž / zo

,

Ž A1a.

where z o is an integration constant, called roughness length. In dimensionless form: uqs

u

ž / u)

1 s

k

ln

z

ž / zo

.

Ž A1b.

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

184

Fig. 4. Two-layer representation of the velocity field and of the concentration field.

When the mass-flow rate of the evaporating liquid is E s yK X EcrE z, with the eddy diffusivity K X describing the concentration field Ž K X s Krs , with the turbulent Schmidt number s ., one obtains in an analogous manner a logarithmic concentration profile Žsee Brutsaert w1x or Schlichting w8x.. Es z Es z zo cs y c s ln s ln q ln , Ž A2a. k u) z ov k u) zo z ov

ž /

ž / ž /

where the integration constant z ov is designated vapour roughness length. In dimensionless form: cs y c s z cqs s ln . Ž A2b. k z ov Ž Eru ) . When the Schmidt number is Žclose to. unity, the scaled velocity and concentration fields are similar. If the upper edge of the viscous sublayer is denoted by z s , one obtains: 1 zs uq ln , Ž A1c. s s k zo and: s zs cq ln , Ž A2c. s s k z ov

ž /

ž /

ž /

where the subscript s denotes the quantities at the upper edge of the viscous sublayer ŽFig. 4.. In this simplified representation of the velocity and concentration fields, the transition layer, where both turbulent transport and molecular diffusion are present, has been neglected Ždashed line in Fig. 4.. The viscous sublayer is ‘patched’ to the fully turbulent logarithmic layer at height z s z s . z s can be eliminated from Eq. ŽA1c. and Eq. ŽA2c. to yield: zo k k q ln s Ž cq b Ž Sc,Re o . . Ž A3. s y s us . s z ov s s

ž /

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

185

Eq. ŽA3. can be substituted into Eq. ŽA2a. with the result: c cs

s

E

s1y

ln

cs u ) k

z

ž / zo

qb .

Ž A4.

For Sc G 1, the molecular diffusion sublayer is contained inside the viscous sublayer of the velocity field. Empirical correlations, analogous to Eq. ŽA3., have been set up for Schmidt numbers deviating from unity. They are universal and depend on the Schmidt number only w1x. Note: Sometimes z is replaced by Ž z y d o ., which is equivalent to a shift in reference level. The d o is called a Žzero-plane. displacement height ŽBrutsaert w1x..

Appendix B The differential equation Eq. Ž23. and the algebraic equation Eq. Ž20. can be solved for the unknown dimensionless mass-flow density flux EX and the boundary-layer thickness dc .With the substitution: 1

dc

ž / z1

1q2 n

sG

,

Ž B1.

Eq. Ž23. can be written: 1

d

C1 d x

w EXG x s EX ,

X

Ž B2.

and Eq. Ž20.: Es

kn

y1

n

s

X

G

1q 2 n

qb

.

Ž B3.

EX of Eq. ŽB3. can be substituted into Eq. ŽB2. to obtain: 1

d

C1 d x X

~°G

¢

kn

G

1q 2 n

¶• ßs

y1

n

s

qb

kn

y1

n

s G

1q2 n

qb

.

Ž B4.

The differentiation on the LHS of Eq. ŽB4. is performed and the dependent variable G and the independent variable xX are separated to give: 1qn 1q2n

n

Ž gG .

1q 2 n

X

d Ž gG . s C1d Ž g x . with g s

n

Ž gG .

1q 2 n

q1

s

ž / k nb

1q2 n n

.

Ž B5.

q1

The equation can be integrated starting at xX s 0 Žthe beginning of the concentration boundary layer. with dc Ž x s 0. s 0, i.e., G Ž xX s 0. s 0, or Žg xX . s 0: ŽgG . s 0.

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

186

An approximate analytic solution can be obtained because the integrand in Eq. ŽB5. is close to unity in all cases. 1. In the limit G ™ 0, the integrand in Eq. ŽB5. is unity and the solution is G s C1 xX . 2. In the limit G ™ `, the integrand in Eq. ŽB5. tends towards Ž1 q n.rŽ1 q 2 n. and the solution is G s C1w1 q nrŽ1 q n.x xX . The solution of interest lies between these two extremes. The average of the two preceding extreme conditions can be considered as an approximation:

G s C1 1 q

1

n

xX ,

2 1qn

Ž B6.

valid for the whole range of G . This solution is confirmed by a perturbation analysis where ŽgG . is expanded with respect to the small parameter ´ s nrŽ1 q 2 n.. The difference between the resulting expression for G and Eq. ŽB6. is small of order n2 . It has been confirmed by numerical integration of Eq. ŽB5., for two representative values of n Ži.e., n s 1r7 and n s 1r3. that the relationship is almost linear for the range of xX of interest Ž0 F xX - Lrz o .. The thickness of the concentration boundary layer can be obtained when Eq. ŽB6. is substituted into Eq. ŽB1.: 1

dc

ž / z1

X

s Ž C1 C 2 x .

1q2 n

.

Ž B7.

Eq. ŽB5. can be written in a form suitable for a numerical integration ŽRunge–Kutta.: n

dG dx

X

s C1

Ž gG . 1qn 1q2n

1q2 n

q1

.

n

Ž gG .

1q 2 n

Ž B8.

q1

The integration is performed with respect to the independent variable xX because its range is well defined, i.e., xX s 0 % LX . By integration, one obtains G Ž xX . and in particular G Ž LX .. Eq. ŽB3. yields EX Ž LX ..

References w1x W. Brutsaert, Evaporation into the Atmosphere, Kluwer Academic Publishers, 1982. w2x D.M. Webber, P.W.M. Brighton, A mathematical model for a spreading vaporizing pool, in: S. Hartwig ŽEd.., Heavy Gas and Risk Assessment, Reidel, 1986. w3x J. Boke, S. Hartwig, Einsatz von Schaumen zur Konsequenzminderung bei Storfallen mit tiefkalten ¨ ¨ ¨ ¨ Flussigkeiten IV. Symposion, Schwere Gase und Sicherheitsanalyse, Bonn 26.r27.09.91, in: S. Hartwig ¨ ŽEd.., Battelle-Institut e.V., Frankfurt a.Main. w4x T. Wren, D. Martin, H. Wilkinson, Predicting evaporation rates of petrol pools in tunnels, 2nd International Conference, Safety in Road and Rail Tunnels, Granada, Spain, 3–6 April, 1995. w5x M.F. Fingas, A literature review of the physics and predictive modelling of oil spill evaporation, J. Hazard. Mater. 42 Ž1995. 157–175. w6x O.G. Sutton, Micrometeorology, A Study of Physical Processes in the Lowest Layers of the Earth’s Atmosphere, McGraw-Hill, 1953.

J.P. Kunschr Journal of Hazardous Materials 59 (1998) 167–187

187

w7x P.W.M. Brighton, Evaporation from a plane liquid surface into a turbulent boundary layer, J. Fluid Mech. 159 Ž1985. 323–345. w8x H. Schlichting, Boundary-Layer Theory, 7th edn., McGraw-Hill, 1987. w9x R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, 3rd edn., Wiley, 1963. w10x J.C.R. Hunt, A.H. Weber, A Lagrangian statistical analysis of diffusion from a ground level source in a turbulent boundary layer, Q.J.R. Met. Soc. 105 Ž1979. 423–443. w11x R. Reijnhart, R. Rose, Vapour cloud dispersion and the evaporation of volatile liquids in atmospheric wind fields: II. Wind tunnel experiments, Atmos. Environ. 14 Ž1980. 759–762. w12x R. Reijnhart, J. Piepers, L.H. Toneman, Vapour cloud dispersion and the evaporation of volatile liquids in atmospheric wind fields: I. Theoretical model, Atmos. Environ. 14 Ž1980. 751–758. w13x H. Grober, S. Erk, U. Grigull, Die Grundgesetze der Warmeubertragung, 3rd edn., Springer-Verlag, ¨ ¨ ¨ Berlin, 1981. w14x D. Mackay, R.S. Matsugu, Evaporation rates of liquid hydrocarbon spills on land and water, Can. J. Chem. Eng., Vol. 51, August, 1973. w15x J.P. Holman, Heat Transfer, 4th edn., McGraw-Hill, 1976.