Application of computational fluid dynamics for LNG vapor dispersion modeling: A study of key parameters

Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Contents lists available at ScienceDirect

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

Application of computational fluid dynamics for LNG vapor dispersion modeling: A study of key parameters Benjamin R. Cormier, Ruifeng Qi, GeunWoong Yun, Yingchun Zhang, M. Sam Mannan* Mary Kay O’Connor Process Safety Center, Artie McFerrin Department of Chemical Engineering, Texas A&M University System, College Station, TX 77843-3122, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 26 May 2008 Received in revised form 25 September 2008 Accepted 9 December 2008

The increased demand for Liquefied Natural Gas (LNG) has led to the construction of several new LNG terminals in the United States (US) and around the world. To ensure the safety of the public, consequence modeling is used to estimate exclusion distances. For LNG industry, the purpose of identifying these exclusion distances is to protect the public from being reached by flammable vapors during a release and they are determined by one-half of the Lower Flammability Limit (half LFL, 2.5% v/v). Since LNG vapors are heavier-than-air when released into atmosphere, it goes through several stages, which are respectively characterized as negative, neutral, and positively buoyant as it dilutes. To address this complex phenomenon, several simple models were developed and tested against large scale experimental data for the past three decades. This paper was derived from the development of design and safety specifications for LNG facilities based on experimental and theoretical research at Mary Kay O’Connor Process Safety Center (MKOPSC). Medium-scale LNG tests were performed at the Brayton Fire Training Field (BFTF), Texas A&M University to provide data for this specific research. Computational fluid dynamics (CFD) was used to perform consequence modeling for LNG release. The CFD code showed good agreement with the data collected during the November 2007 test performed at BFTF. This paper showed the simulation setup and the comparison with data collected for two scenarios: release on water and on dry concrete. Once the model was tuned against experimental data, it was used in a sensitivity analysis on parameters to assess the effects on the LFL distance and the concentration levels. Furthermore, three turbulence models were compared. The source term was composed of turbulence intensity at the source, LNG pool geometry, mass evaporation rate, and LNG pool area. The vapor dispersion parameters were wind velocity, sensible heat flux, and obstacles effects. It was concluded that at low wind velocity, the source term parameters strongly influenced the LFL distance and the concentration level. On the other hand, at high wind velocity, the source term parameter had a slight effect on the LFL distance and the concentration levels. Ó 2008 Elsevier Ltd. All rights reserved.

Keywords: CFD modeling LNG Vapor dispersion Source term Parametric analysis

1. Introduction There is a substantial amount of scientific literature on consequence modeling for gas heavier-than-air (or dense gas). Simple models were developed and were compared to experimental data to obtain a certain level of confidence in the prediction of exclusion distances. Today there is a need to assess complex scenarios (including complex geometry and obstacles). The applicability of the simple models is limited in complex scenarios. Several groups

* Corresponding author. Tel.: þ1 979 862 3985. E-mail address: [email protected] (M. Sam Mannan). 0950-4230/$ – see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jlp.2008.12.004

of research have been focusing on using CFD code, which seems to be an ideal candidate for the current needs of evaluating for complex scenarios. Even though the open literature has been growing on the use of CFD as consequence modeling, there are little details on how to address complex scenarios and how to verify the calculation involved. CFD consequence modeling is mainly used to assess a specific scenario which makes the task of developing standard guidelines difficult. This paper intended to show the capability of CFD modeling as applied to LNG vapor dispersion modeling through complex scenarios. First, the CFD model was calibrated against the experimental data collected at BFTF. Then the effects of key parameters involved in the setup were evaluated. The study of the

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

key parameters will give any engineer an understanding and a way to estimate the uncertainties associated with the simulation results. CFD code is capable of handling obstacles in a three-dimensional environment, which presents itself as an ideal candidate for consequence analysis tools for complex scenarios. It is not only capable of assessing the effectiveness of passive barriers such as dike (that may have an effect on the exclusion zones), but also it may be used as a design and plant layout tool. The setup of initial and boundary conditions is critical for the accuracy of prediction. In the open literature, there are at least four CFD codes that may be used for LNG vapor dispersion modeling where details in the setting up may be found: FEM3 CFD (Luketa-Hanlin, Koopman, & Ermak, 2007), FLACS (Dharmavaram, Hanna, & Hansen, 2005; Hanna, Hansen, & Dharmavaram, 2004), FLUENT (Gavelli, Bullister, & Kytomaa, 2006), and CFX (Sklavounos & Rigas, 2004). CFX and FLUENT are general-purpose CFD codes. Generic CFD codes are not calibrated for a specific problem like FEM3 CFD, or FLACS which were specifically developed for consequence modeling. The use of a generic CFD code was more appropriate for this type of work than specific CFD codes because of the flexibility of setting up. The model was based only on physics. For the vapor dispersion model, CFX uses the Reynolds Averaged Navier–Stokes (RANS) equations and is also based on the finite volume method for the conversion of partial differential equations and auxiliary conditions into a discrete system of equations. CFX has shown good results in the cryogenic release in atmosphere from Sklavounos and Rigas (2005). The authors verified the results of their CFX model against LNG and hydrogen trials. Even through the release trials were performed over flat terrain, the results were in agreement with the historical data. The experimental setup and data analysis details were presented in a separate paper. This paper only uses part of the analyzed data to setup and calibrate the CFD model. Part of the data collected during the November 2007 test was used to improve the source term modeling.

333

behind vapor dispersion is described in Fig. 1. Both box/simple models and CFD modeling for vapor dispersion have the same type of inputs. The level of details that may be captured by CFD modeling is far superior to box models. In Fig. 1, the algorithm may be divided in three sections: (i) problem definition, (ii) solver, and (iii) post processor. Typically, solver and post processor are not a variable in the study of the key parameters. Using a commercial CFD code, which has proven viable in other applications, reduces the uncertainties associated with the computational code calculations. The current best practice to assess the uncertainties of a model is to run a sensitivity analysis on the key parameters. As described in Fig. 1, the parameters are defined at the problem definition. The problem definition is composed of four sets of parameters: (i) domain definition, (ii) source term, (iii) atmospheric conditions, and (iv) ground effects. 2.1. Domain definition CFD modeling requires the determination of the domain, which is then divided into finite elements with mesh technique. In each element, the governing and turbulence equations are solved simultaneously with respect to a set of reference parameters, initial and boundary conditions. 2.1.1. Domain size and mesh The domain represents the volume of interest. The size may be estimated either from a coarse grid or by using box models. Then the domain is refined through a process to determine the numbers of elements required to obtain a grid independent solution. This practice is called grid independence test (Olvera & Choudhuri, 2006) or numerical experiments, i.e., repeating the calculation on a series of successively refined grids. For sufficiently small sizes, the rate of convergence is governed by the order of principal truncation error component (Ferziger & Peric, 1999).

2. Vapor dispersion modeling The approach for this work was to study in detail the parameters involved in vapor dispersion, especially for LNG. The algorithm

2.1.2. Governing equations The governing equations are solved simultaneously by using RANS. The continuity equations are:

Fig. 1. Algorithm for the vapor dispersion.

334

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

vr þ V$ðrUÞ ¼ 0 vt

(1)

and the momentum equation is:

   T vðrUÞ þ V,ðrU5UÞ  V$ meff VU ¼ Vp0 þ V$ meff VU þB vt (2) where

meff ¼ m þ mt

(3)

2 p0 ¼ p þ rk 3

(4)

meff ¼ the effective viscosity accounting for turbulence mt ¼ the turbulence viscosity

vðrhtot Þ vp  þ V$ðrUhtot Þ ¼ V$ðlVTÞ þ V$ðU$sÞ þ U$SM þ SE vt vt (5) SE ¼ the energy source T ¼ the temperature s ¼ viscous stress SM ¼ the momentum source. The energy and momentum source were set to zero since the source is introduced by one of the boundary conditions. htot is the total enthalpy, related to the static enthalpy and h is the function of temperature and pressure by:

1 htot ¼ h þ U 2 2

(6)

2.1.3. Turbulence model The two-equation turbulence model k–3, based on RANS equations, seems to overestimate the concentration for heavier-than-air gas modeling within the low atmospheric boundary layer (Sklavounos & Rigas, 2004). A stress model such as Reynolds Stress model may be used to increase accuracy as suggested by Gavelli et al. (2008). The choice of the turbulence depends on the accuracy wanted and the time constraints. A more accurate turbulence model may increase the running time. The two-equation turbulence model k–3 was chosen because of running time benefit. The values of k and 3 come directly from the differential transport equations for the turbulence kinetic energy (TKE) and turbulence dissipation rate (ANSYS, 2007):

   m vðrkÞ þ V$ðrUkÞ ¼ V$ m þ t Vk þ Pk  r3 sk vt

mþ





