Quantitative Method of the Structural Damage Identification of Gas Explosion Based on Case Study: The Shanxi “11.23” Explosion Investigation

Huanjuan Zhao, Yiran Yan and Xinming Qian

(1.Mine Emergency Technology Research Center(University of Science and Technology Beijing), State Administration of Production Safety Supervision and Administration, School of Civil and Resources Engineering, University of Science and Technology Beijing, Beijing 100083, China; 2.State Key Laboratory of High-Efficient Mining and Safety of Metal Mines (University of Science and Technology Beijing), Ministry of Education, Beijing 100083, China; 3.State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China)

There is a possibility for explosions to occur in building structures containing combustible gases. Explosions may lead to severe financial losses and mass casualties[1]. For example, on August 12, 2015, an accidental explosion occurred in Binhai New District of Tianjin,China, which generated a shock wave with a radius of 3 000 m and resulted in 6 000 people becoming homeless. Another example is “the 11.23 accident” which occurred in Shouyang town, Shanxi Province on November 23, 2012, resulted in mass destruction. Because of such kind of accidents,it is necessary to establish an effective and quantitative method of analyzing explosion accidents and setting up a systematic analytical database that can lower the risk of future explosions. In order to establish such a system, it is of crucial importance to conduct risk assessments[2]to determine the cause of accidents, and provide anti-explosion protection to structures[3]. S.Sarshar[4]illustrates the challenges and proposals for managing major accident risk including investigation. Furthermore, Haojie Xu[5]studies the structural damage identification based on a modified Cuckoo Search algorithm. Such measures would benefit investigations of existing explosion sites and help prevent future accidents. However, the accuracy of conventional analytical methods used to review explosion accidents fails to present quantitative analysis. Moreover, complex factors, such as the characteristics of the structure and the type of gas involved in the explosion, make it difficult to determine the definitive cause of the explosion.

Currently, empirical methods or semi-empirical methods based on assumptions carried out in open areas are frequently applied to compute explosion characteristics in engineering projects[6]; however, such results differ significantly from the distribution of explosion pressure recorded in actual structures and lead to relatively large errors in the analysis of accidental explosions. Moreover, it is unrealistic to perform a completely reproduced experiment of the accident. Importantly,Baker says that PHYSICAL simulation experiments, under known and specified conditions (e.g. wind speed, humidity), can reach more precise conclusions regarding a known event[7].

As such, research groups have conducted in-depth and detailed theoretical and experimental studies concerning the damage effects on various structures under the influence of explosion loads, and have proposed prevention and control theories of explosions for different structure types, as well as corresponding design plans[8]. For instance, Remennikov et al.[9]points out the importance of evaluating explosion effects on nearby buildings when determining the explosion load of a specific building, and Luccioni et al.[10]used the dynamic analysis software AUTODYN to simulate the transmission and reflection of shock waves in an explosion. Meanwhile, a number of domestic organizations and researchers have also made significant contributions to this research field. For example, Wang Haifu et al.[11]applied a numerical simulation approach to study the explosion responses of various structures.

Davis S G uses FLACS simulations to investigate the chain of the explosion events and provide a more complete understanding of the evidence, including near-field blast damage[12].

However,little research regarding the quantitative damage effects of gas explosions on structures has been conducted. While 2,4,6-trinitrotoluene (TNT) conversion is widely used to study gas explosions, research on explosion effects concerning the intrinsic properties of gases is rare[13]. Qian Xinming et al.[14]utilized numerical simulations in order to study the effects of gas explosions, but this method has not been systematized.

In reality, other than ignition, gas species, gas volumes, and the location of explosions are all factors that might influence explosion effects[15-18], and thus all require a quantitative analysis. By analyzing key issues involved in accidental explosions, we developed a semi-quantitative analysis method of accident investigation involving industrial chemicals. Based on investigations of various explosion accidents, and using a specific, complicated accident as an example, we used AUTODYN to quantitatively study the damage effects of a gas explosion on building structures. Using the results from our analysis, safety guidelines are provided for the prevention of combustible gases inside structures. Furthermore, our study considers the inner and outer structures of a building, which impact the transmission of an explosion and the interactions between nearby structures.

1 Data Analysis of Accidents

1.1 Analysis of typical characteristics

