Numerical simulation of evaporation of volatile liquids

Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

Contents lists available at ScienceDirect

Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp

Numerical simulation of evaporation of volatile liquids A.D. Galeev*, A.A. Salin, S.I. Ponikarov Kazan National Research Technological University, Department of Mechanical Engineering, 68 Karl Marx Str., 420015 Kazan, Russian Federation

a r t i c l e i n f o

a b s t r a c t

Article history: Received 21 November 2014 Received in revised form 27 June 2015 Accepted 28 August 2015 Available online 1 September 2015

The paper presents the results of the validation of the developed pool evaporation model using literature and our own experimental data. The proposed model was used to examine the effect of wind velocity and pool sizes on the evaporation rate of volatile liquid (hexane). Contrary to the semi-empirical evaporation model widely used in hazard assessment, stronger dependence of evaporation rate on pool size at low wind speeds is obtained. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Liquid spill Pool evaporation Model validation Numerical simulation

1. Introduction The correct prediction of the consequences of hypothetical accidental scenarios is an essential component for developing and implementing appropriate protective capabilities and plan mitigation measures. This requires the use of the reliable mathematical models describing how toxic chemicals or flammable substances are released and dispersed in the air. A complete model for an event can roughly be divided into three parts: source modeling, dispersion modeling and effect modeling (Vik and Pettersson Reif, 2011). There is a lot of work concerning the modeling of hazardous substances dispersion in the environment, whereas little attention has been given to source term modeling. In emergency spills of stable liquids, the input of hazardous substances into the environment is the result of evaporation from the pool surface. Generally, there are two main types of mathematical models to evaluate the evaporation rate from a pool: semi-empirical and computational fluid dynamics (CFD) models. In the area of hazard analysis of accidental chemical spills, many investigators and hazard reference books have applied Mackay and Matsugu mass transfer equation (Kawamura and Mackay, 1987; Mackay and Matsugu, 1973). In the model (Kapias and Griffiths, 1999), the mass transfer coefficient is calculated using Brighton's theory (Brighton, 1987). This analytical model takes into account the effects of surface roughness, friction velocity of the airflow, and the

* Corresponding author. E-mail address: [email protected] (A.D. Galeev). http://dx.doi.org/10.1016/j.jlp.2015.08.011 0950-4230/© 2015 Elsevier Ltd. All rights reserved.

effect of high vapor pressure on the mass transfer process. Khajehnajafi and Pourdarvish (2011) used Churchill's equation (Churchill, 1976) where the Nusselt number is determined for the complete range of Re and Pr covering laminar, transition and turbulent regimes. This equation is based on experimental data for forced convection from isothermal flat plates. The Churchill's equation is applied to mass transfer calculation substituting the Sh and Sc for the Nu and Pr. The mathematical model, represented by Habib et al. (2009), based on numerical calculation of differential transient equations for boundary layer together with algebraic turbulence model of Cebeci-Chang (Cebeci and Chang, 1978), including correction for surface roughness. The above models have been developed on the assumption that the vapor behaves as a passive contaminant, i.e. it does not affect the flow field. This assumption may be invalid when the molecular weight of evaporating component differs from the molecular weight of environment air (Desoutter et al., 2009). By using the CFD, it is possible to overcome the limitations of existing semi-empirical models. This approach by solving threedimensional conservation equations for mass, momentum and energy makes it possible to take into account the complex interaction between the evaporation process and the vapor dispersion. Many works have used CFD to examine the evaporation rate of water (Raimundo et al., 2014) or dilute solutions of ammonia (Rong et al., 2010; Rong et al., 2011; Saha et al., 2011). In the work (Rong et al., 2010) the effects of airflow and aqueous ammonia solution temperature on ammonia transfer are investigated and the numerical sensitivity studies are performed by Rong et al. (2011) to

40

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

