- Email: [email protected]
Numerical simulation of LNG dispersion under two-phase release conditions
Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
Contents lists available at SciVerse ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
Numerical simulation of LNG dispersion under two-phase release conditions S.G. Giannissi a, b, *, A.G. Venetsanos a, N. Markatos b, J.G. Bartzis c a
Environmental Research Laboratory, National Centre for Scientific Research Demokritos, 15310 Aghia Paraskevi, Attikis, Greece Department of Process Analysis and Plant Design, School of Chemical Engineering, National Technical University of Athens, Heroon Polytechniou 9, 15780 Zografou, Greece c Department of Mechanical Engineering, University of Western Macedonia, Parko Agiou Dimitriou, West Macedonia, 50100 Kozani, Greece b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 9 January 2012 Received in revised form 10 August 2012 Accepted 21 November 2012
The use of LNG (liquefied natural gas) as fuel brings up issues regarding safety and acceptable risk. The potential hazards associated with an accidental LNG spill should be evaluated, and a useful tool in LNG safety assessment is computational fluid dynamics (CFD) simulation. In this paper, the ADREA-HF code has been applied to simulate LNG dispersion in open-obstructed environment based on Falcon Series Experiments. During these experiments LNG was released and dispersed over water surface. The spill area is confined with a billboard upwind of the water pond. FA1 trial was chosen to be simulated, because its release and weather conditions (high total spill volume and release rate, low wind speed) allow the gravitational force to influence the cold, dense vapor cloud and can be considered as a benchmark for LNG dispersion in fenced area. The source was modeled with two different approaches: as vapor pool and as two phase jet and the predicted methane concentration at sensors’ location was compared with the experimental one. It is verified that the source model affect to a great extent the LNG dispersion and the best case was the one modeling the source as two phase jet. However, the numerical results in the case of two phase jet source underestimate the methane concentration for most of the sensors. Finally, the paper discusses the effect of neglecting the 9.3 experimental wind direction, which leads to the symmetry assumption with respect to wind and therefore less computational costs. It was found that this effect is small in case of a jet source but large in the case of a pool source. Ó 2012 Elsevier Ltd. All rights reserved.
Keywords: LNG dispersion Two phase jet Falcon series Stable conditions ADREA-HF code CFD
1. Introduction In the last years the need to depend on other than petroleum energy resources, such as renewable energy sources or natural gas is growing rapidly. Especially, the use of natural gas as fuel has increased over the past few years. That along with the fact that natural gas is a flammable gas and for its storage and transportation usually liquefaction is necessary, led the research and industry community to evaluate the hazards of a potential LNG spill. Therefore, several experiments and computational studies related to LNG release have been conducted in the past. The LNG dispersion experiments that have been performed are the Burro Series (Burro Series Data Report, 1982), Coyote Series
* Corresponding author. Environmental Research Laboratory, National Centre for Scientific Research Demokritos, 15310 Aghia Paraskevi, Attikis, Greece. Tel.: þ30 210 6503416. E-mail addresses: [email protected] (S.G. Giannissi), venets@ ipta.demokritos.gr (A.G. Venetsanos), [email protected] (N. Markatos), bartzis@ uowm.gr (J.G. Bartzis). 0950-4230/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jlp.2012.11.010
(Coyote Series Data Report, Oct. 1983), Falcon Series (Falcon Series Data Report, June 1990), Maplin Sands tests (Colenbrander, Evans, & Puttock, 1984a), Esso tests (Feldbauer, Heigl, McQueen, Whipp, & May, 1972), Shell jettision tests (Kneebone & Prew, 1974), Avocet (Koopman, Bowman, & Ermak, 1979) and BFTF (Cormier, Qi, Yun, Zhang, & Mannan, 2009; Qi, Cormier, & Mannan, 2010). In all these experiments the release is above water surface, except for BFTF experiments, which also include the case of release above concrete. From the above experiments the effect of obstacles in an open environment was present only in the Falcon series and the BFTF experiments. In the present work the Falcon experiments were selected for simulation, because of the presence of obstacles. Specifically, FA1 trial was chosen, because of its release and weather conditions (high total spill volume and release rate, low wind speed) that permitted the gravitational force to influence the dense cold vapor cloud dispersion. Moreover, during that trial no rapid phase transition (RPT) explosions or fire were occurred, so the data of all sensors are available. Previous simulation work related to the Falcon experiments was performed by Chan (1992) using the Fem3A code (Chan, 1992). The
246
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
turbulence submodel that was used is the K-theory local equilibrium model. The value of K-diffusivity within the fence was estimated using some parameters that were adapted in order to obtain results consistent with experiment. The source model that was used in that study assumed that the spilled LNG evaporated as fast as it was spilled and over the entire water pond surface area. Although the experimentally measured wind direction formed an angle of 9.3 with respect to the fence alignment, Chan et al. assumed that the direction of the wind was parallel to the alignment of the fence. The predicted results were consistent with the experiment, however the model overall underpredicted the concentrations. Later on Gavelli, Bullister, and Kytomaa (2008 and 2009) simulated FA1 with the commercial CFD code Fluent. His first work (Gavelli et al., 2008) discusses the effect of vapor barriers on the area of ignitable vapor cloud, concluding that with barriers the hazardous area is reduced and impounded inside the fence. The source was modeled as an evaporating pool from the entire water pond surface. The model underpredicted the gas concentrations at the sensors downwind the fence and Gavelli et al. suggested that additional mixing of the cloud inside the fence would tend to improve the numerical prediction. Having noticed that the turbulence induced by the spill and by the evaporating LNG pool affects the LNG dispersion, in his latter study (Gavelli, Chernovsky, Bullister, & Kytomaa, 2009) he introduced a method to estimate the turbulent intensity by using the experimental video frames, and compared the calculated results with the experimental results. Four independent evaporating LNG pools (one beneath each spill spider arm) were assumed to be the LNG source. The pools were assumed to spread in time. The pools’ horizontal spreading velocities were estimated from the experimental video frames. The vapor from the pool was assumed vertically injected with velocity of approximately 0.07 m/s corresponding to an LNG vaporization flux of 0.127 kg/m2s. Three different turbulent intensity levels at the source were tested, in order to evaluate the sensitivity of the vapor cloud dispersion to the source-level turbulence. His investigation pointed out that turbulence created by the spreading and vaporization of an LNG pool onto water surface is important and a detailed definition of the LNG vapor source term including the turbulence intensity in the simulation process is critical in order to obtain acceptable results. However, the numerical results underestimated the gas concentration at the sensors downwind the fence in all three cases of different source-levels of turbulent intensity. In both studies the direction of the wind formed an angle of 9.3 with respect to the alignment of the fence. Focusing on the last observation that the spreading and vaporization of an LNG pool onto water surface creates an important amount of turbulence along with the fact that when modeling the source as an evaporating pool an assumption about the pool size and the evaporation rate per unit area should be made, the present study used a new approach of modeling the source. In the new approach the source is treated as a two phase jet, enabling simultaneous modeling of pool formation and spreading as well as the subsequent methane vapor dispersion. The computational results were compared with the experimental results and with the results that were obtained with modeling the source as a vapor pool. Additionally the effect of neglecting the 9.3 experimental wind direction, which leads to the symmetry assumption with respect to wind and therefore less computational costs was tested by performing calculations with and without wind direction. The numerical simulations were performed with the ADREA-HF three dimensional, time dependent, fully compressible finite volume CFD code (Bartzis, 1991; Würtz et al., 1996; Venetsanos, Papanikolaou, & Bartzis, 2010), that has been successfully used in the past for prediction of pollutants and hazardous gases dispersion in open and closed environments. In the last years main focus has been given to prediction of hydrogen dispersion for both cryogenic and compressed
releases (Baraldi, Venetsanos, Papanikolaou, Heitsch, & Dallas, 2009; Statharas, Venetsanos, Bartzis, Würtz, & Schmidtchen, 2000; Venetsanos, Baraldi, Adams, Heggem, & Wilkening, 2008; Venetsanos & Bartzis, 2007), concluding that ADREA-HF code is an appropriate CFD tool for risk assessment of hydrogen applications. 2. Description of Falcon series experiments The Falcon Series experiments were conducted in the Frenchman Flat area in Nevada in 1987 by Lawrence Livermore National Laboratory (LLNL) for the Department of Transportation (DOT) and the Gas Research Institute (GRI) as part of a joint government/industry study. A series of five LNG spills were performed in order to evaluate the effectiveness of walls and barriers on reducing the hazard distance associated with a large scale LNG release and to obtain extensive data for computer models validation. The LNG was spilled into a rectangular water pond (60 40 0.76 m) through a multi-exit spill “spider” for uniform distribution above water pond. The “spider” consisted of four arms and each arm was oriented 90 from adjacent arms. The exits had 0.11 m diameter and the LNG was released straight downward above the water surface (Fig.1). There was a circulation system in the pond, so that the water temperature would be stable and the vaporization rate could be considered nearly equal to the spill rate. The domain was a confined area with an 88 m long, 44 m wide and 8.7 m tall fence. A billboard 13.3 m tall and 17.1 m wide was located upwind of the water pond. During the Falcon series, a total of 77 gas concentration sensors were deployed. The sensors were located at 50, 150, and 250 m downwind of the fence and at 1, 5, 11, and 17 m height. The facility is depicted in Fig. 2. The figure is taken from the graphical user interface of ADREA-HF code named Edes. The rectangular blue area is the water pond, and the release points are where the four yellow dots indicate. The spheres show the sensors’ location in the experiment. Along with the gas sensors, several instruments were deployed, in order to measure the air temperature, pressure and humidity, wind speed and direction, turbulence intensity, and heat flux from the ground. The release conditions of Falcon series experiments varies from 8.7 to 30.3 m3/min spill rate, 20.6e66.4 m3 spill volume, and 78e138.8 s spill duration. The wind speed measured at 2 m height ranged from 1.7 to 5.2 m/s and the atmospheric stability conditions were neutral to stable. Table 1 shows analytically the experimental conditions of FA1 test. 3. Mathematical formulation The governing equations that describe the mean flow are the 3-D time-dependent conservation equations for mixture mass, momentum, enthalpy and component (methane) total mass fraction.
Fig. 1. The multi-exit spider configuration for the LNG release (Falcon Series Data Report, June 1990).
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
X
247
qi ¼ 1
(6)
qi ¼ qiv þ qil
(7)
Turbulence was modeled using the standard k-3 model (Launder & Spalding, 1974), in which buoyancy effects were included (Markatos & Pericleous, 1984; Venetsanos et al., 2010). The value of turbulent Schmidt and Prandtl numbers (s) was 0.72. The phase distribution of the binary mixture (air and methane) is performed assuming thermodynamic equilibrium. Methane in liquid phase is assumed to appear instantaneously when the mixture temperature falls below the mixture dew temperature. If this happens then the amount of liquid methane is calculated using Raoult’s law for ideal mixtures. It is important to take under consideration the heat exchange between the ground and the bulk flow. Therefore, the ground heat transfer was modeled by solving a transient one-dimensional temperature equation inside the ground (Statharas et al., 2000). Fig. 2. The experimental layout of the Falcon tests (ADREA-HF view). The spheres represent the sensors’ locations.
Mixture mass (continuity equation):
vr vrui þ ¼ 0 vt vxi
(1)
Mixture NaviereStokes (momentum equations)
vrui vrui uj vP v þ ¼ þ vxi vxj vt vxj
mef
vui vuj þ vxj vxi
!! þ rgi
mef vH s vxj
! þ
4.1. The 3D steady state wind problem
dP dt
(3)
Component total mass fraction:
vrqi vruj qi v þ ¼ vxj vt vxj
mef vqi s vxj
! (4)
The mixture density is related to component densities and mass fractions through:
1
r
¼
X qi
ri
¼
Xqiv
riv
þ
qil
The simulations were performed in two steps. First a 3D steady state wind field problem was solved to obtain the 3D-wind field over the obstructed terrain. The computational domain for this problem covered a wider area around the obstacles. Then the transient dispersion problem was solved in a domain more restricted and refined near the sources. The first problem provided initial and inflow boundary conditions for the second problem. Due to the different grids used, transfer of values from one problem to the next was performed using linear interpolation.
(2)
Mixture Enthalpy:
vrH vruj H v þ ¼ vt vxj vxj
4. Modeling strategy
(5)
ril
The domain dimensions for the 3D steady state wind field simulation were 500 500 50 m in the x, y, z directions respectively. The total number of cells was 132020. Inflow velocity and temperature profiles were obtained from a least squares fit of the available experimental data. The derived profiles and the experimental data are shown in Fig. 3. The figures also shows the wind velocity and temperature profiles generated by the MonineObukhov similarity theory, using the reported experimental data (L ¼ 4.963) and described by the following equations:
z u z * $ ln jm zo L k z T* z $ ln jh TðzÞ ¼ To þ zo L k uðzÞ ¼
Where for stable conditions (L 0):
jm Table 1 FA1 trial release and weather conditions and boundary layer data. Spill rate (m3/min) Spill duration (s) Ambient Pressure (atm) Ambient Temperature (K)/upper-lower Wind Speed (m/s) Atmospheric Stability conditions Friction velocity, u* (m/s) Dynamical temperature, T* (K) Surface temperature, To (K) MonineObukhov length (m) Roughness length (m)
28.7 138.8 0.9 307.25e305.95 1.7(@ 2 m) Stable 0.0605 0.0577 304.5 4.963 0.008
(8)
z L
¼ jh
z L
z ¼ 5$ L
(9)
It can be observed that the velocity and temperature profiles produced by the Similarity Theory differ from the respective experimental. This behavior can be attributed to the fact that FA1 trial was performed in the evening during summer, when the atmospheric conditions were not yet stabilized. The average measured wind direction was at an angle of 9.3 with respect of the alignment of the fence. Therefore, u- and vvelocity components were computed by multiplying the wind velocity by the cosine and by the sine of the angle respectively. Both components were set as input in the 3D steady state problem. Inflow profiles for the turbulent kinetic energy and its dissipation rate were generated using the following relationships:
248
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
(b)
(a)
9
315
8 313
7 311 Temperature(K)
Velocity (m/s)
6 5 4 3
309
307
2 experimental data least squares fit data
1
305
experimental data least squares fit data
similatiry theory
0
similarity theory
303
0
10
20
30
40
50
Height (m)
0
10
20
30
40
50
Height (m)
Fig. 3. Wind velocity (a) and ambient temperature (b) profiles determined from a least squares fit of the available experimental data (boxes) and compared against similarity theory.
1 =
k ¼
1
cm
2
u2* ;
3
¼
u3* kz
(10)
where cm ¼ 0:09 and k ¼ 0:41. 4.2. The 3D transient dispersion problem The computational domain was 376 400 50 m in the x, y, z directions respectively, see Fig. 4. The total number of cells was 605 744. The minimum horizontal cell size was 0.11 m near the sources and expanding away from them with an expansion factor of 1.12. The minimum cell size in the vertical direction was 0.4 m near the ground and expanding away from it with an expansion factor of 1.07. The abovementioned grid provided grid-independent results. This was verified by performing a separate grid dependency study which examined grids consisting of 269 008, 399 840, 605 744 and 621 792 cells. At the inlet boundaries inflow (Dirichlet) boundary conditions were used with values provided by the previous 3D steady state wind field problem. Due to the wind direction the boundaries at the
west and north domain were considered both inlet boundaries. The outlet boundaries are prescribed by applying zero gradient for all variables except methane mass fraction and temperature, for which either a zero gradient boundary condition was applied if outflow occurs or a given value boundary condition (equal to the initial value) if inflow occurs The top boundary is a free boundary. According to Luketa-Hanlin et al. (Luketa-Hanlin, Koopman, & Ermak, 2007) symmetry boundary condition along the top plane is more appropriate. Therefore, the normal component of velocity was set equal to zero and for all other variables zero gradient boundary condition was applied. At the bottom boundary standard rough wall functions were applied and a zero gradient for methane mass fraction. On the ground part of the bottom boundary the heat flux was obtained by solving the one-dimensional energy equation inside the ground. On the water part of the bottom boundary (the pond) a constant temperature was assumed equal to 301.55 K. The source was modeled with two different methods: as a constant area evaporating pool and as a two-phase jet. In the constant area evaporating pool formulation it is assumed that the liquid methane evaporates instantly as soon as it touches the ground. This assumption is valid since the release is above water. Water has large heat capacity, so its temperature can be
Fig. 4. The mesh display on plane y ¼ 0 and the bottom section and an aspect of the grid refinement close to the LNG release points (ADREA-HF view). The yellow arrows show the location of the release points (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
conserved almost constant during the spill process. Furthermore, there was a recirculation system in the pond, in order to keep the water’s temperature constant. The liquid methane pool was assumed to have during the entire release period, the horizontal dimensions of the water pond (2400 m2). This area was assumed to release saturated vapor methane (at 90 890 Pa, 110.25 K) upwards with a vertical velocity of 0.054 m/s. This velocity was calculated from the given mass spill rate (205.71 kg/s), the density of saturated vapor methane (1.59 kg/m3) and the assumed pool area. In the two-phase jet formulation it is assumed that methane was released vertically downwards under two-phase flow conditions exactly from the experimental release points from 4 jets with diameters 0.11 m (as in the experiment). The amount of liquid methane that flashed vaporized at the source was calculated by assuming isenthalpic expansion from a pressure of 450 000 Pa (before opening the valve) to 90 890 Pa (ambient pressure). The temperature upwind and downwind the exit was the saturated temperature at each pressure. The calculated boundary conditions on each jet area were vertical velocity 544.57 m/s, ambient pressure, temperature 110.25 K and void fraction (vapor volume/ total volume) 98.04%. The density of the two-phase methane was 9.937 kg/m3. The two phase jet formulation enables the prediction of both the methane vapor dispersion and the pool formation and spreading. If liquid methane mass fraction has been predicted in the boundary cells on the bottom domain (ground level) then the total horizontal area of these cells is assumed to be the liquid pool area. This procedure is followed for every time step, and thus the prediction of the spreading of the pool is possible.
249
The discretization of the differential equations is performed in ADREA-HF with the control volume approach. For the time integration and for discretization of convective terms the first order fully implicit scheme and the first order upwind scheme were used respectively. The time step was set automatically by the code starting from a minimum of 104 s with a restriction of maximum Courant number equal to 2. 5. Results and discussion 5.1. Qualitative comparison with the experiment Fig. 5 shows the predicted methane concentration (v/v) time histories for both types of source modeling compared against the experimental for three sensors at locations 50, 150 and 250 m downwind of the impoundment on the plane y ¼ 0 and at 1 m above the ground surface. It can be observed that for both source modeling formulations the predictions underestimate the maximum methane concentration with respect to the experiment. With the two phase jet formulation the maximum predicted concentration is in better agreement to the experimental at all sensors, except the one near the release point (at the location 50 m downwind of the fence). Another effect of the source modeling can be observed in the predicted arrival time of the vapor cloud. The arrival time predicted using the two phase jet formulation is shorter and in much better agreement to the experiment than with the vapor pool formulation. In the case of the pool the vapor cloud is impounded inside the fence for longer duration than in the case of two phase jet. The main reasons for this behavior are the high momentum and the
Fig. 5. Methane concentration (by volume) time histories compared against the experimental data for FA1 trial. The sensors are in different distances downwind the release point at plane y ¼ 0 and at 1 m elevation.
250
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
Fig. 6. The predicted LFL (5% v/v) iso-surface at 60 s in the case of vapor pool (left) in comparison with the respective predicted LFL iso-surface in the case of two phase jet (right).
Fig. 7. FA1 predicted vapor concentration contours (left) on the plane x ¼ 150 m at 140 s (top) and 160 s (bottom) compared with the experimental contours (right). Source is modeled as two phase jet.
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
251
evaporating area assumption made for the vapor pool formulation is not totally consistent. After the ceasing of the release (at 138.8 s) the pool area disappears almost instantaneously. 5.2. Quantitative comparison with the experiment For a quantitative comparison of the predicted results compared to the experimental statistical performance measures were applied. The statistical measures used in the present work are the fractional bias (FB), normalized mean square error (NMSE), geometric mean bias (MG) and geometric mean variance (VG).
FB ¼ 2
NMSE ¼
Fig. 8. The liquid pool area of methane versus time for the whole domain.
Table 2 Statistical performance measures for both types of source modeling. Trial
Ideal value FA1 vapor pool FA1 two phase jet
Statistical number FB
NMSE
MG
ln(VG)
0 0.62 0.60
0 2.02 1.15
1 0.002 0.65
0 196 1.0
turbulence associated/generated by the two-phase jet source. Both mechanisms contribute to the fast overfilling of the fence and the subsequent downwind dispersion of the cloud. Fig. 6 presents the predicted LFL (5% v/v) iso-surface at 60 s for both types of source modeling. The trapping of the LFL cloud inside the impoundment is clearly shown in the case of the vapor pool prediction. In Fig. 7, the predicted vapor concentration contours with the case of two phase jet for two different times (140 and 160 s) on the crosswind plane x ¼ 150 m are compared with the respective observed concentration contours. The overall agreement between model prediction and the experiment is fairly good. Fig. 8 represents the time history of the liquid pool area of methane for the whole domain. According to the figure the pool area gradually increases and reaches its highest value (61 m2) at approximately 136 s. This behavior indicates that the constant
Cp Co
!
Co þ Cp 2 Cp Co Cp Co
(11)
(12)
MG ¼ exp ln Cp ln Co
(13)
2 lnðVGÞ ¼ ln Cp ln Co
(14)
where C is for the time average methane concentration and X is the mean value of X, based on all available sensors. Subscripts (p) and (o) denote predicted and observed values respectively. The ideal value of FB, ln(VG) and NMSE is zero, while the ideal value of MG is 1. A negative sign in FB and values of MG lower than 1.0 indicate underestimation of the concentration compared to the experimental values. Table 2 displays the statistical performance measures for both types of source modeling. It can be observed that the case of source modeling as two phase jet gives statistical measures closer to the ideal values and with much less spreading. It can also be observed that both formulations generally underpredict the concentrations. One possible reason for the underestimation of the gas concentration is the fact that the effect of the spider arm structure from which the LNG was released was not taken into account. As it is mentioned in the study of Gavelli et al. (Gavelli et al., 2009) the spider arm hindered the vapor cloud to spread in its direction (upwind the fence). As a result, the cloud spreads more downwind the release. The underprediction of the gas concentration with the vapor pool formulation can mainly be attributed to the fact that the wind drifts the cloud to its direction. Therefore, at the sensors along x-axis and at the sensors in the north area of the computational domain very low gas concentration levels are predicted.
Fig. 9. Scatter plot for observed versus predicted maximum (left) and average (right) gas concentration for 18 sensors from FA1 (source modelling as two phase jet).
252
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
Fig. 10 presents the observed versus the predicted arrival time of the maximum gas concentration. The points represent the arrival time of maximum concentration for 18 sensors. It can be observed that the majority of the points are in the vicinity of the diagonal line and so the simulation is in good agreement with the experiment as far as the arrival time of maximum gas concentration is concerned. The points with smaller times correspond to the sensors nearest to the source and it can be concluded by the scatter graph that the prediction of the arrival time for those sensors is better than the prediction for the sensors far from the source. However, the model tends to predict the arrival time of maximum concentration approximately within a factor of 2. 5.3. Effect of wind direction on the vapor cloud dispersion
Fig. 10. Scatter graph of observed versus predicted arrival time (s) of maximum concentration for 18 sensors from FA1 (source modelling as two phase jet).
Fig. 9 depicts the scatter plot for predictions versus observations, using the maximum and average concentration of 18 sensors for the two phase jet case. According to these scatter plots the simulation underestimates both the maximum and average concentration approximately by a factor of 3.
In previous work (Chan et al. (Chan, 1992)), the average measured wind direction at an angle of 9.3 with respect to the alignment of the fence was neglected. The wind direction was assumed parallel to the fence and the x-z plane was the plane of symmetry, in order to decrease computational time. In the present work the effect of neglecting the wind direction on vapor cloud dispersion was examined. Fig. 11 and Fig. 12 show the predicted methane concentration (v/v) time histories with wind direction parallel to the fence compared against the experimental and the prediction with wind direction at an angle of 9.3 with respect to the alignment of the
Fig. 11. Methane concentration (by volume) time histories compared against the experimental data and the predicted methane concentration with the measured wind direction for the vapor pool formulation. The sensors are in different distances downwind the release point at symmetry axis (y ¼ 0) and at 1 m elevation.
Fig. 12. Methane concentration (by volume) time histories compared against the experimental data and the predicted methane concentration with the measured wind direction for the two phase jet formulation. The sensors are in different distances downwind the release point at symmetry axis (y ¼ 0) and at 1 m elevation.
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
253
Fig. 13. The predicted vapor concentration contours with the vapor pool formulation at 280 sec with wind direction shifted at an angle of 9.3 (left) and wind direction parallel to the x-axis (right).
Fig. 14. The predicted vapor concentration contours with the two phase jet formulation at 100 sec with wind direction shifted at an angle of 9.3 (left) and wind direction parallel to the x-axis (right).
fence for the vapor pool and two phase jet formulation respectively. The sensors are at locations downwind of the impoundment on the symmetry plane (y ¼ 0). The effect of the wind direction on vapor cloud dispersion is greater in the case of vapor pool than in the case of two phase jet. For the vapor pool formulation the prediction with wind direction parallel to the alignment of the fence seems in better agreement with the experiment. With the measured wind direction the predicted methane cloud disperses at the angle of wind direction (see Fig. 13), and consequently lower gas concentration levels reach the sensors along the x-axis. For the two phase jet formulation the effect of the wind direction on vapor cloud dispersion is less significant. The wind drifted the vapor cloud to its direction, but the spread is significant in both lateral sides (see Fig. 14), in contrast to the vapor pool formulation. The aforementioned discrepancy, with respect to the wind direction, between the two source modeling formulations can be attributed to the high jet source momentum. 6. Conclusions In this work, the ADREA-HF code was applied to simulate the FA1 experiment related to atmospheric dispersion of LNG in the presence of obstacles. Two different source modeling formulations were applied: a) a constant area vapor pool formulation without
any account of extra turbulence generated by the spill itself, emitting vapor methane vertically upwards and b) a two-phase jet formulation with 4 jets of diameter and location as in the experiment emitting two-phase methane with a void fraction of 98% vertically downwards. In the second case the turbulence generated due to the spill is predicted. The predictions were compared with the experimental concentrations both qualitatively and quantitatively. The quantitative assessment was based on statistical performance measures. It was found that for the case of the two phase jet the model is in general better agreement with the experiment but with underestimation of maximum concentrations. For the case of constant area vapor pool the model generally underestimates concentrations and overestimates the arrival times. One possible reason for the underestimation of the gas concentration is the fact that the effect of the spider arm structure hindering the vapor cloud to spread upwind the fence was not taken into account. While, the main reason for the underestimation of the gas concentration with the vapor pool formulation is that the cloud is drifted in the direction of the wind, and at the sensors along x-axis and in the north area of the computational domain the model predicts very low gas concentration levels. In the present work, a sensitivity study was also performed on the effect of wind direction (approximately 10 variation) on the vapor cloud dispersion for both source modeling formulations. The
254
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
effect of wind direction was found important for the vapor pool formulation and not significant for the jet formulation. The main reason for this discrepancy was attributed to the high two-phase jet momentum. Nomenclature xi ui P T gi qi t H
s k cm u*
k
T* To L z zo Greek 3
mef r,ri
Cartesian j co-ordinate (m) i component of velocity (ms1) pressure (Pa) temperature (K) gravity acceleration in the i-direction (ms2) mass fraction of i component (dimensionless) time (s) Enthalpy (Jkg1) turbulent Prandtl and Schmidt number (dimensionless) turbulent kinetic energy (m2s2) dimensionless constant friction velocity (ms1) von Karman constant (dimensionless) dynamical temperature (K) surface temperature (K) MonineObukhov length (m) height (m) roughness length (m)
Turbulent energy dissipation rate (m2s3) effective viscosity (kgm-1s1) mixture density, i-component density (kgm3)
Subscripts v vapor l liquid References Baraldi, D., Venetsanos, A. G., Papanikolaou, E., Heitsch, M., & Dallas, V. (2009). Numerical analysis of release, dispersion and combustion of liquid hydrogen in a mock-up hydrogen refuelling station. Journal of Loss Prevention in the Process Industries, 22, 303e315. Bartzis, J. C. (1991). ‘ADREA-HF: A three dimensional finite volume code for vapour cloud dispersion in complex terrain’. Report EUR 13580.
Burro Series Data Report. (1982). LLNL/NWC Report No.UCID-19075, v.1 2. Berkeley, CA: Lawrence Livermore National Laboratory. Chan, S. T. (1992). Numerical simulations of LNG vapor dispersion from a fenced storage area. Journal of Hazardous Materials, 30, 195e224. Colenbrander, G. W., Evans, A., & Puttock, J. S. (1984a). Spill tests of LNG and refrigerated liquid propane on the sea’ Maplin Sands 1980, Dispersion Data digest, Trial 27. Shell Research Ltd, Thornton Research Center, Report TNER.84.028, May 1984. Cormier, B. R., Qi, R., Yun, G. W., Zhang, Y., & Mannan, M. S. (2009). Application of computational fluid dynamics for LNG vapor dispersion modeling: a study of key parameters. Journal of Loss Prevention in the Process Industries, 22, 332e352. Coyote Series Data Report. (Oct. 1983). LLNL/NWC, UCID-19953, Vol. 1 2. Falcon Series Data Report; Gas research Institute, 1987 LNG barrier Verification field trials, GRI Report No. 89/0138, Chicago, IL, June 1990. Feldbauer, G. F., Heigl, J. J., McQueen, W., Whipp, R. H., & May, W. G. (1972). Spills of LNG on water: Vaporization and downwind drift of combustible mixtures. ESSO Research and Engg. Co, Report no EE61E-72. Gavelli, F., Bullister, E., & Kytomaa, H. (2008). Application of CFD(Fluent) to LNG spills into geometrically complex environments. Journal of Hazardous Materials, 159, 158e168. Gavelli, F., Chernovsky, M. K., Bullister, E., & Kytomaa, H. K. (2009). Quantification of source-level turbulence during LNG spills onto water. Journal of Loss Prevention in the Process Industries, 22, 809e819. Kneebone, A., & Prew, L. R. ‘Shipboard jettison test of LNG onto the sea’. In Proceedings on the 4th International Conference on LNG, Algiers, 1974; pp. 1e25. Koopman, R. P., Bowman, B. R., & Ermak, D. L. (1979). Data and calculation of dispersion on 5m3LNG spill tests. Lawrence Livermore Laboratory. Launder, B. E., & Spalding, D. B. (1974). The numerical computation of turbulent flow. Computer Methods in Applied Mechanics and Engineering, 3(2), 269e289. Luketa-Hanlin, A., Koopman, P. R., & Ermak, D. L. (2007). On the application of computational fluid dynamics codes for liquefied natural gas dispersion. Journal of Hazardous Materials, 140(3), 504e517. Markatos, N. C., & Pericleous, K. A. (1984). Laminar and turbulent natural convection in an enclosed cavity. International Journal of Heat and Mass Transfer, 27(5), 755e772. Qi, R., Cormier, B. R., & Mannan, M. S. (2010). Numerical simulations of LNG vapour dispersion in Brayton Fire Training Field tests with ANSYS CFX. Journal of Hazardous Materials, 183, 51e61. Statharas, J. C., Venetsanos, A. G., Bartzis, J. G., Würtz, J., & Schmidtchen, U. (2000). Analysis of data from spilling experiments performed with liquid hydrogen. Journal of Hazardous Materials A, 77(1e3), 57e75. Venetsanos, A. G., & Bartzis, J. G. (2007). CFD modeling of large-scale LH2 spills in open environment. International Journal Hydrogen Energy, 32, 2171e2177. Venetsanos, A. G., Baraldi, D., Adams, P., Heggem, P. S., & Wilkening, H. (2008). CFD modelling of hydrogen release, dispersion and combustion for automotive scenarios. Journal of Loss Prevention in the Process Industries, 21, 162e184. Venetsanos, A. G., Papanikolaou, E., & Bartzis, J. G. (2010). The ADREA-HF CFD code for consequence assessment of hydrogen applications. International Journal of Hydrogen Energy, 35, 3908e3918. Würtz, J., Bartzis, J. G., Venetsanos, A. G., Andronopoulos, S., Statharas, J., & Nijsing, R. (1996). A dense vapour dispersion code package for applications in the chemical and process industry. Journal of Hazardous Materials, 46, 273e284.
Contents lists available at SciVerse ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
Numerical simulation of LNG dispersion under two-phase release conditions S.G. Giannissi a, b, *, A.G. Venetsanos a, N. Markatos b, J.G. Bartzis c a
Environmental Research Laboratory, National Centre for Scientific Research Demokritos, 15310 Aghia Paraskevi, Attikis, Greece Department of Process Analysis and Plant Design, School of Chemical Engineering, National Technical University of Athens, Heroon Polytechniou 9, 15780 Zografou, Greece c Department of Mechanical Engineering, University of Western Macedonia, Parko Agiou Dimitriou, West Macedonia, 50100 Kozani, Greece b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 9 January 2012 Received in revised form 10 August 2012 Accepted 21 November 2012
The use of LNG (liquefied natural gas) as fuel brings up issues regarding safety and acceptable risk. The potential hazards associated with an accidental LNG spill should be evaluated, and a useful tool in LNG safety assessment is computational fluid dynamics (CFD) simulation. In this paper, the ADREA-HF code has been applied to simulate LNG dispersion in open-obstructed environment based on Falcon Series Experiments. During these experiments LNG was released and dispersed over water surface. The spill area is confined with a billboard upwind of the water pond. FA1 trial was chosen to be simulated, because its release and weather conditions (high total spill volume and release rate, low wind speed) allow the gravitational force to influence the cold, dense vapor cloud and can be considered as a benchmark for LNG dispersion in fenced area. The source was modeled with two different approaches: as vapor pool and as two phase jet and the predicted methane concentration at sensors’ location was compared with the experimental one. It is verified that the source model affect to a great extent the LNG dispersion and the best case was the one modeling the source as two phase jet. However, the numerical results in the case of two phase jet source underestimate the methane concentration for most of the sensors. Finally, the paper discusses the effect of neglecting the 9.3 experimental wind direction, which leads to the symmetry assumption with respect to wind and therefore less computational costs. It was found that this effect is small in case of a jet source but large in the case of a pool source. Ó 2012 Elsevier Ltd. All rights reserved.
Keywords: LNG dispersion Two phase jet Falcon series Stable conditions ADREA-HF code CFD
1. Introduction In the last years the need to depend on other than petroleum energy resources, such as renewable energy sources or natural gas is growing rapidly. Especially, the use of natural gas as fuel has increased over the past few years. That along with the fact that natural gas is a flammable gas and for its storage and transportation usually liquefaction is necessary, led the research and industry community to evaluate the hazards of a potential LNG spill. Therefore, several experiments and computational studies related to LNG release have been conducted in the past. The LNG dispersion experiments that have been performed are the Burro Series (Burro Series Data Report, 1982), Coyote Series
* Corresponding author. Environmental Research Laboratory, National Centre for Scientific Research Demokritos, 15310 Aghia Paraskevi, Attikis, Greece. Tel.: þ30 210 6503416. E-mail addresses: [email protected] (S.G. Giannissi), venets@ ipta.demokritos.gr (A.G. Venetsanos), [email protected] (N. Markatos), bartzis@ uowm.gr (J.G. Bartzis). 0950-4230/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jlp.2012.11.010
(Coyote Series Data Report, Oct. 1983), Falcon Series (Falcon Series Data Report, June 1990), Maplin Sands tests (Colenbrander, Evans, & Puttock, 1984a), Esso tests (Feldbauer, Heigl, McQueen, Whipp, & May, 1972), Shell jettision tests (Kneebone & Prew, 1974), Avocet (Koopman, Bowman, & Ermak, 1979) and BFTF (Cormier, Qi, Yun, Zhang, & Mannan, 2009; Qi, Cormier, & Mannan, 2010). In all these experiments the release is above water surface, except for BFTF experiments, which also include the case of release above concrete. From the above experiments the effect of obstacles in an open environment was present only in the Falcon series and the BFTF experiments. In the present work the Falcon experiments were selected for simulation, because of the presence of obstacles. Specifically, FA1 trial was chosen, because of its release and weather conditions (high total spill volume and release rate, low wind speed) that permitted the gravitational force to influence the dense cold vapor cloud dispersion. Moreover, during that trial no rapid phase transition (RPT) explosions or fire were occurred, so the data of all sensors are available. Previous simulation work related to the Falcon experiments was performed by Chan (1992) using the Fem3A code (Chan, 1992). The
246
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
turbulence submodel that was used is the K-theory local equilibrium model. The value of K-diffusivity within the fence was estimated using some parameters that were adapted in order to obtain results consistent with experiment. The source model that was used in that study assumed that the spilled LNG evaporated as fast as it was spilled and over the entire water pond surface area. Although the experimentally measured wind direction formed an angle of 9.3 with respect to the fence alignment, Chan et al. assumed that the direction of the wind was parallel to the alignment of the fence. The predicted results were consistent with the experiment, however the model overall underpredicted the concentrations. Later on Gavelli, Bullister, and Kytomaa (2008 and 2009) simulated FA1 with the commercial CFD code Fluent. His first work (Gavelli et al., 2008) discusses the effect of vapor barriers on the area of ignitable vapor cloud, concluding that with barriers the hazardous area is reduced and impounded inside the fence. The source was modeled as an evaporating pool from the entire water pond surface. The model underpredicted the gas concentrations at the sensors downwind the fence and Gavelli et al. suggested that additional mixing of the cloud inside the fence would tend to improve the numerical prediction. Having noticed that the turbulence induced by the spill and by the evaporating LNG pool affects the LNG dispersion, in his latter study (Gavelli, Chernovsky, Bullister, & Kytomaa, 2009) he introduced a method to estimate the turbulent intensity by using the experimental video frames, and compared the calculated results with the experimental results. Four independent evaporating LNG pools (one beneath each spill spider arm) were assumed to be the LNG source. The pools were assumed to spread in time. The pools’ horizontal spreading velocities were estimated from the experimental video frames. The vapor from the pool was assumed vertically injected with velocity of approximately 0.07 m/s corresponding to an LNG vaporization flux of 0.127 kg/m2s. Three different turbulent intensity levels at the source were tested, in order to evaluate the sensitivity of the vapor cloud dispersion to the source-level turbulence. His investigation pointed out that turbulence created by the spreading and vaporization of an LNG pool onto water surface is important and a detailed definition of the LNG vapor source term including the turbulence intensity in the simulation process is critical in order to obtain acceptable results. However, the numerical results underestimated the gas concentration at the sensors downwind the fence in all three cases of different source-levels of turbulent intensity. In both studies the direction of the wind formed an angle of 9.3 with respect to the alignment of the fence. Focusing on the last observation that the spreading and vaporization of an LNG pool onto water surface creates an important amount of turbulence along with the fact that when modeling the source as an evaporating pool an assumption about the pool size and the evaporation rate per unit area should be made, the present study used a new approach of modeling the source. In the new approach the source is treated as a two phase jet, enabling simultaneous modeling of pool formation and spreading as well as the subsequent methane vapor dispersion. The computational results were compared with the experimental results and with the results that were obtained with modeling the source as a vapor pool. Additionally the effect of neglecting the 9.3 experimental wind direction, which leads to the symmetry assumption with respect to wind and therefore less computational costs was tested by performing calculations with and without wind direction. The numerical simulations were performed with the ADREA-HF three dimensional, time dependent, fully compressible finite volume CFD code (Bartzis, 1991; Würtz et al., 1996; Venetsanos, Papanikolaou, & Bartzis, 2010), that has been successfully used in the past for prediction of pollutants and hazardous gases dispersion in open and closed environments. In the last years main focus has been given to prediction of hydrogen dispersion for both cryogenic and compressed
releases (Baraldi, Venetsanos, Papanikolaou, Heitsch, & Dallas, 2009; Statharas, Venetsanos, Bartzis, Würtz, & Schmidtchen, 2000; Venetsanos, Baraldi, Adams, Heggem, & Wilkening, 2008; Venetsanos & Bartzis, 2007), concluding that ADREA-HF code is an appropriate CFD tool for risk assessment of hydrogen applications. 2. Description of Falcon series experiments The Falcon Series experiments were conducted in the Frenchman Flat area in Nevada in 1987 by Lawrence Livermore National Laboratory (LLNL) for the Department of Transportation (DOT) and the Gas Research Institute (GRI) as part of a joint government/industry study. A series of five LNG spills were performed in order to evaluate the effectiveness of walls and barriers on reducing the hazard distance associated with a large scale LNG release and to obtain extensive data for computer models validation. The LNG was spilled into a rectangular water pond (60 40 0.76 m) through a multi-exit spill “spider” for uniform distribution above water pond. The “spider” consisted of four arms and each arm was oriented 90 from adjacent arms. The exits had 0.11 m diameter and the LNG was released straight downward above the water surface (Fig.1). There was a circulation system in the pond, so that the water temperature would be stable and the vaporization rate could be considered nearly equal to the spill rate. The domain was a confined area with an 88 m long, 44 m wide and 8.7 m tall fence. A billboard 13.3 m tall and 17.1 m wide was located upwind of the water pond. During the Falcon series, a total of 77 gas concentration sensors were deployed. The sensors were located at 50, 150, and 250 m downwind of the fence and at 1, 5, 11, and 17 m height. The facility is depicted in Fig. 2. The figure is taken from the graphical user interface of ADREA-HF code named Edes. The rectangular blue area is the water pond, and the release points are where the four yellow dots indicate. The spheres show the sensors’ location in the experiment. Along with the gas sensors, several instruments were deployed, in order to measure the air temperature, pressure and humidity, wind speed and direction, turbulence intensity, and heat flux from the ground. The release conditions of Falcon series experiments varies from 8.7 to 30.3 m3/min spill rate, 20.6e66.4 m3 spill volume, and 78e138.8 s spill duration. The wind speed measured at 2 m height ranged from 1.7 to 5.2 m/s and the atmospheric stability conditions were neutral to stable. Table 1 shows analytically the experimental conditions of FA1 test. 3. Mathematical formulation The governing equations that describe the mean flow are the 3-D time-dependent conservation equations for mixture mass, momentum, enthalpy and component (methane) total mass fraction.
Fig. 1. The multi-exit spider configuration for the LNG release (Falcon Series Data Report, June 1990).
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
X
247
qi ¼ 1
(6)
qi ¼ qiv þ qil
(7)
Turbulence was modeled using the standard k-3 model (Launder & Spalding, 1974), in which buoyancy effects were included (Markatos & Pericleous, 1984; Venetsanos et al., 2010). The value of turbulent Schmidt and Prandtl numbers (s) was 0.72. The phase distribution of the binary mixture (air and methane) is performed assuming thermodynamic equilibrium. Methane in liquid phase is assumed to appear instantaneously when the mixture temperature falls below the mixture dew temperature. If this happens then the amount of liquid methane is calculated using Raoult’s law for ideal mixtures. It is important to take under consideration the heat exchange between the ground and the bulk flow. Therefore, the ground heat transfer was modeled by solving a transient one-dimensional temperature equation inside the ground (Statharas et al., 2000). Fig. 2. The experimental layout of the Falcon tests (ADREA-HF view). The spheres represent the sensors’ locations.
Mixture mass (continuity equation):
vr vrui þ ¼ 0 vt vxi
(1)
Mixture NaviereStokes (momentum equations)
vrui vrui uj vP v þ ¼ þ vxi vxj vt vxj
mef
vui vuj þ vxj vxi
!! þ rgi
mef vH s vxj
! þ
4.1. The 3D steady state wind problem
dP dt
(3)
Component total mass fraction:
vrqi vruj qi v þ ¼ vxj vt vxj
mef vqi s vxj
! (4)
The mixture density is related to component densities and mass fractions through:
1
r
¼
X qi
ri
¼
Xqiv
riv
þ
qil
The simulations were performed in two steps. First a 3D steady state wind field problem was solved to obtain the 3D-wind field over the obstructed terrain. The computational domain for this problem covered a wider area around the obstacles. Then the transient dispersion problem was solved in a domain more restricted and refined near the sources. The first problem provided initial and inflow boundary conditions for the second problem. Due to the different grids used, transfer of values from one problem to the next was performed using linear interpolation.
(2)
Mixture Enthalpy:
vrH vruj H v þ ¼ vt vxj vxj
4. Modeling strategy
(5)
ril
The domain dimensions for the 3D steady state wind field simulation were 500 500 50 m in the x, y, z directions respectively. The total number of cells was 132020. Inflow velocity and temperature profiles were obtained from a least squares fit of the available experimental data. The derived profiles and the experimental data are shown in Fig. 3. The figures also shows the wind velocity and temperature profiles generated by the MonineObukhov similarity theory, using the reported experimental data (L ¼ 4.963) and described by the following equations:
z u z * $ ln jm zo L k z T* z $ ln jh TðzÞ ¼ To þ zo L k uðzÞ ¼
Where for stable conditions (L 0):
jm Table 1 FA1 trial release and weather conditions and boundary layer data. Spill rate (m3/min) Spill duration (s) Ambient Pressure (atm) Ambient Temperature (K)/upper-lower Wind Speed (m/s) Atmospheric Stability conditions Friction velocity, u* (m/s) Dynamical temperature, T* (K) Surface temperature, To (K) MonineObukhov length (m) Roughness length (m)
28.7 138.8 0.9 307.25e305.95 1.7(@ 2 m) Stable 0.0605 0.0577 304.5 4.963 0.008
(8)
z L
¼ jh
z L
z ¼ 5$ L
(9)
It can be observed that the velocity and temperature profiles produced by the Similarity Theory differ from the respective experimental. This behavior can be attributed to the fact that FA1 trial was performed in the evening during summer, when the atmospheric conditions were not yet stabilized. The average measured wind direction was at an angle of 9.3 with respect of the alignment of the fence. Therefore, u- and vvelocity components were computed by multiplying the wind velocity by the cosine and by the sine of the angle respectively. Both components were set as input in the 3D steady state problem. Inflow profiles for the turbulent kinetic energy and its dissipation rate were generated using the following relationships:
248
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
(b)
(a)
9
315
8 313
7 311 Temperature(K)
Velocity (m/s)
6 5 4 3
309
307
2 experimental data least squares fit data
1
305
experimental data least squares fit data
similatiry theory
0
similarity theory
303
0
10
20
30
40
50
Height (m)
0
10
20
30
40
50
Height (m)
Fig. 3. Wind velocity (a) and ambient temperature (b) profiles determined from a least squares fit of the available experimental data (boxes) and compared against similarity theory.
1 =
k ¼
1
cm
2
u2* ;
3
¼
u3* kz
(10)
where cm ¼ 0:09 and k ¼ 0:41. 4.2. The 3D transient dispersion problem The computational domain was 376 400 50 m in the x, y, z directions respectively, see Fig. 4. The total number of cells was 605 744. The minimum horizontal cell size was 0.11 m near the sources and expanding away from them with an expansion factor of 1.12. The minimum cell size in the vertical direction was 0.4 m near the ground and expanding away from it with an expansion factor of 1.07. The abovementioned grid provided grid-independent results. This was verified by performing a separate grid dependency study which examined grids consisting of 269 008, 399 840, 605 744 and 621 792 cells. At the inlet boundaries inflow (Dirichlet) boundary conditions were used with values provided by the previous 3D steady state wind field problem. Due to the wind direction the boundaries at the
west and north domain were considered both inlet boundaries. The outlet boundaries are prescribed by applying zero gradient for all variables except methane mass fraction and temperature, for which either a zero gradient boundary condition was applied if outflow occurs or a given value boundary condition (equal to the initial value) if inflow occurs The top boundary is a free boundary. According to Luketa-Hanlin et al. (Luketa-Hanlin, Koopman, & Ermak, 2007) symmetry boundary condition along the top plane is more appropriate. Therefore, the normal component of velocity was set equal to zero and for all other variables zero gradient boundary condition was applied. At the bottom boundary standard rough wall functions were applied and a zero gradient for methane mass fraction. On the ground part of the bottom boundary the heat flux was obtained by solving the one-dimensional energy equation inside the ground. On the water part of the bottom boundary (the pond) a constant temperature was assumed equal to 301.55 K. The source was modeled with two different methods: as a constant area evaporating pool and as a two-phase jet. In the constant area evaporating pool formulation it is assumed that the liquid methane evaporates instantly as soon as it touches the ground. This assumption is valid since the release is above water. Water has large heat capacity, so its temperature can be
Fig. 4. The mesh display on plane y ¼ 0 and the bottom section and an aspect of the grid refinement close to the LNG release points (ADREA-HF view). The yellow arrows show the location of the release points (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
conserved almost constant during the spill process. Furthermore, there was a recirculation system in the pond, in order to keep the water’s temperature constant. The liquid methane pool was assumed to have during the entire release period, the horizontal dimensions of the water pond (2400 m2). This area was assumed to release saturated vapor methane (at 90 890 Pa, 110.25 K) upwards with a vertical velocity of 0.054 m/s. This velocity was calculated from the given mass spill rate (205.71 kg/s), the density of saturated vapor methane (1.59 kg/m3) and the assumed pool area. In the two-phase jet formulation it is assumed that methane was released vertically downwards under two-phase flow conditions exactly from the experimental release points from 4 jets with diameters 0.11 m (as in the experiment). The amount of liquid methane that flashed vaporized at the source was calculated by assuming isenthalpic expansion from a pressure of 450 000 Pa (before opening the valve) to 90 890 Pa (ambient pressure). The temperature upwind and downwind the exit was the saturated temperature at each pressure. The calculated boundary conditions on each jet area were vertical velocity 544.57 m/s, ambient pressure, temperature 110.25 K and void fraction (vapor volume/ total volume) 98.04%. The density of the two-phase methane was 9.937 kg/m3. The two phase jet formulation enables the prediction of both the methane vapor dispersion and the pool formation and spreading. If liquid methane mass fraction has been predicted in the boundary cells on the bottom domain (ground level) then the total horizontal area of these cells is assumed to be the liquid pool area. This procedure is followed for every time step, and thus the prediction of the spreading of the pool is possible.
249
The discretization of the differential equations is performed in ADREA-HF with the control volume approach. For the time integration and for discretization of convective terms the first order fully implicit scheme and the first order upwind scheme were used respectively. The time step was set automatically by the code starting from a minimum of 104 s with a restriction of maximum Courant number equal to 2. 5. Results and discussion 5.1. Qualitative comparison with the experiment Fig. 5 shows the predicted methane concentration (v/v) time histories for both types of source modeling compared against the experimental for three sensors at locations 50, 150 and 250 m downwind of the impoundment on the plane y ¼ 0 and at 1 m above the ground surface. It can be observed that for both source modeling formulations the predictions underestimate the maximum methane concentration with respect to the experiment. With the two phase jet formulation the maximum predicted concentration is in better agreement to the experimental at all sensors, except the one near the release point (at the location 50 m downwind of the fence). Another effect of the source modeling can be observed in the predicted arrival time of the vapor cloud. The arrival time predicted using the two phase jet formulation is shorter and in much better agreement to the experiment than with the vapor pool formulation. In the case of the pool the vapor cloud is impounded inside the fence for longer duration than in the case of two phase jet. The main reasons for this behavior are the high momentum and the
Fig. 5. Methane concentration (by volume) time histories compared against the experimental data for FA1 trial. The sensors are in different distances downwind the release point at plane y ¼ 0 and at 1 m elevation.
250
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
Fig. 6. The predicted LFL (5% v/v) iso-surface at 60 s in the case of vapor pool (left) in comparison with the respective predicted LFL iso-surface in the case of two phase jet (right).
Fig. 7. FA1 predicted vapor concentration contours (left) on the plane x ¼ 150 m at 140 s (top) and 160 s (bottom) compared with the experimental contours (right). Source is modeled as two phase jet.
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
251
evaporating area assumption made for the vapor pool formulation is not totally consistent. After the ceasing of the release (at 138.8 s) the pool area disappears almost instantaneously. 5.2. Quantitative comparison with the experiment For a quantitative comparison of the predicted results compared to the experimental statistical performance measures were applied. The statistical measures used in the present work are the fractional bias (FB), normalized mean square error (NMSE), geometric mean bias (MG) and geometric mean variance (VG).
FB ¼ 2
NMSE ¼
Fig. 8. The liquid pool area of methane versus time for the whole domain.
Table 2 Statistical performance measures for both types of source modeling. Trial
Ideal value FA1 vapor pool FA1 two phase jet
Statistical number FB
NMSE
MG
ln(VG)
0 0.62 0.60
0 2.02 1.15
1 0.002 0.65
0 196 1.0
turbulence associated/generated by the two-phase jet source. Both mechanisms contribute to the fast overfilling of the fence and the subsequent downwind dispersion of the cloud. Fig. 6 presents the predicted LFL (5% v/v) iso-surface at 60 s for both types of source modeling. The trapping of the LFL cloud inside the impoundment is clearly shown in the case of the vapor pool prediction. In Fig. 7, the predicted vapor concentration contours with the case of two phase jet for two different times (140 and 160 s) on the crosswind plane x ¼ 150 m are compared with the respective observed concentration contours. The overall agreement between model prediction and the experiment is fairly good. Fig. 8 represents the time history of the liquid pool area of methane for the whole domain. According to the figure the pool area gradually increases and reaches its highest value (61 m2) at approximately 136 s. This behavior indicates that the constant
Cp Co
!
Co þ Cp 2 Cp Co Cp Co
(11)
(12)
MG ¼ exp ln Cp ln Co
(13)
2 lnðVGÞ ¼ ln Cp ln Co
(14)
where C is for the time average methane concentration and X is the mean value of X, based on all available sensors. Subscripts (p) and (o) denote predicted and observed values respectively. The ideal value of FB, ln(VG) and NMSE is zero, while the ideal value of MG is 1. A negative sign in FB and values of MG lower than 1.0 indicate underestimation of the concentration compared to the experimental values. Table 2 displays the statistical performance measures for both types of source modeling. It can be observed that the case of source modeling as two phase jet gives statistical measures closer to the ideal values and with much less spreading. It can also be observed that both formulations generally underpredict the concentrations. One possible reason for the underestimation of the gas concentration is the fact that the effect of the spider arm structure from which the LNG was released was not taken into account. As it is mentioned in the study of Gavelli et al. (Gavelli et al., 2009) the spider arm hindered the vapor cloud to spread in its direction (upwind the fence). As a result, the cloud spreads more downwind the release. The underprediction of the gas concentration with the vapor pool formulation can mainly be attributed to the fact that the wind drifts the cloud to its direction. Therefore, at the sensors along x-axis and at the sensors in the north area of the computational domain very low gas concentration levels are predicted.
Fig. 9. Scatter plot for observed versus predicted maximum (left) and average (right) gas concentration for 18 sensors from FA1 (source modelling as two phase jet).
252
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
Fig. 10 presents the observed versus the predicted arrival time of the maximum gas concentration. The points represent the arrival time of maximum concentration for 18 sensors. It can be observed that the majority of the points are in the vicinity of the diagonal line and so the simulation is in good agreement with the experiment as far as the arrival time of maximum gas concentration is concerned. The points with smaller times correspond to the sensors nearest to the source and it can be concluded by the scatter graph that the prediction of the arrival time for those sensors is better than the prediction for the sensors far from the source. However, the model tends to predict the arrival time of maximum concentration approximately within a factor of 2. 5.3. Effect of wind direction on the vapor cloud dispersion
Fig. 10. Scatter graph of observed versus predicted arrival time (s) of maximum concentration for 18 sensors from FA1 (source modelling as two phase jet).
Fig. 9 depicts the scatter plot for predictions versus observations, using the maximum and average concentration of 18 sensors for the two phase jet case. According to these scatter plots the simulation underestimates both the maximum and average concentration approximately by a factor of 3.
In previous work (Chan et al. (Chan, 1992)), the average measured wind direction at an angle of 9.3 with respect to the alignment of the fence was neglected. The wind direction was assumed parallel to the fence and the x-z plane was the plane of symmetry, in order to decrease computational time. In the present work the effect of neglecting the wind direction on vapor cloud dispersion was examined. Fig. 11 and Fig. 12 show the predicted methane concentration (v/v) time histories with wind direction parallel to the fence compared against the experimental and the prediction with wind direction at an angle of 9.3 with respect to the alignment of the
Fig. 11. Methane concentration (by volume) time histories compared against the experimental data and the predicted methane concentration with the measured wind direction for the vapor pool formulation. The sensors are in different distances downwind the release point at symmetry axis (y ¼ 0) and at 1 m elevation.
Fig. 12. Methane concentration (by volume) time histories compared against the experimental data and the predicted methane concentration with the measured wind direction for the two phase jet formulation. The sensors are in different distances downwind the release point at symmetry axis (y ¼ 0) and at 1 m elevation.
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
253
Fig. 13. The predicted vapor concentration contours with the vapor pool formulation at 280 sec with wind direction shifted at an angle of 9.3 (left) and wind direction parallel to the x-axis (right).
Fig. 14. The predicted vapor concentration contours with the two phase jet formulation at 100 sec with wind direction shifted at an angle of 9.3 (left) and wind direction parallel to the x-axis (right).
fence for the vapor pool and two phase jet formulation respectively. The sensors are at locations downwind of the impoundment on the symmetry plane (y ¼ 0). The effect of the wind direction on vapor cloud dispersion is greater in the case of vapor pool than in the case of two phase jet. For the vapor pool formulation the prediction with wind direction parallel to the alignment of the fence seems in better agreement with the experiment. With the measured wind direction the predicted methane cloud disperses at the angle of wind direction (see Fig. 13), and consequently lower gas concentration levels reach the sensors along the x-axis. For the two phase jet formulation the effect of the wind direction on vapor cloud dispersion is less significant. The wind drifted the vapor cloud to its direction, but the spread is significant in both lateral sides (see Fig. 14), in contrast to the vapor pool formulation. The aforementioned discrepancy, with respect to the wind direction, between the two source modeling formulations can be attributed to the high jet source momentum. 6. Conclusions In this work, the ADREA-HF code was applied to simulate the FA1 experiment related to atmospheric dispersion of LNG in the presence of obstacles. Two different source modeling formulations were applied: a) a constant area vapor pool formulation without
any account of extra turbulence generated by the spill itself, emitting vapor methane vertically upwards and b) a two-phase jet formulation with 4 jets of diameter and location as in the experiment emitting two-phase methane with a void fraction of 98% vertically downwards. In the second case the turbulence generated due to the spill is predicted. The predictions were compared with the experimental concentrations both qualitatively and quantitatively. The quantitative assessment was based on statistical performance measures. It was found that for the case of the two phase jet the model is in general better agreement with the experiment but with underestimation of maximum concentrations. For the case of constant area vapor pool the model generally underestimates concentrations and overestimates the arrival times. One possible reason for the underestimation of the gas concentration is the fact that the effect of the spider arm structure hindering the vapor cloud to spread upwind the fence was not taken into account. While, the main reason for the underestimation of the gas concentration with the vapor pool formulation is that the cloud is drifted in the direction of the wind, and at the sensors along x-axis and in the north area of the computational domain the model predicts very low gas concentration levels. In the present work, a sensitivity study was also performed on the effect of wind direction (approximately 10 variation) on the vapor cloud dispersion for both source modeling formulations. The
254
S.G. Giannissi et al. / Journal of Loss Prevention in the Process Industries 26 (2013) 245e254
effect of wind direction was found important for the vapor pool formulation and not significant for the jet formulation. The main reason for this discrepancy was attributed to the high two-phase jet momentum. Nomenclature xi ui P T gi qi t H
s k cm u*
k
T* To L z zo Greek 3
mef r,ri
Cartesian j co-ordinate (m) i component of velocity (ms1) pressure (Pa) temperature (K) gravity acceleration in the i-direction (ms2) mass fraction of i component (dimensionless) time (s) Enthalpy (Jkg1) turbulent Prandtl and Schmidt number (dimensionless) turbulent kinetic energy (m2s2) dimensionless constant friction velocity (ms1) von Karman constant (dimensionless) dynamical temperature (K) surface temperature (K) MonineObukhov length (m) height (m) roughness length (m)
Turbulent energy dissipation rate (m2s3) effective viscosity (kgm-1s1) mixture density, i-component density (kgm3)
Subscripts v vapor l liquid References Baraldi, D., Venetsanos, A. G., Papanikolaou, E., Heitsch, M., & Dallas, V. (2009). Numerical analysis of release, dispersion and combustion of liquid hydrogen in a mock-up hydrogen refuelling station. Journal of Loss Prevention in the Process Industries, 22, 303e315. Bartzis, J. C. (1991). ‘ADREA-HF: A three dimensional finite volume code for vapour cloud dispersion in complex terrain’. Report EUR 13580.
Burro Series Data Report. (1982). LLNL/NWC Report No.UCID-19075, v.1 2. Berkeley, CA: Lawrence Livermore National Laboratory. Chan, S. T. (1992). Numerical simulations of LNG vapor dispersion from a fenced storage area. Journal of Hazardous Materials, 30, 195e224. Colenbrander, G. W., Evans, A., & Puttock, J. S. (1984a). Spill tests of LNG and refrigerated liquid propane on the sea’ Maplin Sands 1980, Dispersion Data digest, Trial 27. Shell Research Ltd, Thornton Research Center, Report TNER.84.028, May 1984. Cormier, B. R., Qi, R., Yun, G. W., Zhang, Y., & Mannan, M. S. (2009). Application of computational fluid dynamics for LNG vapor dispersion modeling: a study of key parameters. Journal of Loss Prevention in the Process Industries, 22, 332e352. Coyote Series Data Report. (Oct. 1983). LLNL/NWC, UCID-19953, Vol. 1 2. Falcon Series Data Report; Gas research Institute, 1987 LNG barrier Verification field trials, GRI Report No. 89/0138, Chicago, IL, June 1990. Feldbauer, G. F., Heigl, J. J., McQueen, W., Whipp, R. H., & May, W. G. (1972). Spills of LNG on water: Vaporization and downwind drift of combustible mixtures. ESSO Research and Engg. Co, Report no EE61E-72. Gavelli, F., Bullister, E., & Kytomaa, H. (2008). Application of CFD(Fluent) to LNG spills into geometrically complex environments. Journal of Hazardous Materials, 159, 158e168. Gavelli, F., Chernovsky, M. K., Bullister, E., & Kytomaa, H. K. (2009). Quantification of source-level turbulence during LNG spills onto water. Journal of Loss Prevention in the Process Industries, 22, 809e819. Kneebone, A., & Prew, L. R. ‘Shipboard jettison test of LNG onto the sea’. In Proceedings on the 4th International Conference on LNG, Algiers, 1974; pp. 1e25. Koopman, R. P., Bowman, B. R., & Ermak, D. L. (1979). Data and calculation of dispersion on 5m3LNG spill tests. Lawrence Livermore Laboratory. Launder, B. E., & Spalding, D. B. (1974). The numerical computation of turbulent flow. Computer Methods in Applied Mechanics and Engineering, 3(2), 269e289. Luketa-Hanlin, A., Koopman, P. R., & Ermak, D. L. (2007). On the application of computational fluid dynamics codes for liquefied natural gas dispersion. Journal of Hazardous Materials, 140(3), 504e517. Markatos, N. C., & Pericleous, K. A. (1984). Laminar and turbulent natural convection in an enclosed cavity. International Journal of Heat and Mass Transfer, 27(5), 755e772. Qi, R., Cormier, B. R., & Mannan, M. S. (2010). Numerical simulations of LNG vapour dispersion in Brayton Fire Training Field tests with ANSYS CFX. Journal of Hazardous Materials, 183, 51e61. Statharas, J. C., Venetsanos, A. G., Bartzis, J. G., Würtz, J., & Schmidtchen, U. (2000). Analysis of data from spilling experiments performed with liquid hydrogen. Journal of Hazardous Materials A, 77(1e3), 57e75. Venetsanos, A. G., & Bartzis, J. G. (2007). CFD modeling of large-scale LH2 spills in open environment. International Journal Hydrogen Energy, 32, 2171e2177. Venetsanos, A. G., Baraldi, D., Adams, P., Heggem, P. S., & Wilkening, H. (2008). CFD modelling of hydrogen release, dispersion and combustion for automotive scenarios. Journal of Loss Prevention in the Process Industries, 21, 162e184. Venetsanos, A. G., Papanikolaou, E., & Bartzis, J. G. (2010). The ADREA-HF CFD code for consequence assessment of hydrogen applications. International Journal of Hydrogen Energy, 35, 3908e3918. Würtz, J., Bartzis, J. G., Venetsanos, A. G., Andronopoulos, S., Statharas, J., & Nijsing, R. (1996). A dense vapour dispersion code package for applications in the chemical and process industry. Journal of Hazardous Materials, 46, 273e284.