Using information gathered from the investigation of the original accidental explosion site, including measurement data, the original structures at the site of the accident were restored and simplified, and the overall dimensions of each structure were reserved, as illustrated in Fig.1A. The original structure contained windows in the northern wall, a pot on a windowsill, and 3 cabinets in the northwestern corner against the wall. Additionally, there was a staircase against the western wall, 3 partitions that formed 4 cubicles against the eastern wall, and 2 freezers against the southern partition (situated north to south). The accident site is shown in Fig. 1B, where glass from the windows along the northern wall and glass from the light inlet window in the well platform are all broken and thrown outwards. The damage to the wooden partition in the north part of the structure is severe, with only part of the southern most partition remaining, as shown in Fig. 1C. From Fig. 1D, it can be seen that the damage to the wooden partition below the staircase located in the southwest portion of the structure is trivial. Fig.1E shows no explosion damage to the doors of the disinfection cabinet and no damage to the electric water heater and the sink, whereas in Fig. 1F, there are apparent dents on the top surface of both freezers. Between the third and fourth grid in the window well, counting from the west side and between the steel frame and the wall of the window well, the pot (now distorted and damaged) can be seen.

A—the windows; B—the south-east part; C—the north-east part; D—the sterilizer; E—the north freezer; F—the south freezer; G—the pot after the occurance of the accident Fig.1 Damage of the combustible gas explosion accident

The pot was impacted from below and flew outward in an oblique upward direction via the window, Fig. 1G. Its trajectory can be resolved to displacements in 3 directions: horizontal displacementL1; lateral displacementL2; and vertical displacementL3.L1andL2were measured to be 1.22 m and 2.66 m. The pot’s trajectory has an aspect ratio of 2.66/1.22, given in Fig. 2. The top of the building structure is divided into 4 zones with the use of girders; lacing wires in each zone show bend deformations(of differing extents), as shown in Fig. 3. Overall, it can be seen from the bending degree of roof wires in Fig. 3 that the western side of the structure is heavier than the eastern side, the northern side weighs more than the southern side, and the northwestern portion of the structure is the heaviest. Factors are labeled asC1,C2…Cn, with ∮(c) being the ultimate cause of the accident.

Fig.2 Pot’s flight trajectory

This explosion accident is unique and complex, and has various potential causes. Structures that were damaged in the building include 2 freezers, the eastern partition, the pot (location changes), windows (glass), and lacing wires inside the ceiling. Compared to other areas in the roof, the northwestern portion of the roof experienced more damages. The cross-sectional area of lacing wires inside the ceiling is small enough so that it has limited influence on the transmission of shock waves; however, since it still affects the computation of step size, it must still be simplified and eliminated. The destructive conditions of the roof and ceiling are thus represented by pressure values.

A—the southwest corner; B—the southeast corner; C—the northeast corner; D—the northwest cornerFig.3 Damage details of the roof

1.2 Analysis of explosion causes

1.2.1Suspicious gas species

A natural gas pipeline is buried approximately 3.25-3.45 m away from the northern exterior wall of the structure. The natural gas pipeline is buried in the east-west direction at a depth of 0.80-0.85 m (Fig. 4A). In the outer wall of this pipeline, there is a hole with a diameter of 8 cm that faces the fifth steel girder of the window well (counting from the west), which is shown in Fig.4B and Fig.4C. An empty liquefied petroleum gas tank of 15 kg is located inside the northwestern portion of the building, see Fig. 4D. The gas that caused the explosion was determined to be liquid petroleum gas (LPG) or natural gas.

Fig.4 Natural gas pipe and the liquid petroleum gas tank found in the accident scene

Once the LPG is inside a tank escapes, it immediately diffuses and burns or explodes when it meets an ignition source. The density of LPG is higher than the density of air, and has the volume fraction of around 2%-10%. The main constituent of natural gas is methane, which is the major cause of natural gas explosions. Natural gas flows upwards in the air and its explosion limit is around 5%-15% (volume fraction).

1.2.2Gas parameters of the explosion source

① Gas contents and distributions

The height of the LPG tank in the structure is 0.9 m. In general,LPG diffuses from the bottom of a tank upwards. The vertex of the natural gas pipeline outside the basement is nearly the same height as the roof. As such, the natural gas from the pipeline diffuses upward first. By reviewing the characteristics of the accident, we determined that the source of the explosion could be located in the western area inside the building (the northwestern corner or anywhere along the western side). Locations and gas species of the explosion source are arranged in pairs, and, after preliminary computations, scenarios concerning the distributions of 5 gases are provided in Tab. 1. Hereafter, scenario 1, 2, 3, 4 and 5 are respectively known as HGAS I, HGAS II, HLPG I, HLPG IIand HLPG II. In the same way, HGASmeans scenario 1 and 2 while HLPGmeans scenario 3, 4 and 5.

