Liquefied natural gas (LNG) leaks often lead to cascading accident disasters, including vapor cloud release, explosion, fire, and toxic gas release. The origin and evolution of each accidental disaster must be considered when assessing safety. This paper discusses the safety assessment project of an LNG gas storage station in Xuzhou, China. Multiple conceivable disasters due to the leakage of LNG storage tanks are simulated and analyzed using the computational fluid dynamic software FLACS. We studied different wind speeds interacting with the flammable vapor cloud area and creating frostbite in areas of low temperature. Diffusion simulations of vapor cloud explosion (VCE), thermal radiation, and the distribution of toxic substances were performed. The overpressure-impulse criterion was used to calculate the influence range of VCE. Heat flux, heat dose, and heat flux-heat dose criteria were used to calculate the safe distance for personnel in the event of fire. Based on the calculation results of the three latter criteria, this paper recommends using the heat flux criterion to evaluate fire accidents. The danger zone of each accident was compared. VCE accidents yielded the largest area.

1. Introduction

Liquefied natural gas (LNG) is an efficient, clean, and relatively inexpensive energy source and widely used in every daily life, industrial production, and other fields [1]. It has become one of the three pillars of the world’s energy with the volatile international oil price and the rapid growth of the demand for clean energy [2]. Ensuring people’s safety is a paramount concern due to the flammable and explosive nature of LNG, with leakage causing many serious accidents and enormous loss. For example, in 1944, an LNG leakage explosion occurred at an LNG station in Cleveland, USA. The explosion spread over 14 blocks, killed 136 people, and injured 225. The damage extended to 0.75 square miles [3]. In 1997, a spherical tank at an oil refinery in India leaked and exploded, resulting in 25 storage tanks and 19 buildings burned down. More than 60 people were killed, and the damage was nearly $150 million [4]. In 2004, an LNG receiving terminal in Algeria leaked and caused a fire, resulting in 27 deaths, 74 injuries, and over 1 billion U.S. dollars in damages [5]. In 2011, an LPG cylinder leakage accident at a restaurant in Xi’an China resulted in 11 deaths and 31 injuries and also affected 12 shops and severely damaged 53 vehicles with a direct economic loss of 9.9 million Yuan [6].

The diffusion and explosion models and the combustion fireball from LNG leakage are studied by many researchers. Kim et al. [7] proposed a new semianalytical method to estimate the evaporation rate of the LNG liquid pool based on the heat conduction equation and diffusion data. Wang et al. [8] introduced the atmospheric transmission rate τ in the original TNO dynamic model. The optimization model of theoretical predictions fully agrees with the experimental data. Yi et al. [9] established a fire spread model based on the fire dynamics simulator (FDS) and analyzed the fire spread results for a fire scenario in different functional places under the same conditions and the spread of fire outside of the room. The results showed that visibility was a more influential factor in determining the critical time required for fire to become a hazard comparing with CO mass fraction and temperature. Dadashzadeh et al. [10] proposed a methodology for toxicity risk assessment during an LNG fire. It is announced that high concentrations of combustion products and longer exposure time at the process facility increase risks. However, the direct consequence of fire was not considered.

Baalisampang et al. [11] developed an evaluation method for cascade accidents caused by LNG leakage. This method uses computational fluid dynamics (CFD) software FLACS simulation of LNG leakage on a series of disasters and considers the evolving scenario, namely, from a liquid pool after LNG evaporation followed by combustion, explosion, and yielding combustion products. Lv et al. [12] established a reduced model of 160,000 m3 LNG storage tank according to the ratio of 1 : 100, carried out LNG explosion tests, and simulated the explosion process with FLACS. They then established the correlation between the test model and the full-size model explosion overpressure to evaluate the maximum explosion overpressure of large LNG storage tanks. Li et al. [13] used FLACS to simulate the development of a flame, overpressure after leakage, and an explosion of a natural gas pipeline. It was found that the length of natural gas chambers and different ignition positions had little influence on overpressure. This finding supported stronger structural designs and a division of the fire prevention zone. Ren [14] studied the changing law of temperature in the process of large-area LNG leakage and explosion using numerical methods. Moreover, the numerical model considered the influence of wind speed on the temperature of the escaped liquefied natural gas. Hu et al. [15] proposed a topological network model and inversion method for LNG tank fire accidents. They established a topological model for LNG tank fire inversion and found the shortest path, according to the weighted edge topological network structure and to determine the fire location. The location of the fire source obtained by this method was identical to the actual accident location. Bao et al. [16] carried out experimental research on the deflagration of unrestricted methane-air mixture. The impacts of methane concentration and volumetric on the deflagration loads were demonstrated. They also proposed a new model for the prediction of deflagration loads of borderless methane-air mixture. Li et al. [1720] studied the performances of infilled masonry walls under vented gas explosion by filed blast tests and proposed corresponding retrofitting techniques, such as bonding with CFRP, BFRP, and GFRP materials.

