Hydrogen is widely applied in the different industries. Particular danger of its application is conditioned by the large energy of combustion, fast transition from burning to detonation, and, as a consequence, powerful explosion of the cloud in the atmosphere after release. Hydrogen differs essentially from other explosive compressed gases and liquids, firstly, by very small density, and, secondly, by very low liquid stage temperature. The quantity of hydrogen participating in explosion is defined by conditions of its evaporation (in case of release and spillage in liquefied form), dispersion and mixing with air. These processes will be influenced by specified above properties of liquefied and compressed hydrogen [1]. Spilled liquids evaporate from spillage surface mixing with fresh air and forming dangerous mixtures. These physical processes substantially depend on shape of the spill spot [2], environment conditions, especially, wind speed [3], earth surface relief [4].
Usually, accidental release of dangerous liquids takes place after malfunction of storage or transportation equipment (Figure 1Figure 1Figure 1Figure 1). Released admixture mixes with fresh air and creates air-gas cloud with some gas concentration, and the mixture could explode. Accidental gas releases and could explosions cause significant damage to the environment and create such dangerous factors as shock-impulse loads on humans and building constructions [5] in hazard zones, toxic inhalation dozes (if admixtures are toxic) [2], and thermal radiation dozes from high-temperature combustion products [6] and fire flame fluxes [7]. Safety experts need all the information about distribution of dangerous factors around accident epicenter in order to assess the consequences for environment and provide protection measures to mitigate them [8] reaching acceptable risk levels and check if protection construction could withstand explosion loads without destruction [9]. Mathematical modeling of physical processes during accidental release and explosion of combustible gases (Figure 2), instead of full experiment, could significantly reduce the costs of risk assessment work. Atmospheric dispersion of gaseous admixtures is traditionally modeled under the assumption of Gaussian distribution of the admixture concentration and on the basis of corresponding analytical functional dependences [10]. But for neutral, heavy and light gases only the partial solutions are obtained, and such important factors as relief and gas compressibility under explosion condition can not be taken into account. The most adequate description of the physical processes of dispersion of chemically reacting gases is possible only using the Navier-Stokes system of non-stationary equations for compressible gas [11]. Currently, numerical simulation of turbulent flows is carried out by solving the Reynolds-Favreaveraged Navier-Stokes equations, supplemented by a model of turbulence [12]. However, most turbulence models do not describe with an equal degree of adequacy the various types of flows that can appear [13]. This is especially true for flows with intense flow breaks and/or large pressure and temperature gradients. That is why, more simple than Navier-Stokes equations model but sophisticated enough in order to take into account relief, flow compressibility, and to avoid difficulties to select adequate turbulence model is needed aiming to extract mass concentration admixture and overpressure distributions as hazardous risk factors.
In presented work an attempt to solve numerically a three-dimensional problem of the admixture dispersion in the atmosphere and explosion equally applied for different density gases (including hydrogen) have been made.
A mathematical model of the physical process of the emergency release and dispersion of a gas impurity into the atmosphere, its dispersion in the surface layer of the atmosphere with the formation of a gas-air mixture and its explosion allows us to obtain non-stationary fields of the mass concentration of the impurity and pressure in the accident zone. These distributions of hazardous parameters can be further used to calculate the fields of hazardous factors (excess pressure and impulse of the compression phase in the shock wave front) to assess the consequences of their impact on the environment (operating personnel and man-made facility structures) in order to develop recommendations for reducing the operation risks of high danger enterprises.
To evaluate the spacious distribution of hazardous admixture mass concentration during the release and dispersion processes the basic equations system of air-gas mixture movement in the near-Earth atmosphere layer is used based on Euler approach supplemented with admixture conservation equation. An instant explosion model is used to define overpressure fields and asses explosion consequences. Admixture mass concentration is controlled during the calculation process to localize the area where admixture concentration is within flammability limits, and in the moment of explosion air-gas mixture flow parameters are replaced on detonation parameters. After that explosion products concentration and pressure distributions are controlled to evaluate overpressure and impulse at the explosion wave front. This useful information can be used in future versions of the model to calculate human damage probability using probit analysis model.
The computational domain is a parallelepiped located in the right Cartesian coordinate system (Figure 3). It is divided into spatial cells which dimensions are determined by the scale of the characteristic features of the area (surface roughness, objects dimensions).
An adequate description of physical processes of nonreactive gas mixture with air and further admixture dissipation in the atmosphere (or enclosed ventilated premise) is possible only with the use of time-dependent Navier-Stocks equations for compressible gas. The limited resources of modern computers do not allow obtaining effectively the direct numerical solution of these equations. As a rule, numerical modeling of large space turbulent flows is carried out by solving of Reynolds-averaged Navier-Stocks equations [14] added by turbulence model [15]. However the majority of turbulence models do not describe with an identical degree of adequacy all various types of flows. It especially relates to flows with intensive stream separation and large-scale gradients of pressure and temperature. As a result of the flow structural analysis and decomposition of the full gas-dynamic mathematical model it is assumed that the convective mass exchange mainly influences on the process considered. Thus, for the description of processes of two-component gas mixture and dissipation it is enough to use the simplified Navier-Stocks equations received by the truncation of viscous members (Euler approach with source members).
The calculated space Ω is a parallelepiped located in the right-hand Cartesian system of coordinates (X, Y, Z) with the basis in plane XOZ (axis Y is directed opposite to the Earth's gravity). The calculated space is broken into spatial cells. The full system of the time-dependent equations describing three-dimensional two-component gas mixture flow looks like [16], [17], [18]:
where a, b, c, d, f represent the following vector-columns:
where t represents time, u, v, w are the components of air speed vector, P, -pressure and density, and E is full energy of a volume unit of gas mixture:
where e -internal energy of gas mass unit; components of the vector 𝑓 ⃗ -projections of the distributed volumetric sources; g -gravitational acceleration; 𝑒 𝑠 -intensity of a thermal emission in gas volume unit as a result of chemical reaction.
If consider leaked and evaporated explosive admixture [19] and combustion products as separate gases [20], the law of each admixture component transfer, taking into account a diffusion speed, looks like [21]:
where Q -relative mass density of an admixture (the ratio of gaseous admixture substance density to the mixture density); 𝜌 𝑄 𝑡 -an admixture density change rate as a result of diffusion (according to Fick law [22], 𝜌 𝑄 𝑡 = 𝑑𝑖𝑣(𝜌𝑄 𝐷 𝑔𝑟𝑎𝑑𝑄), and the factor of diffusion 𝑄 𝐷 is defined according to Berljand [23]); 𝜌 𝑄 𝑠 -an admixture density change rate as a result of chemical reaction.
The system of the equations (1, 8) is completed by the mixture component heat-transfer
properties equations [24]. For ideal gas the value of e is related to the values of P and of the mixture by the following dependence 𝑒 = 𝑃 ((𝑘 − 1)𝜌) ⁄ .
It is assumed that the air flow quantity component of velocity does not surpass sound speed. Entry boundary conditions are set on the finite-difference cells surfaces through which atmospheric air enters. The approaching flow is defined by values of total enthalpy
entropy function
flow velocity vector (angles 𝛼 𝑥 , 𝛼 𝑧 ), and relative admixture mass density Q (Q 1 if the gaseous admixture flows in). The entry flow parameters are defined by equations (3, 4) (if angles 𝛼 𝑥 , 𝛼 𝑧 are set) using "left" Riemannian invariant correlation [24]. On impermeable computational cells' surfaces the "no passing" conditions are satisfied: 𝑞 𝑛 = 0 where 𝑛 ⃗⃗ is a vector of normal to considered surface. Exit boundary conditions are set on the computational cells surfaces through which the mixture flows out (except for the atmospheric pressure 𝑃 𝐴 , the "right" Riemannian invariant correlation [24] is used).
At start time in all "gaseous" cells of the computational space the parameters of an ambient air are assigned. In cells, where the admixture cloud takes place, relative mass concentration of an admixture equals 1 (100%). In cells with hydrogen evaporation (or outflow) the law of admixture consumption variation is set.
The vector equation ( 1) is a consequence of the mass, impulse and energy conservation laws which can be presented in the integrated form for each calculated cell
where V is a volume of an elementary computational cell, 𝜎 ⃗-a limiting surface of the given cell which has an external normal vector 𝑛 ⃗⃗ (𝜎 ⃗ = 𝜎 𝑛 ⃗⃗)
, 𝐴 ̂ -a tensor of the flow density of conservative variables 𝑎 ⃗ which columns are vectors b, c and d, accordingly.
The mixture component transfer law ( 8) can be presented also in the integrated form for each computational cell:
The equations (11,12) suppose occurrence and existence of parameters break surfaces of two types: shock waves and tangential breaks. The functions, satisfying to the equations (11,12), can be considered as the gas dynamics equations generalized solutions. The use of integrated conservation laws as initial for construction of finite-difference equations provides the formation of breakup solutions without isolation of breaks.
The set of gas-dynamic parameters in all computational cells at the moment of time n t represents the known solution. Gas-dynamic parameters at the moment of time
are calculated by means of explicit finite-difference approximations for equations (5) according to Godunov method [24]. At first stage continuous distribution of parameters is replaced with piecewise constant integral-averaged values in each computational cell. At the same time the borders of a cell represent unstable surfaces of arbitrary breakup which disintegrate to steady wave elements: a shock wave, a contact surface and a wave of rarefaction. For each such breakup the streams of mass, impulse and energy through sides of gas cells are defined. The stability of the finite-difference scheme is provided due to a choice of time step size.
On the basis of mathematical model the computer system of the engineering analysis of the gas-dynamic processes of release, mixture and dispersion was developed which is used in the research computer information system. It allows predicting an admixture concentration change in time and space in the calculated area and computing shock-wave parameters formed after detonation in the atmosphere during fuel-air mixture dispersion process with the use of personal computers for practically reasonable time.
Since we did not possess hydrogen experimental data, we have used propane and methane experimental data to verify mathematical model and developed computer system.
In order to verify the developed model for gaseous admixture dispersion in the atmosphere the comparison of the computation results with experimental data [25] was conducted (Figure 4Figure 4Figure 4).
For validation of calculation model of atmospheric gas cloud explosions a comparison of calculated results with experimental data [26] іs conducted (Figure 5).
An explosion of a stoichiometric propane-air mix cloud was calculated at conditions of experiment: volume of the gas mixture cloud -1495 m 3 ; energy of explosion -4940 МJ. On the basis of these data initial conditions for calculation, pressure and temperature of combustion products in the cloud, have been defined:
where V -volume of the gas mixture cloud; E -energy of explosion; k -adiabatic coefficient; PA -atmospheric pressure; mix -molar mass of the mixture; mix -mixture density; R -universal gas constant.
For obtaining parameters of gas detonation, more sophisticated models [27] can be used to assess explosion hazards [28]. In work [26] on the basis of experiments of explosion of clouds of acetylene, propane and methane with air and propane and methane with oxygen the correlation dependence was received that allows to define the overpressure [26] and analitical correlation dependance [26] In Figure 5 the results of numerical calculations are presented which well enough coincide with experimental data and with correlation dependence [26].
Processes of gaseous hydrogen release (spilled liquid hydrogen evaporation), hydrogen-air mixture formation, explosion and further dispersion of hydrogen-air cloud in the atmosphere (taking into account the movement of the air, gravity, presence of buildings, and thermodynamic gas properties) are considered. Usually, the evaluation of safety measures of a hydrogen fueling station is provided using physical modeling [30], and then quantitative risk assessment can be done [31] to satisfy safety regulation documents [32].
Some possible scenarios at hydrogen dispensing station are modeled [29]. The station has the large cryogenic hydrogen tank (5.7 m 3 ) that feeds three 12-cylinder packages in total amount 799.2 m 3 in which gaseous hydrogen is stored at ambient temperature. The hydrogen under pressure is dispensed to vehicles' tanks. The station is separated from residential area by three zones.
• The region in the immediate vicinity of the equipment to protect the personnel from small leaks.
An exclusion zone in the immediate area of the facility (23 m) to provide protection against unplanned minor releases of hydrogen.
An additional margin, as large as 53 m is necessary to protect against large potentially catastrophic release of hydrogen. The most dangerous scenarios from the potential catastrophe point of view have been numerically simulated using the created mathematical model and developed computer system.
1. Spillage from the tank of all volume of liquid hydrogen, consecutive evaporation of cold liquid, mixture of gaseous hydrogen with moving air and the further dispersion of a mix towards residential area.
2. Release of all volume of the gaseous hydrogen compressed under a high pressure from cylinders of distribution with formation of a cloud and its further dispersion in a stream of air.
3. Explosion of a hydrogen-air cloud that is formed as a result of evaporation of spilled liquid hydrogen or instant gaseous hydrogen release.
According to the scenario 1, numerical computations of an instant spillage of all volume of liquid hydrogen from the cryogenic tank of dispensing station with formation of a spillage pond with total area of 65.7 m 2 were conducted. Liquid hydrogen evaporation productivity was evaluated of 11.56 kg/sec [29]. The wind blew with speeds of 1, 3 and 10 m/sec, and its direction was towards residential area. The calculation space has the following dimensions: the length along OZ axis -70 m, width (OX) -22 m, and height (OY) -20 m.
Hydrogen mass concentration distributions near the surface of the ground are presented on Figure 6. With an increase in wind speed the rate of dispersion increased too (Figure 7). According to comparison of hydrogen time-space concentration distributions at different wind conditions in the case of spilled liquid hydrogen its buoyancy is insignificant because of low vapor temperature and, as a consequence, small difference between densities of hydrogen and air. Mass concentration fields are very important data for safety experts to carry out following risk assessment. They can evaluate the total amount of hydrogen mass that potentially could explode (mass concentration has to be between flammability limits). This parameter characterizes the power of the explosion and the hazardous area where overpressure exceeds acceptable safe levels.
An instant release of all gaseous hydrogen from a package of dispensing high-pressure cylinders with the formation of a cloud in volume 799.2 m 3 have been numerically simulated according to the scenario 2. According to results of numerical experiments, buoyancy of hydrogen is essential only at low wind speeds (Figure 8). It should be noted that the fields of mass concentration of flammable impurities significantly depend not only on wind speed, but also on the presence of obstacles in the form of buildings in the actual space, which distort the flow velocity fields, which will subsequently affect the distribution of pressures in the zone of emergency release and With an increase of wind speed the air flow prevents the cloud movement up pressing it to the ground. Calculation results are similar to the results of the physical experiment [33] provided in the photographs of the hydrogen cloud plume, which arises as a result of evaporation from the liquid hydrogen spill spot under different wind conditions (Figure 9). Obviously, if wind speed is greater then hydrogen buoyancy, hydrogen behaves as a neutral gas. Both in the case of spilled hydrogen dispersion and dispersion of instantly released gaseous hydrogen, the presence of construction facilities and residential buildings essentially affect flow symmetry decreasing transfer and mixing processes in the space before the buildings and accelerating these processes between constructions.
In computations on dispersion the dependence of hydrogen mass in the cloud between top and bottom concentration limits of flammability was determined (Figure 10). The largest mass of hydrogen and, consequently, maximal explosion energy will be under low wind speed conditions. But, as evident from Figure 10, a, an explosive cloud with less hydrogen mass and larger energy of explosion can be located closer to residential constructions. In case of instant release of compressed gaseous hydrogen, the cloud is compact enough (Fig. 10, b) and it can be considered as spherical cloud. But, in case of presence of constant source of assigned productivity (hydrogen spillage evaporation), the hydrogen mixture that explodes according to the mathematical model occupies the space of complex shape (Figure 11) in the calculated area. The air-gas mixture with mass concentration that is lower than bottom flammability limit does not explode because of lack of fuel component. If mass concentration is greater than top flammability limit, the mixture will not explode too because of lack of oxidant. An increase of pressure takes place between residential constructions, as it would be expected. Overpressure affecting the walls of buildings (on the right side of the Figure 12, a) is about 28 kPa that allows suggesting the possibility of serious destructions of residential buildings and industrial facilities threatening to health or even life of industrial object personnel and residential area inhabitants.
The three-dimensional model of the release and explosion of gaseous admixtures in the atmosphere has been developed. The validation of the numerical results shows acceptable accuracy in comparison with known experimental data.
Numerical computations for dispersion of evaporated spilled liquid and instantly released gaseous hydrogen and for explosion of the hydrogen-air cloud generated during dispersion have been carried out.
Obtained results can be usually predicted at qualitative level but presented numerical computations have allowed making quantitative forecasting with no contradictory physical picture.
Numerical modelling of the harmful admixture dispersion in the atmosphere generated as a result of air-fuel mixture distribution in three-dimensional space with the use of developed computer information system is applicable to engineering calculations for different technological systems, including ones which work on liquefied and gaseous hydrogen.
Presented mathematical model can be used in more complex information system as a predictor of distribution of such a hazardous flow parameter as pressure during an accidental gas explosions. Calculated pressure field is an origin of such hazardous factors values as maximum overpressure and impulse of pressure phase at the front of the explosion wave which moves away from an accident epicenter and influences harmfully on humans and buildings. This information about the explosion wave can be used to evaluate distribution of damage probability and to assess the risks of dangerous industrial objects by safety experts.