Model for Cryogenic Flashing LNG Leak

Abstract

The growth of liquefied natural gas (LNG)’s importance for curbing greenhouse gas emissions has increased the interest in understanding LNG’s risks, particularly regarding small-diameter leaks (<25 mm). Recently, INERIS performed a series of pressurized LNG release experiments for orifice sizes of up to 9 mm. Based on INERIS findings and the Isenthalpic Homogeneous Equilibrium Model (HEM), this paper created a leak model for flashing LNG leak. The leak model consists of nine equations and quantifies leak parameters for risk assessment. One potential use of this leak model is providing an equivalent leak boundary condition for computational fluid dynamic (CFD) simulation to predict gas dispersion. Using the leak model as input, FLACS gas dispersion simulation was carried out for one INERIS experiment leak case. Compared to TR56, the Singapore safety guideline for LNG bunkering, the dispersion result does not contradict the expected plume reach. Further validation with risk analysis tools and actual experiments is needed to confirm that the leak model is fit for use.

1. Introduction

According to the DNV GL energy outlook 2019, natural gas will become one of the world’s primary energy sources by 2025 [1]. The main reason for natural gas’ growth is greenhouse gas emission-curbing regulations [2]. Natural gas is efficiently stored at cryogenic temperature (approximately −162 °C) with 600 times more density than its gaseous state. Usually used for power generation, liquefied natural gas (LNG) is now seen as the leading alternative source of fuel in the shipping industry [3]. Due to its cryogenic storage temperature, LNG leaks will continuously vaporize into natural gas (NG), leading to a potential fire that can result in fatality and asset loss. A recent liquefied gas leak from a tanker truck in China led to 19 deaths and more than 170 injured [4].

The US Health and Safety Laboratory published a report on the work done thus far to describe LNG spillage (also known as source term) for risk analysis [5]. The LNG source term is categorized broadly as liquid jet, two-phase jet, evaporating pools, and sub-sea release. The factors affecting the source term modeling include the LNG storage condition (pressure and temperature), orifice size and shape, ambient temperature, leak trajectory, atomization of the liquid spray, and the entrainment rate of fresh air. When determining the source term, it is important to deploy techniques to capture the uncertainty of the factors because it can affect the extent of potential consequences. Dorota and colleagues have demonstrated the use of the fuzzy method and Monte Carlo method to improve the accuracy of predicting LNG dispersion [6]. Additionally, their sensitivity analysis indicates that LNG release rate is a major source of uncertainty when defining the source term.

Walter and colleagues summarized the approaches developed to model LNG dispersion from the source term [7]. The developed dispersion models are categorized as the integral model and the Navier–Stokes model. Between the two models, the Navier–Stokes model is preferred because it provides the most complete representation of the fundamental fluid dynamics. However, Navier–Stokes modeling is time consuming and the integral model is used for preliminary risk assessment.

Today LNG accident concern focuses on two-phase pressurized leakage through connection cracks. The specification for the bunkering of liquefied natural gas-fueled vessels (ISO 20519) lists the release of LNG through a broken instrument connection (25 mm hole) as the likely event [8]. A full-bore rupture requires a vessel or road tanker to pull away unintentionally, which is gross negligence. The United Kingdom and Singapore risk guidance regulations treat a crack leak as a major concern for risk rather than a bore rupture [9,10].

When a pressurized cryogenic liquid leaks out into the atmosphere from a crack, it will be in two phases, liquid and vapor. The depressurization process leads to mechanical fragmentation of the liquid and forms vapor [11]. Additionally, the liquid will undergo the rapid vaporization of liquid due to the sudden pressure drop. This rapid vaporization is known as flashing. Conventionally, the non-vaporized liquid, which falls and forms a puddle on the ground, is called rainout [12].