The calculated value by the theoretical model often differs significantly from the actual value of the project due to the influences of obstacle distribution, block rate, type of environment, and other site-specific factors. LNG leakage accidents are often accompanied by a series of cascading accidents such as gas leakage, explosion, fire, and flash fire. The accidents that may happen during the LNG storage tank leakage and the corresponding conditions are shown in Table 1. It is necessary to analyze the accident types accurately and establish a refined model for specific cases to obtain more accurate calculation results. Taking an LNG storage station in Xuzhou, China, as the engineering background, this paper studied the disastrous consequences of LNG storage tank leakage and explosion by using the CFD software FLACS and divided the safe distance of personnel by the probability method, providing guidance and reference for engineering construction and disaster prevention.

2. Evaluation Methods

The evaluation method in this study was developed based on the evaluation process reported by Baalisampan et al. [11]. The original process can be basically divided into five steps:(1)The LNG leakage formed a liquid pool. The commercially available software FLACS was used to model and simulate the evaporation of LNG.(2)A series of potential transition phenomena, such as LNG vapor diffusion, pool fires, and steam cloud explosions, were considered. The phenomenological changes in LNG leakage, as well as the effects of atmospheric and thermal radiation, were used to assess potential cascading events. The heat load obtained from the pool fire had to be used as the ignition source of the VCE to meet the minimum ignition energy for the transition to occur.(3)The heat released by the pool fire was assumed as a potential ignition source to ignite and consume the vapor cloud at that location. It is because that VCE [21, 22] might be produced if the ignition delay is between 5 and 10 minutes.(4)The probability function was used to analyze the consequences of fire, explosion, or combustion products and the impacts of thermal radiation on property and personnel.(5)Comprehensive assessment with disasters of fire, explosion, and combustion.

This original process above introduced the progression of LNG leakage accidents, however, with some deficiencies yet. Four of these deficiencies were improved in this study.(1)The storage temperature of LNG is around −160°C, and vaporization absorbs a large amount of heat, with a sharp thermal drop in the surrounding space. Frostbite can be caused [14] when the temperature of human skin tissue drops below the freezing point (273 K). Rescue personnel must take rescue measures within a safe temperature range. Therefore, the impacts of LNG leakage on the ambient temperature were considered in this study.(2)The wind speed can directly affect the distribution of LNG vapor clouds and combustion products, as well as the thermal load of the pool fire, which in turn indirectly affects the explosion load distribution of the vapor cloud. Wind speed is, therefore, an important factor in LNG leakage accidents.(3)The impact of explosion overpressure is insufficient. Impulse is also an important evaluation factor. A comprehensive consideration of overpressure and impulse can more accurately assess the impact of an explosion. The improved evaluation process is shown in Figure 1.(4)We evaluated fire accidents from three aspects of heat flux, heat measurement, and their comprehensive influences, respectively, and compare the calculation results of the three and finally found out the most optimal evaluation index.

3. Mathematical Model

3.1. Basic Assumptions

