Energy & Fuels 2005, 19, 2535-2544
2535
Fuel Gas Dispersion under Cryogenic Release Conditions Spyros Sklavounos and Fotis Rigas* National Technical University of Athens, School of Chemical Engineering, 15700 Athens, Greece Received February 8, 2005. Revised Manuscript Received July 21, 2005
In this paper, a computational approach to the atmospheric dispersion that is associated with cryogenic releases of the gaseous fuels hydrogen and natural gas was attempted. A computational fluid dynamics (CFD) method was followed, using the code CFX 5.7 for the simulation of largescale liquefied hydrogen and liquefied natural gas spill trials, which provided successful prediction of the concentration histories, against experimental records. The main outcome that was determined through the results was that fuel gases, such as natural gas and hydrogen, disperse as heavy gases rather than light gases during a cryogenic release, despite their positive buoyancy under normal conditions. This conclusion is of major importance from the standpoint of safety, in regard to chemical installations that handle liquefied gases, given that the gravity-driven dispersion of flammable gases increases the risk of ignition and, hence, accidental fire or explosion. Thus, the behavior of the released gas should be known when selecting the appropriate dispersion model that will be used in industrial safety studies. Furthermore, it was concluded that CFD may be effectively used in risk analysis procedures, providing reliable results even for the complex case of nonisothermal atmospheric gas release.
1. Introduction Recently, a big effort has been dedicated to the introduction of new energy vectors in the world energy market for the purpose of partial substitution of conventional fuels.1,2 Besides electrical energy, major interest has been focused on natural gas and renewable energy sources such as hydrogen, huge amounts of which are used in the chemical and power industries. Handling, storage, and transportation of the aforementioned gases are processes that usually include their liquefaction for economical transportation and storage, because of their very low gas densities under ambient conditions.3 The liquefaction is accomplished at significantly low temperatures, together with high-pressure application. However, a major safety issue arises regarding the dispersion behavior of such gases while being released under those extreme conditions, because of the possibility for them to be involved in vapor cloud explosions and fires.4 The behavior of the released gas determines the appropriateness (or not) of several empirical dispersion models commonly used in plant safety studies and land use planning procedures,5 although it has an important role in the establishment * Author to whom correspondence should be addressed. Tel.: +302107723267. Fax: +302107723163. E-mail:
[email protected]. (1) Green Paper: Towards a European Strategy for the Security of Energy Supply; Report COM(2002) 321, Commission of the European Communities, Brussels, Belgium, 2002. (2) Momirlan, M.; Veziroglu, T. N. Renew. Sust. Energy Rev. 2002, 6, 141-179. (3) Carson, P.; Mumford, C. Hazardous Chemicals Handbook; Butterworth: Oxford, U.K., 1994. (4) Deaves, D. M. J. Loss Prev. Process Ind. 1992, 5, 219-227. (5) Lees, F. P. Loss Prevention in the Process Industries; Butterworth: Oxford, U.K., 1996.
of accident prevention and mitigation measures in industrial sites. The interest for the assessment of hazards associated with accidental fuel gas cryogenic spills has led famous research institutes in the past to sponsor experiments that would focus on studying these incidents. Thus, cryogenic gas releases have been studied experimentally in both laboratory-scale6 and large-scale7-9 tests, providing useful findings for the way that gas dispersion proceeds and a unique series of experimental data that could be used for validation purposes of computer modeling procedures. A comprehensive review of data on cryogenic liquid spills has been presented by Thyer in 2003.10 What dominates after the spill is the violent vaporization of the liquefied substance, which leads to the formation of a cold pseudo-dense vapor cloud. The atmospheric dispersion is affected by heat-transfer processes that define the temperature and, hence, the density of the cloud. Although routine box models are capable of estimating gas dispersion from cryogenic spills, they are subject to certain deficiencies (for (6) Verfondern, K.; Dienhart, B. Int. J. Hydrogen Energy 1997, 22, 649-660. (7) Burro Series Data Report; LLNL/NWC Report No. UCID-19075, V. 1-2, Lawrence Livermore National Laboratory, Berkeley, CA, 1982. (8) Chirivella, J. E.; Witcofski, R. D. Experimental Results from Fast 1500-Gallon LH2 Spills. In Cryogenic Properties: Processes and Applications; Kidnay, A. J., Hiza, M. J., Frederking, T. H. K., Kerney, P. J., Wenzel, L. A., Eds.; AIChE Symposium Series, No. 251; American Institute of Chemical Engineers (AIChE): New York, 1986; pp 120140. (9) Puttock, J. S.; Blackmore, D. R.; Colenbrander, G. W. J. Hazard. Mater. 1982, 6, 13-41. (10) Thyer, A. M. J. Hazard. Mater. 2003, A99, 31-40.
10.1021/ef0500383 CCC: $30.25 © 2005 American Chemical Society Published on Web 08/25/2005
2536
Energy & Fuels, Vol. 19, No. 6, 2005
instance, one-dimensional modeling).11 In the meantime, the rapid increase in computational power allows for the utilization of more-powerful tools, such as the computational fluid dynamics (CFD) codes. Hence, in recent papers.12,13 CFD codes have been considered to be useful approximation tools for particular gas dispersion cases, with encouraging results. In this paper, a CFD approach of cryogenic gas release and dispersion is presented, in which the dispersion simulation of liquefied natural gas (LNG) and liquefied hydrogen (LH2) spill trials was attempted. The former were performed in 1980 at the Naval Weapons Center (NWC) in China Lake, CA, and are widely known as the Burro series tests. The experiments involved eight spills of LNG, each covering an area of ∼40 m3, onto an artificial basin filled with water. The LH2 spill trials were performed by NASA in 1980 at the White Sands test facility in New Mexico. They included seven spills of LH2, each covering an area of ∼5.1 m3, onto a slightly elevated solid platform. Both trials are extensively described, and relevant experimental data are provided in published material.14,15 For the computations, the CFD code CFX-5.7 was utilized; it performed three-dimensional calculations for the mass, energy, and momentum balances, as well as turbulent atmospheric mixing. CFX is a general-use code and is not calibrated for the solution of specific problems. However, it has been validated so far against isothermal heavy gas dispersion trials.16 In all runs, the thermodynamic properties of the gases were defined as a function of temperature using bibliographic data,17,18 to take into account the effect of cloud temperature variation during gas cloud dispersion. The computational results were compared with the experimental measurements, revealing a reasonably good agreement with the records from the trials. Calculations of concentration integrals with time at certain points within the domain showed that the code predicted the total dose (total fuel mass flow) within an error of 10%. The calculation of statistical performance measures confirmed that the code predicted the experimental concentration measurements within a factor of 2, whereas the corresponding scatter plot of predicted versus observed values showed the agreement between the computational and experimental concentration histories. Visualizations of the results showed that both gases disperse as heavy gases rather than light gases for considerable distances before they lift off, despite the fact that they are much lighter than air under normal (11) Duijm, N. J.; Carissimo, B. Evaluation Methodologies for Dense Gas Dispersion Models. In The Handbook of Hazardous Materials Spills Technology; McGraw-Hill: Fingas, M., Ed.; New York, 2002; Chapter 19. (12) Swain, M. R.; Shriber, J.; Swain, M. N. Energy Fuels 1998, 12, 83-89. (13) Hanna, S. R.; Hansen, O. R.; Dharmavaram, S. Atmos. Environ. 2004, 38, 4675-4687. (14) Description and Analysis of Burro Series 40-m3 LNG Spill Experiments; LLNL/NWC Report No. UCRL-53186, Lawrence Livermore National Laboratory, Berkeley, CA, 1981. (15) Witcofski, R. D.; Chirivella, J. E. Int. J. Hydrogen Energy 1984, 9, 425-435. (16) Sklavounos, S.; Rigas, F. J. Hazard. Mater. 2004, A108, 9-20. (17) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids; McGraw-Hill: Singapore, 1988. (18) Jacobsen, R. T.; Penoncello, S. G.; Lemmon, E. W. Thermodynamic Properties of Cryogenic Fluids; Plenum: New York, 1997.
Sklavounos and Rigas
conditions. Consequently, it is confirmed that the risk for cloud ignition is much higher, because the flammable gases travel along the ground, remaining at low altitude. Moreover, a clear proportional relationship between gas concentration and temperature was observed. The relationship was different for the hydrogen and the natural gas, but kept constant with time for each substance and indicated the strong influence of the cold mixed gas on the temperature. 2. Problem Statement In contrast to isothermal releases, in which the gas is released in the open air under normal conditions and its atmospheric dispersion is dependent on the molar weight of the released material, cryogenic releases are accompanied by a rapid phase change, strong temperature variation, and, eventually, substantial density and buoyancy changes in the vapor cloud that is formed. The extremely low temperature of the vapor cloud (which initially equals the boiling temperature of the released material) leads to the formation of a cold, denser-thanair cloud that may travel significant distances before it lifts off. As a result, there is an increased risk of ignition (for a flammable gas, such as methane and hydrogen) or a toxic effect (for a toxic gas, such as ammonia), because the cloud remains at low height. The cold cloud may lose part of its negative buoyancy, because of the continuous heat input from the ground to the cloud, whereas the significance of this effect is reduced as the ground temperature decreases. Eventually, the initially heavy plume is heated, reducing its density, and it becomes buoyant, expanding with liftoff. Evidently, an effective consequence analysis of cryogenic release scenarios would demand the accurate modeling of heat-transfer processes and variability of the thermodynamic properties of the released substance. The box models being applied in routine risk analysis studies are capable of making such calculations, however, with certain deficiencies. They assume certain horizontal and vertical profile functions for the cloud distribution, which are based on the cloud height, width, and length. Thus, they account for the variability of properties only in the wind blow direction, whereas they lack the possibility to analytically represent a specific geometry that may have a significant effect on dispersion.19 A different approach could be sought by means of CFD. CFD codes are powerful tools capable of accurately solving complex flow phenomena with a minimum degree of empiricism. In the past, they have been applied in gas dispersion problems to provide numerical concentration predictions with valuable results.20,21 The CFX code used in this work solves the three-dimensional Reynolds-averaged Navier-Stokes equations for mass, momentum, and energy balances,22 and, among others, it incorporates the Shear Stress Transport (SST) model (19) Kunsch, J. P.; Fannelop, T. K. J. Hazard. Mater. 1995, 43, 169193. (20) Swain, M. R.; Grilliot, E. S.; Swain, M. N. Chem. Health Safety 1999, 6, 28-32. (21) Pereira, J. C. F.; Chen, X. Q. J. Hazard. Mater. 1996, 46, 253272. (22) Chung, T. J. Computational Fluid Dynamics; Cambridge University: Cambridge, U.K., 2002.
Fuel Gas Dispersion under Cryogenic Release
Energy & Fuels, Vol. 19, No. 6, 2005 2537
Table 1. Computational Domain Dimensions and Related Mesh Construction Parameters parameter
global element scale factor, m
refinement factor on spill area
smooth transition factor
domain dimensions (X × Y × Z)
LNG domain LH2 domain
2 8
0.25 0.12
1.15 1.1
600 m × 300 m × 50 m 200 m × 100 m × 45 m
for turbulence modeling.23 The numerical solution of these equations is based on the Finite Volume Method (FVM)24-26 and proceeds to the discretization of the computational domain into several volume elements that constitute an unstructured tetrahedral mesh. The second-order backward Euler transient scheme27 was used in the numerical solution to increase the accuracy of the results and reduce numerical diffusion. Analytical information about the incorporated models may be found in a previous paper;16 therefore, they are not discussed in the present paper. CFX code utilizes two approaches for modeling the logarithmic layer profile: the “wall function” method and the “low-Reynolds-number” method. The former uses empirical formulas that impose suitable conditions near the wall without resolving the boundary layer, thus saving computational cost, whereas the latter resolves the details of the boundary layer profile, thus leading to higher run-time requirements. The wall function method is used by default when the k- turbulence model is used in a simulation, indeed utilizing “scalable” wall functions, which have the advantage of being able to be applied on arbitrarily fine meshes. On the other hand, when a turbulence model based on the ω-equation applies (such as the SST model in the current work), the scalable wall function approach is used in conjunction with the low-Reynolds-number method as the mesh is refined, so that the accuracy of the solution of boundary layer profile is enhanced. 3. Description of the Trials 3.1. Liquefied Natural Gas Spills. The Burro Series of LNG spill trials were sponsored by the United States Department of Energy (U. S. DOE) and conducted by Lawrence Livermore National Laboratory (LLNL) and the Naval Weapons Center (NWC) at the China Lake site in California during the summer of 1980 to determine the transport and dispersion of vapor from spills of LNG on water. The eight tests that were performed were all nominally 40-m3 spills onto a water surface at spill rates fluctuating between 11 m3/min and 18 m3/ min. The water test basin was circular, having an average diameter of 58 m, whereas the water surface was lying ∼1.5 m below ground level. The LNG was fed from a cryogenic liquid storage tank through a 0.25-mdiameter stainless-steel pipeline and released straight downward onto the water surface at the center of the pond. A splash plate had been intentionally installed at the release point to limit penetration of the LNG into the water. The water is an efficient heat supply, so the LNG boiled instantly. Gas concentrations were mea(23) Park, S. H.; Kwon, J. H. AIAA J. 2004, 42, 1348-1357. (24) Versteeg, H. K.; Malalasekera, W. An Introduction to Computational Fluid DynamicssThe Finite Volume Method; Longman: New York, 1995. (25) Kallinderis, Y. Appl. Numer. Math. 1996, 20, 387-406. (26) Stynes, M. J. Comput. Appl. Math. 1995, 63, 83-90. (27) ANSYS CFX-5.7 Solver Theory Manual; CFX, Ltd.: Oxfordshire, U.K., 2003.
sured downwind at several distances from the release point, using infrared gas sensors. Sufficient data for a simulation setup, accompanied by experimental records, for comparison with the computations, were observed for trial 8,14 which was considered further in the computational procedure. Proceeding in the simulation, the geometry of interest was built in the computer-aided design (CAD) interface of the code, the spatial dimensions of which are given in Table 1. The entire domain was subsequently subdivided into a total number of 779 000 volume elements, forming the mesh shown in Figure 1. As one can see, special care was taken for the refinement of the computational mesh onto and around the spill area, where the mixing of released gas with atmospheric air initializes. 3.2. Liquefied Hydrogen Spills. LH2 spills were conducted by NASA Langley Research Center at the White Sands test facility in 1980. The objective of the program was to study the phenomena associated with the rapid release of large amounts of LH2, with an emphasis on the generation and dispersion of the flammable hydrogen cloud. The trials included short releases of ∼5.1 m3 of LH2 in the open air. The LH2 was fed through a 0.15-m-diameter pipeline and spilled on a steel plate located in the center of a slightly elevated 9.1-m-diameter pond with compacted sand as the bottom.15 Concentration measurements were obtained using sample bottle clusters that were positioned downwind at predetermined time intervals. In addition, concentration histories were obtained by relating temperature measurements with a hydrogen concentration, assuming adiabatic mixing of the atmospheric air with hydrogen at its normal boiling point. Details of the data obtained during the experiments are presented for trial 6, which was of primary interest, because the short spill time (and, hence, higher spill rate), which would result in the formation of a largesized gas cloud. Therefore, trial 6 was judged suitable for further analysis. An aspect of the totally 532 000volume-element computational mesh constructed for the LH2 spill trial simulation is presented in Figure 2. Again, the mesh was selectively refined on the upper surface of the pond, along with its perimeter. The spatial dimensions of the domain are displayed in Table 1, whereas detailed meteorological data and other information for the experimental conditions of both LNG and LH2 trials are given in Table 2. 4. Simulation Setup The first step in the simulation preparation was to build an appropriate domain for each trial in the CAD system of the program and, afterward, create an optimal mesh that would allow a good approximating solution without prohibiting demand on computational time. The computational grids of the domains are partially presented in Figures 1 and 2. Near the gas inflow surface, as well as the sidewalls of the dike and the platform,
2538
Energy & Fuels, Vol. 19, No. 6, 2005
Sklavounos and Rigas
Figure 1. General aspect of liquefied natural gas (LNG) spill computational domain and spatial cross-section mesh display.
some inflation layers were constructed, in addition to the initial tetrahedral mesh. Inflation layers consist of prismatic elements that contribute to better modeling of physics in regions where sharp changes of fluid properties occur. In this case, the mixing in the vicinity of sources geometries causes rapid variations in air and gas stream velocities and temperatures, as well as in mixture composition. Data on the meshes construction parameters, as well as the dimensions of the corresponding domains, are given in Table 1. 4.1. Boundary Conditions. The equations relating to the flow may be numerically defined by the specification of certain conditions on the external boundaries of the domain. The following choices of boundary conditions are available in CFX code: (a) Inflow boundary, where the fluid or component(s) of a fluid mixture enters the domain of interest. It is strictly a one-way flow boundary condition and is defined by the temperature and the mass flow rate (or the velocity, or the static pressure) of the inlet. (b) Outflow boundary, where the fluid (or a fluid mixture) exits the domain. It is also a one-way flow boundary condition and it is defined by the mass flow rate (or the velocity, or the static pressure) of the outlet.
(c) Opening boundary, where the direction of the fluid flow is unknown a priori, for example, because of the generation of recirculation eddies. Fluid temperature and relative pressure should be known for the definition of this type of boundary. (d) Wall boundary, which can be fixed or slipping, in relation to the motion of the fluid. Walls are impermeable to fluid flow and allow the permeation of heat into and out of the domain. After a boundary is denoted as a wall, the relative velocity of the fluid on it is set to zero, according to the nonslip condition. The boundary conditions set for the present simulations are mentioned: the lower surface (ground) was defined as a wall with a temperature equal to the atmospheric temperature, to take into account ground heating effects; air inflow was set according to the data of Table 2 by denoting the experimentally recorded average wind velocity vectors (1.8, 0, 0) in the X-, Y-, and Z-directions. Data on ground surface roughness in the experiments were not available, so a uniform velocity profile was assumed, according to the average wind speed values of Table 2. The dikes were denoted as the inlets of gas into the domain, whereas the remaining outer surfaces were set as openings that
Fuel Gas Dispersion under Cryogenic Release
Energy & Fuels, Vol. 19, No. 6, 2005 2539
Figure 2. General aspect of liquefied hydrogen (LH2) spill computational domain and spatial cross-section mesh display. Table 2. Weather and Release Conditions during LNG Burro No. 8 and LH2 No. 6 Trials parameter
spill rate (m3/s)
liquid spill volume (m3)
spill duration (s)
average wind speed (m/s)
atmospheric temperature (K)
barometric pressure (kPa)
release temperature (K)
LNG trial 8 LH2 trial 6
0.266 0.135
28.4 5.11
107 38
1.80 1.75
306.1 288
94.1 78.6
111.5 20.48
represented the atmosphere surrounding the release area. Instead of using subroutines, CFX carries an automatic math expression editor, in which simple expressions can be made, so that can be involved later in the computations. Thus, in regard to the gas inlet, a properly adapted step-function was made to provide the mass inflow rate in the transient simulation, according to the realistic spill data displayed in Table 2:
m ˘ ) s˘ × step
( ) t1 - t0 tc2
(1)
where m ˘ is the mass inflow rate (given in units of kg/ s), s˘ the spilling rate (given in units of kg/s), t0 the spill initiation (t0 ) 0 s), t1 the spill duration (in seconds), and tc is a constant (tc ) 1 s). The release temperature was set equal to 20.48 and 111.5 K for LH2 and LNG release, respectively.
In problem physics definition, both hydrogen-air and methane-air fluid pairs were considered as real gases of variable composition mixture. Thermodynamic properties of each gas were individually set in a process that properties variability with temperature were taken into account through experimental data found in the literature.17,18 5. Results and Discussion The simulation procedure was subdivided into two parts. First, the problem was solved in steady state with zero mass inflow for the gas, to obtain an initial values file for the flow within the domain. After the problem had been converged in steady state with reasonable accuracy, the simulation was set-up again, but in transient form. In Figure 3, the flow velocity profile within the computational domain of the LH2 release simulation is
2540
Energy & Fuels, Vol. 19, No. 6, 2005
Sklavounos and Rigas
Figure 3. Flow velocity profile in the steady-state run for LH2 release simulation.
demonstrated as an example. Generally, the proximity between the downstream boundary and portions of the domain should be avoided where the flow direction varies or even turbulent eddies are formed. Otherwise, the solution process may become unstable and flow modeling may be insufficient. In the current case, the velocity profile in the outlet boundary (which has been intentionally situated at a relatively large distance away from the release source) shows that the flow pattern has been smoothed at the edge of the domain, in contrast with the complex flow profile near the release source. Obviously, the large changes in flow velocity and direction observed near the release source appear because of the presence of the platform. In both cases (steady-state and transient mode), the threshold value for the root-mean-square (RMS) of residuals set equal to 10-4. The RMS is a measure of the local imbalance of each conservative control volume equation and relates directly to whether the equations have been accurately solved. For the LH2 trial, in which the overall phenomenon had a comparatively short duration (∼60 s, including cloud passage from the points of interest), small time steps were given: 20 × 0.2 s + 112 × 0.5 s. Regarding the LNG trial, which had considerably larger duration (∼600 s), time steps were equal to 10 × 0.5 s + 595 × 1 s. Transient simulation times for the two trials were 15 and 54 h, respectively, on a personal computer (PC) platform with a 2400 MHz AMD processor and 1500 MB RAM.
Figure 4. Computational and experimental methane concentration histories at point (57, 0, 1) for the LNG spill trial.
5.1. Comparison with Experimental Records. Figures 4-7 show the experimental gas concentration curves versus the computational curves for the LNG and LH2 trials. The relative errors presented in Table 3 show that the predicted values lie within a factor of 2 of the observed values; therefore, it may be inferred that the values of the discrete magnitudes of cloud passage duration, maximal concentration, and total dose predicted via the CFX code are consistent with the data from the field tests. It may be concluded through the relative error estimations that the code underestimated the maximal concentrations and calculated shorter duration times for
Fuel Gas Dispersion under Cryogenic Release
Energy & Fuels, Vol. 19, No. 6, 2005 2541 Table 3. Estimation of Relative Errors between Computational and Experimental Maximal Concentrations, Arrival Times, and Total Doses at Measurement Points Relative Error (%) cloud passage maximal duration concentration
quantitya [1] LNG trial (57, 0, 1) [2] LNG trial (140, 0, 1) [3] LH2 trial (9.1, 0, 1) [4] LH2 trial (18.3, 0, 9.4)
-9.3 -14.5 -2.7 +1.7
-4.4 -23.2 -41 -40
+9.3 -9.6 -8.2 +14.3
7.1
27.2
10.4
mean absolute error a
total dose
Coordinates with reference to the center of the dikes.
Table 4. Evaluation of Statistical Performance Measures for LNG and LH2 Trials
Figure 5. Computational and experimental methane concentration histories at point (140, 0, 1) for the LNG spill trial.
statistical parameter
FB
MG
MRSE
NMSE
VG
ideal value [1] LNG trial (57, 0, 1) [2] LNG trial (140, 0, 1) [3] LH2 trial (9.1, 0, 1) [4] LH2 trial (18.3, 0, 9.4) combined data
0 -0.149 0.0215 -0.089 -0.0826 -0.075
1.0 0.852 1.019 0.851 0.732 0.8635
0 0.0497 0.032 0.1458 0.526 0.188
0 0.109 0.000073 0.078 0.44 0.157
1.0 1.052 1.033 1.177 2.559 1.45
eters that must be calculated are the fractional bias (FB), the geometric mean bias (MG), the geometric mean variance (VG), the mean relative square error (MRSE), and the normalized mean square error (NMSE):
Fractional Bias:
( )
FB ) 2
Geometric Mean Bias: Figure 6. Computational and experimental hydrogen concentration histories at point (9.1, 0, 1) for the LH2 spill trial.
X0 - Xp
(2)
X0 + Xp
MG ) exp(ln X0 - ln Xp) (3)
Geometric Mean Variance: VG ) exp[(ln X0 - ln Xp)2] (4) Mean Relative Square Error: MRSE ) 4
(
)
Xp - X0 Xp + X0
2
(5)
Normalized Mean Square Error: NMSE ) Figure 7. Computational and experimental hydrogen concentration histories at point (18.3, 0, 9.4) for the LH2 spill trial.
the cloud; however, minor differences are found among the computational and experimental total mass flows, namely, the dose (concentration × time) of each gas. Differences in arrival times might be due to instant changes in wind speed (wind gust) during the experiment execution that would affect the speed of entrained gas. The wind speed changes, as well as random changes in the blow direction, may affect the turbulent mixing of gas with air and, hence, cause rapid wide gas concentration variation.28 To obtain a statistically valid conclusion for the reliability of the computational results, as far as concentration histories are concerned, several statistical performance measures may be applied.29-31 The param(28) Nielsen, M. J. Loss Prev. Process Ind. 1991, 4, 29-34.
(Xp - X0)2 Xp × X0
(6)
where X0 is an observed (experimental) quantity, Xp is the corresponding predicted (modeled) quantity, and the overbar indicates an average. According to the calculated values of these parameters, which are displayed in Table 4, and the corresponding ideal values, it can be inferred that simulation results are observed to be in good agreement with the experimental results. A perfect model compared against observations would be placed at the MG ) 1 and VG ) 1 point on a MG vs VG diagram, such as that shown in Figure 8. A model that has no random scatter but suffers a mean bias would be placed somewhere along the parabolic curve, (29) Mohan, M.; Panwar, T. S.; Singh, M. P. Atmos. Environ. 1995, 29, 2075-2087. (30) Hanna, S. R.; Chang, J. C. Atmos. Environ. 2001, 35, 22312242. (31) Duijm, N. J.; Ott, S.; Nielsen, M. J. Loss Prev. Process Ind. 1996, 9, 323-338.
2542
Energy & Fuels, Vol. 19, No. 6, 2005
Sklavounos and Rigas
Figure 8. Plot of geometric mean (MG) versus geometric variance (VG) for the points of interest of LNG and LH2 trials comparing gas concentration histories. The points are numbered according to the sorting in Tables 3 and 4. Table 5. Coefficients of the Linear Regression Model C ) a - bT between Gas Concentration (C) and Temperature (T) Coefficient Value
Figure 9. Scatter plot for observed versus predicted concentration histories for the data from the LNG and LH2 trials.
which represents the minimum possible value of VG that corresponds to a particular MG. The dotted vertical lines on the same figure represent the MG values that correspond to a factor of 2 difference in the means. Obviously, code predictions fall between these lines, and, therefore, simulation results are within a factor of 2 of the observed values. Moreover, it could be said that, generally, the code slightly overpredicts gas concentrations (regardless of the underprediction of concentration maxima), because the majority of statistical test points are to the left of the centerline. In Figure 9, an overall graphical representation of experimentally observed versus computationally predicted cryogenic gas concentrations is depicted, taking into account the four points of interest of the two trials. The scatter plot shows clearly that the predicted versus observed concentration values are posed in the vicinity of the diagonal line, so it can be concluded that the simulations provided good approximations of gas concentration histories for both experimental trials. The data in Figures 8 and 9 were obtained by comparing gas concentration histories, whereas their averaging over all measurement positions is displayed in the last row of Table 4. In Figure 10, calculated concentration isopleths on the Z-X plane of symmetry is illustrated for hydrogen spill dispersion and compared with the corresponding experimental contours presented in the paper of Witcofski
trial
a
b
correlation coefficient, R2
LNG LH2
1.254 1.108
0.0042 0.0038
0.98 0.99
and Chirivella.15 The experimental contours were derived from correlating temperature data obtained during the trial with concentration and assuming adiabatic gasair mixing. Throughout the comparison, it is evident that simulation results are found in proximity with those of the trial, predicting the gas cloud dimensions well, in addition to the concentration distribution within it. The LFL is calculated to be extended 45 m away from the source (versus 38 m in the experiment), whereas the flammable cloud height is calculated to be equal to 22 m (versus 19 m in the experiment). However, the LFL distance at ground level has been overestimated (36 m versus 19 m in the experiment). 5.2. Concentration-Temperature Correlation. Elaboration of the results obtained through the simulations revealed a strong relationship between the computed gas concentration and temperature. The relationship was affirmed for a series of random points within each computational domain. Initially, concentration and temperature histories in those points were plotted giving graphs such that shown in Figure 11. As observed in this figure, concentration changes with time occur in the opposite way of the temperature changes. Thus, plotting gas concentration (C) against temperature (T) (Figure 12), a linear regression model of the form C ) a - bT is obtained. Model coefficient values for each simulation, as well as the regression correlation coefficients (R2), are displayed in Table 5, showing a statistically significant approximation. Nevertheless, a correlation between gas concentration and temperature should be expected, because the mixing of cold gas with air causes a decrease in temperature and, indeed, the higher the gas concentration, the higher the drop of temperature. If we keep in mind the assumption of adiabatic mixing between the dispersing gas and the entrained air (Section 3.2), then the gas
Fuel Gas Dispersion under Cryogenic Release
Energy & Fuels, Vol. 19, No. 6, 2005 2543
Figure 10. Comparison between experimental and computational gas concentration isopleths on the Z-X symmetry plane for the LH2 trial (time: 21 s).
Figure 12. Change of concentration versus temperature described by linear regression (LH2 spill, point 1).
concentration C is analogous to the total enthalpy change ∆H:32
where β is a coefficient, whereas HR and HT are the enthalpy budgets at atmospheric temperature and a time-varying total temperature at a random point of interest, respectively. Apparently, the first term on the right-hand side of eq 8 remains constant, while the second varies with temperature. Generally, when a cold substance is mixed with ambient air (at ambient temperature), neglecting surface heat transfer, then the negative buoyancy is also conserved, provided that
C ) β × ∆H S C ) β × HR - β × HT
(32) Nielsen, M. Spreading of Cold Dense Clouds. In The Handbook of Hazardous Materials Spills Technology; Fingas, M., Ed.; McGrawHill: New York, 2002; Chapter 18.
Figure 11. Concentration versus temperature histories for the LH2 spill (point 1), as predicted via the CFX code.
(7) (8)
2544
Energy & Fuels, Vol. 19, No. 6, 2005
differences in molar specific heat are negligible. This also provides the direct proportional relation between concentration and temperature in the cloud. As a result, it may be inferred that the code performed reasonable calculations. Furthermore, it arises that, in field-scale trials, plausible gas concentration predictions may be obtained in a more convenient way by relating them with temperature measurements. 6. Conclusions In this paper, a fuel gases safety-related issue was investigated. In particular, a computational approach to the behavior of hydrogen and natural gas during their atmospheric dispersion following cryogenic liquid spills was attempted. The simulation of two field-scale trials of liquefied hydrogen (LH2) and liquefied natural gas (LNG), performed via the computational fluid dynamics (CFD) code CFX 5.7, showed that CFD may be effectively used in risk analysis procedures in which liquefied fuel spill scenarios are involved. Simulation results, compared to the experimental records, were determined to be in reasonably good agreement, whereas the evaluation of statistical performance measures proved the simulation to be valid. The scatter plot of computational versus
Sklavounos and Rigas
experimental concentrations showed that they are found in proximity with the diagonal line, leading to the conclusion that the code yielded good approximations of gas concentration histories observed in both LNG and LH2 trials. The comparison of temperature and concentration histories at certain points of interest revealed a linear relationship between them. In addition, visualization of gas concentration distribution within the domains demonstrated that the cloud resulting from a cryogenic spill is driven by negative-buoyancy forces that remain at a low height for substantial distances away from the spill area and, indeed, in flammable concentrations. Consequently, even lighter-than-air gases (such as hydrogen and natural gas) should be considered to be heavier than air when they are released under cryogenic conditions. Acknowledgment. This work was supported by the Ministry of National Education and Religious Affairs (Community Support Framework 2000-2006), under the HERACLITUS Research Program. The authors wish also to thank Dr. Aubrey Thyer (Health and Safety Laboratory Consultant) for his considerate provision of refs 8 and 15. EF0500383