mt 3 V3 þ ðC31 Pk  C32 r3Þ s3 k

(7)

The k–3 two-equation model uses the gradient diffusion hypothesis to relate the Reynolds stresses to the mean velocity gradients and the turbulent viscosity. The turbulence viscosity is linked to TKE and dissipation using the relation:

mt ¼ Cm r

k2

(11)

3

 2  Pk ¼ mt VU$ VU þ VU T  V$Uð3mt V$U þ rkÞ þ Pkb 3

The k–3 turbulence model presented in Equations (7) and (8) has default parameters for Cm, C31, and C32 as (0.09, 1.44, and 1.92). Further study showed how to tune these parameters with respect to atmospheric data. Alinot and Masson (2005) recommend using (0.033, 1.17, and 1.92) instead of the default parameters for simulations less than one hundred meters above the ground. These parameters differ from the default setup because they were proven to be more accurate for the atmospheric boundary layer under various thermal stratifications (stable, neutral, and unstable). 2.1.4. Domain reference setup The reference parameters are defined by the acceleration due to gravity, reference pressure, and reference density. A reference density value was set at 1.225 kg/m3 at ground level (Olvera & Choudhuri, 2006). A high precision of reference density is required given the importance of significantly small density variations. The buoyant effect is driven by the gravity and the density of the fluid. Mixing cold gas with ambient air contains important thermal effects that need to be captured by the model. For LNG vapor dispersion, air and methane are usually modeled as gas phase, and their physical properties must be temperature dependent. Methane gas density must be around 1.78 kg/m3 at 110 K to 0.45 kg/m3 at 300 K. The humidity in the air may be added either using vapor water or by changing the thermodynamics properties of air to account for the humidity effect in the energy balance. The assumption of using methane gas instead of LNG composition (e.g., methane and trace heavier hydrocarbons) seems reasonable since methane will vaporize first. CFX uses the Boussinesq model (ANSYS, 2007) to predict the turbulence inside the cloud as function of the thermal expansivity, set to 0.00366 K1:



r  rref ¼ rref b T  Tref



(12)

where

(8)

where

with

C31, C32 ¼ turbulence model constant sk ¼ turbulence model constant for the k equation s3 ¼ turbulence model constant for the 3 equation sr ¼ 0.9 for Boussinesq buoyancy.

Cm ¼ viscosity turbulence constant respect k ¼ turbulence kinetic energy 3 ¼ turbulence dissipation rate.

The total energy balance is:

vðr3Þ þ V$ðrU 3Þ ¼ V$ vt

(10)

with

B ¼ the sum of body forces p0 ¼ the modified pressure U ¼ the velocity vector r ¼ the density.



mt g$Vr rsr

Pkb ¼ 

 1 vr  r vT 

b¼ 

(13) P

with

(9)

b ¼ the thermal expansivity, Tref ¼ the reference temperature rref ¼ the reference density.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

335

2.1.5. Other considerations In the case of a small to medium release where the height and the plume length of the cloud do not exceed the first hundred meters, the Coriolis acceleration may be negligible (Alinot & Masson, 2005). Due to medium-scale release shown in this paper, the Coriolis acceleration was negligible. In a case of a larger release, the Coriolis acceleration should be investigated to estimate its effects. 2.2. Source term The details of the source term are presented in Fig. 1 under boundary condition one. First, the scenario is determined to develop the source term. The size and geometry are setup as a surface boundary of the domain, and the gas phase velocity is specified as an inlet boundary. The gas phase velocity is calculated from the mass evaporation rate (kg/m2 s) by using the gas density. The literature and simple consequence modeling software packages use a mean mass evaporation rate from previous validated tests. The mean value is valid for simple and box models since they do not provide an accurate three-dimensional representation. On the other hand, these models only provide a conservative estimate of the exclusion zones. Mean mass evaporation rate may not represent the evaporation phenomenon correctly. For instance, releases on to open water have proven to be more energetic resulting in local increase of the mass evaporation (e.g., increase of gas phase velocity). Hissong (2007) reported an increase of the heat transfer by a factor between 8 and 10 times during Esso test release on water. During large releases on water, hydrate/ice formation was observed. Hydrate/ice formation is known to influence the mass evaporation rate locally. The simple consequence modeling does not account for these changes in mass evaporation rate. The first scenario developed in this paper describes a release on water. For the release in an impoundment (with physical boundary for the pool spread), as the LNG spread reaches the warm concrete, it will vaporize to cool the concrete ground. High evaporation will occur at constant mass evaporation rate until the concrete is cooled down resulting in a reduction of the heat transfer between the LNG and the concrete. Then, LNG will start to accumulate at the bottom of the impoundment. The mass evaporation rate will be a timedependent function. Depending on the quantity of LNG release, the mass released may be sufficient enough to sustain an LNG pool and the mass evaporation may transit directly into the time-dependent function. The case with a low flow rate was used as the second scenario in this paper. The initial vapor cloud height depends on the combination of source term and wind entrainment velocity. Fig. 2 shows the balance of the forces around the gas phase velocity as it gains altitude. In the case where the wind velocity is dominant, the upward momentum effect is diminished resulting in a low vapor cloud height (e.g., flat cloud). When the wind velocity is low, the upward momentum is more pronounced, resulting in a high cloud height and higher gas concentration downwind. This phenomenon was observed during the November 2007 test at BFTF. As the wind velocity drops, the upward momentum became dominant and the vapor cloud rises upward, increasing significantly the height of the cloud. The turbulence created at the source by evaporation phenomenon may be modeled through the turbulent kinetic energy and dissipation energy by using the equation (Luketa-Hanlin et al., 2007):

k0 ¼ and

2 3 vg Ti 2

(14)

Fig. 2. Balance of forces for LNG release in the estimation of the cloud height.

30 ¼ Cm3=4

k3=2 0:07D

(15)

Where vg is the mean gas phase velocity, Ti is the turbulence intensity parameter, and D is the equivalent LNG pool diameter. The recommended values of turbulence intensity are between 1 and 10% to account for the evaporation effect only. This range is valid for an evaporation effect but not to account the external turbulence created by fluid velocity of LNG landing on the water or wave movement of the water. This factor may be increase above the 10% up bound to account for external effects. During the test performed at BFTF, it was obvious that the combination of evaporation phenomenon and the external turbulence had significant effect on the upward momentum generated inside the fence, which affected the shape of the vapor cloud. Most of the consequence modeling software use a mean mass evaporation rate as their source term and they do not include turbulence parameters for the source. 2.3. Atmospheric condition In consequence modeling, the atmospheric condition is classified by the Pasquill–Gifford stability classes (A–F) with A corresponding to the most unstable and F to moderately stable conditions. The stability conditions are a function of temperature, sensitive heat flux, roughness, wind velocity, and wind direction (negligible in the case of flat terrain) (DeVaull et al., 1995). Unstable conditions occur during the day as the ground is heated by the sun. The warm ground leads to unstable density stratification. This thermal effect leads to high vertical mixing rates, i.e., rapid dilution of the chemical in the atmosphere. Neutral atmospheric condition occurs in the evening when the heat transfer from the ground is very small, nearly zero. During stable stratification, the convection and conduction cool the atmospheric boundary layer leading to limited vertical mixing rates in the atmosphere, i.e., the dilution is insignificant and usually happens during the night. The Monin–Obukhov length demarcates the height below which mechanically generated turbulence dominates in the mixing process. The Monin–Obukhov length is negative for unstable conditions, positive for stable conditions and infinite for neutral conditions (Pasquill, 1974). The atmospheric conditions may be modeled by using an atmosphere boundary layer model related to Monin–Obukhov length. This model predicts temperature, wind velocity, atmospheric stability, and turbulence for the domain. To express the atmospheric conditions based on Monin–Obukhov length (Brown

336

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

et al., 1990), the wind velocity, u0(z) and temperature, T0(z), may be derived as function of vertical z-axis:

u0 ðzÞ ¼

    z  u* z loge  jm zo L ka

T0 ðzÞ ¼ Tw þ

   z T* z loge  jh ka z0 L

(16)

(17)

qw ¼ the sensible heat flux CPair ¼ the heat capacity of air rair ¼ the density of air. The momentum diffusivity and the heat diffusivity are estimated from the equations:

Km ðzÞ ¼

when L > 0

jm fm

z L z L

¼ jh ¼ fh

z  L  z L

¼ 5

z L

¼ 1þ5

z L

jm

jh

  1 þ x2 ¼ 2 loge 2

z  L

fm fh

L

    p 1þx 1 þ x2  2 tan1 ðxÞ þ ¼ 2 loge þ loge 2 2 2

z L

z

x ¼

The turbulence kinetic and energy dissipation term may be expressed as a function of the Monin–Obukhov length as well (Panofsky & Dutton, 1984):

30 ðzÞ ¼ (20)

(21)

1 x2