Obstacles on the ground make the wind field structure nearby irregular, resulting in extraordinarily complex atmospheric turbulences and high uncertainty regarding leakage and diffusion of combustible gasses. We simplified the calculation through the following assumptions [23]:(1)The movement of air-born particles can be regarded as the movement of incompressible fluid when the wind speed in the lower atmosphere is about 10 m/s.(2)Turbulence is random in space and time but can be described statistically. This can be used as the starting point to model turbulence motion. Variables such as wind speed are the sum of their average and disturbance quantities according to the Reynolds hypothesis. This applies to the horizontal component of wind speed .(3)The turbulence viscosity coefficient characterizes the properties of atmospheric turbulence according to the Boussinesq hypothesis. The two-equation turbulence model, namely, k-ε turbulence model, was used to solve the value of .

3.2. Governing Equation

The three-dimensional equations of turbulent diffusion of gas in the atmosphere can be derived as follows.

Mass conservation equation (continuity equation):

Momentum conservation equation:

Energy conservation equation:where T is the temperature (K), p is the gas pressure (Pa), u is the velocity (m/s), ρ is the air density (kg/m3), μ is the turbulent fluid viscosity, and keff is the effective conduction coefficient. For turbulence conditions, the widely used k-ε turbulence model was adopted.

The concentration distribution of continuous leakage can be calculated according to one of the following three formulae based on wind speed u, the relationship between the height of atmospheric mixing layer H0, and the diffusion coefficient of vertical wind direction [2427]:(1)When the wind speed is 1 m/s < u ≤ 1.6 H0, the gas diffusion concentration at a point (x, y, z) in the spatial coordinate system with the leakage source as the origin and the wind direction as the X-axis iswhere C(x, y, z) is the average concentration of the leakage gas at a certain point downwind (x, y, z) (kg/m3); Q is the leakage source strength (kg/s); u is the average wind speed (m/s); H is the effective height of the leakage source (m); , , and are the diffusion parameters in the direction of X, Y, and Z axis, respectively (m); x is the distance downwind (m); y is the wind distance (m); and z is the height above the ground (m).(2)When the wind speed is u < 0.5 m/s, assuming that the steam is evenly distributed around the leakage source in all directions, the concentration c(r) from the leakage source is as follows:where c(r) is the concentration at the distance r from the leakage source (kg/m3); a and b are diffusion coefficients (m); and mΔ is the duration of static wind, , and m is 1, 2, 3, …(3)When the wind speed u is 0.5 m/s < u < 1 m/s, continuous leakage is described by the superposition of instantaneous air mass leakage quantity QΔt in time Δt. The concentration at point (x, y, z) with the leakage source as the origin of the coordinate system and the downwind direction as the X-axis is

4. Numerical Modeling

4.1. Geometric Model

The size and distribution of objects in the actual scene were determined based on the satellite topographic map, as shown in Figure 2(a). Figure 2(c) is the physical map of the LNG storage tank. The numerical model of the LNG storage tank and its surrounding buildings was established by using FLACS, as shown in Figure 2(c), where the wind direction, leakage position, and ignition position in the developed model were also indicated. The monitor point was set every 3 m near the ignition position, designated from 1 to 30.

4.2. Grid Sensitivity Analysis

The mesh size affects the accuracy of the calculation results although FLACS was usually verified in the calculation of microscale model. The results corresponding to mesh sizes can differ more than 10-fold. The sensitivity of the mesh size must be analyzed before the numerical calculation. As the number of grids increases, the calculation accuracy improves, but the calculation scale also increases. Five grid size models of 2.0 m, 1.5 m, 1.25 m, 1.0 m, and 0.5 m were used to calculate the explosion scenario. The results are shown in Figure 3. Overpressure value increases gradually with the decrease in the grid size. The difference between the overpressure values corresponding to the 1.0 m and 0.5 m grids does not exceed 4.3%; it is much smaller than other grid sizes. When the grid size is under 1.0 m, the calculation results barely change. The core diffusion area is set to 1.0 m, and the extended area is set to 2.0 m. The entire area is divided into 1.0 m uniform grids to calculate gas explosion, as shown in Figure 4.

4.3. Parameter Determination

