Study on Hazardous Areas of Hydrogen Fluoride Diffusion Based on CFD Simulation

Hydrogen fluoride (HF) is a typical dangerous gas in the fluorine chemical industry. Its leakage is one of the most common types of accidents in this industry, and it poses a serious threat to personnel safety and health, environmental sanitation, and social stability. In this paper, the process and consequences of an HF leakage accident in a fluorine chemical plant were simulated by using computational fluid dynamics (CFD) simulation software, and hazardous areas (lethal area, severe injury area, light injury area, and maximum allowable concentration area) of HF diffusion were determined according to the HF concentration corresponding to the degree of personal injury. Moreover, the effects of wind speed and height on hazardous areas were analyzed. The research results of this paper provide model support for similar enterprises to predict the consequences of harmful gas leakage accidents, and give suggestions on emergency evacuation and rescue work, which have practical application significance.


Introduction
Hazardous gas leakage and diffusion accidents have a wide range and can cause large losses [1][2][3][4]. Hydrogen fluoride (HF) is a typical dangerous gas in the fluorine chemical industry. There have been many HF leakage accidents. In 2009, in the United States, thousands of people were evacuated urgently because of a hydrofluoric acid leakage in a fluorine chemical plant. In 2012, a HF leakage accident occurred in Gumi, South Korea, which caused necrosis of vegetation within 1 km and toxic gas staying for more than ten days [5]. In 2018, three deaths in China occurred due to HF poisoning.
Considering the cost and uncontrollability of harmful gas leakage and diffusion experiments, the current main research method in this field is computational fluid dynamics (CFD) simulation [20][21][22][23]. In 2011, Tauseef et al. [24] used data from the 26th experiment of Thorney Island in the United States to study the accuracy of various CFD turbulence models. The research results showed that the gas diffusion calculated by the realizable kε model was the closest to the actual situation. The realizable k-ε model can not only capture the gravity drop related to the diffusion of heavy gases, but also clearly describe the gas movement in the presence of obstacles. In 2016, Fiates et al. [25] conducted a leakage and diffusion experiment of carbon dioxide. Then, they used CFD software to reproduce the experiment, and the simulation results were in good agreement with the experimental results. In 2017, Yang et al. [5] used CFD software to recover an HF leakage accident in South Korea in 2012 and studied the environmental impact of the accident. The above studies show that the CFD simulation method is reliable and can accurately describe the phenomenon of harmful gas diffusion in various scenarios. In order to analyze the serious consequences and influencing factors of harmful gas diffusion, in 2013, Wang et al. [26] explored the diffusion law of oil gas after leakage through experiments and simulations, and set the area where the oil gas concentration exceeded the minimum concentration of strong explosion as a hazardous area. The effect of the leakage flow rate on the hazardous area was also discussed in the research. In 2014, Dadashzadeh et al. [27] simulated the diffusion of carbon monoxide, nitrogen dioxide, and methane in a specific area and divided the risk area for a single person according to the highest concentration of different gases that do not cause harm to the human body. In 2019, Souza et al. [28] proposed a method for quickly estimating the extension volume of hazardous areas based on CFD software and studied the effects of transport properties, orifice size, gas molar mass, temperature, pressure, and lower explosion limit on the extension of hazardous areas. In 2020, Barros et al. [29] discussed the effects of wind direction and speed on the volume of gas clouds based on the results of CFD calculations, and predicted the hazard level of the area according to relevant international standards. Overall, the harmful gas diffusion phenomenon has been widely studied, but little attention has been paid to the specific harmful effects of the gases on the human body. Moreover, the changes in the gas hazardous area are affected by many factors. The distribution and changing laws of various hazardous areas at different heights are still unclear.
To improve understanding of how the harmful effects of HF on the human body are reflected in the hazardous area, in this paper, the method of dividing the hazardous area of the fluorine chemical plant according to the maximum allowable concentration in the HF occupational exposure scenarios and the concentrations that cause light, severe, and lethal human injuries is proposed. In addition, in order to explore the influencing factors more comprehensively, the changing law of various hazardous areas at different heights is investigated through CFD software calculation and visualization methods. The results can provide a technical basis for the research and judgment of the risk situation of harmful gases leakage accidents in fluorine chemical enterprises. Safety workers can refer to the prediction model to deploy emergency resources, implement emergency rescue measures, and plan escape routes and safe areas.

Modeling and Verification
To verify the reliability of the CFD simulation method, Fluent was used in this study to restore the Goldfish HF leakage and diffusion experiment.
where ρ is the density, t is the time, u is the velocity, p is the pressure, τ is the shear stress, g is the gravity acceleration, v c and p c are the specific heats, T is the temperature, and T k is the thermal conductivity.
The gas diffusion process is complicated. A suitable turbulence model can more accurately calculate the actual situation. Fluent provides many turbulence models, such as the standard k-ε model [31], realizable k-ε model [32], and Reynolds stress model [33]. Some studies have shown that the standard k-ε model is more suitable for gas flow scenarios outside complex buildings [30,34]. Therefore, this paper chooses the standard k-ε model, which is shown in Equations (4) and (5) where k is the turbulent kinetic energy, ε is the rate of dissipation, 1 C ε , 2 C ε and C μ are constants, 1.44, 1.92, and 0.09, respectively, and k σ and ε σ are the turbulent Prandtl numbers for k and ε, 1.0 and 1.3, respectively.