② Determination of pressure and temperature of mixed gases

A simulation method of high-pressure gas is applied with the consideration of the special properties of combustible gases. Chemical reactions take place instantaneously, and the gas explosion process inside this structure is approximated to be an isometric process where the pressure increases dramatically (State 1: before explosion;State 2: at the moment of explosion).In this study, the main parameters considered are: gas density (ρ), volume (V), explosion location, pressure (P), and temperature (T). According to the conservation of internal mass,ρis constant. Parameters of the initial state and the final state can be connected by 3 principles: a. the principle of conservation of mass; b. the principle of conservation of momentum; and c. the principle of conservation of energy. Through equation of states for a gas at a constant volume, we calculated heat of reaction, enthalpy change, and the pressure and temperature of natural gas at the moment of explosion.

Tab.1 Scenarios methods

2 Numerical Simulation of the Gas Explosion Process Inside the Basement

The impacts of burning and explosion products on the structure were the results of a fliud-structure interaction problem.The failure of critical components is represented by the factorCn(the constant value in the accident results) and constitutes the failure of a critical component.Itsvalue was set to “1” when the failure was consistent with the actual conditions of the accident, and was set to “0” when the failure was inconsistent with the actual conditions of the accident.

2.1 Computational method

Walls were considered rigid. Shell grids were adopted for structures that contain glass, freezers, staircases, and partitions. Euler grids were utilized to represent air and high-pressure products. In this paper the Shell/Euler mixed computational modeling method is adopted and a HP-Z800 high-performance workstation is chosen as our computational machine.

2.2 Discrete model

Fig. 5 shows a discrete model of the building structure and a 3D grid model. A Shell grid model with a minimum grid length of 15 mm and a total grid number of 22 180 was used for the walls and windows. The air domain of the Euler grids were filled with products that contained “high-pressure explosive products” (shown in Fig. 6). Furthermore, the boundary conditions were based on the working conditions of actual structures, where the location of windows along the boundary was set as an outlet boundary to ensure that the explosion products could flow out smoothly.

Fig.5 Model of the basement and inner objects

2.3 Material model

We selected either the dynamic material database in the AUTODYN software or obtained the data experimentally. We chose different material models for different components, as shown in Tab. 2.

For the material models shown in Tab. 2, the equation of state for an ideal gas: the ideal gas constant (γ) is defined as

Fig.6 Computational domain of the 5 scenarios

ObjectMaterialStatefunctionConstitutivemodelFailuremodelWallCONCRETE-LRigidNoneNoneAirAirIdealgasNoneNoneFillinggasHE-AirPolynomialNoneNoneGlassGLASS-EPXYLinearVonMisesHydro(Pmin)PartitionCONCRETE-LLinearDrucker-PragerNoneRoofbeamCONCRETE-LRigidDrucker-PragerNonePotAL2024LinearBilinearhardeningPrincipalstrainFreezerAL2024LinearBilinearhardeningPrincipalstrainSterilizerAL2024LinearBilinearhardeningPrincipalstrainStaircaseSteel1006LinearJohnsoncookPrincipalstrain

p=(γ-1)ρe+pshift

(1)

whereρis the density,eis the internal energy, andpshiftis the initial pressure.

Polynomial Equation of State: the generalized polynomial function of pressure is the function of compressed density, which is a form of the Mie-Gruneisen equation of state.It applies different methods to analyze pressure and tensile force.

Linear Equation of State: provides the definition of the bulk modulus and reference density:

p=Kμ

(2)

whereμ=(ρ/ρ0)-1;Kis the bulk modulus of the material.

Johnson Cook: a model of strain hardening. Strain rate is related to the temperature:

(3)

THis the corresponding temperature which can be expressed asTH=(T-Tinitial)/(Tmelt-Tinitial) whereA,B,C,n, andmare constants of the material.

Bilinear Hardening: in this model the equation of stress is

(4)

whereσis the stress of the material;εis the effective plastic strain;σsis the yield limit;Eis Young’s modulus; andETis the shear modulus.