Studies of a non-cryogenic orifice pressurized leak nature have been performed. The Center for Chemical Process Safety (CCPS) studied the release of pressurized water and trichlorofluoromethane (CFC-11) through a small orifice (6.35 mm) [13]. The CCPS study found that significant mechanical fragmentation of the liquid droplet is possible, even when the stagnant liquid temperature is above its ambient boiling temperature. In another study, Hervieu and Veneau measured the droplet size and velocity distribution of liquefied propane leaks through an 8 mm orifice [14]. A Phase Doppler Particle Analyzer (PDPA) was used to measure liquid drop size and velocity. The study concluded that vaporization takes place along the jet in the flow direction. Additionally, higher pressure release is associated with more intense flashing, leading to stronger fragmentation and a higher vaporization rate. The French National Institute for Industrial Environment and Risks (INERIS) performed semi-outdoor large-scale experiments to study propane and butane release through orifices of up to 25 mm [15]. The purpose of the dispersion experiment is to develop models of flashing releases as encountered in realistic industrial environments. Both non-impinged and impinged jet releases were studied by INERIS. Rainouts were collected and measured. The study concluded that mechanical fragmentation and not flashing are the main mechanisms for vapor generation for propane and butane. For non-impinged leaks, there is hardly any rainout for propane due to its low boiling point. Up to 15% of the total mass release was recorded as rainout for the propane impinged leaks. However, the collected rainout is suspected to be a mixture of ice and propane. Hence, the exact propane rainout mass could not be determined. Recently, INERIS performed an experiment named SPARCLING, which releases subcooled LNG up to 8 bar through various orifice sizes (up to 9 mm). The study measures the size distribution and velocities of LNG droplets using a dual Phase Doppler Anemometer [16]. There are two major findings from SPARCLING. The first observation is that an orifice leak <3 mm will always be in two phases up to 9 bar. The second finding is that LNG release above 1.5 bar will not have any rainout.

The work conducted thus far to model flashing leak is for non-cryogenic sources. It would be worthwhile to come up with a model specifically for cryogenic liquids such as LNG. The leak model can be used as the source term for risk assessment of accidental cryogenic liquid leakage. One potential application would be using the model to provide the boundary condition for computational fluid dynamic (CFD) simulation to predict flammable gas dispersion. The current approach to simulate a flashing leak directly is computationally expensive and impractical to use commercially [17].

This paper used SPARCLING findings and the Isenthalpic Homogeneous Equilibrium Model (HEM) to create a leak model for cryogenic flashing leak. Section 2 explains the leak model formulation. Section 3 examines the model applicability by comparing with established LNG safety guidelines. To gauge the leak model adequacy, it was applied in a SPARCLING leak scenario and simulated with CFD software, FLame ACceleration Simulator (FLACS). The dispersion result is then compared with the Singapore safety guideline for LNG bunkering. The last section provides a conclusion and discusses future work to improve the leak model.

2. Materials and Methods

The proposed model tracks the different stages of an LNG leak. Figure 1 shows that the leak starts at the breach, undergoes de-pressurization in the expansion zone, undergoes air mixing in the entrainment zone, and ends at the dispersion zone [18]. In the entrainment zone, LNG that does not vaporize will fall onto the ground to form a rainout. The leak model consists of 9 equations to quantify the leak scenario. Equations (1) and (2) are for the breach orifice zone. Equations (3)–(7) cover the expansion zone leak development. Equations (8) and (9) quantify the leak parameter at the entrainment zone.

The mass flow rate across the orifice for subcooled stagnation condition is provided, as in [19]:

Q = A₁ ∗ C_d √2 × [P₀ − P_V(T₀)] × ρ₀

(1)

Q is the mass flow rate. A₁ is the orifice cross-sectional area. C_d is the discharge coefficient set as 0.62, which is typical for fluid flow across a sharp orifice. P₀ is the stagnation pressure. P_V(T₀) is vapour pressure at stagnation temperature T₀. ρ₀ is the LNG liquid density at stagnation pressure.

The leak velocity, V₁, at the orifice is calculated by:

V₁ = Q/(A₁ × ρ₀)

(2)

In the expansion zone, the stagnation pressure will be depressurized to atmospheric pressure and will undergo volume expansion. If the stagnation temperature is higher than the LNG ambient boiling temperature, a portion of the liquid will be vaporized to reduce the remaining liquid temperature to boiling temperature. The depressurization process occurs rapidly, and heat from nearby air will not be absorbed by the cold LNG. Hence, both the liquid and gas LNG will be at the same temperature (LNG ambient boiling temperature).

To calculate the two-phase composition from the depressurization, HEM was used to determine the vapour mass fraction at the end of the expansion zone [18]. There are three assumptions made in HEM.