Model Validation
Considering the size of the computational domain of the simulation in this study, the HF concentration data of the Goldfish 2, 3 experiments at 300 m from the leak source [5] were selected in this paper for reliability verification.
In the Goldfish 2 and 3 experiments, HF was stored in a tank in liquid form. The density of the fluid was 956.50 kg/m 3 , and the gauge pressure of the tank was 6.67 × 10 5 Pa. The liquefied HF transformed into a gaseous state at the moment of leakage. The research by Yang et al. [5] showed that there was not much difference between the results of single-phase and two-phase simulations. Therefore, the single-phase simulation method was adopted in this paper. The leakage direction was consistent with the wind direction, the atmospheric stability grade was D, and the ground surface roughness was 0.003. Other relevant information and data are shown in Table 1 [19]. The final calculation domain had a length of 360 m, a width of 50 m, and a height of 30 m. The leaking storage tank was simplified as a cylinder with a diameter of 2 m and a height of 2 m. After the grids were divided in the ICEM software, the total number of grids obtained was 445,304. Fluent calculation results showed that the corresponding concentrations of GF2 and 3 were 2.3587 × 10 −5 kg/m 3 and 2.2021 × 10 −5 kg/m 3 , respectively. The comparison with the actual experimental measurement values is shown in Table 2. The relative error of the data is within the acceptable range. Therefore, the simulation method can be applied to divide and analyze the hazardous areas of a certain plant area after an HF leakage accident.

Simulation Application Background
In this study, a fluorine chemical plant in Quzhou City, Zhejiang Province is selected for numerical simulation of the HF leakage accident. According to the wind direction and wind speed data of Quzhou in 2019 and 2020 published by the China Meteorological Administration, the most common wind direction in this area was the northeast wind. In 2019, it accounted for 65.48% of the annual days, and in 2020, it accounted for 39.08% of the annual days. The most common wind speed was 1.6-5.4 m/s. In 2019, it accounted for 35.62% of the annual days, and in 2020, it accounted for 60% of the annual days. Therefore, the simulated wind direction in this study is the northeast wind, and the simulated wind speed is 0, 1, 2, 4, and 6 m/s, respectively. The storage conditions of HF are the same as the Goldfish 2, 3 experiments, and the simulation is still single-phase. According to the actual situation of the local environment, the model parameters are revised. This study establishes a HF leakage model with a length of 570 m, a width of 400 m, and a height of 100 m. The total number of grids is 823,302.
Based on the human injuries caused by HF [36,37], the diffusion area is divided into four categories: lethal area, severe injury area, light injury area, and maximum allowable concentration area (MAC area). The description of these areas is listed in Table 3.

Diffusion Law of HF under Different Wind Speeds
This study calculates the HF leakage diffusion process when the wind speed is 0, 1, 2, 4, and 6 m/s. In the calculation results, the phenomena at wind speeds of 1 m/s and 4 m/s are very typical. Figures 1 and 2 show the three-dimensional concentration distribution changes of the HF diffusion process under these two conditions. As shown in Figure  1, when the wind speed is 1 m/s, HF diffuses in all directions after leaking, forming a planar influence zone. The closer to the leak hole, the higher the HF concentration. The closer to the bottom of the buildings, the higher the HF concentration. When the wind speed is 4 m/s, the three-dimensional concentration distribution changes of HF are shown in Figure 2. After HF leaks, it spreads to the downward wind direction and forms a bandshaped influence zone at 420 s. The closer to the leak hole, the higher the HF concentration. However, compared with the wind speed of 1 m/s, HF no longer accumulates at the bottom of buildings near the leak hole, but spreads far away. To accurately understand the harmful effects of HF on humans, this study simulates the diffusion of HF at a height of 1.8 m (1.8 m is the most common height for workers in the fluorine chemical plant) under different wind speeds. The results are shown in Figure  3. In a no-wind environment, the leaked HF spreads more evenly in all directions. The denser the buildings, the slower the dispersion rate. After 240 s, the distance that HF disperses in the z-axis direction is greater than the distance that it diffuses in the x-axis direction. When the wind speed is 1-2 m/s, the HF diffusion law is less affected by the wind, which is similar to the outcome of the no-wind environment. Due to the blocking effect of dense buildings, HF tends to disperse in the upward wind direction after 240 s. When the wind speed is 4-6 m/s, the wind speed has a greater effect on the diffusion of HF. HF forms a band-shaped area in the downwind dispersion process, and the greater the wind speed, the narrower the band area.