The LNG storage tank volume was 157.9 m3, the initial inside pressure 0.6 MPa, and the LNG temperature was −162°C. The atmospheric pressure was assumed as 1.01 × 105 Pa, the atmospheric temperature was set to 20°C, and the wind speed was 0 m/s, 4 m/s, and 8 m/s, respectively. This allowed the discussion of the influence of wind speed on the gas diffusion. The relative turbulence intensity was 0.1, and the turbulence length scale was 0.01 m. The leakage hole diameter of LNG storage tank was 100 mm, and the leakage hole was assumed circular. The leak rate was calculated using the leak module of FLACS with an initial leak rate of 38.25 kg/s.

We used the preference wind to set up the wind outflow boundary and the others by using nozzle, as provided in Table 2.

The gas diffusion velocity was related to the surface roughness Z0. Higher surface roughness will yield slower gas diffusion. The ground roughness was set to 1.0 according to the suggested standard value provided in Table 3 [28].

The atmospheric stability was determined by factors like wind speed, radiation amount, and cloud cover. According to the Pasquill–Gifford classification standard of atmospheric stability in Table 4, the atmospheric stability corresponding to the four wind speeds is B, C, and D, respectively. The initial condition parameters are listed in Table 5.

5. Results and Discussion

5.1. Analysis of Diffusion Results

LNG leakage and diffusion were simulated under three wind conditions with a wind speed of 0 m/s, 4 m/s, and 8 m/s. Figure 5 shows the flammable zone of LNG under each wind condition at the time of 190 s, where the methane flammable concentration was between 5.0% and 15.0%. The combustible area was 12409 m2, 5097 m2, and 870 m2. Clouds travel faster and spread more widely with the wind. However, wind also hastens the dilution of fuel in the air. Explosions are prevented when the methane concentration falls below 5%. Therefore, the presence of actually wind reduces the LNG diffusion risk; the higher the wind speed is, the smaller the combustible area is.

The combustible area reaches a dynamic equilibrium and remains unchanged with time when the fuel vapor cloud speed equals the lost speed. The combustible product has reached its dynamic equilibrium at wind speeds of 4 m/s and 8 m/s at the time of 190 s, while when  = 0 m/s, the combustible zone is still progressing even after 500 s, as shown in Figure 6. It is found that the wind speed is a key factor in the dilution rate, diffusion rate, and influence range of the combustible region. Wind turbines can be added to speed up dilution of the steam cloud and reduce the risk of explosion in addition to controlling the leakage source to prevent LNG leakage-related disasters.

The evaporation of LNG causes the temperature of the surrounding space to drop rapidly, so rescuers must pay attention to the potential of frostbite during emergency repairs. Although the influence range of low temperature is smaller than that of diffusion, its impact lasts longer. Once a leakage accident occurs, it is accompanied by low-temperature injury. The occurrence of fire and VCE requires certain ignition conditions to be presented. The risk of frostbite is higher than that of fire and VCE. Figure 7 shows a diagram of the temperature distribution on-site under different wind speeds at 190 s of leakage. In the absence of wind, diffusion of LNG is slow, and the fuel is rather concentrated around the storage tank, which decreases the temperature in the whole site rapidly. The area of frostbite (T ≤ 273 K) is significantly larger than that at wind speeds of 4 m/s and 8 m/s. The temperature in the lower air inlet decreases in the whole station when the wind speed is 8 m/s, leading to an area of frostbite up to 7 m around the leaking tank. Wind can effectively accelerate the diffusion of LNG and reduce the frostbite area and the duration of low temperature.

5.2. Risk Assessment of Fire Accident

The commonly used thermal radiation assessment criteria include the thermal dose, the heat flux, and the heat flux-heat dose criteria. The heat flux criterion stipulates that the critical heat flux value of possible burns for personnel is 4 kW/m2. The critical heat dose for a one-time burn affecting personnel is 125 kJ/m2 [23]. The criterion holds that whether the target can be harmed is not determined by either heat dose or heat flux alone but in conjunction. If heat dose and heat flux are used as abscissa and ordinate, respectively, the critical state of target damage corresponds to a critical curve, as shown in Figure 8. The lower left of the curve is the safe area; the upper left is the harmful zone. The asymptotes are q = qcr and Q=Qcr, respectively, corresponding to critical heat flux and critical heat dose.