• Both the liquid and gas phase temperature will be LNG ambient boiling temperature at the end of the expansion zone;
• The heat exchange from vaporization is much greater than heat gain from surrounding so the warming effect is ignored;
• The vapour mass fraction is much lower than 1.

The vapour mass fraction at the end of the expansion zone is calculated as follows:

X = C_pl × (T₀ − T₁)/L_v

(3)

X is the vapour mass fraction. C_pl is LNG liquid heat capacity at P₀ and T₀ conditions. T₁ is LNG ambient boiling temperature set at 112 K. L_V is LNG latent heat of vaporization.

The two-phase leak density at the end of the expansion zone is given as follows:

ρ₂ = (1 − X) × ρ_L + X × ρ_g

(4)

ρ₂ is the leak density at end of the expansion zone. ρ_L and ρ_g are the respective LNG liquid and gas densities at boiling point temperature in atmospheric pressure.

Based on the conservation of mass and momentum, the condition at the beginning and end of the expansion zone can be expressed, respectively, as follows:

ρ₁ ∗ V₁ ∗ A₁ = ρ₂ × V₂ × A₂

(5)

ρ₁ × V₁² × A₁ + A₁ × P₁ = ρ₂ × V₂² × A₂ + A₂ × P₂

(6)

ρ₁ is the LNG density at the start of the expansion zone. V₁ is the leak velocity calculated from Equation (2). A₁ is the orifice cross-sectional area. ρ₂ is the density of LNG calculated from Equation (4). V₂ is the depressurized velocity at the end of the expansion zone. A₂ is the expanded cross-sectional leak area at the end of the expansion zone. P₂ is atmospheric pressure.

Combing Equations (5) and (6) will form quadratic Equation (7). Equation (7) was solved to find V₂, which was back substituted to Equation (5) to find A₂.

ρ₁ × V₁ × A₁ × ρ₂ × V₂² − (ρ₁ × V₁² × A₁ + A₁ × P₁) × ρ₂ × V₂ + ρ₁ × V₁ × A₁ × P₂ = 0

(7)

According to the SPARCLING experiment conducted by INERIS, pressurized LNG leak beyond 1.5 bar has zero rainout, implying that all liquid is vaporized in the entrainment zone [16]. Given that the vaporization occurs rapidly, it is assumed that the LNG gas temperature remains constant. The leak velocity decreases as it propagates away from the breach. By the conservation of mass, the reduced leak velocity must lead to a proportional increase in the cross-sectional area. The entrainment zone leak area, A₃, is represented as follows:

ρ₂ × V₂ × A₂/ρ_g × V₃ = A₃

(8)

V₃ is the reduced velocity when vaporization is completed. There is insufficient information to determine when vaporization is completed; hence, V₃ is arbitrarily set as ¼ of V₂, and Equation (8) becomes:

4 × ρ₂ × A₂/ρ_g = A₃

(9)

Equations (1)–(9) form the leak model and the calculated parameters provide an equivalent leak boundary condition for CFD simulation. By eliminating the expansion and entrainment zone simulation, which is computationally expensive, the model reduces the CFD simulation time [20].

3. Results

3.1. Leak Scenario

The simulated scenario was chosen from SPARCLING, a 9 mm-diameter leak with an 8 bar initial stagnation pressure. The LNG composition is based on Australian LNG, which has a methane concentration of 99.8%. Given the extremely high methane concentration, LNG was considered as pure methane in the simulation. The assumption of pure methane composition provides a conservative gas dispersion result as heavier alkanes diffuse slower. The leak inventory volume is based on a cryogenic hose of 10 m, 8 inch in length [21]. A 5 s delay to cut off the upstream supply from the emergency shutdown valve was considered, and this increased the total leak volume [22]. The proposed model used the initial pressurized leak parameter at the breach orifice and mapped out the de-pressurized and flashed leak parameter at the entrainment zone.

Table 1. Leak parameter before and after leak model application.

Jet Region	Parameter	Value
Breach orifice	A1	0.0000636 m2
	V1	35.03 m/s
	P1	8.0 bar
	ρ1	418.3 kg/m3