Von Mises: the yield surface and shear modulus are both defined as constant. Hydro: a hydrostatic tensile stress. If a negative pressure is reached, failure occurs. Principal strain: if the maximum main strain or shear strain exceeds their corresponding failure strain, then failure occurs. After the date was validated, we withdrew the results.

3 Discussion

3.1 High-pressure gas pressure contours

In order to reflect the pressure responses of a high-pressure gas, a cross-sectional screenshot is provided from the north-south direction. Fig.7 shows the pressure contours in the Euler computational domain at several representative moments, where shock waves overlap after rebounding inside the room.

Fig.7 Computational domain pressure in typical time

This incurs relatively high pressure in certain areas, and could possibly destroy weaker structures. The peak pressure values in HLPGand Hgas, respectively, both reached an order of 104Pa; however, the peak pressure in HLPGwas slightly higher than that in Hgasdue to the greater heat combustion. Furthermore, pressures in HLPGand Hgasboth decreased to a low value along with the explosion propagation.

Shock waves inside the room interact with each other. Pressure responses can be relatively complex, and can be roughly reflected as follows: in Hgas I, there was a trend in which peak pressure values first occurred in the northwestern corner of the structure and then in the southwestern corner, with the northwestern corner bearing a higher pressure; in HGAS II, pressure changed from the northwestern corner of the structure to the southwestern corner, where the northwestern corner and the area above the southwestern corner had relatively high pressure (inconsistent with the characteristics of the accident); in HLPG I, pressure changed between the northwestern corner and the southwestern corner from the bottom of the structure to the top, where the northwestern corner had a higher pressure; HLPG II, pressure changed between the northwestern corner and the southwestern corner from the bottom of the structure to the top, where the northwestern corner and the middle area had higher pressures (inconsistent with the characteristics of the accident); in HLPG III, the pressure changed between the northwestern corner and the southwestern corner from the bottom of the structure to the top, where the northwestern corner had a higher pressure.

3.2 Analysis of material status in each structure

Fig. 8 shows material status in different structures at several representative moments where there is an apparent stress concentration zone. The stress of HGASand HLPGboth reached an order of 105Pa.

Fig.8 Material state diagram in typical time

① In each scenario,the glass (window panes) and partitions in the eastern portion of the structure were among the first items to be badly damaged or destroyed. The damage speed of the partitions in the scenarios where the northwestern corner and southwestern corner were filled with explosive gas was significantly faster than the one in the scheme where merely the northwestern corner was filled with the explosive gas. Shock wave overpressure (9-16 kPa) impacted the top surface of the freezers, and caused local dent displacement. ② In HGAS, instead of flying outwards in an oblique upward direction, the pot moved along the windowsill due to pressure being applied on its top surface. This trajectory is inconsistent with the accident site. In HLPG, the pot was impacted from below and flew outward in an oblique upward direction via the window. Its trajectory can be resolved to displacements in 3 directions: horizontal displacementL1; lateral displacementL2; and vertical displacementL3. When the horizontal flight distanceL1=1.22 m, the flight trajectory of the pot must satisfy the requirement ofL2∶L1=2.66∶1.22, as shown in Fig. 9. In HLPG1, the trajectory error of the pot was extremely low (0.458 7%), as indicated in Fig. 9 where the dotted line represents the actual data.

Ultimately, as demonstrated by this study, explosion gas species, volumes of gas(es), and explosion locations have a significant influence on the damages incurred. Gas species is the most significant factor as it has more influences on the resulting explosion than gas volumes and explosion locations. Results from HLPG Ibest match the actual destruction observed in the original accident, and it is consistent with information provided by witnesses. Thus, the results of HLPG Ivalidate the computational analysis methods outlined above. Furthermore, forces applied to several typical monitoring points on the roof are analyzed in details, shown as below.

Fig.9 Comparison of the pot’s flight trajectory

3.3 Pressure analysis of typical monitoring points on the roof

The pressure-time curves of typical monitoring points are shown in Fig.10. In the northwestern portion of the roof the initial pressure was: PHLPG I0=PHLPG II0=PHLPG III0=0. The sample point in HGASwas located inside the filling gas and the values of PHGAS I0and PHGAS II0were relatively big. The pressure in HLPGincreased rapidly, with its subsequent peak pressure being higher than that in HGAS.