(23)

 16z 1=4 L

(24)

1

L ¼ the Monin–Obukhov length z ¼ the vertical coordinate (z ¼ 0, at the ground level) z0 ¼ the surface roughness g ¼ the acceleration of gravity ka ¼ the von Karman constant, 0.42

u2* Tw ka gT*

(25)

and

Tw ¼ Ta



 Po m P

(26)

u* ¼ the turbulence friction velocity Tw ¼ the zero level potential temperature T* ¼ the scaling potential temperature Ta ¼ the ambient temperature The sensible heat flux is defined from the estimated variables as well:

q_ w ¼ T* rair CPair u* with

fh ðz=LÞ

u3* z f ka z h L

k0 ðzÞ ¼ 5:48u2*



fh ðLzÞ fm ðLzÞ

(30) 12

(31)

where the constant 5.48 has been experimentally determined for neutral atmospheric layer. A recent test in March 2008 showed an increase of humidity from 51% (outside of the cloud) to 63% (inside the cloud). From this observation, it was concluded that the humidity in a vapor cloud comes from two sources, the humidity already present in the air and the humidity picked up by the release on water. These two sources of water condensate at the contact of cold LNG vapor into droplets. The mechanism of the droplet formation and droplet sizes involved with the humidity from the water surface is not known specifically for LNG release on water. 2.4. Ground effects

The parameters in the wind and temperature equation (i.e., equations (16) and (17)) have the following parameters derived from the historical data:

L ¼

(29)

(19)

(22)

with 

u* ka z

Kh ðzÞ ¼

1 ¼ x ¼

L

(28)

(18)

when L < 0

z

u* ka z

fm ðz=LÞ

(27)

The ground effects are characterized by the obstacles geometry that may affect the simulation results. Each surface may have a different heat flux (or temperature) with respect to the sun exposure and surface roughness. The ground effect also may be significant for free obstacles simulation. For more complex scenarios, obstacles are represented by geometry shape where a surface roughness or/and temperature may be set. Peterson (1990) showed that the effect of large homogeneous roughness may significantly reduce the concentrations downwind. Heterogeneous roughness, when the surface roughness changes as the features of the terrain changes (i.e., from concrete to grass), may also have an effect on the concentration downwind. 3. Parameters’ effect on vapor dispersion One objective of this work was to identify not only key parameters but also to estimate the effect of the parameters on the vapor dispersion. As presented in the previous section, the three sets of parameters were regrouped into models to represent the boundary conditions of the domain. Each parameter had an effect on the cloud concentration, length, or shape. To understand or better control the hazard of a flammable vapor cloud, it is essential to know the effects of various parameters on a vapor cloud. It may be possible to control some of these parameters to diminish the hazard through design. Table 1 regroups the parameters discussed in the previous section with their effects on the vapor dispersion if known.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

337

Table 1 Effects of parameter on the vapor dispersion. Parameters

Properties/range/estimation

Unknown or known effects

Release rate Pool shape Mass evaporation rate

Quantity Rectangular, circular, physical boundary

Relatively known, increase the size of the cloud as the rate increase Unknown

0.029–0.185 kg/m2 s recommended mean values Experiment may provide an evaporation profile

Unknown effects from low or high mean values Unknown effects of mass evaporation rate profile inside a mean value

Pool area

Estimated from the mass evaporation rate, pool shape, and release rate Referred as turbulence factor, from 1 to 10

Unknown Unknown

1–10%

Unknown

In the range of two hundred positive or negative

Known, thermal effect to the cloud:

Turbulence factor due to the release scenario Turbulence intensity due to evaporation phenomena Sensible heat flux from surrounding

Negative sensible heat flux increases the downwind distance by limiting the vertical mixing rate. Positive increases the vertical mixing rate inside the cloud, which reduce the downwind distance Wind velocity

Atmospheric boundary layer, wind profile respect to height

Temperature ambient

Temperature profile respect to height and sensitive heat flux

Known for large wind velocity Unknown for low wind velocity Forced thermal mixing with cold LNG vapor Reduction of downwind distance

Surface roughness Obstacles Wind direction

Values known i.e. passive barrier Obstacles related

Known effects for homogeneous roughness terrain Unknown Unknown

4. Simulation setup and calibration The numerical simulation was performed for two types of release, on water and on concrete. Both cases were performed during November 2007 at the BFTF. The data collected was used to calibrate the model with respect to a set of conditions, and then a sensitivity analysis was performed on selected parameters. The tests carried out were designed with two objectives: (i) refine the source term (ii) and collect data with obstacles to compare with CFD modeling result. Thereby gas concentrations and temperature measurements were collected outside and inside the release area. A vapor fence was used around the release area to create obstacles effect. For the release onto water, thermocouples were placed underneath water to capture the heat transfer to the water. 4.1. Release on water with vapor fence 4.1.1. Site description and domain definition The larger pit (10 m  6.7 m  1.3 m) at the BFTF was filled with 1.3 m of water and a 1.3 m vapor fence was installed at the edge and around the water surface, shown in Fig. 3. Two foam generators were installed at 1.3 m above the ground to apply high expansion foam (500:1). Foam generators were not utilized during this test, but they were used for mitigation as necessary. LNG was released on water through a ramp to ease the vertical velocity, but the horizontal fluid velocity was conserved. The main obstacles of the site were identified as the vapor fence and the foam generators. Thereby, they were modeled in the domain. Gas sensors were placed in the downwind direction of the release area. The gas poles are noticeable in Fig. 3 in the back of the release area (white and orange bands). The strip on the pole may be used as a scale, each of which was 0.3 m. A large domain was developed to enclose the sensors and also to avoid undesirable effects from the boundary conditions on the vapor dispersion. A factor of three was used to estimate the downwind distance of the domain instead of the half LFL downwind distance. This factor was used to reduce effect from the boundary conditions even though

open boundary conditions were used. Based on Cartesian coordinates (x, y, z), the size of the domain was a rectangular parallelepiped from (66, 56, 0) to (55, 88, 15) m with a reference point (0, 0, 0), shown in Fig. 4. The maximum LFL downwind distance was reported around 25 m in the y-axis. CFX codes have the open boundary conditions which let the fluids freely come in and out of the domain, which reduces considerably the effect of the boundary condition on the results compared to a conventional input boundary condition that will force the fluid inlet. Positive y-axis and positive x-axis represented South East (SE) and South West (SW) directions respectively. Positive z-axis was the height of the site with zero being at the ground level. Wall boundary conditions were used to simulate the vapor fence, the water surface that was not covered by LNG, the foam generators and the concrete ground. The mesh was refined up to 0.024 m in the area where high gradient was expected (i.e., at the source term). The refinement was high at the source term (i.e., inside the fence) and at the water area, but the refinement was mild around the foam generators and the vapor fence. A hexahedral mesh was used as an inflation layer above the source term and then a tetrahedral mesh for the rest of the domain, as shown in Fig. 4. The foam generators were also included in this refinement procedure due to the circular geometry and their location above the source term. The ground of the domain had a minor refinement to enhance the buoyant effect in the cloud. Several mesh refinements were performed to find the number of elements necessary to obtain an independent grid solution. A simulation without the source term (i.e., no methane inlet) was used to determine the optimal mesh and the effect of boundary conditions by monitoring the wind profile with respect to the z-axis and TKE along the x-axis. These profiles were compared through different mesh sizes. Then once the wind and TKE profiles reached an independent solution compared to the mesh size, it was concluded that the solution is independent of the mesh size. The optimal mesh was defined to be 730,000 elements for the release on water. Fig. 5 shows the difference of the wind profile from one side of the domain to the

338

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 3. Test representation.

other side with respect to the z-axis direction. The wind had less than 6% change ‘‘effect’’ form the boundary conditions, which corresponded to 0.04 m/s change in the wind velocity. 4.1.2. Site description and domain definition The boundary conditions setup was based on the data collected during the test. The release onto water test was performed twice at different weather conditions. The characteristics of the tests are presented in Table 2. The atmospheric boundary condition was set to opening boundary condition. The wind and temperature profile was developed using equations (16)–(24) as a function of the Monin– Obukhov length. The friction velocity, roughness, and potential temperature were the only inputs to this model. Other variables presented in Table 2, Richardson number, heat and momentum diffusivity, were used to verify the accuracy of the wind profile

using equations (25), (28) and (29). The turbulence for the atmospheric boundary condition was set through the k–3 two-equation turbulence model, where k and 3 followed the equations (30) and (31) also related to the Monin–Obukhov length. The values presented in Table 2 only gave a general description of the conditions for the overall tests. During the TEEX1, the wind velocity at 2.3 m varied from 1 m/s to 3 m/s and the wind direction varied from SE to S. For a fair comparison with the collected data, the wind velocity and direction were set according to the particular time and then compared. For different wind profile, the friction velocity was adjusted which governed the rest of the variables with respect to the time selected. The humidity already present in the air (about 32  1% relative humidity for TEEX1 and TEEX2) was accounted for in the modeling by changing the thermodynamic properties of the air. Because of the lack of information, the condensation effect of the water from