Entrainment zone	A3	0.0239 m2
	V3	21.70 m/s
	P2	1.0 bar
	ρ3	1.797 kg/m3

summarizes the initial leak parameters at the breach orifice and entrainment zone after leak model application.

The breach occurs at the hose joint connection point with the manifold. From a risk perspective, knowing how far a flammable leak source can disperse is critical for safety. According to Singapore risk analysis technical guidance, leak direction, wind direction, and wind speed should be considered when assessing the risk [10].

Three probable leak directions, vertical upward, vertical downward, and horizontal, were considered. The vertical leaks only consider one wind direction, since the two directions are perpendicular to each other. For the horizontal leak, three different wind directions, east (tailwind), west (headwind), and northeast (crosswind), were considered to cover probable wind interaction with the leak. Figure 2 shows the three leak directions and the three wind directions.

According to Nubli and co-workers, a lower wind speed has a slower dilution rate and leads to greater gas dispersion [23]. The wind speed was set at 1 m/s (reference height of 10 m), which belongs to roughly 1 percentile of Singapore’s average wind speed range [24].

3.2. CFD Initial Condition and Output Control

The entrainment zone leak parameter shown in Table 1 served as the leak boundary condition for the gas dispersion simulation. The CFD software utilized was FLACS version 10.7, which has been validated by numerous LNG experiment data [25]. In FLACS simulation, Favre-averaged transport equations for mass, momentum, enthalpy, turbulent kinetic energy, the rate of dissipation of turbulent kinetic energy, mass-fraction of fuel, and mixture-fraction were solved on a structured Cartesian grid using a finite volume method. The standard k-ε model for turbulence was used in FLACS simulation. For more information on the governing equations for fluid flow, please refer to the FLACS-CFD 20.1 User Manual, 2020 [25]. The CFD initial condition was set at an average temperature of 25 °C and 1 bar ambient pressure. Ground roughness was set at 0.03 m, suitable for terrain with few obstacles. The chosen Pasquill class was stable (F), considering the low wind speed. The Courant number (CFLC) based on sound velocity was set at 10, and Courant number (CFLV) based on flow velocity was set at 1.

A grid convergence test was performed at three grid sizes (1 m, 0.5 m, and 0.25 m), and the gas stoichiometric concentration volume (Q8) was tracked to gauge the discretization error. The Q8 difference between consecutive grid sizes was less than 5%, and the 0.25 m grid was chosen for the simulation. The grid around the leak was locally refined as per FLACS best practice guidelines to 0.219 m. The neighbouring cells increased gradually at approximately a 10% rate till it reached 0.25 m. Figure 3 shows the simulation grid setup. The number of controlled volumes for the simulation ranges from roughly 800,000 to 1,100,000.

3.3. Simulation Finding

Natural gas (NG) is flammable between 5% and 15% methane concentration in the air. The 5% concentration is called the lower flammable limit (LFL), and the 15% concentration is called the upper flammable limit (UFL). The LFL distance covered by each simulated scenario is presented in

Table 2. LFL distance of leaks.

Leak Direction	Wind Direction	Furthest LFL Distance (m)
Vertical upward	North	16.72
Vertical downward	North	35.70
Horizontal	East (tailwind)	49.83
Horizontal	West (headwind)	32.53
Horizontal	North-east (crosswind)	47.04

. The LFL distance is the furthest distance from the leak source where the methane concentration has been diluted to 5% in air.

Among the five simulated scenarios, the horizontal leak with tailwind had the furthest LFL distance at 49.83 m, which can be explained by three reasons. The first reason is that it was non-impinged, unlike the vertical downward leak. An impinged leak has the help of obstruction to redirect dispersion flow, leading to a smaller LFL distance. Figure 4a shows the plane view of the downward leak overall spread across the ground. Figure 4b shows the section view of the downward leak at the leak source. The colour depicts the methane concentration at the location. Deep red represents methane concentration at UFL or higher, and deep blue represents the methane concentration at LFL. Concentrations of lower than LFL are not shown in the figure. The gas height measured 0.65 m from the ground, reflecting a low laying gas dispersion, which is a phenomenon observed in past LNG experiments [26,27,28].