In the northeastern portion of the roof the initial pressure was: PHGAS I0=PHGAS II0=PHLPG I0=PHLPG II0=PHLPG III0=0. The pressure values in HGASand HLPGwere both lower than the corresponding peak pressures in the northwestern corner.

In the southwestern portion of the roof the initial pressure was: PHGAS I0=PHLPGI0=PHLPG II0=PHLPG III0=0. The pressure values of HGASand HLPGwere both lower than the corresponding peak pressure in the northwestern corner.

In the southeastern portion of the roof the initial pressure was: PHGAS I0=PHGAS II0=PHLPGI0=PHLPG II0=PHLPG III0=0, followed by the main results of PHGAS II>PHLPG III>PHLPG II>PHLPG I>PHGASI. In HLPG, the pressure in the northwestern corner of the roof was higher than that of its surroundings; however, in HLPGII, damage to the floors in the southwestern corner was relatively severe. In HGAS, there was no significant difference between the pressure measured in the center of the northwestern portion of the roof and that of its surroundings.

The pressure analysis outlined above demonstrates that HLPGIbest matches the characteristics of the accident. It also illustrates that the width/height ratio of accumulated gases influences the generated pressure and impact direction.

3.4 Impact analysis of personnel locations

Fig. 11 shows monitoring points near personnel locations, where point 28 is the location monitoring point of Person No. 1 (behind the freezer) at a height of 50 cm above the ground; point 29 is the location monitoring point of Person No. 2 (near the sink) at a height of 90 cm above the ground. Corresponding pressure-time curves are shown in Fig. 11.

Fig.11 Pressure-time curves of critical personnel locations

The pressure and the ascending velocity of HGASat the location of Person No. 1 were lower than those of HLPG. Among the 3 HLPGscenarios, the peak value of HLPGⅡwas slightly higher than the ones in the other 2 scenarios.

The pressure and the ascending velocity of HGASat the location of Person No. 2 were lower than HLPG. The initial pressures in the 3 HLPGscenarios were all fairly high, but decreased immediately and rapidly. The peak value of HLPGⅡwas slightly higher than that in the other 2 scenarios.

All the pressures were lower than 0.05 MPa, which means that shock waves could have caused injuries to personnel (such injuries would not have been fatal).The influence of gas properties on explosion results are summarized in Tab.3 (allowing data based on a quantitative analysis to be easily compared). Tab. 4 provides the constant values of key components in Tab.3. The valueCis 1 if the result agrees well with the homologous accident factor, otherwise the valueCis 0.

Finally,constant rates for the 5 scenarios can be calculated as Additionally, potential factors and effects relevant to the analysis are illustrated in Tab.5.

Tab.3 Statistic table of typical simulation results

Note: Data in the table above indicates the height above the ground; slight damage is less than some damage; some damage is less than damage

Tab.4 Constant value of key components

Tab.5 Influence of gas properties on explosions

4 Conclusions

① In order to recreate the actual structures involved in the explosion,the original accident was examined. The fracture strengths of different locations were inspected, while the damage characteristics of each structure (freezers, partitions, the pot, window glass, and lacing wires of the ceiling) are determined, and the explosion gas was confirmed as liquefied petroleum gas or natural gas. Based on the gas properties and device conditions, it was determined that the source of the explosion was either in the northern corner of the structure or in the internal portion of the western side of the building. After preliminary computations, the gas distribution of the 5 scenarios were quantified.

② A physical model consistent with the actual conditions was used in order to simulate the damage effects of several gas explosion scenarios on the structures. Significant differences in each scenario are obtained, then we analyzed the pressure contours and material status of different structures, and compared pressure characteristics of typical monitoring points on the roof. Computational results from scenario 1 (HLPG I) were consistent with the primary characteristics of the original accident, which validates the use of explosion gas species, gas concentrations, and leaking locations as analysis parameters. We confirmed the results and demonstrated the feasibility of a high-pressure gas simulation method. Ultimately, our method accurately evaluates gas properties, displays quantitative advantages in analyzing accidents, and directly describes the evolution process of gas explosions inside each structure. Furthermore, our method provides an additional approach to study explosion prevention, which identifies weak structures, and generates data that could improve the explosion impact resistance of structures.