Fig. 4. Mesh details and obstacles.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

339

Fig. 6. Source term characteristics from the water temperature profile during TEEX1. Fig. 5. Wind velocity change from (4.5, 50, z) to (4.5, 80, z) with respect to z-axis.

the pool was not incorporated in the modeling (i.e., size droplet formation). The source term was defined as inlet boundary condition in the domain. The source term is composed of five elements: geometry, area, gas phase velocity profile, turbulence, and location. The test was designed to create a free-obstacle pool formation (i.e., no interaction with the boundary of the pit). The gas phase velocity profile was derived from the mass evaporation rate data collected, shown in Fig. 6. This contour was obtained by Kriging interpolation (ASTM D5523, 1996) from the twelve measured points. The figure represents the vertical heat transfer through the water located 45 cm underneath the water surface. The mass evaporation rate was found in the range between 0.1 and 1 kg/m2 s. The pool area was optimized with three conditions: (i) liquid mass flow rate of the discharge hose (5.2 kg/s with liquid density of 425 kg/m3), (ii) mean value of mass evaporation (also estimated from 0.15 to 0.18 kg/m2 s) and (iii) the mass evaporation rate profile (Fig. 6). Due to uncertainties associated with the data collected and the interpolation methodology used, several iterations of the following method were needed to find the optimal area: 1. The initial area was estimated using the mean mass evaporation rate and the liquid mass rate of the discharge hose.

Table 2 Summary of the collected data for release on water. Test Name

TEEX1

TEEX2

Mean flow rate, m3/min Valve open time, s Mean wind velocity @ 2.3 m Mean wind direction @ 2.3 m Wind profile At 2.3 m At 10 m u*, m/s Temperature Profile At 2.3 m At 10 m T0(z)  Tw,  C Absolute air pressure, bar Relative Humidity, % Dew Point,  C Stability Class Solar Radiation, W/m2 Sensible heat flux, W/m2 Heat diffusivity @ 2.3 m Momentum diffusivity @ 2.3 m Richardson Number @ 2.3 m Monin–Obukhov length, m Roughness length, m

0.75  0.09 619 1.2  0.38 160  18

0.63a 763 2.3  0.23 140  34

1.2  0.38 1.9  0.58 0.053

2.3  0.23 3.3  0.78 0.175

15.9  0.26 15.9  0.19 0.024 1.013  0.001 32.6  1.6 30.8  0.98 B 593.7  48.9 1.54 0.024 0.021 0.24 8.3 0.00085

16.5  0.61 16.8  0.71 0.10 1.012  0.001 31.2  1.2 31.8  0.97 C 704.4  10.4 44.8 0.30 0.30 0.026 85.8 0.00085

a

No data, estimated from total volume.

2. The total liquid evaporation rate was estimated by integrating the mass evaporation profile over the initial area (from step 1) 3. If the integrated mass evaporated was equal to the liquid mass rate of the discharge hose, the LNG pool area was optimized. If the integrated mass evaporated was too large, the area of the pool was reduced and if the integrated mass evaporated was too small compared to liquid discharge hose, then the pool area was increased. For TEEX1 and TEEX2, the optimal pool area was found to be 15 and 13 m2 with the mass evaporation profile, respectively. By using the methane gas density, the velocity profile was derived and integrated into the inlet boundary. Two pool geometries were used during this study, circular and rectangular. A 3-inch diameter pipe was used to release LNG onto water. The horizontal fluid velocity of the discharge hose created turbulence inside the vapor fence. This turbulence was simulated by increasing the values of turbulence intensity of the source term or by incorporating an inlet jet representing the effect of the fluid velocity. An inlet jet velocity of 4.3 m/s was added to simulate the steady state fluid flow (no vapor formation out of the discharge hose). The ground was set as a wall boundary. Wall boundary acts a solid (impermeable) boundary to fluid flow. It does not let any flow to go through but allows heat flux. For this model, the ground was set with a surface roughness and sensible heat flux corresponding to the sensible heat flux reported in Table 2. The water surface was set to wall boundary conditions with a measured temperature of 293 K and a surface roughness for water at 4  104 m from the literature. The vapor fence and foam generators were set with a smooth surface and an adiabatic heat source rather than a surface roughness or/and a material temperature. 4.1.3. Simulation results The data collected during the November 2007 test provided information to help in the evaluation of the visible behavior/ features of cloud. The visual data as well as sensor data were useful in the verification of the output from the simulation results. Fig. 7 shows the side view of the release area with two types of camera at the same time frame. Methane gas at room temperature is invisible to the naked eyes. When LNG was released, the cryogenic temperature condensed the moisture in the air making a white cloud visible to the naked eyes. To be able to monitor the movement of the methane, two types of cameras were synchronized together. A normal and an infrared hydrocarbon camera were used. The infrared hydrocarbon camera shows the carbon atoms. In the hydrocarbon industry, an infrared camera is commonly used to detect leaks/vapor of hydrocarbon invisible to the naked eye. The black color represents the cold carbon atoms; the white color represents the warm carbon atoms. This picture shows the height

340

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 7. Video pictures of release on water, infrared hydrocarbon camera (top) and normal camera (bottom), TEEX1.

of the white condensate cloud and the angle of the cloud traveling at the edge of the barrier. TEEX1 and TEEX2 release on water presented interesting features during the generation of the vapors. High gas phase velocity was noted leaving the release area during the tests which indicates high turbulence inside the vapor fence. High evaporation of LNG and fluid velocity combined with low wind velocity created a large vapor cloud above the release area (concept illustrated in Fig. 2). To be able to perform the data analysis using CFD code, three cases were defined and then simulated as a continuous simulation. Continuous simulations represent the steady state for a set of parameters. For the TEEX1 tests, the wind velocity was set at 3 m/s at 2.3 m and wind direction as SE. For the TEEX2 tests, the wind velocity was set at 7.4 and 6.8 m/s with a wind direction to SE and S respectively, shown in Fig. 8. 4.1.3.1. TEEX1 simulation results. The friction velocity was modified to accommodate the wind velocity of 3 m/s, but the other parameters were kept constant as shown in Table 2. The gas phase velocity profile was developed following the three steps presented previously to optimize the LNG pool area. Two LNG pool geometries were also tested. The LFL was used as a monitoring value for the downwind distance at 0.5 m high (in the z-axis) and at the discharge hose location (with 4.5 m value in the x-axis) for each run and is reported in Table 3. Several runs were performed to assess the effect of TKE located at the source term and LNG pool shape parameters, reported in Table 3. TKE and LNG pool geometry were varied because they were not directly measured during the test. TKE parameter was increased up to 5.8 m2/s2 to account not only for evaporation turbulence but also for the turbulence created by the fluid velocity from the discharge hose. The turbulence intensity had a small effect on the LFL distance, but the LNG pool geometry had a greater effect.

The rectangular LNG pool geometry gave a mean LFL distance of 19.5  1.5 m. The circular LNG pool geometry gave a mean LFL distance of 15.5  0.5 m, which was reported as being the LFL distance during the November 2007 test with wind velocity of 3 m/ s. The simulation results with a circular LNG pool geometry were in agreement with the collected data. The results reported in Table 3

Fig. 8. LFL distance estimated from the data collected during TEEX1 and TEEX2.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352 Table 3 TEEX1 simulation summary for a wind velocity of 3 m/s. Run

LNG pool geometry shape

Turbulent kinetic energy at the source term, m2/s2

Turbulence intensity, with Eq. (13), Ti, %

LFL distance SE (4.5, y, 0.5), m

TEEX1_1 TEEX1_2 TEEX1_3 TEEX1_4 TEEX1_5 TEEX1_6 TEEX1_7 TEEX1_8

Rectangular Rectangular Rectangular Rectangular Circular Circular Circular Circular

5.8  106 5.8  104 5.8  102 5.8 5.8  106 5.8  104 5.8  102 5.8

1 10 100 1000 1 10 100 1000

20 18 20 21 15 16 15 16