The Changing Law of Hazardous Areas at Different Heights
To analyze the effect of height on hazardous areas, the changes in four kinds of hazardous areas at heights of 0 m, 1.8 m, 5 m, 10 m, 15 m, and 20 m when the wind speed is 2 m/s are calculated, respectively. The calculation results are shown in Figure 4.  When the height is less than 5 m, the change in the severe and light injury areas is like that of the lethal areas. When the height is more than 5 m, their areas fluctuate greatly, but the changing trend also increases first and then decreases with an increase in the height. Figure 4d shows that when the time is 600 s, the MAC areas at heights of 0 m, 1.8 m, 5 m, 10 m, 15 m, and 20 m are 7343 m 2 , 7751 m 2 , 8219 m 2 , 8668 m 2 , 11,252 m 2 , and 22,806 m 2 . When the height is less than 5 m, the changing trend of the MAC areas is like that of the other kinds of hazardous areas. When the height is more than 5 m, the MAC areas will increase significantly as the height increases.

The Changing Law of Hazardous Areas under Different Wind Speeds
To study the effect of wind speed on hazardous areas, the changes of four kinds of hazardous areas at the height of 1.8 m when the wind speed is 0 m/s, 1 m/s, 2 m/s, 4 m/s, and 6 m/s are, respectively, analyzed. The results are shown in Figure 5. When the wind speed is 0-2 m/s, the growth rate of the four hazardous areas is relatively flat, which is due to the uniform diffusion of HF in all directions after leakage at low wind speeds. The reason for the small fluctuation is that the dense buildings on the dispersion path block the HF. When the wind speed is 4-6 m/s, four hazardous areas fluctuate greatly. This is because HF spreads in the downward wind direction in a band shape at high wind speeds, which can be greatly affected by the wind. Figure 5a shows that in 0-150 s, the greater the wind speed, the larger the lethal area. This is because, at the initial stage of diffusion, the wind transmits a large amount of highconcentration HF to the fluorine chemical plant area. The greater the wind speed, the faster the transmission speed. In 150-500 s, the greater the wind speed, the smaller the lethal area. When the wind speed is 0-2 m/s, the largest lethal area is 887.13% larger than that when the wind speed is 4-6 m/s. The reason for this difference is that in the middle and late periods of dispersion, the greater the wind speed, the greater the dilution effect on HF. After 500 s, due to the blocking effect of dense buildings on the dispersion path, the size relationship of the lethal area under different wind speeds is unstable.
According to Figure 5b-d, the greater the wind speed, the larger the severe injury area, the light injury area, and the MAC area. This is because when the wind speed is 0-2 m/s, the HF spreads uniformly in all directions in a planar shape, and the main body of the planar shape is the lethal area. The severe injury, light injury, and MAC areas are ring-shaped surrounding the lethal area. When the wind speed is 4-6 m/s, the lethal concentration of HF forms a band-shaped area. The greater the wind speed, the narrower the band-shaped area, the smaller the lethal area, and the larger the severe injury area, light injury area, and MAC area.
When the wind speed is 4-6 m/s, the changes in the severe injury area, the light injury area, and the MAC area are extremely unstable. This is because greater wind speeds have greater disturbance to HF diffusion, and dense buildings on the dispersion path have a blocking effect on it.

Conclusions
In this paper, CFD simulation software was used to predict the process and consequences of an HF leakage accident in the fluorine chemical plant. The diffusion law of HF under different wind speeds was studied. After fully considering the injury effect of HF on the human body, the affected area was divided into lethal area, severe injury area, light injury area, and MAC area. The effects of wind speed and height on hazardous areas were discussed. The main conclusions can be drawn as follows: • When the wind speed is less than 2 m/s, HF diffuses in all directions in a planar shape. When the wind speed is greater than 4 m/s, HF spreads in a downward wind direction in a band. The greater the wind speed, the greater the effect on HF diffusion.

•
When the wind speed is 2 m/s and the height is less than 5 m, the change in hazardous areas is almost independent of the height. At 600 s, the lethal area fluctuates within 43,751-45,211 m 2 , the severe injury area fluctuates within 5088-5458 m 2 , the light injury area fluctuates within 2065-2514 m 2 , and the MAC area fluctuates within 7343-8219 m 2 . When the height is greater than 5 m, as the height increases, the lethal area gradually decreases, the MAC area gradually increases, and the severe and light injury areas first increase and then decrease.

•
When the height is 1.8 m and the wind speed is 0-2 m/s, the growth rate of hazardous areas is flat. When the wind speed is 4-6 m/s, the hazardous areas fluctuate greatly. The reason for this difference is that at low wind speeds, HF spreads uniformly in all directions. However, under high wind speeds, HF diffuses in a band shape. The severe injury area, light injury area, and MAC area increase with the increase in wind speed. For the lethal area, the area increases with the increase in wind speed in the range of 0-150 s and decreases with the increase in wind speed in the range of 150-500 s. Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.