The statement of the problem and algorithm of computational modeling of the processes of formation of the hydrogen-air mixture in the atmosphere, its explosion (taking into account chemical interaction) and dispersion of the combustion materials in the open space with complex relief are presented. The finitedifference scheme was developed for the case of the three-dimensional system of gas dynamics equations complemented by the mass conservation laws of the gas admixture and combustion materials. The algorithm of the computation of thermal and physical parameters of the gas mixture appearing as a result of the instantaneous explosion taking into account chemical interaction was developed. The algorithm of computational solution of the difference scheme obtained on the basis of Godunov method was considered. The verification of the mathematical model showed its acceptable accuracy in comparison with known experimental data. It allows using the developed model for the modeling of pressure and thermal consequences of possible failures at the industrial enterprises which store and use hydrogen. The computational modeling of an explosion of the gas hydrogen cloud appearing as a result of instantaneous destruction of high pressure containers at the fueling station was carried out. The analysis of different ways of protection of the surrounding buildings from destructive effects of the shock wave was conducted. The recommendations considering the choice of dimensions of the protection area around the fuelling station were worked out.
Hydrogen is widely used as an alternative transport fuel [1] which is harmless to the environment [2]. However, such characteristics of hydrogen [3] as low density, high combustion energy and rapid transition from burning to detonation cause the problems of safe storage [4] and delivery of the hydrogen fuel [5] and right allocation of fuelling stations with regard to residential areas [6]. The equipment seal failures, destruction of storage volumes of high compressed hydrogen, cause its release into the atmosphere and formation of an explosive hydrogen-air mixture [7]. It creates a real threat of the hydrogen inflammation, detonation explosions [8] and, as a result, considerable material damage [9].
A physical experiment modeling these gas-dynamics phenomena is very expensive [10]. Its results cannot be easily used under real conditions of industrial enterprises. For this reason, it is more advantageous to use computational methods, which can have varying degrees of complexity and accuracy [11]–[26]. As a rule, overpressure and shock wave impulse are determined using semi-empiric equations of regression to predict construction loadings after an explosion [17]. However, experimental data are usually obtained in the open space, not taking into account complex relief [27]. Modern computation methods of burning-to-detonation transition (e.g. [28]– [32]) are developed for model problems. Therefore, the problem of the creation of the mathematical model describing adequately time-dependent processes of formation of explosive gas mixtures in the three-dimensional space (due to evaporation from liquid hydrogen spill [33] taking into account ambient conditions [34]), their explosion (taking into account chemical interaction of the mixture components) and further dispersion of the combustion materials in the atmosphere is very important. A computer system realizing such a mathematical model will allow analyzing and forecasting three-dimensional fields of the explosive admixture concentration, mixture thermodynamics parameters in time (before and after explosion) and space, evaluating possible after-explosion infrastructure destruction, suggesting rationally designed protection devices which will provide a sufficient level of safety for people in the accident zone [35] and their sustainable strength so as not to collapse themselves [36]. Moreover, an adequate mathematical model of the reacting gas mixture flow will allow not only to evaluate the dimensions of the protective device, but also to select the material for its manufacture for an explosion of different power [37]. It is also necessary to take into account the model’s ability to provide the necessary data for deterministic [38] or probabilistic [12] approaches to assess the consequences of accidents for the surrounding area.
An adequate description of physical processes of hydrogen release and further mixture dispersion in the atmosphere is possible only using the system of Navier-Stokes time-dependent equations (NSE) for compressible gas. The limited capabilities of the computing technology available to engineers do not allow carrying out direct computational solution of these equations effectively. In many cases, the computation modeling of turbulent flows is realized by solving Reynolds-averaged NSE, complemented by a turbulence model [39]–[41]. But a great majority of the turbulence models do not describe different types of flows adequately. It especially relates to flows with an intensive boundary flow separation and/or heavy changes of such flow parameters as temperature or pressure. Therefore, it is necessary to create new mathematical models and finite-difference schemes for a computation modeling of such flows.
The main purposes of this work are to develop a simplified mathematical model describing adequately time-dependent processes of explosive gas mixtures formation in the three-dimensional space, their explosion and dispersion of combustion materials in the atmosphere as well as to create a computational modeling algorithm of these processes.
It is assumed that the mass, impulse and energy transfer impacts for the most part on actual physical processes. Thus, it is reasonable to apply roughened NSE without taking in account viscosity (Euler equations with sources [42]).
The computation area Ω is a parallelepiped located in the right-hand coordinate system (fig. 1) and it is partitioned to spatial cells, the scale of which depends on characteristic sizes of the computation area (roughness of streamlined surfaces, dimensions of streamlined objects). The system of gas mixture motion differential equations can be represented as follows [43]: ⃗ + (⃗ + ⃗ + ⃗ = ⃗, where vectors ⃗, (⃗, ⃗, ⃗, ⃗ look like this:
⃗ = [, , , , ]", (⃗ = [, + #, , , ( + )]", ⃗ = [, , + #, , ( + )]", ⃗ = [, , , + #, ( + )]",
⃗ = [0, 0, −, 0, −]",
= @ + 1D2 (# + # + #)E,
(
The equation of motion of every admixture (combustible gas, combustion materials) is as follows [43]: () + () + () + () = $!, where Q represents the admixture mass concentration, $! – mixture component diffusion source member (according to Fick dependence [44], [45] looks like this $! = (%), and the diffusion coefficient % can be evaluated using the theory of Berlyand [46]–[48])
The set of equations (
It is assumed that any component of air flow velocity is subsonic. The approaching flow is
defined by the values of total enthalpy & = D − 1 D + 1D2 (# + # + #), entropy function
& = ⁄', flow velocity vector (angles (, )), and relative admixture mass concentration (£1
if the gas admixture flows in). The entry flow parameters are defined by equations (
At start time in all "gaseous" cells of the computational space, the ambient parameters are assigned. In the cells occupied by an admixture cloud, the relative mass concentration of the admixture equals = 1 (100%).
According to the suggested combustion model, it is assumed that the global instantaneous chemical reaction takes place in the volume where the concentration of the admixture is in the inflammability range (a control volume). It means that the values of the parameters of the gas mixture in the control volume are instantly changed by the values of the corresponding parameters of combustion materials and remnants of the one of the mixture components (combustible gas in the case of thin mixture or air in the case of rich mixture). It is assumed that the flame front propagates with infinitely large speed.
The main purposes of this work are to develop a simplified mathematical model describing adequately time-dependent processes of explosive gas mixtures formation in the three-dimensional space, their explosion and dispersion of combustion materials in the atmosphere as well as to create a computational modeling algorithm of these processes.
Mass of combustible participating in burning is determined for computational cells where the admixture concentration is in the range between minimal and maximal concentration limits of inflammability *+, ≤ ≤ *-(: .. = R(∆). &.. = R(∆).
= R(∆).
. = − .. − &...
The mass of the combustible not participating in the burning process is determined only for the computational cells where the admixture concentration > *-(:
The total mixture mass in the volume where the burning process occurs is determined for computational cells where the admixture concentration *+, ≤ ≤ *-(:
On the other hand, the total mixture mass m includes the masses of an oxidant m¢ , burning combustible m¢¢ and not burning combustible m0¢¢. Hence the oxidant mass in the mixture is:
The mass concentration (averaged on computation space volume) of mixture components are determined as follows:
.. = ..⁄; &.. = &..⁄; . = .⁄ = (1 − .. − &..).
The excessive air factor in the mixture a is as follows:
= .⁄(&..) = (1 − .. − &..)⁄(&..), where & = /.0⁄..is a stoichiometric number, /.0 – air mass which is necessary in theory for the complete combustion of 1 kg of the fuel.
The lower combustion value of the admixture Hu is set from the tables of thermophysical properties of the matters. The molar mass µc and the adiabatic coefficient kc of the combustion materials are determined on the basis of reversibility hypothesis of the realized chemical reactions.
In the case when ∝≥ 1 the thermophysical properties of the gas mixture after an explosion are determined as follows: = 1D[ 1 − (& + 1).. − &.. + (& + 1)..
. 1
+ .&...\,
(
The vector equation (
d + e g = d ,
2 = [1 − (& + 1).. − &..]2. + (& + 1)..21 + &..2.., (
1 + & − (& + 1).
..
\, 2 = _1 − 1⬚`2.. + 1⬚21, = 5/..0( − 1) + - = = 2⁄3.
5(1 − .. − &..)⬚.. ( − 1) &..
+ -, =
.
3 = _1 − 1⬚`3.. + 1⬚31, (20)
In any case such properties of the gas mixture as adiabatic coefficient, pressure, temperature and density are as follows:
5, = ⁄. (24)
From here, it is assumed that an explosion is instantaneous; combustion takes place in the
permanent volume occupied by an explosive mixture with the combustible concentration within
the limits of inflammability. After the explosion in the localized volume the fluid dynamics
parameters of two-reactant mixture (air and combustible gas) are modified to the fluid dynamics
parameters of the three-component mixture (air, products of combustion and fuel remnants).
(
A release of the gas inflammables into the atmosphere and their explosion cause the formation and propagation of the shock waves, personnel affection and damage of the vitally important objects. An overpressure in the front of the shock-wave generated by an explosion is usually used to evaluate building surface loadings.
The validation of the developed model against experimental results was carried out. An explosion of the hemispheric homogeneous stoichiometric hydrogen-air mixture cloud was modeled (experiments at Fraunhofer ICT [54]) under the following conditions: the total volume of the cloud 2094 m3; the initial pressure 98.9 kPa; the initial temperature 283 K; the diameter of the hemispheric cloud 20 m. The pressure development at the distances of 35 m (the point B on the fig. 2) and 80 m (the point C) from the epicenter of the explosion (the point A) was examined during the computation.
The computation space had the following dimensions: the length of 200 m; the width of 100 m, the height of 30 m. The dimensions of the computational grid were 200x100x30 cells. The computer and code had the following characteristics: processor Intel® Core(TM) i7–3630QM CPU @ 2.40 GHz 2.40 GHz, 16 Gb RAM, Windows 7, 0.2 h CPU time.
The flow pressure distribution near the ground at the moment when the shock-wave passes the point B is presented on the fig. 2. It is obvious that a zone of decreased pressure has been formed around the epicenter of the explosion. The further propagation of the explosion shock-wave along the computation area to the control point C is accompanied by decreasing of the intensity of overpressure (fig. 3).
The overpressure dynamics at the control milestones B and C is presented on the fig. 4-5 against the experimental results and the computational results obtained using other different codes [54]– [62]. More sharp form of the calculated curve can be explained by particular features of the accepted combustion model (global instantaneous chemical reaction). More intensive decrease of the overpressure as the shock-wave propagates from the point B to the point C may be obviously referred to the first order scheme of the Godunov method.
As a whole the computational results have quite a good agreement with the experimental data. It allows using the developed mathematical model and code for a simulation of the large-scale explosions of the hydrogen-air mixtures in the open atmosphere and prediction of pressure consequences of such explosions on surrounding buildings.
An explosion of stoichiometric propane-air mixture (the volume of combustible cloud – 1495 m3, energy of explosion – 4640 MJ) was also carried out in order to validate the developed model and code.
The computational results (fig. 6) show acceptable accuracy comparing with experimental data and regressive dependence [17] (in the fig. 6: & = ⁄!⁄8. is a dynamic radius, where R – distance from the epicenter of explosion, m; E – energy of explosion, J).
The hydrogen dispensing station [63] and local area around it is considered. The station has a cryogenic tank with liquid hydrogen (5,7 m3). The tank feeds 36 high-pressure cylinders (18,4 m3 in overall) with gas hydrogen. During operation of the fuel station different emergencies caused by gas and liquid hydrogen leaks from the defective equipment (or as a result of its destruction) can happen. It results in formation of the explosive cloud of the hydrogen-air mixture and its dispersion in the atmosphere. One of the most dangerous scenarios (taking into account potential catastrophic consequences threatening to the equipment of the station, its personnel and near-by residential area buildings) is an explosion of hydrogen-air mixture as a result of large-scale instantaneous release into the atmosphere of all the volume of compressed gas hydrogen from the all high pressure dispensing cylinders [63]–[65].
This scenario was modeled using the developed mathematical model. The hydrogen cloud had such initial parameters: volume 798 m3; hydrogen mass 687.4 kg; hydrogen mass concentration 100%; pressure 1031371.3 Pa; temperature 288 K. It was assumed that the cloud had been dispersing during 0.06 s after instantaneous destruction of high pressure cylinders (a physical explosion). Then a detonation explosion of the hydrogen-air mixture (a chemical explosion) took place. The chemical explosion caused appearance of high temperature combustion materials and shock wave which had a pressure impact on the station equipment and residential buildings. The location plan is presented in the fig. 7.
The height of the buildings was 5–12 m. The location of the station (positions A and B in the fig. 7) was varied in relation to the buildings location. The different types of the shield installations (banks, bumper wall) around the station against destructive influence of the shock wave overpressure were also considered. The pressure history was analyzed at the explosion epicenter and point on the building (A, B and C on the fig. 7, accordingly).
6.1. Hydrogen cloud explosion nearby infrastructure objects
The case when the fuel station location is close to the residential area is considered (position B in the fig. 7). Any shield installation from destructive influence of an explosion shock wave is used. The distribution of the hydrogen volume concentration in the just before an explosion is presented in the fig. 8. Obviously, the radius of a hemispheric zone of the detonation burning makes approximately 25 m and surpasses the height of the building.
The pressure distribution near the ground and in YOZ-plane when the pressure at point C reaches the maximal is presented in fig. 9. The analysis shows that pressure in the wave front on the building near the ground is approximately twice as higher as in open space.
The pressure history in the control points B and C is presented in the fig. 10. Two peaks of pressure in the point B (fig. 10 a) correspond to the moments of time when physical and chemical explosions occur. As the explosion epicenter is moved away from the residential area, the value of the shock wave pressure amplitude decreases quickly (fig. 10 b).
The case when the fuel station location is distant from the residential area is also considered (position A in the fig. 7). The distance from the fuel station to the nearest buildings is selected from the recommendations [63]–[65]. Any shield installation against the destructive influence of an explosion shock wave is used. The analysis of the hydrogen volume concentration distribution shows that the size and the form of the detonation burning zone are similar to the case with the near-by location of the station (fig. 8). The pressure distribution in the planes XOZ (near the ground) and YOZ when overpressure in the control point C reaches the maximum is presented in the fig. 11.
The pressure history in points A and C is presented in the fig. 12. Two peaks of pressure in the point A (fig. 12 a) correspond to the moments when explosion occur. The distancing of the explosion epicenter from the residential area significantly decreases (approximately five times) the maximal pressure loading on the building walls (fig. 12 b, 10 b).
A similar to 5.2 case of the distant location of the fuel station from the residential area is considered (position A in the fig. 7). To shield the buildings against the destructive impact of the explosion the banks (7 m high) surrounding the station are installed. The analysis of the hydrogen volume concentration distribution shows that the overall dimensions and the form of the detonation zone have been changed under the influence of complex relief of banks (fig. 13).
The pressure distribution in the cross-planes XOZ (near the ground) and YOZ (through the point C on the building) when the overpressure in the control point C reaches maximum (fig. 14) has changed insignificantly comparing to the remote explosion 5.2. 6.4. Distant partly banked explosion of the hydrogen cloud
A similar to 5.3 case of the distant location of the fuel station from the residential area is considered (position A in the fig. 7). To shield the buildings from the impact of an explosion the banks (7 m high) partly surrounding (in the north-east) the station are installed. Overall dimensions and a form of the detonation zone have been changed under the influence of the banks relief (fig. 15).
The pressure distribution in the planes XOZ (near the ground) and YOZ when the overpressure at the control point C reaches the maximum is presented in the fig. 16. The partly banking of an explosion site has insignificantly changed the overall pressure field comparing to the explosion 5.3. 6.5. Distant explosion partly surrounded with higher banks
A similar to 5.4 case of the distant location of the fuel station from the residential area is considered (position A in the fig. 7). To protect the buildings from an explosion impact the higher banks (13 m high) which partly (in the north-east) and more distantly surrounding the station are installed (fig. 17, 18 a). The pressure distribution in the planes XOZ (near the ground) and YOZ when the overpressure at the control point C reaches the maximum is presented in the fig. 18. The data analysis shows that the higher partly banking of an explosion site allows decreasing slightly the pressure loading on the buildings (fig. 18 b). 6.6. Distant hydrogen explosion with the use of bumper walls t, s
A similar to 5.2 case of distant location of the fuel station from the residential area is considered. To protect the buildings from an explosion impact the bumper walls (8 m high, 3 m thick) are installed immediately in front of the buildings (fig. 19 a). The pressure distribution in the planes XOZ (near the ground) and YOZ when the pressure at the control point C reaches the maximum is presented in fig. 19. The data analysis makes it clear that an installation of the bumper wall near the buildings causes a decrease of pressure by approximately 10 percent comparing with a case without protection (fig. 20).
It should be noted that all the considered types of the protective installations do not allow bringing maximal overpressure down in the control point on the wall of the building to the safe level.
The mathematical model of the gas-dynamics processes of the two-agent explosive gas mixture formation, its explosion and dispersion of the combustion materials in the open atmosphere was developed. The finite-difference approximation was developed for the case of three-dimensional system of the gas dynamics equations complemented by the mass conservation laws of the gas admixture and combustion materials. The algorithm of the computation of the thermo-physical parameters of the gas mixture resulting after instantaneous explosion taking into account the chemical interaction was developed. The smoothing effect of the first-order finite-difference scheme can be compensated by a denser computational grid. The verification of the mathematical model showed an acceptable accuracy in comparison with the known experimental data that allowed using it for the modeling of consequences of the possible failures at industrial objects which store and use hydrogen.
The computational modeling of the gas hydrogen explosion at the fuel station was carried out.
The analysis of the different ways of protecting the surrounding buildings from the shock wave
destructive impact was conducted. It was revealed that the considered types of the protective
installations (partial or complete banking, bumper walls) had an influence on the pressure
distribution in the computation area but did not allow bringing the maximal overpressure down to
the safe level. It was concluded that a bumper wall immediately in front of the protected object was
one of the most effective protective installation. It is necessary to take into account a
threedimensional character of the shock wave in order to select safe dimensions of the protection zone
around the hydrogen storage facilities.
[20] A. J. Laderman, Detonability of hydrogen-oxygen mixtures in large vessels at low initial
pressures, AIAA Journal4(
Combustion and Flame 101(
CRC Press, Boca Raton, 2012. doi: 10.1201/b12884.
[42] Y. Skob, V. Men’shikov, M. Ugryumov, Solution of the three-dimensional turbomachinery
blade row flow field problem with allowance for viscosity effects, Fluid Dynamics, 26(
doi:10.1038/1821134a0.
[48] M. E. Miller, G. C. Holzworth, An Atmospheric Diffusion Model for Metropolitan Areas,
Journal of the Air Pollution Control Association 17(
doi:10.1007/978-1-4615-0663-8. [50] V. Guinot, An outline of Godunov-type schemes, in: V. Guinot (Ed.), Godunov-type Schemes,
Elsevier, 2003, 93–116. doi:10.1016/B978-044451155-3/50005-0. [51] P.K. Sweby, Godunov Methods. In: E.F. Toro, (Ed.) Godunov Methods. Springer, New York,
NY, 2001. https://doi.org/10.1007/978-1-4615-0663-8_85. [52] B. Després, Lagrangian Godunov Schemes, in: G. Demidenko, E. Romenski, E. Toro, M.
Dumbser (Eds) Continuum Mechanics, Applied Mathematics and Scientific Computing: Godunov's Legacy, Springer, Cham, 2020. doi:10.1007/978-3-030-38870-6_16. [53] D. D. Knight, Elements of Numerical Methods for Compressible Flows, Cambridge University
Press, Cambridge, England, 2006. doi:10.1017/CBO9780511617447. [54] J. J. García, D. Baraldi, E. Gallego, A. Beccantini, A. Crespo, O.R. Hansen, S. Høiset, A.
Kotchourko, D. Makarov, E. Migoya, V. Molkov, M.M. Voort, J. Yanez, An intercomparison
exercise on the capabilities of CFD models to reproduce a large-scale hydrogen deflagration in
open atmosphere, International Journal of Hydrogen Energy 35(
Papanikolaou, S. Kumar, E. Studer, Y. Dagba, T. Jordan, W. Jahn, S. Høiset, D. Makarov, J.
Piechna, An intercomparison exercise on the capabilities of CFD models to predict distribution
and mixing of H2 in a closed vessel, International Journal of Hydrogen Energy 32(
Huser, W. Jahn, T. Jordan, J.-M. Lacome, H.S. Ledin, D. Makarov, P. Middha, E. Studer, A.V. Tchouvelev, A. Teodorczyk, F. Verbecke, M.M. Van der Voort, An inter-comparison exercise on the capabilities of CFD models to predict the short and long term distribution and mixing