showed an overestimation in the LFL distance. Fig. 9 shows the comparison of concentration profile between the simulation results and the data collected during the November 2007 tests. Even though the LFL distance matched the experimental data, the concentration levels inside the cloud were off as shown in Fig. 9. The gas sensors were installed within the vapor cloud centerline in a southerly direction. As the test proceeded, the wind shifted letting only half the cloud within the array of sensors, presented in the right side of Fig. 9. For the same location, the highest concentration measured was 27% v/v during the test, and the simulation showed a concentration level reaching 32% v/v. During the TEEX1 tests, the temperature profile of the pit was measured at 36 locations inside the fence limits. Fig. 10 shows the concentration derived from these temperature measurements to provide an idea of the cloud shape inside the fence. These data indicated the pattern of the LNG pool inside the wall. LNG cloud was not covering the entire water surface area and the concentration was not uniform inside the fence. High concentrations of methane were detected against the upwind wall. The concentrations shown in Fig. 10 were not computed with the humidity present in the cloud, which resulted in a lower concentration level compared to the real concentration. Nevertheless, the pattern still represented the whereabouts of the vapor cloud and the shape. This shape indicated the high turbulence created inside the fence. If there were no turbulence, the concentration would accumulate evenly over the entire surface and would move gradually upward creating a layer of concentration like a ‘‘plug flow.’’ The shape of the inside cloud and the LFL distance were replicated by using the CFD codes. The concentration levels were accumulating on the upwind wall as shown in the simulation results and the data collected during the test. These agreements also showed the applicability of CFD codes and their usefulness in the estimation of the concentration levels as well as the effect of obstacles. Simple models cannot reproduce the behavior of the cloud inside the fence.

341

4.1.3.2. TEEX2 simulation results. The model was also tested against the data collected for TEEX2. The friction velocity was set to 0.383 and 0.353 m/s for a wind velocity of 7.4 and 6.8 m/s at 2.3 m above the ground. The air temperature was kept constant as the value reported in Table 2. The gas phase profile was developed using the same procedure shown previously. The results are presented in Table 4. The collected data showed an LFL downwind distance of 11.7  0.5 and 15.9  0.5 m for the high wind velocity of 7.4 m/s for SE and 6.8 m/s for S direction respectively. The numerical results showed good agreement with the collected data again. For the SE direction, the computed LFL downwind distance averaged to 12.5  0.5 and 13.5 m for different pool geometry, which slightly overestimated the reported value. For the S direction, the computed LFL downwind distance averaged to 15.6  0.5 and 17 m for different pool geometry. The reported values were close enough to agree on. The turbulence intensity did not significantly affect the LFL distances in the TEEX1 and TEEX2 tests. The plume shape is illustrated in Fig. 11 for TEEX2. The concentration of the cloud outside the fence was recorded up to 14% v/v during the test, while the simulation results showed a concentration up to 15% v/v. The simulation showed an agreement in matter of LFL distance, concentration level, and the shape. The same concentration profile inside the fence was noticed during TEEX2 shown in Fig. 11 where higher concentration levels were accumulating on the upwind wall. The model was able to reproduce the concentration levels, the shape of the cloud, and the LFL distance. At higher wind speed, it was noticed that the model seemed to get a better agreement, which underlined the need to understand the turbulence effect inside the fence. The turbulence effect on the concentration level in the downwind distance was strongly influenced by the wind velocity. At high wind velocity, the turbulence effects were diminished when compared to low wind velocity. 4.2. LNG release in concrete pit 4.2.1. Site description and domain definition In the scenario of LNG release on concrete, TEEX7, the pit was filled with LNG to carry out a fire test. The larger pit (10 m  6.7 m  1.3 m) at the BFTF was filled with 15 cm of LNG and a 1.3 m vapor fence was installed around the edge of the pit. The pit was 2.6 m deep with a fence of 1.3 m above the ground. Fig. 12 shows the mesh definition for this test. The foam generators are the same design as presented in the previous section The grid was refined around the pit area, inside, and around the obstacles where a large gradient is expected. At the bottom of the pit, the mesh was also refined using a hexahedral mesh of 6 cm. The

Fig. 9. Comparison of concentration profiles at 0.5 m above the ground level, simulation (left) and collected data (right), TEEX1.

342

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 10. Iso-surface concentration inside the pit based on temperature estimation, TEEX1.

domain size was (43, 35, 1.3) and (34, 60, 16) in Cartesian coordinates. The grid independent test was also performed and it was found that the optimal grid size was 750,000 elements. To setup the source term, the mass evaporation rate had to be determined first. Different from LNG release on water, the heat sink in this scenario is not infinite. The estimation of the mass evaporation rate profile was calculated through theoretical analysis. The study from Briscoe and Shaw (1980) was modified to estimate the mass evaporation rate. Since the study is applicable for high flow (i.e., formation of pool is instantaneous), these equations were modified for low flow release. The heat balance around the vaporizing liquid pool on land includes the solar radiation, long wave radiation, convection from wind, latent heat loss, and the conduction from concrete. The dominant heat is from the concrete. But as the concrete started cooling due to the contact with cryogenic liquid, the heat flux decreased, resulting in reduction of the mass evaporation rate. The contribution of long wave radiation, solar radiation and convection from the wind is less than 4%. Thereby the radiation and convection heat were ignored. The mass evaporation rate due to the concrete heat transfer was given by

mcond ¼

kðTA  TB Þ

(32)

lðpatÞ1=2

With a low flow rate, the initial flow rate did not cool the pit enough to sustain an immediate pool formation. The pit had to be cooled to

TEEX2_1 TEEX2_2 TEEX2_3 TEEX2_4 TEEX2_5 TEEX2_6 TEEX2_7 TEEX2_8 TEEX2_9 TEEX2_10 TEEX2_11

Wind Turbulent LFL distance velocity at SE kinetic (4.5, y, 0.5) S pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi energy, m2/s2 2.3 m, m/s ð x2 þ y2 ; 0:5Þ, m

SE SE SE SE SE SE S S S S S

5.8  106 5.8  104 5.8  102 5.8 5.8  104 5.8  102 5.8  106 5.8  102 5.8 5.8  104 6.1  102

7.4 7.4 7.4 7.4 7.4 7.4 6.8 6.8 6.8 6.8 6.8

t2 ¼

 kðTground  TLNG ÞA 2 pa mLNG l 1



1 mLNG l

Z

t1

0

(33)

  k Tground  TLNG ðpatÞ1=2

A dt

(34)

where A ¼ the surface area of the pit k ¼ the thermal conductivity of the ground a ¼ the thermal diffusivity of the ground mLNG ¼ the mass rate of LNG released l ¼ the latent heat for LNG Tground ¼ the temperature of the surface TLNG ¼ the temperature of the liquid LNG t1 ¼ the relative time required to equalize heat flux of the concrete to the heat flux of the mass released for a continuous release t2 ¼ the actual time to reach the equilibrium for the mass released.