The second reason is that horizontal leak experienced lower wind speed as compared to vertical upward leak. The FLACS wind boundary was modelled based on the Monin and Obukhov theory, and the wind velocity increased logarithmically with elevation [29]. Figure 5a shows the section view of wind velocity (north direction) distribution. It can be observed that the wind speed above 10 m elevation is approximately 2 m/s, 10 times more than the wind speed of 0.2 m/s found below 2 m elevation. For the horizontal leak, the gas dispersion barely surpassed 1 m elevation. In the case of the vertical upward leak, the leak reached up to a height of 11 m with an LFL distance of 16.72 m, as shown in Figure 5b. The higher wind speed experienced by the vertical upward leak led to faster gas dilution, resulting in a shorter LFL distance.

The third reason is that a tailwind leak will not suffer “push back” as a headwind would. The opposing headwind reduced the leak velocity, leading to a smaller LFL distance. Figure 6a,b show the section view of tailwind and headwind leak, respectively. The LFL distance was 49.83 m for the tailwind leak and 32.53 m for the headwind.

The LFL distance difference between the horizontal leak with tailwind and crosswind was 2.79 m. The slight difference is likely due to the manifold and hose geometry influencing the wind flow around the leak source. Overall, the five LFL distances differences are explained, and the results do not contradict each other.

3.4. Comparison with Established LNG Safety Guideline

To gauge the adequacy of the proposed leak model to use for risk assessment, the gas dispersion result was compared with the LNG bunkering technical reference, TR56. This technical reference was produced by Enterprise Singapore and recognized by Singapore maritime authorities [30]. TR56 provides an LFL distance vs. stagnation pressure graph, shown in Figure 7, to advise how far flammable leak NG can disperse. The graph is based on a continuous release through a 25 mm-diameter hole. No further information on the surrounding geometry and wind condition was provided. Therefore, it is assumed the TR56 wind condition and geometry are identical to those in Section 3.1 leak scenario.

For an 8 bar 25 mm leak, the worst LFL distance was 60 m, roughly 10 m further than the proposed model result. The LFL distance for a 25 mm hole was expected to be further than a 9 mm hole since the leak area is bigger. However, considering that the cross-sectional area of the 25 mm hole was approximately 7.7 times that of a 9 mm hole, the difference of 10 m LFL distance seems to have been underestimated. One explanation for the minor LFL distance difference could be that a portion of the 25 mm leaked LNG did not vaporize immediately and form a rainout. This means that not all the leak LNG contributed to the immediate gas dispersion. In the case of the 9 mm leak, the SPARCLING experiment showed that no rainout was formed for this leak size. Therefore, the NG dispersion rate of the 25 mm hole will be greater than that of the 9 mm hole but not necessarily proportional to their cross-sectional area difference. At first glance, the proposed leak model prediction is probable.

4. Discussion and Conclusions

The existing modelling work performed for flashing leak is for non-cryogenic pressurized sources, such as propane and butane. These models are not compatible with cryogenic liquids, such as LNG. This paper is an attempt to produce a source term model for cryogenic flashing leak that will be useful for risk assessment, such as LNG gas dispersion. One potential use of the leak model is providing an equivalent leak boundary condition for CFD simulation and helping to reduce computational time. The leak model formulation is based on the recent LNG experiment SPARCLING findings and the Isenthalpic Homogeneous Equilibrium Model (HEM). A comparison of the LNG bunkering technical reference, TR56, with the FLACS simulation using the leak model did not show contradicting results.

To further verify the leak model fidelity, a comparison can be made with commercial software, such as Process Hazard Analysis Software (PHAST) [31]. However, commercial software uses an integral model, which is unable to analyse impinged leaks. Considering the rise in LNG usage in the marine industry, reserving 60 m (refer to Figure 7) for LNG bunkering will be a challenge for port facilities. The quickest way to reduce the LFL distance would be to use a barrier to redirect LNG dispersion. CFD software coupled with the leak model could potentially be used to evaluate the barrier effectiveness to reduce gas dispersion. Hence, it would be worthwhile to conduct an LNG leak experiment to validate the leak model and CFD software. Additionally, the experiment can confirm the model assumption of when vaporization is completed. Finally, the experiment should explore leak orifices of between 10 mm and 25 mm, which is currently lacking in the modern literature.