Based on the burns of exposed skin, Lees [29] obtained the critical heat flux-time curve of exposed skin under the action of thermal radiation, as shown in Figure 9. Unbearable pain is experienced when the temperature below the skin at a depth of 0.1 mm exceeds 44.80°C. Pain increases sharply and then gradually subsides until it completely disappears with increasing temperature, indicating that the exposed skin has been completely burned out. At a less than 40 s exposure time, the curve is an inclined line with a constant heat dose between 60 and 70 kJ/m2, which is consistent with the recommended value of 65 kJ/m2 in the literature [30]. Beyond an exposure time of 120 s, the curve is basically a horizontal line, and the corresponding heat flux is constant at 1.4 to 1.5 kW/m2. When the heat flux is less than or equal to 1.4 kW/m2, there is no tissue damage regardless of the exposure time. Normal blood flow reduces the possibility for the local temperature to reach or exceed 44.80°C.

The fire accident started at 300 s after the leakage of LNG, and it can be known from the simulation results that the influence range of heat radiation reached the maximum after 25 s of flame combustion, as shown in Figure 10, and the radius of the danger is 43.82 m. In addition, the heat dose value was measured through the monitor points, and Figure 11 shows the heat dose values of points 21 to 30.

Buettner provided the empirical formula, assuming that 20% of the skin is exposed. Pietersen proposed the current common probability equation of thermal radiation injury in 1990 [3133].

The probability of death when the skin is exposed is

Probability of death with clothing protection (20% skin exposed) is

Probability of second-degree burn (serious injury) with clothing protection (20% skin exposed) is

Probability of first-degree burns (minor injuries) with clothing protection (20% skin exposed) iswhere I is the heat flux received by the human body (w/m2); t is the exposure time (s); and Pr is the damage probability variable. When Pr = 5, the corresponding damage percentage is 50%.

We used the principle of exposure time of 10 s and 50% probability and determined the radius R1 of death, the radius R2 of serious injury, and the radius R3 of minor injury, as listed in Table 6. The influence radius of the heat flux-heat dose criterion is 45 m, the influence radius of the heat flux criterion is similar, and the influence radius of the heat dose criterion is only 30.4 m, as shown in Figure 12. The calculation of heat flux-heat dose criterion is more complicated than the heat flux criterion. Therefore, we propose to use the heat flux criterion to evaluate the heat of combustion in the fire and explosion scenarios.

5.3. VCE Analysis and Risk Assessment

The range of combustible zone can be determined with the results of leakage. The simulated 3D cloud image of the VCE is shown in Figure 13.

In the 1970s, the United States ballistic research laboratory (BRL) and the naval weapons laboratory (NOL) proposed a set of damage criteria of blast wave, namely, the overpressure-impulse criteria [31]. They state that the overpressure and impulse of the detonation wave together determine whether and to what extent the target is damaged. If different combinations of overpressure and impulse satisfy certain conditions, they can cause the same degree of damage to the target. In this paper, the damage range of VCE was evaluated by using this criterion. The ordinate in Figure 14 [23] is the peak overpressure of the propagated blast wave, and the abscissa is the respective impulse. The curve is the equal-damage curve of the target, where the lines and are the asymptote of the equal-damage curve and and represent the critical impulse and critical overpressure of the target damage, respectively. The damage-free zone is to the left of the vertical asymptote and below the horizontal asymptote, while the damage zone is to the right of the vertical asymptote and above the horizontal asymptote. Based on the damage criterion, Yan [34] fitted the curves of severe injury, moderate injury, and minor injury according to the biological test results, and the corresponding calculation formula is as follows:where C is a constant value and C values corresponding to serious injury, injury, and minor injury are 4.9 × 104, 2.3 × 104, and 1.0 × 104, respectively. Critical overpressure is , and critical impulse is .