③ Pressure-time curves for monitoring points in the same place show the influence of width-height ratio of accumulated gases on the pressure and impact directions. Explosion effects are affected by gas species (relative density influences the impact direction and heat of combustion influences the pressure), volumes of gases, and explosion locations. Gas species influences explosion effects the most, followed by gas volumes and explosion locations.

Ultimately, an effective, feasible, and quantitative method was established for analyzing accidental explosions. In our method, accident characteristics with recognizable impact direction and impact strength are determined first. Then combustion gas species and leakage are analyzed and quantitative computation is performed. By comparing the influence of gas properties on explosion consequences, it is capable of us to proposepossible strategies for explosion prevention. Finally, the dangers and dynamic processes of combustible gases in limited spaces were demonstrated.

[1] Henrych Josef. The dynamics of explosion and its use[M]. New York: Elsevier Scientific Publishers, 1979.

[2] Zhou Jianfeng. SPA-fuzzy method based real-time risk assessment for major hazard installations storing flammable gas [J]. Safety Science,2010, 48(6):819-822.

[3] Jeom Kee Paik, Jerzy Czujko, Bong Ju Kim, et al. Quantitative assessment of hydrocarbon explosion and fire risks in offshore installations [J]. Marine Structures, 2011 24(2): 73-96.

[4] Sizarta Sarshar, Stein Haugen, Ann Britt Skjerve. Challenges and proposals for managing major accident risk through the planning process [J]. Journal of Loss Prevention in the Process Industries, 2016, 39: 93-105.

[5] Xu H J, Liu J K, Lu Z R. Structural damage identification based on modified Cuckoo search algorithm [J]. Structural Engineering and Mechanics, 2016,58(1): 163-174.

[6] Mohammad Dadashzadeh, Rouzbeh Abbassi, Faisal Khan. Explosion modeling and analysis of BP Deepwater Horizon accident [J]. Safety Science, 2013, 57: 150-160.

[7] Baker A R. The role of physical simulation experiments in the investigation of accidents [J]. Journal of Occupational Accidents, 1982, 4(1): 33-45.

[8] Sun Jianhua, Fan Hongyu, Wei Chunrong, et al. Numerical simulation on foam ceramic blasting block device under the action of explosion transform [J]. Safety Science, 2012, 50(4): 588-592.

[9] Remennikov A M, Rosa T A. Modeling blast loads on buildings in complex city geometries [J].Computers & Structures, 2005, 83(27): 2197-2205.

[10] Luccioni B, Ambrosini D, Danesi R. Blast load assessment using hydro codes [J]. Engineering Structures,2006, 28(12): 1736-1744.

[11] Wang Haifu, Zhang Yuanfeng, Yu Qingbo, et al. Study on initial mechanism of reactive fragment to covered explosion[J]. Transaction of Beijing Institute of Technology, 2012,32(8): 786-789, 823. (in Chinese)

[12] Davis S G, Engel D, Gavelli F, et al. Advanced methods for determining the origin of vapor cloud explosions case study: the 2006 Danvers explosion investigation [J]. Fire Technology, 2014, 50: 823.

[13] Yoshitomo Inaba, Tetsuo Nishihara, Mark A Groeth, et al. Study on explosion characteristics of natural gas and methane in semi-open space for the HTTR hydrogen production system [J]. Nuclear Engineering and Design, 2004, 232(1): 111-119.

[15] Qian Xinming, Zhao Huanjuan, Nan Yuxiang, et al. Simulation analysis of destructive effect for gas explosion happened in thin-shell construction [J]. Journal of Vibration and Shock,2015, 34(9): 82-87.

[16] Radulescu M I, Sharpr G J, Lee J H S, et al. The ignition mechanism in irregular structure gaseous detonations [J]. Proceedings of the Combustion Institute, 2005, 30: 1859-1867

[17] Raphael Moura, Michael Beer, Edoardo Patelli, et al. Learning from major accidents to improve system design [J]. Safety Science, 2016, 84: 37-45.

[17] Li Weijun,Zhang Laibin, Liang Wei. An accident causation analysis and taxonomy (ACAT) model of complex industrial system from both system safety and control theory perspectives[J]. Safety Science,2017,92:94-103.

[18] Duan Xiaoyu, Cui Qingzhong, Guo Xueyong, et al. Experimental investigation on shock wave characteristics of aluminized explosives in air blast[J]. Journal of Beijing Institute of Technology, 2017, 26(2): 165-173.