Abstract
Currently, there is very limited understanding of a gas explosion process inside residential buildings. In this study, a numerical model of gas explosion in a residential building was developed using Computational Fluid Dynamics (CFD). The numerical simulations were performed for different gas cloud filling regions and equivalence ratios to identify the initial scenario, and the simulation results were compared with the real consequences of gas explosion. Additionally, the temporal and spatial evolvement characteristics of explosion overpressure and indoor temperature were analyzed. Furthermore, the effects of vent area ratio and the activation pressure of vent panels in the kitchen were investigated to propose effective mitigation measures for the gas explosions inside residential buildings. The results show that the simulation results reproduced by the CFD model are in good agreement with the real accident consequences. During the explosion process, the overpressure distribution in a room is almost uniform at the same moment and there exists little spatial difference. The maximum temperature can reach up to 1953°C, which can cause secondary fire accidents easily. The maximum flame speed is in the range of 34.3 m/s and 230.9 m/s. It indicates that gas explosion inside residential buildings is a typical deflagration process. When the vent area ratio is less than 0.3, the overpressure peaks decrease rapidly with the increase of the vent area ratio. However, when the vent area ratio is larger than 0.3, the overpressure peaks are almost independent on the vent area ratio. There is a proportional relationship between the overpressure peaks and the activation pressure of vent panels. These achievements provide reliable reference data for the accident investigation of gas explosion and subsequent treatment.
1. Introduction
Natural gas is becoming an important energy source option and its world production capacity keeps increasing rapidly [1]. Due to the insufficiency of effective safety management for gas companies and safety awareness for residential users, gas leakage and explosion accidents have occurred frequently inside residential buildings in recent years. It has been reported that there were 531 gas explosion accidents that occurred within residential buildings in China in 2016 [2], causing significant property loss and human injuries. It is of great importance to identify the initial accident scenarios, analyze the dynamic characteristics of gas explosion, and propose effective mitigation measures for gas explosion and subsequent treatment inside residential buildings.
Currently, many scholars have carried out some theoretical and experimental studies on building damage under explosive load [3, 4]. There are mainly three theoretical explosion models, including TNT equivalent, TNO multienergy method, and Baker-Strehlow method [5]. However, these three explosion models mentioned above are only suitable for evaluating the explosion effects in the far field and cannot be used to predict the temporal and spatial evolvement characteristics of explosion overpressure and temperature in the near field [6]. Thus, these three methods cannot be directly applied to investigate gas explosion accidents within residential buildings and evaluate explosion consequences.
With the development of computer techniques and finite element analysis software, numerical simulation based on Computational Fluid Dynamics (CFD) has become one of the major methods to study gas explosion process in confined spaces. For instance, Wang et al. [7] performed CFD simulations to investigate the effects of ignition location and vent size on the gas explosion process in a vented room. Huser et al. [8] established CFD models to analyze the effects of gas cloud size, gas cloud location and shape, ignition location and vent area fraction on explosion overpressure in process plants. Ferrara et al. [9] developed a two-dimensional axisymmetric CFD model to analyze the effects of ignition position, duct diameter, and length. Cen et al. [10] performed CFD simulations to investigate the effects of ignition location, vent area ratio, vent activation pressure, the mass per unit area of panels, and volume blockage ratio on the peak pressures in compressor compartments. Karnesky et al. [11] studied the effects of vent area and obstacles on explosion overpressure in vented enclosures based on a three-dimensional gas dynamic model with constant burning velocity. Pedersen et al. [12] established CFD models to investigate the effects of vent size and configuration and flow interactions between two enclosures on the gas explosion venting process in a twin-compartment enclosure. These CFD models can provide more detailed and accurate results of gas explosion process in the near field than theoretical and experimental studies alone. Unfortunately, there is a paucity of literature for applying CFD models to analyze the dynamic characteristics of gas explosion inside residential buildings and its effective mitigation measures.
The objective of this research is to study the dynamic characteristics of gas explosion and propose corresponding proper mitigation measures inside residential buildings based on CFD models. The numerical simulations are performed for different gas cloud filling regions and equivalence ratios to identify the initial scenario. Meanwhile, the simulation results are compared with the explosion consequences associated with a real accident. In addition, the temporal and spatial evolvement characteristics of explosion overpressure and indoor temperature are analyzed. Furthermore, the effects of vent area ratio and the activation pressure of vent panels in the kitchen are investigated to propose effective mitigation measures for gas explosions inside residential buildings.
2. Accident Description
At about 11:10 am on March , 2016, a natural gas explosion accident occurred in Apartment 101 on the first floor inside a residential building in Xindu District in Chengdu, Sichuan, China. The location and floor layout of Apartment 101 are shown in Figure 1. In this accident, one person got just slightly injured, while the doors, windows, ceilings, and walls were damaged seriously, causing a direct economic loss of more than 1 million Chinese yuan.
As shown in Figure 2, there were different damage degrees for different rooms inside Apartment 101. The study room and balcony were damaged most seriously, in which ceilings were blown away and walls were severely cracked. However, the kitchen and bathroom were slightly damaged, in which walls and ceilings still remained intact, only the equipment hanging on interior walls fell, and the window frames were bent out of shape. The ceilings in the living room and the second bedroom were also blown away and appeared cracks in both rooms. Furthermore, the secondary fire accident was observed in the master bedroom and there were visible black smoke traces on interior walls.
(a) Study room
(b) Balcony
(c) Kitchen
(d) Bathroom
(e) Living room
(f) Second bedroom
(g) Dining room
(h) Master bedroom
According to the accident investigation conducted later, the investigating organization generally concluded that natural gas had leaked gradually from the hose pipe in the kitchen during the whole night prior to ignition and subsequently diffused to other rooms. When the householder came in from the master bedroom to the living room and turned on the lighting switch in the next morning, a natural gas explosion occurred immediately. Thus, the ignition source that caused this explosion accident is confirmed to be the lighting switch in the living room.
3. Numerical Model
3.1. Governing Equations
The explosion venting of combustible gas in confined spaces is a complex and transient combustion process, which should satisfy the conservation equations for mass, momentum, enthalpy, and mass fraction of chemical species [13, 14].
Mass conservation equation:where is the volume porosity; is the density, kg/m3; is the area porosity in the jth direction; is the velocity component in the jth direction, m/s.
Momentum conservation equation:where, is the pressure, Pa; is the stress tensor, N/m2; is the gravitational acceleration in the ith direction, m/s2; is the flow resistance due to walls, N; and is the flow resistance in the xi direction due to subgrid obstructions, N.
The flow resistance in the direction due to subgrid obstructions in (2) can be expressed aswhere the friction factor is a nondimensional constant depending on the type and orientation of obstruction; is the velocity component in the ith direction, m/s; is the area porosity in the ith direction; is the coordinate in the ith direction.
Enthalpy conservation equation:where is the enthalpy, J/kg; is the effective viscosity, Pa·s; is the Prandtl-Schmidt number; typically = 0.7; is the heat transfer rate between the fluid and the internal obstructions in the volume, J/s.
Mass fraction of chemical species conservation equation: where is the mass fraction of a chemical species; is the Prandtl-Schmidt number; typically = 0.7; and is the fuel reaction rate inside the control volume, kg/(m3·s).
In this work, the model is used to describe turbulence properties during the explosion process. The conservation equations for the kinetic energy of turbulence () and its rate of dissipation () are expressed as follows [13, 14].
Turbulent kinetic energy conservation equation:
Dissipation rate of turbulent kinetic energy conservation equation:where is the turbulent kinetic energy, m2/s2; is the production of dissipation, m2/s3; and are the Prandtl-Schmidt numbers; typically = 1.0 and =1.3; C2 is a constant in (7); typically C2=1.92.
The stress tensor in the above equations is given bywhere is Kronecker delta function. =1 if , and =0 if .
The effective viscosity is defined as follows:where is a constant; typically =0.09. The second term of the right hand side of the equation above is known as the turbulent viscosity or eddy viscosity.
3.2. Physical and Computational Model
Gas explosion accident occurred in Apartment 101 on the first floor, which seriously damaged the ceilings and caused cracks, as shown in Figures 2(a), 2(e), and 2(f). Under this situation, the shock wave would spread to Apartment 201 on the second floor through ceiling cracks. However, other apartments around Apartment 101 were barely damaged based on onsite investigation. Thus, Apartment 101 and 201, which were damaged seriously in this accident, are selected to perform numerical simulation analysis in this work. According to the accident data from the explosion scene, a physical model that is consistent with actual situation developed based on the commercial CFD code FLACS, as shown in Figure 3(a). The indoor furniture and major equipment are also established in this model after onsite investigation, as shown in Figure 3(b). When the accident was triggered, the windows and doors that were entirely closed are specified as vent panels, as illustrated in Figure 4.
(a) Three-dimension model of Apartment 101 and 201
(b) Indoor layout
(a) Vent panel of windows
(b) Vent panel of doors
(c) Vent panel of cracks
Considering the reflection effect of pressure wave (illustrated in Figure 1), the entire computational domain has the dimensions of 15 m long, 24 m wide, and 9 m high. All boundaries of the whole computational domain are set as the EULER boundary, except along the balcony direction which is set as the PLANE_WAVE boundary. The computational domain is initially defined with the standard atmospheric pressure and a temperature of 20°C. Due to the fact that the leakage pressure and rate are very small in this study, the initial turbulence intensity prior to ignition could be ignored approximately. To make the vent panels fully resolved on the grid which follows the FLACS recommended guidelines for explosions [14], a uniform 0.15 m grid cell size is adopted to ensure calculation accuracy in this study. To compare and analyze the variations of explosion overpressure and temperature indoors, eight monitoring points (i.e., M1-M8) are arranged on the plane =1.5 m. To evaluate the damage degree of the explosion overpressure on the people’s head and ears, a typical monitoring point (i.e., M9) is arranged on the plane =1.6 m, as shown in Figure 5. The ignition source is located at the lighting switch in the living room, as presented in Figure 5.
3.3. Model Validation
In order to validate the present numerical model, the comparisons of the experimental data and numerical results for the 12.5% methane-air explosions are performed in a 12 m3 concrete chamber [15]. The chamber has the dimensions of 2 m long, 2 m wide and 3 m high, and has a square vent of 0.64 m2 located on the front wall. Six pressure transducers are mounted on the chamber to record the explosion overpressure. The ignition source is located in the center of the chamber. Two different explosion venting scenarios, i.e., the vent activation pressure = 21 kPa and = 55 kPa, are taken into account in this paper.
The present numerical results are compared with the experimental results recorded by the pressure transducer P6 [15], as presented in Figure 6. It can be found that the numerical results are in good agreement with the experimental results for both two different explosion scenarios.
4. Dynamic Characteristics of Indoor Gas Explosion
4.1. Identification of Initial Scenarios
In this paper, the numerical simulations are performed for different gas filling schemes to identify the initial accident scenario. In total, six kinds of gas filling schemes are chosen, as listed in Table 1. The variables mainly include gas cloud filling region (L×W×H) and equivalence ratio (Eq).
For scheme No.1 and No.2, the gas filling region has an overall size of 9.7 m × 9.95 m × 3 m. Figure 7 shows the pressure contours on the plane =1.5 m for these two different schemes. As shown in Figure 7, the maximum overpressure generated in the explosion process for Eq=1.0 and Eq=0.9 is 0.64 bar and 0.40 bar, respectively. According to the damage criterion of explosion overpressure on buildings [16], a maximum overpressure of 0.64 bar can cause the room-pillar to be broken off and the room frame loosened. For comparison, a maximum overpressure of 0.40 bar can cause a large crack in the walls, which is consistent with the real accident consequences. Figure 8 presents the pressure-time curves obtained at M9 for these two different schemes. As displayed in Figure 8, the maximum overpressure generated at M9 for Eq=1.0 and Eq=0.9 is 0.47 bar and 0.31 bar, respectively. Based on the injury criterion of explosion overpressure on personnel [17], a maximum overpressure of 0.47 bar and 0.31 bar can both cause the auditory organ injury, visceral hemorrhage, and fracture, which is not in accordance with the fact that the householder was just slightly injured. Thus, scheme No.1 and No.2 are inconsistent with the real accident consequences.
(a) Eq=0.9
(b) Eq=1.0
For scheme No.3 and No.4, the gas filling region has an overall size of 9.7 m × 8.5 m × 3 m. Figure 9 shows the pressure contours on the plane =1.5 m for these two different schemes. As shown in Figure 9, the maximum overpressure generated in the explosion process for Eq=1.0 and Eq=0.9 is 0.63 bar and 0.39 bar, respectively. According to the damage criterion of explosion overpressure on buildings [16], a maximum overpressure of 0.63 bar can cause the room-pillar to be broken off and the room frame loosened. For comparison, a maximum overpressure of 0.39 bar can cause a large crack in the walls, which is consistent with the real accident consequences. Figure 10 presents the pressure-time curves obtained at M9 for these two different schemes. As displayed in Figure 10, the maximum overpressure generated at M9 for Eq=1.0 and Eq=0.9 is 0.46 bar and 0.30 bar, respectively. Based on the injury criterion of explosion overpressure on personnel [17], a maximum overpressure of 0.46 bar and 0.30 bar can both cause the auditory organ injury, visceral hemorrhage, and fracture, which is not in accordance with the fact that the householder was just slightly injured. Thus, scheme No.3 and No.4 are also inconsistent with the real accident consequences.
(a) Eq=0.9
(b) Eq=1.0
For scheme No.5 and No.6, the gas filling region has an overall size of 9.7 m ×7.5 m × 3 m. Figure 11 shows the pressure contours on the plane z=1.5 m for these two different schemes. As shown in Figure 11, the maximum overpressure generated in the explosion process for Eq=1.0 and Eq=0.9 is 0.60 bar and 0.38 bar, respectively. According to the damage criterion of explosion overpressure on buildings [16], a maximum overpressure of 0.60 bar can cause the room-pillar to be broken off and the room frame loosened. For comparison, a maximum overpressure of 0.38 bar can cause a large crack in the walls, which is consistent with the real accident consequences. Figure 12 presents the pressure-time curves obtained at M9 for these two different schemes. As displayed in Figure 12, the maximum overpressure generated at M9 for Eq=1.0 and Eq=0.9 is 0.40 bar and 0.24 bar, respectively. Based on the injury criterion of explosion overpressure on personnel [17], a maximum overpressure of 0.40 bar can cause the auditory organ injury, visceral hemorrhage,s and fracture, which is not in accordance with the fact that the householder just got minor injury. In contrast, a maximum overpressure of 0.24 bar would make person injured slightly. Therefore, the initial scenario prior to ignition is identified as that the entire fuel region is filled with flammable gas cloud of a volume of 9.7 m × 7.5 m × 3 m and a gas concentration of 8.6%.
(a) Eq=0.9
(b) Eq=1.0
To further validate the reliability of CFD model, the simulation results are compared with the actual accident consequences shown in Figure 2. Figure 13 shows the pressure-time profiles obtained at different monitoring points. The comparison of the destructive effect of explosion overpressure on buildings between simulation and actual accident is given in Table 2. As can be seen, all the overpressure peaks simulated for the initial scenario identified earlier are within the overpressure peaks estimated based upon the actual accident scene. It means that the simulation results reproduced by the CFD model are in good agreement with the actual explosion consequences. Therefore, it is feasible and reliable to investigate indoor gas explosion accident based on the CFD model.
4.2. Pressure Evolvement and Distribution
Figure 14 shows the evolvement and distribution of pressure along with time on the plane =1.5 m. As can be seen, after the ignition in the initial accident room, the shock wave spreads quickly to other rooms through the inner doors, as presented in Figure 14(b). The pressure wave gradually covers the entire apartment at nearly 0.41 s. Meanwhile, the high-pressure gas rushes out of the apartment and spreads to the outside through the windows. Subsequently, a maximum explosion overpressure will be generated in the balcony at 0.47 s. The reason of this phenomenon lies in the fact that the window was closed in the balcony prior to ignition and the farthest flame propagation distance that would be produced is from ignition source to the balcony. Under these situations, the unburned gas compressibility inside the balcony and the flow resistance caused by interior walls would be enhanced during the propagation process. Thus, the higher turbulent intensity is generated, which contributes to speeding up the flame combustion velocity and creating a higher overpressure peak. After all combustion energy is released, the strength of the shock wave gradually decays, as illustrated in Figure 14(f).
(a)
(b)
(c)
(d)
(e)
(f)
In addition, the pressure contours also indicate that the overpressure distribution in a room is almost uniform at the same moment and there exists little spatial difference. This is mainly due to the high propagation speed of pressure wave in confined spaces.
4.3. Temperature Evolvement and Distribution
Figure 15 shows the evolvement and distribution of temperature along with time on the plane z=1.5 m. As can be seen, after igniting the premixed gas, burned gas expands rapidly and pushes unburned gas ahead of the flame front. The temperature rises rapidly around the ignition source, and the flame front diffuses from the ignition point to the surrounding. The flame front has a hemispherical shape and a smooth surface, as shown in Figure 15(a), which indicates that the flame propagation starts in a laminar manner. Subsequently, the propagation direction of flame front surface would be changed due to the constraint and barrier effects of interior walls. When the flame front enters each room though the inner doors, the flame begins to transform into irregular motion and the combustion process becomes unstable. The reason of this phenomenon lies in the fact that the cross section of gas flow becomes smaller, causing the increment of the airflow speed. Under these situations, the turbulence would be generated, which leads to the wrinkle and deformation of flame front surface, as displayed in Figure 15(b). Consequently, a larger flame surface area is produced, and hence the flame speed is increased. The flame quickly enters each room and then rushes out of the residential building through the windows, as presented in Figure 15(d). Eventually, the whole apartment would involve in an instantaneous high temperature environment.
(a)
(b)
(c)
(d)
Figure 16 shows the variations of flame speed along with the time at different monitoring points. After igniting the premixed gas, the flame speed rises slowly and the explosion pressure increases gradually. When the explosion pressure exceeds the activation pressure of vent panels, the vent panels are activated and the unburned gas is discharged, resulting in a reduction in flame speed. When the flame front interacts with indoor furniture and equipment, the turbulence effect is greatly enhanced and a larger flame surface area would be produced. Under these situations, the flame speed increases rapidly. Subsequently, the flame speed slows down with the emissions of the combustion products.
As displayed in Figure 16, the maximum flame speed is in the range of 34.3 m/s and 230.9 m/s, which is less than the sound speed under the corresponding pressure. It means that gas explosion inside residential buildings is a typical deflagration process [5], and thus a higher flame burning velocity would be generated in this process. Therefore, the explosive products would not deposit and leave black smoke traces on interior walls.
Figure 17 presents the variations of temperature along with the time at different monitoring points. As can be seen, the maximum temperature obtained at different monitoring points can reach up to 1953°C, which can cause secondary fire accidents easily inside residential buildings. In this accident, the cotton and cloth in the master bedroom were ignited by the instantaneous high temperature, which results in incomplete combustion. Consequently, a significant amount of black smoke traces was left on the walls and ground in the master bedroom. However, other rooms only experienced natural gas deflagration process, and there were no secondary fire accidents. Thus, the walls in these rooms still remained white and clean.
5. Mitigation Measures
Gas leakage would be easily caused by the high-frequency failure of gas facilities (i.e., hose pipes, valves, and cookers) in the kitchen. Once there is a poor ventilation condition, the accumulation of natural gas released might result in violent gas explosions. Thus, the leakage alarm instrument and automatic cut off valve are usually installed in the kitchen to prevent gas explosions. According to Code for Design of City Gas Engineering [18], adequate ventilation is ensured in the kitchen to mitigate gas explosions. However, the doors and windows within the kitchen are usually closed, especially in cold regions in northern China. Therefore, there is a need to study the effects of vent area ratio and the activation pressure of vent panels to propose effective mitigation measures for gas explosions inside residential buildings. In this section, numerical simulations are performed for different explosion venting scenarios by varying one parameter to analyze the variations of overpressure peak in the kitchen. This kitchen has the dimensions of 3.3 m long, 1.6 m wide, and 3 m high, as shown in Figure 18. Two vent panels are configured within the kitchen. Each vent panel has the dimensions of 0.9 m wide and 2.3 m high. A summary of the initial data for these simulations is shown in Table 3, including vent area ratio and the activation pressure of vent panels. In this paper, the vent area ratio is defined as , where is the vent area and is the volume of the kitchen [19].
5.1. Vent Area Ratio
In this section, numerical simulations are performed by varying the width of vent panels in the kitchen. The relevant parameters listed in Table 3 are used to analyze the effect of vent area ratio () on the overpressure peaks. Figure 19 shows the typical overpressure peak variations for the activation pressure of vent panels of 0.02 bar, with varying the vent area ratio from 0.1 to 0.7 (i.e., the width of vent panels ranging from 0.14 m to 0.96 m).
As can be seen, the overpressure peak decreases with the increase of vent area ratio. This is due to the fact that the larger vent area results in more combustible gas expelling from the kitchen in the unit time. In addition, the overpressure peak decreases rapidly with the increase of vent area ratio when , while it is almost independent on vent area ratio when . Therefore, the vent area ratio is recommended to be used as the minimum vent area ratio in the venting design of the kitchen from the point of view of mitigating gas explosions. Furthermore, Figure 19 also illustrates that the overpressure peak in a more enclosed space is much greater than that in an open space. It indicates that natural gas should not be used in places with poor ventilation conditions, such as a basement or a room without windows. If natural gas must be used in these places, the forced ventilation facilities need to be equipped.
5.2. Activation Pressure of Vent Panels
It is in common that the broken pressure of glasses is above 0.02 bar [20, 21]. So the activation pressure of vent panels varying from 0.02 bar to 0.14 bar is considered to analyze the effect of activation pressure on the overpressure peaks. All other parameters remain unchanged, as given in Table 3. Figure 20 shows the typical overpressure peak variations at , with varying the activation pressure of vent panels from 0.02 bar to 0.14 bar.
As can be seen, there is a proportional relationship between the overpressure peak and the activation pressure of vent panels. The reason of this phenomenon lies in the fact that the increase of activation pressure will induce a delay in the opening time of vent panels, which directly contributes to enhancing the unburned gas compressibility and turbulent intensity and thus causes an increment in overpressures. However, the activation pressure is closely related to the strength and thickness of vent panels [15]. The larger the strength and thickness of vent panels are, the greater the activation pressure is. Thus, the low-intensity materials should be selected as the vent panels in the venting design of the kitchen from the point of view of mitigating gas explosions, such as the single glazed windows or doors.
6. Conclusions
The numerical model of gas explosion inside a residential building is developed based on the CFD code FLACS. The numerical simulations are performed for different gas filling schemes to identify the initial scenario. Meanwhile, the simulation results are compared with the explosion consequences associated with a real accident. In addition, the dynamic characteristics of explosion overpressure and indoor temperature are analyzed. Furthermore, the effects of vent area ratio and the activation pressure of vent panels in the kitchen are investigated to propose effective mitigation measures for gas explosions inside residential buildings. The conclusions are mainly provided as follows.
(1) The initial scenario prior to ignition is identified as that the entire fuel region is filled with flammable gas cloud of a volume of 9.7 m × 7.5 m × 3 m and a gas concentration of 8.6%.
(2) All the overpressure peaks simulated for initial scenario are within the overpressure peaks estimated for the actual accident scene. It indicates that the simulation results are in good agreement with the real accident consequences. Thus, it is feasible and reliable to investigate indoor gas explosion accidents based on the CFD model.
(3) The overpressure distribution in a room is almost uniform at the same moment and there exists little spatial difference.
(4) The maximum flame speed is in the range of 34.3 m/s and 230.9 m/s, which means that the gas explosion inside residential buildings is a typical deflagration process. Thus, the explosive products would not deposit and leave black smoke traces on interior walls.
(5) The maximum temperature can reach up to 1953°C in each room, which can cause secondary fire accidents easily inside residential buildings. In this accident, the cotton and cloth in the master bedroom were ignited by the instantaneous high temperature, which led to incomplete combustion. Thus, a significant amount of black smoke traces was left on the walls and ground in the master bedroom.
(6) The overpressure peak decreases rapidly with the increase of vent area ratio when , while it is almost independent on vent area ratio when . There is a proportional relationship between the overpressure peaks and the activation pressure of vent panels.
(7) Natural gas should not be used in places with poor ventilation conditions, such as a basement or a room without windows. If natural gas must be used in these places, the forced ventilation facilities need to be equipped. The single glazed windows or doors are usually installed in the kitchen to mitigate gas explosions.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (No. 51474184) and the Natural Science Foundation of the State Administration of Work Safety in China (No. 2012-387, Sichuan-0021-2016AQ).
Supplementary Materials
The supplementary materials files mainly include a computational model and the corresponding data analysis used in our manuscript to make readers verify the results of our article and replicate the analysis. (Supplementary Materials)