The most vulnerable nonfatal organ of human body impacted by a blast wave is the ear, and the most vulnerable fatal organ is the lung [35]. The consequences of ear and lung injury can be analyzed based on the blast wave. For the probability variables of death from eardrum rupture and pulmonary hemorrhage, the formula which takes impulse into account was used [36]:where is the peak overpressure, ; is the damage probability variable; I is the peak impulse, ; and is the effective overpressure, which is determined jointly with the body posture . The calculation formula is as follows [31]:

Possible casualties can be projected for a steam cloud explosion by dividing the injury area into dead, serious, and minor area. The classification of explosion-related injuries is shown in Table 7.

The radius of eardrum rupture of 1% of the personnel is used as the dangerous distance of VCE to guide the reinforcement of key parts and the emergency evacuation of personnel. The calculated personnel damage radius is shown in Table 8. With the explosion point as the center, the dead areas, severely injured areas, and lightly injured areas were indicated, respectively, in the storage station, as shown in Figure 15. It can be found that almost the entire site is in the danger zone, and its influence range is greater than steam cloud diffusion and fire accidents.

5.4. Combustion Product Impact

The main component of LNG is methane, and the main toxic gas produced in the event of VCE and fire accident is carbon monoxide, which will kill people within 2-3 minutes if its concentration reaches 1.28% [37]. Compared with VCE and thermal radiation accident, the carbon monoxide impact lasts longer and is not easily detected by human body, which is a more dangerous accident. The distribution of carbon monoxide in VCE accident was simulated. The distribution is shown in Figure 16. It can be seen that carbon monoxide can affect the area outside the station, the office building in the downwind is vulnerable, and the exit 2 is covered by the toxic. The evacuation route should be determined according to the wind direction when organizing personnel to evacuate.

5.5. Scope of Each Disaster

Using comprehensive comparison VCE, toxic mole fraction, vapor cloud diffusion, heat radiation, and frostbite, we can see from Figure 17 that the VCE’s sphere of influence is the largest of these five harms in this type of accident, followed by toxic mole fraction and vapor cloud diffusion, and the frostbite area the smallest. The impact of the VCE is the primary factor in formulating preventive measures and emergency plans.

6. Summary

LNG leakage will lead to a series of cascading accidents, such as cryogenic frostbite caused by factors such as LNG evaporation, steam cloud diffusion, explosion, fire, and toxic gas. The cascading accidents often have greater impact than the single disaster, so the assessment of LNG leakage requires a comprehensive analysis of all types of harmful events. We further improved the LNG leak assessment method proposed by Baalisampan et al. Considering the wind speed, cryogenic frostbite, and the new personnel injury assessment method, the safety assessment method was applied to an LNG storage station in Xuzhou, China. Some conclusions are drawn as follows:(1)Wind can effectively accelerate the spread of LNG and reduce the areas of combustion and frostbite. A higher wind speed leads to a more obvious reduction. When the wind speed is 4 m/s, the combustible area reduces by 58.9% compared with that of 0 m/s, while when the wind speed reaches 8 m/s, the combustible area reduces by 93.0%, which shows that the wind speed is an important factor affecting diffusion.(2)Three criteria, heat flux, heat dose, and heat flux-heat dose, were employed to calculate the distances for personal injury in a fire. The range determined by the heat dose criterion is dangerous. The heat flux and the heat flux-heat dose criterion yield similar distances. Given the tedious calculation of the latter and the equivalent result, the heat flux method was recommended to use for the evaluation.(3)The overpressure-impulse criterion was employed and combined with the probability assessments. The radius in which 1% of the personnel experience eardrum rupture is denoted as the VCE dangerous distance.(4)The influence range of VCE, toxic mole fraction, vapor cloud diffusion, heat radiation, and frostbite were comprehensively compared. The influence range of VCE was the largest of all the cascade accidents.

Data Availability

No additional data are available.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors acknowledge the financial support from the National Natural Science Foundation of China (Grant no. 51978166 and 51738011), the Fundamental Research Funds for the Central Universities, and Natural Science Foundation of Jiangsu Province of China (Grant no. BK20180081).