examine the effects of computational geometry and inlet turbulence intensity on ammonia emissions. In these studies, the ammonia was a passive contaminant and has no impacts on airflow. Raimundo et al. (2014) used CFD model to investigate the mass transfer at the free surface of the water tank. In the abovementioned works, the liquideair interface is assumed to be isothermal due to the small liquid volatility. When volatile liquids, such as hexane and acetone, evaporate, the liquid temperature may significantly change due to heat loss. Furthermore, the evaporation rate of volatile liquids may be significantly influenced by Stefan flux, which is not taken into account in above CFD models. Liu et al. (2011) developed a methodology using CFD to simulate the boiling process of liquid nitrogen. The vaporization rate of a boiling liquid is governed by the heat transfer phenomena including conduction, convection and thermal radiation mechachot et al., 2013), while the evaporation rate of nonnisms (Ve boiling liquids is controlled by the removal of vapor from the pool surface by airflow (van den Bosch, 1997). Thus, for boiling and non-boiling liquids the using of CFD code is reasonable to improve the prediction of convective heat and mass transfer, respectively. The present paper aims to validate the developed pool evaporation CFD model using literature and our own experimental data and to study the effect of the pool sizes on the evaporation rate from volatile liquid spill using numerical simulation. In the previous work (Galeev et al., 2013a) the pool evaporation model was validated using only experimental data on liquefied butane evaporation. However, before to use a mathematical model in risk assessment, the model must be carefully tested against experimental data from different sources. Therefore, in the present work we also compared simulation results both with the experimental data on evaporation of ethanol and cyclohexane (Habib et al., 2009) and with our own experimental data on evaporation of acetone.

8 > 0 > > > > > ! > > o n  <1 Ksþ  2:25 þ ln þ Cs $Ks $sin 0:4258$ ln Ksþ  0:811 DB ¼ k 87:75 > > > > >  >1  > > : ln 1 þ Cs $K þ s k

The pool evaporation rate sensitivity due to pool length has not been fully addressed in the literature. In the work (Khajehnajafi and Pourdarvish, 2011) the pool size determines the different flow regimes over the flat surface. Depending on the intensity of turbulence in the main stream, the flow near the surface becomes turbulent at some downstream distance (Schlichting, 1968). The point of transition from laminar to turbulent boundary layer flow may be estimated in terms of the length Reynolds number, ReL ¼ ULn1, where U is the main stream air velocity (ms1), L is the pool length in the direction of the airflow (m) and n is the kinematic viscosity of air (Nielsen et al., 1995). For the length Reynolds number, the critical value of 5  105 is often assumed (Incropera et al., 2006). The critical values for the length Reynolds numbers from 1  105 to 3  105 have been reported (Pasquill, 1943; Bird et al., 1960; Schlichting, 1968; Coulson and Richardson, 1993; Nielsen et al., 1995). When considering accidental spills, the pool sizes are large and ReL exceeds the critical value and the prevalence of turbulent mechanism of the transport of vapor away from the evaporating pool is expected. In the model (Mackay and Matsugu,

1973; Kawamura and Mackay, 1987) it is assumed, that evaporation rate Jg,s depends on pool size L as Jg,s ~ L0.11. The exponent on the pool dimension L was obtained by evaporating water from pools of different sizes (Mackay and Matsugu, 1973). This exponent does not reflect the buoyancy effects because the evaporation rate of water is low. The buoyancy effects are caused by the density difference between the evaporating component and ambient air and they may be significant when a pool of volatile liquid has large sizes and wind speed is low. It is evident that further studies are necessary to bring a better understanding to the role of pool dimensions and buoyancy effects in evaporation of volatile liquid. 2. Pool evaporation model and its validation The mass flow of vapor from the pool surface Jg,s due to evaporation was calculated on the basis of the standard wall functions (Fluent, 2006) taking into account the correction for Stefan flow:




Jg;s

 Yg;s  Yg;P rCm0:25 k0:5 P ¼K Yþ (

Y

þ

¼

Sc$yþ   Sct uþ þ PC

  yþ < yþ C  ; yþ > yþ C

(1)

(2)

yþ ¼

rCm0:25 k0:5 P yP ; m

(3)

uþ ¼

1  þ ln Ey  DB; k

(4)



Ksþ < 2:25



ð2:25 < Ksþ < 90Þ ; 

(5)

 Ksþ > 90

Ksþ ¼

r$Ks $Cm0:25 k0:5 P ; m "

PC ¼ 9:24

K¼

Sc Sct

3=4

(6) #

 1 ½1 þ 0:28expð  0:007Sc=Sct Þ;




 ln 1  Yg;P 1  Yg;s : Yg;s  Yg;P

(7)

(8)

It is assumed that the transition between the fully turbulent region and the viscous layer near the wall occurs at a value yþ C of 11, independent of the concentration of the evaporating component. In Eqs. (1)e(6) the coefficient Cm is assumed to be a constant: Cm ¼ 0.09. In the gas dispersion model, this coefficient is

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

determined by relating Cm to the mean flow deformation as adopted in the Realizable keε turbulence model (Shih et al., 1995). The change of local liquid temperature is calculated from the equation

vTliq qa þ qgrd þ qs  qp þ qar  Jg;s $DHg ¼ ; vt CP;liq $mliq

(9)

where qa is the heat flux from the atmosphere, W/m2; qgrd is the heat flux from the ground to the liquid, W/m2; qs is the heat flux from solar radiation, W/m2; qp is the heat flux emitted by the pool due to long-wave radiation, W/m2; qar is the heat flux absorbed by the pool due to long-wave radiation from the atmosphere, W/m2 The heat flux from the ground to the liquid phase, qgrd ¼ lgrd(vTgrd/vy)y¼0 is found from the numerical solution of three-dimensional nonstationary heat conduction equation for the substrate as was done in the work (Galeev et al., 2013a, b). The heat flux from the atmosphere qa is calculated using wall functions (Fluent, 2006). The qs, qp, qar were determined from the equations given in the paper (Kawamura and Mackay, 1987). The change in liquid mass is calculated from the equation:

vmliq ¼ Jg;s : vt

(10)

The Equations (9) and (10) were incorporated into the CFD code Fluent as user-defined scalar (UDS) transport equations without convection and diffusion terms.

41

The mass flux Jg,s is used as a boundary condition for the dispersion model in the pool area. The atmospheric dispersion of vapor is calculated by solving the three-dimensional Reynoldsaveraged NaviereStokes equations and the energy and species equations. The governing equations are closed using the Realizable keε equations for turbulence. To calculate the dispersion of vapor in the atmosphere, it is important to define the correct profiles for the mean velocity, turbulence kinetic energy and its dissipation rate on the inlet to the computational domain. These inlet profiles should be maintained throughout a calculation domain. The fully developed profiles provided by Richards and Hoxey (1993) are mathematically consistent, i.e. they are a solution of the mathematical models describing a homogeneous atmospheric boundary layer (Parente et al., 2011). However, the assumption of constant kinetic energy, k is not consistent with wind-tunnel measurements (Robins et al., 2001; Yang et al., 2009) where a variation of k with height is generally observed. Yang et al. (2009) proposed a new set of inlet conditions where the k profile is a function of height. Gorle, et al. (2009) proposed formulations for the turbulence model coefficients to ensure stream-wise homogeneity while using the k profile proposed by Yang et al. (2009). In the present work the inlet profiles were determined through repeated calculations when the velocity and turbulence parameters profiles deduced at the outlet section of the computational domain were set at its inlet boundary and the computation was repeated until the profiles for outlet and inlet boundaries were close. The plots in Figs. 1e3 show the comparison between calculated profiles and profiles by Yang et al.

Fig. 1. Profiles of (a) velocity U, (b) turbulence kinetic energy k, (c) turbulence kinetic energy dissipation rate ε and (d) turbulent viscosity mt from numerical calculation (solid line) and from Yang et al. (2009) (dashed line) at the wind velocity of 1 m/s at the height of 10 m.

42

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

Fig. 2. Profiles of (a) velocity U, (b) turbulence kinetic energy k, (c) turbulence kinetic energy dissipation rate ε and (d) turbulent viscosity mt from numerical calculation (solid line) and from Yang et al. (2009) (dashed line) at the wind velocity of 2.5 m/s at the height of 10 m.

(2009). These profiles were obtained at wind velocities of 1, 2.5 and 5 m/s and under neutral atmospheric stability. The plots show that there is no significant difference in the velocity and ε at the wall while there is the difference in k. The similar discrepancy between experimental and simulated profiles of k is observed in the work (Mokhtarzadeh-Dehghan et al., 2012). At close velocity profiles, the impact of using different turbulence parameters profiles on the results of vapor dispersion calculation will be governed by the parameter mt/Sct. Due to the fact that coefficient Cm ¼ 0.028 in correlations of Yang et al. (2009) is much lower than the same coefficient in the used Realizable keε model (is equal to about 0.09e0.14 depending on height from wall), the turbulent viscosities at the wall practically coincide. By optimizing the value of turbulent Schmidt number, it is possible to reach a good agreement between the predicted and measured parameters for both approaches. In the work (Galeev et al., 2013b), the results of calculation of neutral gas dispersion with Sct ¼ 0.5 and calculated profiles, agreed better with experiment data than with Sct ¼ 0.7, whereas in the work (Galeev and Ponikarov, 2014) the good agreement between simulation and experiment on heavy gas dispersion was obtained at Sct ¼ 0.7. In the present work numerical simulations are performed with calculated profiles for mean velocity and turbulence parameters and with Sct ¼ 0.7. To validate the CFD evaporation model, the predicted results were compared with the experimental data obtained by Habib et al. (2009). They carried out experiments in the open air for evaporation of ethanol and cyclohexane at constant temperatures. One facility comprised a flat terrain without buildings whilst the other was situated in a very rough area. The pool consisted of a circular basin with a diameter of 0.74 m which was insulated from the

ground. The experimental data for cyclohexane are presented in the work (Khajehnajafi and Pourdarvish, 2011). The difference between the pool evaporation model used in the present study and the one used in our previous works (Galeev et al., 2013a; Galeev et al., 2013b; Galeev and Ponikarov, 2014) is that the model in the present study takes into account roughness for very rough area when the roughness height larger than the height of the centre point of the wall-adjacent cell. The two parameters of the rough wall model, i.e. the roughness height Ks and roughness constant Cs, were determined according to the relationship between Ks, Cs and z0 (Blocken et al., 2007):

KS ¼ 9:793$z0 =Cs :

(11)

In the calculation, it is assumed that the aerodynamic roughness of flat terrain z0 is equal to pool surface roughness z0 ¼ 0.0002 m (Ks ¼ 0.00195 m) and Cs ¼ 1. For the very rough area the aerodynamic roughness is set equal to 0.04 m (Ks ¼ 0.39 m). Although in the work (Troen and Petersen, 1989) for area with many trees and buildings it is recommended to use the value of z0 ¼ 0.1 m we assumed a lower value to avoid the divergence problems during iterative process. Ks should be smaller than the height of the centre point of the wall-adjacent cell. The height of wall-adjacent cells was 0.1 m, therefore, in the simulation of a very rough area, the roughness height, Ks, was set to the value of 0.05 m and from Eq. (11), the corresponding value of constant Cs was 7.83. Fluent Manual states that constant Cs should be smaller or equal to 1 and it is not possible to set a higher value via a user interface. Therefore, a higher value is defined through a User Defined

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

43

Fig. 3. Profiles of (a) velocity U, (b) turbulence kinetic energy k, (c) turbulence kinetic energy dissipation rate ε and (d) turbulent viscosity mt from numerical calculation (solid line) and from Yang et al. (2009) (dashed line) at the wind velocity of 5 m/s at the height of 10 m.

Function. This approach to set roughness parameters in Fluent was used in the gas dispersion calculation by Labovsky and Jelemensky (2011). The computational domain with dimensions 100 m  500 m  100 m in the x, y, z directions was used. A structured non-uniform grid consisting of parallelepiped cells was generated. The grid cells were clustered in the region of the pool. The number of grid cells was approximately 96,000. The second-order scheme was used for the discretization of the advective terms of the governing equations. The PRESTO! algorithm (Fluent, 2006) to calculate the pressure and the SIMPLE algorithm (Patankar, 1980) to correlate the velocity field and pressure were used. The convergence criterion for dispersed material concentration was 1  104, for other variables was 1  103. The maximum time step of 2 s was used. Accuracy of the models was estimated by the value of percent error (%):

. d ¼ 100$ Jg;s;calc  Jg;s;expt Jg;s;expt :

(12)

Tables 1 and 2 show experimental data and prediction values for evaporation of ethanol and cyclohexane at different temperatures and wind speeds. The computed values of the Mackay-Matsugu equation are obtained with a difference d to the measured ones between 27.5 and 61.1%, with an average value of 47.4%. The computed values of the numerical model are obtained with a difference d to the measured ones between 0.6 and 47.4%, with an average value of 12.7%.

To validate the evaporation model, our own experimental data for acetone evaporation were also used. The liquid was poured into the pan with sizes of 0.6 m  0.4 m  0.035 m. The air velocity was measured at the height of 2 m by an acoustical anemometer with an error ± 0.2 m/s. An averaging time of the wind speed was 2 min. To measure the change in mass of the pan with the liquid, the MettlerToledo XP8002S balance was used with an accuracy of 0.01 g. To measure the temperature of the liquid the thermocouple chromelcopel with an accuracy of 0.1  C was used. The measurements were carried out at the air temperature of 23  C and humidity of 40%. The initial temperature of acetone was equal to 23.5  C. Wind speed was varied during the experiment (Fig. 4). In the calculation the average wind speed for the considered time period was used, which was equal to 1.05 m/s. The solar flux was determined using the equation given in the paper (Kawamura and Mackay, 1987), with taking into account the latitude and longitude of terrain and cloud cover, and its value was 430 W/m2. The heat fluxes due to the long-wave radiation of the atmosphere and the pool surface have been taken into account through the use of formulae from the same source. The heat flux from the ground is assumed to be equal to zero because the pan with liquid was placed on the balance. The aerodynamic roughness of the area around the pool is set equal to 0.01 m (Ks ¼ 0.098 m), that corresponds to the area with short grass (Troen and Petersen, 1989). Fig. 5 shows that the model somewhat underestimates the evaporated mass. The difference between the measured and calculated evaporated mass is 17.5%. The difference between

44

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

3. Sensitivity study and discussion

Table 1 Experimental data and predictions for ethanol evaporation. Pool temperature, К Wind velocity at Mass flow, g/s 2 m height, m/s Experiment Mackay e CFD model Matsugu Eq. Flat terrain 310.15 309.65 325.15 Very rough terrain 305.65 310.15 324.65

1.4 1.9 1.7

0.456 0.533 1.02

0.681 0.855 1.59

0.393 0.523 1.06

1.4 1.8 1.8

0.488 0.679 1.32

0.663 0.99 1.95

0.536 0.9 1.96

Table 2 Experimental data and predictions for cyclohexane evaporation. Pool temperature, К Wind velocity at Mass flow, g/s 2 m height, m/s Experiment Mackay e CFD model Matsugu Eq. Flat terrain 303 310 317 Very rough terrain 303 310 317

2.71 3.05 3.49

1.08 1.63 2.33

1.74 2.49 3.63

1.11 1.66 2.55

1.71 1.62 1.4

0.967 1.35 1.67

1.39 1.80 2.13

1.15 1.46 1.68

experimental and predicted dependencies of liquid temperature is explained by the difference in evaporation rates. It should be noted, that up to about 10 min, the experimental and calculated values of the mass and temperature are close, although the wind velocity in calculation (1.05 m/s) was about 2.5 times greater than the one in the experiment. This could be explained by the fact that at very low wind speeds and at the small pool size the molecular mechanism of vapor removal from the evaporation surface could dominate, when the sensitivity of evaporation to wind velocity is not as high as in the turbulent regime. Conventionally, the Sherwood number, Sh, (the dimensionless parameter describing the mass transfer) is proportional to Re0.5 for laminar flow over a flat surface, while Sh ~ Re0.8 in turbulent flow (Incropera et al., 2006).

This section presents the numerical analysis of the effect of pool sizes on evaporation of hexane. Hexane has high volatility (at a liquid temperature of 303 К the hexane vapor pressure is equal to 25 kPa) and hexane vapor is 3 times heavier than the ambient air. These factors affect pool evaporation and vapor cloud dispersion processes. The pools with sizes of 10 m  10 m, 30 m  30 m, 60 m  60 m and 100 m  100 m are considered. The initial liquid layer height was equal to 0.05 m. The initial temperature of ambient air and pool substrate was 303 K. The numerical study included a series of calculations at wind speeds 1, 2.5 and 5 m/s at the height of 10 m. The solar flux had a range of 0e900 W/m2 (Khajehnajafi and Pourdarvish, 2011), therefore in calculation the average value qs ¼ 450 W/m2 was used. The geometry of the computational domain is shown in Fig. 6. The calculation was performed for half of the pool because of symmetry of the problem. The computational domain with dimensions 2000 m  500 m  500 m in the x-, y-, z-directions was used. A structured nonuniform grid consisting of parallelepiped cells was generated. The number of grid cells was approximately 370,000. The height of wall-adjacent cells was equal to 0.2 m. In downwind (x-axis) and crosswind (z-axis) directions, the grid cell size was 1 m over the pool surface. Additionally to the above domain, a domain under the pool containing the ground of 2 m depth was also included for in the calculations. The profiles of wind velocity, k and ε on the inlet to the computational domain (plane ABCD) were obtained through calculation and are shown in Figs. 1e3. On the lateral and upper boundaries (ABFE, BCGF and DCGH) of the domain the symmetry condition was specified, i.e. velocity component normal to the boundary and normal derivatives of the other variables assumed to be equal to zero. The bottom boundary (ADHE) was defined as a no-slip wall. At the outlet boundary (EFGH) the gauge pressure was fixed at zero. The view of used computational grid near the pool is shown in Fig. 7. The roughness coefficient value Cs ¼ 1 was applied in this study. The roughness height Ks for the pool surface was given as 0.00195 m and for the adjacent territory was 0.098 m. The liquid temperature was imposed as a thermal boundary condition on the wall under the pool while the wall around the pool was considered adiabatic. In the simulation, it was assumed that physical properties of liquid, gas and solid layer did not depend on temperature. The following values have been used: Dm,g ¼ 6.9$106 m2/sec; mg ¼ 6.6$106 kg/(m sec); ma ¼ 1.78$105 kg/(m sec); CP,g ¼ 1680 J/ (kg К); CP,a ¼ 1006 J/(kg К); CP,liq ¼ 2300 J/(kg К); lg ¼ 0.013 W/ (m К); la ¼ 0.0242 W/(m К); DН g ¼ 370,000 J/kg; lgrd ¼ 1.28 W/ (m К); С P,grd ¼ 1130 J/(kg К); rgrd ¼ 2300 kg/m3. Here, the Dm, m, CP, l, r, denote molecular diffusion coefficient, molecular viscosity, specific heat, thermal conductivity and density, respectively. The a, g and grd subscripts indicate air, evaporating gas and ground value, respectively. The vapor pressure Pg of hexane as a function of liquid temperature Tliq was determined using the formula (Green and Perry, 2008):

   .    E ; Pg Tliq ¼ exp A þ B Tliq þ C$ln Tliq þ D$Tliq

Fig. 4. Measured wind velocity.

(13)

where A ¼ 104.65, B ¼ 6995.5, C ¼ 12.702, D ¼ 0.000012381, E ¼ 2. In Fig. 8 the plots of specific evaporated mass mvap and volumeaveraged liquid temperature Tliq against time are presented. The plots show that the specific evaporated mass is higher when the pool sizes are larger, in other conditions being equal. The influence

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

45

Fig. 5. Measurements (dashed line) and predictions (solid line) of liquid mass and temperature.

Fig. 6. Schematic representation of the computational domain.

Fig. 7. The computational grid near the pool.

of pool sizes on the evaporation decreases when the wind velocity increases. As the pool sizes increase the vapor concentration above the surface of the evaporating liquid and, consequently, the vaporair cloud density increases. The negative vertical density gradient appearing above pool surface causes the suppression of turbulence (stable stratification) that creates a large resistance to the transport of vapor away from the evaporating pool and results in a slow

evaporation rate. In addition to the buoyancy effect, the increasing of vapor concentration above the pool surface reduces the driving force of the evaporation. The graphs of temperature change are consistent with the time dependence of mvap. A higher evaporation rate results in a lower pool temperature, since the heat losses increase with the increasing of evaporation rate. At the wind velocity of 1 m/s, the temperature rises throughout the considered period of time because heat flux from the sun to pool prevails over the heat losses due to evaporation. The regression analysis of numerical results has shown that at the wind speed of 1 m/s the average specific evaporated mass depends on a pool size as mvap ~ L0.553 with accuracy of approximation of R2 ¼ 0.9981, at wind speed of 2.5 m/s as mvap ~ L0.318 (R2 ¼ 0.9153) and at wind speed of 5 m/s as mvap ~ L0.0578 (R2 ¼ 0.9755). Thus, at low wind velocities the stronger dependence of evaporation rate on pool size is observed, than the one (mvap ~ L0.11) adopted in the Mackay-Matsugu equation commonly used in hazard assessment (Kawamura and Mackay, 1987; Mackay and Matsugu, 1973). Table 3 presents the values of the physical properties of gases and liquids for the upper and lower limits of the variation range of liquid temperature. At these temperatures the physical properties do not differ significantly, that confirms the validity of the assumption of constant physical properties. The molecular diffusion coefficient was determined using Gilliland's correlation (Gilliland, 1934) and other values were taken from the ChemCAD program database. In Fig. 9 the time dependences of the relative turbulent viscosity and the driving force (DYg ¼ Yg,s-Yg,P) are shown. The relative turbulent viscosity is defined as the dynamic coefficient of average turbulent viscosity above the pool surface (y ¼ 0.1 m) divided by the turbulent viscosity in the absence of the disturbing effect of the source. The relative turbulent viscosity for an undisturbed flow (when there is no hexane vapor) is equal to 1.0 for each of the wind velocities. In the initial period of evaporation, a sharp decay of turbulent viscosity occurs due to the formation of the negative density gradient above pool surface. The importance of buoyancy effect rises as the air velocity decreases and pool sizes increase. The slight growth of the eddy viscosity at wind speeds of 2.5 m/s and 5 m/s is attributed to the reduction of liquid temperature, and, as consequence, the reduction of mass flux of vapor from the pool surface. At low wind velocities, the driving force of evaporation significantly drops as pool sizes increase due to the increasing of the concentration above the pool surface. Since in the model the driving force is both in the numerator and the denominator of the

46

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

Fig. 8. Dependencies of specific evaporated mass (left) and average liquid temperature (right) against time at wind velocities of 1 m/s (a), 2.5 m/s (b) and 5 m/s (c).

Table 3 Physical properties of gases and liquid. Physical property

Dm,g, m2/sec mg, kg/(m sec) ma, kg/(m sec) CP,g, J/(kg К) CP,a, J/(kg К) CP,liq, J/(kg К) lg, W/(m К) la, W/(m К) DН g, J/kg

Temperature, K

Variation, %

286

308

6.48$106 6.5$106 1.78$105 1590 1005 2215 0.0115 0.024 380,000

7.24E-06 6.7$106 1.88$105 1690 1008 2310 0.0135 0.026 362,500

11.7 3.1 5.6 6.3 0.3 4.3 17.4 8.3 4.6

formula for the diffusion flux, its effect on the evaporation rate will be determined by the factor ln((1Yg,P)/(1Yg,s)). At the wind speed of 1 m/s when the pool sizes increase from 10 m  10 m to 100 m  100 m the average factor reduces by 11%, at the wind speed of 2.5 m/s e by 6.4% and at the wind speed of 5 М/с e by 3.5%, while average specific evaporated mass mvap reduces by 71.5%, 53% and 12.8%, respectively. Thus, changes in the evaporation rate with increasing of pool sizes are determined mainly by the effect of turbulence suppression over the pool surface, and not by changes in the driving force. In the paper (Britter, 1989), a criterion that indicates when the plume may be considered passive, i.e. when the influence of the density difference on the inertia small and may be neglected, is given. For continuous source this criterion is

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

47

Fig. 9. Dependencies of average relative turbulent viscosity (left) and driving force (right) at wind speeds of 1 m/s (a), 2.5 m/s (b) and 5 М/с (c).


 1=3 . U  0:15; B ¼ g00 $q0 L

(14)

where g00 ¼ gðrg  ra Þ=ra ; q0 is the volume flow rate; L is the source dimension and U is the ambient velocity. For determining the relative importance of the flow regimes (buoyancy-dominated, stably stratified and passive dispersion) the Richardson number, Ric, is also used (Havens, 1992):

. Ric ¼ g00 $H u2* ;

(15)

where H ¼ q0/UL is the characteristic cloud height and u* is the friction velocity If Ric  30 the flow is negative buoyancy dominated, if 1Ric  30, the shear flow is stably stratified and if Ric  1 the passive dispersion occurs (Havens, 1992).

In Table 4 the values of B и Ric for each of considered cases are given. In the calculation of these criteria, the average values of the gas flow rate have been substituted in expressions (14) and (15). The values of Britter's criterion are less than 0.15 only when the pool sizes are 10 m  10 m and the wind speed is 5 m/s. This result is confirmed by the graphs (see Fig. 9), which show that the change of turbulence due to buoyancy effects is small at wind speed of 5 m/ s and at the pool sizes of 10 m  10 m, while in the rest cases the appreciable suppression of turbulence (stable stratification) above pool surface takes place. Generally, at the wind speed of 5 m/s the Britter's criterion is less, or slightly greater than the critical value, whereas at low wind speeds and extended pools it significantly exceeds the critical value. The values of Richardson number are also consistent with the numerical calculations. The Ric in most cases is in the range from 1 to 30 indicating the prevalence of stably stratified shear flow. The Ric < 1 is obtained only for the case where

48

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49

yP

Table 4 Values of B and Ric depending on pool sizes and wind velocity. Pool sizes

10 m  10 m 30 m  30 m 60 m  60 m 100 m  100 m

z0

Britter's criterion B (Richardson number Ric) 1 m/s

2.5 m/s

5 m/s

0.35 0.42 0.46 0.51

0.20 0.27 0.31 0.33

0.12 0.17 0.21 0.25

(12.4) (20.3) (26.7) (35.9)

(2.2) (5.5) (8.4) (10.3)

(0.5) (1.3) (2.6) (4.1)

the pool sizes are 10 m  10 m and the wind speed is 5 m/s and Ric > 30 is obtained when the pool sizes are 100 m  100 m and the wind speed is 1 m/s.

4. Conclusion The validation of the developed pool evaporation CFD model against both literature and our own experimental data is performed. A good agreement is obtained between the measured and the predicted evaporation rates. Based on the proposed model a sensitivity analysis is conducted to determine the effect of pool sizes on evaporation of a volatile liquid (hexane). The numerical analysis has shown that the evaporation rate falls as the pool sizes increase due to the enhancement of buoyancy effects above the pool surface. The importance of the pool sizes effect decreases as the air velocity increases. It was found that in the dependence of evaporation rate on pool size Jg,s ~ Ln the exponent on pool size varies between 0.553 and 0.0578 for wind speeds between 1 and 5 m/s, while in risk analysis handbooks it is generally accepted that Jg,s is proportional to L0.11 regardless of wind speed.

Nomenclature

DB CP Cs Cm Dm E DH Jg,s

roughness function specific heat, J/(kg K) roughness constant coefficient in the turbulence model molecular diffusion coefficient, m2/s empirical constant equal to 9.1 heat of vaporization, J/kg mass flux of vapor from the pool surface due to evaporation, kg/(m2 s) K correction factor for Stefan flow Ks roughness height, m Kþ non-dimensional roughness height s k turbulent kinetic energy, m2/s2 L characteristic pool size, m mliq mass of the liquid per unit surface area of the pool, kg/m2 Pg(Tliq) vapor pressure of evaporating component at liquid temperature Tliq PC function which takes into account the resistance of the diffusion sublayer to mass transfer Sc and Sct molecular and turbulent Schmidt numbers, respectively T absolute temperature, K Tþ non-dimensional temperature t time, s U wind velocity, m/s uþ non-dimensional velocity Yg mass fraction of the evaporating component in the gas phase, kg/kg Yþ non-dimensional mass-fraction yþ non-dimensional distance yþ non-dimensional diffusion sublayer thickness C

normal distance from the pool surface to the neighboring node of the computational grid aerodynamic roughness, m

Greek symbols ε turbulence kinetic energy dissipation rate, m2/s3 k von Karman constant equal to 0.41 l coefficient of thermal conductivity, W/(m K) m coefficient of molecular viscosity, kg/(m s) mt coefficient of turbulent viscosity, kg/(m s) r density, kg/m3 Subscripts no index vapor-air mixture a air g evaporating component grd ground liq pool P centroid of the wall-adjacent cell s surface of pool References Bird, R.B., Stewart, W.E., Lightfoot, E.N., 1960. Transport Phenomena. Wiley & Sons, New York. Blocken, B., Stathopoulos, T., Carmeliet, J., 2007. CFD simulation of the atmospheric boundary layer: wall function problems. Atmos. Environ. 41, 238e252. Brighton, P.W.M., 1987. Evaporation from a Plane Liquid Surface into a Turbulent Boundary Layer. UK Atomic Energy Authority, Safety and Reliability Directorate. SRD/HSE-report R375. Britter, R.E., 1989. Atmospheric dispersion of dense gases. Annu. Rev. Fluid Mech. 21, 317e344. Cebeci, T., Chang, K.C., 1978. Calculation of incompressible rough-wall boundarylayer flows. AIAA J. 16, 730. Churchill, S.W., 1976. A comprehensive correlating equation for forced convection from flat plates. AIChE J. 22, 264e268. Coulson, J.M., Richardson, J.F., 1993. Chemical Engineering (Fluid Flow, Heat Transfer and Mass Transjk), fourth ed. Pergamon Press, Oxford. Desoutter, G., Habchi, C., Cuenot, B., Poinsot, T., 2009. DNS and modeling of the turbulent boundary layer over an evaporating liquid film. Int. J. Heat Mass Transf. 52, 6028e6041. Fluent, 2006. Fluent 6.3 User's Guide. Fluent Inc, Lebanon, New Hampshire, USA. Galeev, A.D., Starovoytova, E.V., Ponikarov, S.I., 2013a. Numerical simulation of the consequences of liquefied ammonia instantaneous release using FLUENT software. Process Saf. Environ. Prot. 91, 191e201. Galeev, A.D., Salin, A.A., Ponikarov, S.I., 2013b. Consequence analysis of aqueous ammonia spill using computational fluid dynamics. J. Loss Prev. Process Ind. 26, 628e638. Galeev, A.D., Ponikarov, S.I., 2014. Numerical analysis of toxic cloud generation and dispersion: a case study of the ethylene oxide spill. Process Saf. Environ. Prot. 92, 702e713. Gilliland, E.R., 1934. Diffusion coefficients in gaseous systems. Ind. Eng. Chem. 26, 681e685. Gorle, C., van Beeck, J., Rambaud, P., van Tendeloo, G., 2009. CFD modeling of small particle dispersion: the influence of the turbulence kinetic energy in the atmospheric boundary layer. Atmos. Environ. 43, 673e681. Green, D.W., Perry, R.H., 2008. Chemical Engineer's Handbook, eighth ed. McGrawHill. Habib, A., Schalau, B., Acikalin, A., Steinbach, J., 2009. Transient calculation of the boundary layer flow over spills. Chem. Eng. Technol. 32, 306e311. Havens, J.A., 1992. Review of dense gas dispersion field experiments. J. Loss Prev. Process Ind. 5, 28e41. Incropera, F.P., DeWitt, D.P., Bergman, T.L., Lavine, A.S., 2006. Fundamentals of Heat and Mass Transfer. Wiley. Kapias, T., Griffiths, R.F., 1999. Modelling the Behaviour of Spillages of Sulfur Trioxide and Oleum. Health and Safety Executive. Contract Research Report 217/ 1999. Kawamura, P.I., Mackay, D., 1987. The evaporation of volatile liquids. J. Hazard. Mater. 15, 343e364. Khajehnajafi, S., Pourdarvish, R., 2011. Correlations for mass transfer from a liquid spill: comparisons and recommendations. Process Saf. Prog. 30, 178e184. Labovsky, J., Jelemensky, L., 2011. Verification of CFD pollution dispersion modelling based on experimental data. J. Loss Prev. Process Ind. 24, 166e177. Liu, Y., Olewski, T., Vechot, L., Gao, X., Mannan, S., 2011. Modelling of a cryogenic liquid pool boiling using CFD code. In: Proc. 14th Annual Symposium, Mary Kay O'Connor Process Safety Center “Beyond Regulatory Compliance: Making Safety Second Nature”. Texas A&M University, College Station, Texas, USA, pp. 25e27. October 2011.

A.D. Galeev et al. / Journal of Loss Prevention in the Process Industries 38 (2015) 39e49 Mackay, D., Matsugu, R.S., 1973. Evaporation rates of liquid hydrocarbon spills on land and water. Can. J. Chem. Eng. 51, 434e439. Mokhtarzadeh-Dehghan, M.R., Akcayoglu, A., Robins, A.G., 2012. Numerical study and comparison with experiment of dispersion of a heavier-than-air gas in a simulated neutral atmospheric boundary layer. J. Wind Eng. Ind. Aerodyn. 110, 10e24. Nielsen, F., Olsen, E., Fredenslund, A., 1995. Prediction of isothermal evaporation rates of pure volatile organic compounds in occupational environments d A theoretical approach based on laminar boundary layer theory. Ann. Occup. Hyg. 39, 497e511. , C., van Beeck, J., Benocci, C., 2011. Improved keε model and wall Parente, A., Gorle function formulation for the RANS simulation of ABL flows. J. Wind Eng. Ind. Aerodyn. 99, 267e278. Pasquill, F., 1943. Evaporation from a plane, free liquid surface into a turbulent air stream. Proc. R. Soc. Lond. A 182, 75e95. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington, DC. Raimundo, A.M., Gaspar, A.R., Oliveira, A.V.M., Quintela, D.A., 2014. Wind tunnel measurements and numerical simulations of water evaporation in forced convection airflow. Int. J. Therm. Sci. 86, 28e40. Richards, P., Hoxey, R., 1993. Appropriate boundary conditions for computational wind engineering models using the kee turbulence model. J. Wind Eng. Ind. Aerodyn. 46e47, 145e153. Robins, A., Castro, I., Hayden, P., Steggel, N., Contini, D., Hesit, D., 2001. .A wind tunnel study of dense gas dispersion in a neutral boundary layer over a rough surface. Atmos. Environ. 35, 2243e2252. Rong, L., Nielsen, P.V., Zhang, G.Q., 2010. Experimental and numerical study on effects of airflow and aqueous ammonium solution temperature on ammonia

49

mass transfer coefficient. J. Air Waste Manag. Assoc. 60, 419e428. Rong, L., Elhadidi, B., Khalifa, H.E., Nielsen, P.V., Zhang, G.Q., 2011. Validation of CFD simulation for ammonia emissions from an aqueous solution. Comput. Electron. Agric. 75, 261e271. Saha, C.K., Wu, W., Zhang, G., Bjerg, B., 2011. Assessing effect of wind tunnel sizes on air velocity and concentration boundary layers and on ammonia emission estimation using computational fluid dynamics (CFD). Comput. Electron. Agric. 78, 49e60. Schlichting, H., 1968. Boundary-layer Theory, sixth ed. McGraw-Hill, New York. Shih, T.H., Liou, W.W., Shabbir, A., Zhu, J., 1995. A new keε eddy-viscosity model for high Reynolds number turbulent flows e model development and validation. Comput. Fluids 24, 227e238. Troen, I., Petersen, E.L., 1989. European Wind Atlas. Risø National Laboratory, Roskilde, Denmark. van den Bosch, C.J.H., 1997. Pool evaporation. In: van den Bosch, C.J.H., Weterings, R.A.P.M. (Eds.), CPR 14E, Methods for the Calculation of Physical Effects Due to Releases of Hazardous Materials (Liquids and Gases), Yellow Book. Committee for the Prevention of Disasters, The Hague, 3.1e3.128. chot, L., Olewski, T., Osorio, C., Basha, O., Liu, Y., Mannan, M.S., 2013. Laboratory Ve scale analysis of the influence of different heat transfer mechanisms on liquid nitrogen vaporization rate. J. Loss Prev. Process Ind. 26, 398e409. Vik, T., Pettersson Reif, B.A., 2011. Implementation of a New and Improved Evaporation Model in Fluent. Norwegian Defence Research Establishment (FFI). FFIrapport 2011/00116. Yang, Y., Gu, M., Chen, S., Jin, X., 2009. New inflow boundary conditions for modeling the neutral equilibrium atmospheric boundary layer in computational wind engineering. J. Wind Eng. Ind. Aerodyn. 97, 88e95.