(

Wind Geometry direction pool shape Rectangular Rectangular Rectangular Rectangular Circular Circular Rectangular Rectangular Rectangular Circular Circular

t1 ¼

Once the actual time to reach the equilibrium is calculated, the mass evaporation rate may be derived from the following:

Table 4 TEEX2 simulation results for SE at 7.4 m/s and S at 6.8 m/s. Run

the point that the pool would start forming. The pool starts forming when the heat flux of the concrete is equivalent to the heat flux provided by cryogenic liquid discharge and may be formulated with the following equation:

12 12.5 12 13 13.5 13.5 15 15 17 17 17

mevap ðtÞ ¼

mLNG kðTground TLNG Þ

1=2 A

lðpaðt1 ðtt2 ÞÞÞ

t  t2 t > t2

(35)

For the LNG flow rate used in this test, it requires 750 gallons per minute (gpm) (20 kg/s) to obtain an immediate pool formation. With a flow rate of 90 gpm (2.4 kg/s) and assuming 15% flashed (not contact with the concrete), it would take one hour to fill the pit up to 15 cm. This time estimated was used for a total pit surface area of 67.5 m2. The discharge pipe was placed almost at the bottom with a ‘‘sock’’ to prevent any splash. At the beginning of the release, a large cloud was noticed probably due to high evaporation rate from the cooling of the pipe and the initial cooling of concrete. As

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

343

Fig. 11. Concentration contour for S and SE wind direction with different turbulence intensity for TEEX2.

the cryogenic liquid started building up at the bottom of the pit, the visible cloud became smaller. Fig. 13 shows the cloud through the normal camera and the infrared hydrocarbon camera at synchronized time and illustrates the size of the cloud as it came out of the vapor fence. The green line represents the hydrocarbon vapor

captured by the infrared camera. The red line represents the edge of the fence and the foam generator for reference (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.). Then the lines were superposed with the regular camera. The size of the cloud was

Fig. 12. Concrete pit mesh definition.

344

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 13. Video picture of release on concrete with normal and infrared hydrocarbon camera.

larger with the hydrocarbon camera than with the normal camera for the same picture scale. The height of the white cloud (i.e. regular camera) was lower than the black cloud (i.e. infrared camera). This was due to the fact that methane gas is not visible at room temperature. The difference between the two cameras was more pronounced here due to lower concentration of water in the cloud. At the beginning of the release, the turbulence intensity should be around its maximum and once the LNG starts building up at the bottom of the pit and reaches the steady state, the turbulence intensity should be at its minimum. The release on concrete took one hour, while the wind velocity varied from 1 to 0.5 m/s with a constant wind direction SE. Due to extreme low wind velocity, the wind model presented in Equation (16) was not utilized. Therefore, a simple wind profile, called power law correlation, was used:

u ¼ u10

 z p 10

(36)

where u10 is the wind velocity at 10 m, p the dimensionless parameter, whose value depends upon atmospheric stability category and surface roughness (CCPS, 1999). The power law wind velocity may also be used when the wind profile is not known. With the low wind velocity and the stable conditions, the atmospheric stability was estimated to be an F on the Pasquill–Gifford stability class, therefore the dimensionless parameter, p, was set to 0.55. Table 5 shows the parameters for release on concrete. For transient simulation, the gas phase velocity, the wind velocity, and the turbulence intensity were changing during the time of the test. The CFX’s expression editor was used. The mass

Table 5 TEEX7 release on concrete parameters. Time

u10, m/s

T, K

Wind direction

Turbulence intensity, %

0–1800 1800–3000

1–0.6 0.6

295–293 293–292

SE SE

10–1 1

evaporation rate presented in equation (37) was setup in the expression editor where it was converted to gas phase velocity:

vg ðtÞ ¼

mevap ðtÞ ArLNG

(37)

where vg is the gas phase velocity in m/s, and rLNG is the gas density of methane at 110 K. In both cases (continuous and transient mode), the threshold value for the root-mean-square (RMS) of residuals was set equal to 104. 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. 4.2.2. Simulation results The temperature profile of the pit was converted into a concentration profile shown in Fig. 14. The concentration profile showed a layer formation from the evaporating pool at the bottom of the pit. The 15 cm at the bottom provided a constant and uniform gas phase velocity. The main observation was that the concentration layers inside the fence were different compared to the release on water. This indicated the absence of turbulence at the source term. Based on this observation, the turbulence intensity was set to one and the simulation results are presented in Fig. 15. The concentration profile was in agreement with the data collected during the TEEX7 tests. The 50% v/v layer was at the same level in both sets of data. The outside visible cloud had a flat shape, illustrated in Figs. 13 and 15. At the edge of the fence, it was noticeable that the cold LNG vapors were falling down to the ground with a wind velocity close to 0.6 m/s at 10 m, which also meant that wind velocity at the ground level was even lower. This effect increased the presence of white cloud around the outside fence. The simulation results are presented in Fig. 15. The concentration data was plotted against the simulation results from the edge of the fence. During TEEX7 and because of low wind speed, the heavy vapor cloud fell around the fence creating a cloud around the fence. The simulation was able to reproduce this

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

345

Fig. 14. Concentration profile for TEEX7 inside the fence and pit. Collected data (top) and simulation results (bottom).

behavior as shown in Fig. 15. The same angle was given to the pictures for comparison. The data collected over the hour was compared to the simulation results presented in Fig. 16. The preliminary results showed that the simulation was under predicting the concentration at the beginning and the end of the

profile. The delay in the beginning of the curve was due to the mass evaporation profile, which provides an instantaneous mass evaporation rate in the simulation, but in reality, the flow rate at the liquid discharge rate at the end of the discharge hose was progressively increasing until the full rate was reached. The end of

Fig. 15. Illustration of the vapor falling down from the edge of the fence. Picture from TEEX7 (right) and simulation results with low turbulence intensity (right).

346

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 16. Concentration profile at 1.3 m at the centerline of the cloud at the edge of the fence.

the profile may be investigated by increasing the turbulence intensity. 5. Parametric analysis The parametric analysis was conducted in three areas: turbulence models, source term, and the atmospheric condition. Different turbulence models were compared to evaluate their effect on the concentration profile and shape of the cloud. The source term is composed of several parameters: each one was varied as the other ones were kept constant. The source term parametric analysis focused on the turbulence effect and geometric configuration. The atmospheric condition parameters were varied using the same procedure as the source term. 5.1. Turbulence equation effects RNG k–3 and SSG Reynolds stress turbulence models were compared against the k–3 turbulence model. The RNG k–3 turbulence model is based on renormalization group analysis of the Navier–Stokes equations. The transport equations for turbulence generation and dissipation are the same as those for the standard k–3 model, but the model constants differ. SSG had shown superior predictive performance compared to eddy-viscosity models in cases, such as buoyant flows or free shear flows with strong anisotropy. SSG model solves differential transport equations individually for each Reynolds stress component. Compared to the k–3 turbulence model, the Reynolds Stress model has six additional transport equations that are solved for each time step or outer coefficient loop in the flow solver. The source terms in the Reynolds Stress equations are also more complex than those of the k–3 turbulence model. As a result of these factors, outer loop convergence may be slower for the Reynolds Stress model than for the k–3 turbulence model. To assess the effect of each turbulence model, the following scenario was changed from the original calibration. The released mass was set to 6.7 kg/s. The released mass was modified (i) to obtain a large cloud in order to increase the parameter effects and (ii) to obtain a point of comparison between two mass released rates. The friction velocity and the potential temperature were set to 0.053 m/s and 0.023 K, which corresponded to low wind velocity, around 1.2 m/s at 2.3 m. At this condition, the vapor cloud behavior was gravity driven, which increased the effect of the turbulence in the vapor cloud concentration profile. The mass evaporation profile shown before was kept as well with an evaporating pool of 15 m2 and the simulation results are shown in Fig. 17.

Fig. 17 shows not only the concentration at the centerline of the cloud but also the concentration contour at 0.5 m above ground for each turbulence model. The fraction of methane was separated into four sections. The concentration line and concentration contours were needed for efficient comparison. The preliminary observation showed that the turbulence models had a minor effect on the concentration profile of cloud and the spreading of the cloud. The computed LFL downwind distances were 31.6, 30.5 and 26.8 m at the centerline for k–3, RNG k–3, and SSG Reynolds stress turbulence model respectively. The concentration contours showed the difference in the concentration level and each model had a slightly different cloud shape. The concentration contours for RNG k–3 and SSG were wider than with k–3 turbulence model. With the RNG k–3 and SSG model, the vapor cloud was coming out on the corner rather than from the downwind fence compared to k–3 model results. The shape of the contour was slightly different for each turbulence model; the SSG Reynolds stress turbulence model had a wider shape. It was noticed that the LFL distance was reduced as the cloud became wider. Meanwhile, the surface area for each simulation was observed to be the same for the fraction greater than 5% v/v. The three turbulence models studied here did not demonstrate a significant difference (about 10% on the average value) in matter of LFL distance. 5.2. Source term effects Four parameters for the source term were selected and their effects on the LFL distance as well as on the downwind concentration were analyzed. The turbulence intensity at the source, gas phase velocity, pool shape, and area varied at a time while other parameters were kept constant. The same scenario as described in the above was kept. The mass evaporation rate was in the recommended range, between 0.029 and 0.185 kg/m2 s. For this recommendation, it was possible to build two cases: (i) with a maximum LNG pool possible of 67.5 m2, i.e., total water surface area, the minimal possible mass evaporation rate was 0.1 kg/m2 s., and (ii) with a mass evaporation rate of 0.185 kg/m2 s, the LNG pool area was set to 35 m2. The geometry of the pool was also modified from a rectangular to circular shape. The scenario developed with the study on turbulence models was kept the same and the results are presented in Fig. 18 by using concentration contours.

Fig. 17. Simulation results from the three turbulence models. The concentration along the centerline of the cloud and the contours.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

The different LNG pool areas and geometries were simulated using a TKE at 5.8  106 m2/s2 at the source term. At low turbulence, LFL distance, shape and methane concentration inside the vapor cloud were influenced by the different source term configuration. The initial formation of the cloud above the evaporating pool influenced the concentration downwind. The LNG pool geometry had an effect on the upper flammable limit (UFL) and LFL as well as the LNG pool area. The source term area had an influence in the width of the cloud. As observed in the previous turbulence study, the wider the cloud was, the shorter the LFL distance. While the use of several combinations of source term was a benefit, a certain level of uncertainty had to be estimated as a result. As shown in Fig. 18, the methane fraction was fairly in the same range of LFL downwind distance. The averaged LFL downwind distance was estimated to be 52.5  4.1 m. There was a deviation of about 7.8% between each scenario, which was for low turbulence at the source term. The second set of runs was kept with the same setup but with a TKE at 5.8  102 m2/s2 shown in Fig. 19. This TKE represented an increase of turbulence intensity by a factor of 100 (from 1% to 100%). This increase was justified to simulate the effect of external turbulence that was created by the release of LNG onto water. The same pattern was observed for the rectangular shape with an increase of downwind distance as the LNG surface area decreased. The circular pool showed a completely different result. The cloud was wider which resulted in a shorter LFL distance. TKE increased the turbulence at the source term resulting in a wider cloud around the fence when Figs. 18 and 19 are compared. As TKE increased at the source term, the effects in the downwind direction were more pronounced for low wind velocity. The averaged LFL downwind distance was estimated to be 47.2  4.7 m with a deviation of 9.9% between each scenario. The variance was larger with a T at 100% than at 1%. The significant variance in the LFL downwind distance with different T strongly indicates that TKE values at the source term were a key parameter. The size and the geometry were also considered as being of medium importance. The change in the LFL distances was further analyzed by investigating the vertical concentration profile inside the fence using a 67.5 m2 LNG pool where the TKE was set to 5.8  106 and 5.8  102 m2/s2, shown in Fig. 20. The simulation results showed a significant change in the vertical concentration contour showing the effects of the turbulence at the source term. The maximum height for each simulation

Fig. 18. Downwind distance for different LNG pool geometry with k0 at 5.8  106 m2/s2 (Ti at 1%) with a flow rate of 6.7 kg/s.

347

Fig. 19. Downwind distance for different LNG pool geometry with 5.8  102 m2/s2 (Ti at 100%) with a flow rate of 6.7 kg/s.

was measured at 2.8 and 3.9 m at the 5% v/v fraction for 5.8  106 and 5.8  102 m2/s2 respectively. The reduction of LFL distance may be explained through the change in initial height of the cloud. High turbulence intensity increased the mixing at the source term resulting in a high upward momentum for the vapors, thereby reducing the concentration downwind. At the edge of the fence, the gas sensors indicated higher concentration at 1.3 m than that at 0.5 m. For the infrared camera, the height of the cloud was estimated at 4.75 m. When comparing with the simulation results, the data indicated a source term with high TKE (above 5.8  102 m2/ s2) during the experimental test. The increase of height and width of the cloud distributed the concentration through a threedimensional space which resulted in LFL downwind distance. The effect of the LNG pool geometry was also analyzed by looking at the vapors mixing inside the fence. Fig. 21 shows the velocity profile above the source for a 67.5 and 35 m2 LNG pool with a TKE at 5.8  106 m2/s2. The figure was combined with two types of views for each scenario to get a three-dimensional analysis, ‘‘bird view’’ and a vertical view. ‘‘Bird view’’ is a combination of side and top view. This view was used to illustrate the effect of the pool geometry on the air circulation above the source term and in between the fence. The change of velocity was set to the horizontal contour to localize the upward and downward air movement with respect to height. The vertical view showed the air circulation by using arrows, which indicated the magnitude and the direction of the air velocity. The arrows helped to locate the formation of eddy inside and outside the fence. The velocity profile near the source term revealed the formation of eddies. Each eddy was marked with a number to identify its location. For the 67.5 m2 LNG pool, eddies only formed near the top of the fence and at the upwind fence, illustrated in Fig. 21. This resulted in producing a low height for the vapor cloud as shown in Fig. 20. The source term with a 67.5 m2 LNG pool covered the entire area inside the fence, leaving no room for turbulence inside the fence and only localized eddies at the upwind fence due to the wind velocity. Several eddies (up to four in the vertical plan) were noticeable from the analysis of the simulation results for velocity profile for a 35 m2 circular LNG pool. The circular LNG pool of 35 m2 had more free space between the source term and the fence, which permitted the formation of several eddies around the LNG pool, shown as three and four in the top of Fig. 21. Eddy increased the mixing effects of the methane vapors with air before the vapors were entrained downwind by air. The difference in the number of eddies formed inside the fence was related to

348

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 20. Methane fraction contour for 67.5 m2 LNG pool with a TKE 5.8  106 and 5.8  102 m2/s2 at the source term.

Fig. 21. Velocity profile above the source term for a 67.5 (bottom) and 35 m2 (top) LNG pool with Ti at 1%.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

the size of the LNG pool. Several large air pockets inside the fence were noticeable on the video from the test, which confirmed the shape of the pool and eddy formation observed in the analysis of the simulation. Either with low or high TKE, it was noticeable that the concentration levels were forming uniformly above the water surface for a 67.5 m2 LNG pool. This effect was due to the uniform mass evaporation rate used over the entire surface area, which created a plug flow effect as the vapors rose inside the fence. TKE affected only the upper part of the concentration profile in the pit. TKE had more effect in the concentration profile as well as in height with a smaller LNG pool area. The space between the fence and the LNG pool helped in the formation of eddies; thereby, it increased the mixing effect inside the fence. A low wind velocity was set up according to the data collected for the TEEX1 tests in Table 2 to study the effect of the TKE on the concentration and the shape of the vapor cloud. The optimized LNG pool size was used with the gas phase velocity profile. The circular pool was used to allow room for the eddies to form. Fig. 22 shows the concentrations along the downwind distance at 0.5 m above the ground level along with their respective contours. The simulation results showed a decrease of the concentration downwind as the TKE increased up to 5.8 m2/s2. The change of LFL distance was significant for turbulence intensity from 5.8  102 to 5.8 m2/s2. The decrease of the LFL distance was mainly due to turbulence at the source term. As the eddies increased at the source term and inside the fence due to large TKE value, the escaped the fence on the upwind corner rather than the downwind corner with a low TKE value. This effect is illustrated in Fig. 22 with the contour where it is noticeable the higher concentration was appearing on the side of the upwind fence rather than on the downwind fence. The cloud was also becoming wider as the TKE increased, resulting in lower LFL downwind distance.

5.3. Atmospheric conditions’ effect The sensible heat flux and wind velocity influence the thermal mechanism inside the cloud as it travels downwind. The heat flux changes the mixing rate inside the cloud, while the wind entrains and mixes the vapor with air. Both mechanisms affect the shape of the cloud and the length. Low wind velocity will cause the cloud to be influenced by gravity, resulting in a wide and large cloud

349

traveling at ground level. High wind velocity increases the entrainment of the cloud and forces the cloud to dilute quickly. The shape of the cloud will be elongated and narrow. Positive sensible heat flux will increase the mechanical turbulence inside the cloud, which will result in dilution of the cloud. Negative sensible heat flux will suppress the mechanical turbulence inside the cloud, which will lead to a stable density profile. The objective of this section was to demonstrate the capability of the CFD code to reproduce these described effects by using the basic physics behind the code (i.e. using mass, energy, and momentum balance). Several runs were performed to study the effect of positive and negative sensible heat flux as well as the wind velocity while keeping other parameters constant. Since the temperature and wind profiles differ when the heat flux is positive or negative, the potential temperature and the velocity friction were modified to keep the same wind profile. The run QwB-001 was used as a reference case to compare the effect of the sensible heat flux and wind velocity on the 10% v/v and LFL distance at 0.5 m above the ground. The release rate was set to 5 kg/s. The sensible heat flux was increased up to 250 W/ m2 (positive and negative). This heat flux value seems to be reasonable compared to the heat flux from a sunny day, around 800– 900 W/m2. The effect of the scenario is reported in Table 6. The results showed a small reduction of the LFL distance when the heat flux from the ground was positive. On the other hand, the negative heat flux created larger LFL distance. High concentration levels were observed in the downwind direction for the negative heat flux compared to the base case. Fig. 23 shows the effect of sensible heat flux on the vapor cloud. A negative sensible heat flux increased not only the LFL distance (i.e., 5% v/v contour) but also the concentration up to 25% v/v outside the fence at 0.5 m above the ground level. The positive sensible heat flux decreases both the concentration and the LFL distance. The wind velocity was increased from 3 m/s to 9 m/s. Fig. 24 shows the concentration profile at 0.5 m above the ground for a constant negative sensible heat flux of 250 W/m2. The results in the figure, as well as in Table 6, show the increase of LFL distance as the wind velocity is reduced. The negative sensible heat flux showed a large gradient change compare to positive sensible heat flux. High wind velocity also showed that in both cases, the plume length was diminished. Fig. 24 also shows that with a low wind velocity and negative heat sensible flux, the cloud had high level of concentration and a pronounced wide cloud. 6. Discussion The wind velocity, obstacles, sensible heat flux, and the released mass affected the LFL distance and concentration levels in numerical simulation that were performed in this research. The increase of wind velocity augmented the mixing effect with vapors and ambient air, which resulted in a reduction of the LFL distance. The

Table 6 Summary of test performed on the sensible heat flux.

Fig. 22. Concentration in the downwind direction for 15 m2 LNG pool with a TKE value from 5.8  106 to 5.8 m2/s2.

Run

Sensible heat flux, W/m2

Added Heat Flux, W/m2

Wind velocity at 2.5 m, m/s

LFL distance (4.5, y, 0.5), m

10% v/v distance (4.5, y, 0.5), m

QwB_001 QwP_001 QwP_002 QwP_003 QwP_004 QwN_001 QwN_002 QwN_003 QwN_004

45 0 0 0 0 0 0 0 0

0 150 250 250 250 150 250 250 250

6.8 6.8 6.8 3 9 6.8 6.8 3 9

15 14 12.5 14 13.5 22 23 26.5 18

1 1 0.5 1 None 9 12.5 15 1

350

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

Fig. 23. Concentration contour with constant wind velocity of 6.8 m/s, sensible heat flux, – 250 (bottom right), 50 (top, reference case), and 250 W/m2 (bottom left).

LFL distances were changed by the effects of the obstacles for different wind direction. Negative heat flux increased LFL distance and concentration levels, while the positive heat flux reduced LFL distance and concentration levels. At any level of wind velocity, these parameters influenced LFL distance and concentration levels, but the effects were not linear. The turbulence at the source, pool geometry, and mass evaporation rate parameters seem to be a function of the wind velocity. At higher wind velocities (above 2.5 m/s), these parameters had slight effects either on LFL distance or concentration levels in the downwind direction. At lower wind velocities, the turbulence, pool geometry and mass evaporation rate had a noticeable effect on the downwind distance. LNG pool area was a key parameter because of change of mass evaporation rate and the turbulence effect inside the fence. With open water surface area, the formation of eddies was limited. In general, the circular LNG pool geometric simulation produced a slightly lower LFL distance, but there was no strong evidence to explain these results except that a circular pool will provide more space on the corner of the fence to create larger eddies than with a rectangular shape. The mass evaporation rate had a different effect on the LFL distance. An accurate mass evaporation profile was useful in the estimation of high methane concentration and helped to determine early strategies to dilute the vapor cloud by water curtains. The difference in LFL distance for each evaporation profile was relatively small, but it may be significant for larger releases. TKE had the main effect on the vapor dispersion compared to any other variable in the source term modeling. This parameter

increases the size of the cloud upward, letting the vapors mix with ambient air. This may result in an increase of concentration near the source but decrease in the LFL distance. A Ti value above 100% seems to agree with the height of the cloud from collected data. It was also concluded that large value of TKE increased the width of the cloud, thereby reducing the downwind concentration. A turbulence intensity value of 1000% seems large; however, during the analysis of the data from CFD modeling it was considered possible when the external turbulence is included in the source model. With a Ti of 1000%, the estimated k0 and 30 values were 5.8 and 3.6 m2/s2 at the source for the results presented in Fig. 22. These values were in the same range as the recent CFD simulation results from Gavelli. Gavelli reported a value of 9 m2/s2 for k0 value at the source (calculated from the video feed) for a Falcon series with a wind velocity of 1.7 m/s (Gavelli, Bullister, & Kytomaa, 2008). These two results may justify the use of high value for TKE in the case of LNG release on water. During the course of these simulations, it was obvious that the wind velocity affected the parameters differently. The CFD model was capable of modeling these phenomena correctly. Table 7 shows the classification of importance for each key parameter discussed here. It also provides a range of validity with respect to the wind velocity. For this work, the low wind velocity range was from 0.5 to 2.5 m/s and the high wind velocity range was from 2.5 to 7 m/s where the observations from Table 7 were made. LNG release on concrete showed a lower turbulence intensity compared to LNG release on water. The source term developed in the release on concrete seems to agree with the data collected, which confirmed the mass evaporation profile.

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

351

Fig. 24. Concentration contour with negative sensible heat flux of 250 W/m2, wind velocity from 3 to 9 m/s.

The humidity was modeled as part of the air by changing the thermodynamic properties. Further investigation of the droplet formation is needed. As the infrared camera showed the misconception between the white condensate cloud and the hydrocarbon vapors, it is important to understand the formation of these droplets to be able to predict (i) the formation of the white condensate cloud and (ii) its relationship with the hydrocarbon vapors. Furthermore, the humidity should be studied to determine its effect on the dense behavior.

Table 7 Estimation of the importance for the key parameters involved in LNG vapor dispersion modeling, release on water. Parameters

Importance

Validity range

Release rate

High

Pool shape, geometry Pool area Gas phase velocity profile (derived from the mass evaporation rate) Turbulence at the source Sensible heat flux from ground

Medium High High

At large and low wind velocity At low wind velocity At low wind velocity At low wind velocity

High High

Wind velocity

High

Temperature ambient Surface roughness Wind direction

Medium Medium Medium

At low wind velocity At large and low wind velocity At large and low wind velocity Not tested here Not tested here At large and low wind velocity

7. Conclusions This paper identified the sensitivity of different parameters for LNG vapor dispersion. It was found that some of these parameters depend on the level of wind velocity to be relevant to the simulation. For simulation with a low wind velocity, where the gravity effect was dominant, the parameters involving the source term had a significant influence on the downwind concentration. The parameters which influenced the momentum balance between the wind and the gas phase were considered key parameters. This paper also provides a methodology to assess LNG releases on water or on concrete (i.e., an impoundment). The CFD model used in this paper gave good agreement with the experimental data for wind velocity higher than 2.5 m/s (at 2.3 m above the ground level). This paper provided a qualitative assessment of the parameters on vapor dispersion. The quantitative assessment of the parameters should be studied further. Most of the work presented was focused on LNG releases on water, especially the sensitivity analysis. Further work should be conducted to perform the same type of analysis for LNG releases on concrete and to assess the effect on the exclusion distances. Acknowledgment The authors of this paper would like to acknowledge the Brayton Fire Training Field of Texas A&M University System who provided

352

B.R. Cormier et al. / Journal of Loss Prevention in the Process Industries 22 (2009) 332–352

guidance in the realization and the execution of these tests in a safe manner. The authors also would like to thank BP Global SPU Gas for their financial support and guidance in the realization of these tests.

References Alinot, C., & Masson, C. (2005). k–3 model for the atmospheric boundary layer under various thermal stratifications. Journal of Solar Energy Engineering, 127(2005), 438–443. ANSYS. (2007). CFX-11 solver theory manual. Oxfordshire: CFX Ltd. ASTM D-5523. (1996). Standard guide for selection of Kriging methods in geo-statistical site investigations. West Conshohocken, PA. Briscoe, F., & Shaw, P. (1980). Spread and evaporation of liquid. Progress in Energy and Combustion Science, 6(2), 127–140. Brown, T. C., Cederwall, R. T., Chan, S. T., Ermak, D. L., Koopman, R. P., Lamson, K. C., et al. (1990). Falcon series data report 1987 LNG vapor barrier verification field trials. CCPS. (1999). Consequence analysis of chemical releases. New York: AIChE. DeVaull, G. D., King, J. A., Lantzy, R. J., & Fontaine, D. J. (1995). Understanding atmospheric dispersion of accidental releases. New York. Dharmavaram, S., Hanna, S. R., & Hansen, O. R. (2005). Consequence analysis – using a CFD model for industrial sites. Process Safety Progress, 24(4), 316–327. Ferziger, J. H., & Peric, M. (1999). Computational methods for fluid dynamics. (2nd ed.). Springer.

Gavelli, F., Bullister, E., & Kytomaa, H. (2006). Validation of fluent predictions for LNG spills into an impoundment. In Mary Kay O’Connor Process Safety Center 2006 international symposium. Gavelli, F., Bullister, E., & Kytomaa, H. (2008). Application of CFD (Fluent) to LNG spills into geometrically complex environments. Journal of Hazardous Materials, 159(2008), 158–168. Hanna, S. R., Hansen, O. R., & Dharmavaram, S. (2004). FLACS CFD air quality model performance evaluation with Kit Fox, MUST, Prairie Grass, and EMU observations. Atmospheric Environment, 38(2004), 4675–4687. Hissong, D. W. (2007). Keys to modeling LNG spills on water. Journal of Hazardous Materials, 140(3), 465–477. Luketa-Hanlin, A., Koopman, R. P., & Ermak, D. L. (2007). On the application of computational fluid dynamics codes for liquefied natural gas dispersion. Journal of Hazardous Materials, 140(2007), 504–517. Olvera, H. A., & Choudhuri, A. R. (2006). Numerical simulation of hydrogen dispersion in the vicinity of a cubical building in stable stratified atmospheres. International Journal of Hydrogen Energy, 31(2006), 2356–2369. Panofsky, H. A., & Dutton, J. A. (1984). Atmospheric turbulence. New York: John Wiley & Sons. Pasquill, F. (1974). Atmospheric diffusion. (2nd ed.). New York, NY: John Wiley & Sons. Peterson, R. L. (1990). Effect of homogeneous surface roughness on heavier-than-air gas dispersion. Journal of Wind Engineering and Industrial Aerodynamics, 36(1990), 643–653. Sklavounos, S., & Rigas, F. (2004). Validation of turbulence models in heavy gas dispersion over obstacles. Journal of Hazardous Materials, 108(1–2), 9–20. Sklavounos, S., & Rigas, F. (2005). Fuel gas dispersion under cryogenic release conditions. Energy & Fuels, 2005(19), 2535–2544.