On the Features of Numerical Simulation of Hydrogen Self-Ignition under High-Pressure Release

: The paper is devoted to the comparative analysis of different CFD techniques used to solve the problem of high-pressure hydrogen release into the air. Three variations of a contemporary low-dissipation numerical technique (CABARET) are compared with each other and a conventional first-order numerical scheme. It is shown that low dissipation of the numerical scheme defines better resolution of the contact surface between released hydrogen and ambient air. As a result, the spatial structures of the jet and the reaction wave that arise during self-ignition are better resolved, which is useful for predicting the local effects of high-pressure hydrogen release. At the same time, the dissipation has little effect on the induction delay, so critical conditions of self-ignition can be reliably reproduced even via conventional numerical schemes. The test problem setups formulated in the paper can be used as benchmarks for compressible CFD solvers.


Introduction
Nowadays, gaseous combustibles are widely used in many applications, especially those related to power engineering.Combustible gases can also be emitted in the course of various technological processes.For example, hydrogen can be produced as a result of the operation of renewable and nuclear energy systems.In view of this, gas explosions represent severe hazards, and approaches for risk assessments and explosion prevention are in demand.
Computational physics is a useful tool that helps to analyze possible emergency situations related to gas explosions and their consequences.Conventional detailed mathematical models include both gas dynamics and chemical kinetics that allow numerical simulation of all the peculiarities of explosion development [1].At the same time, among the variety of numerical techniques, one should choose the most appropriate one to resolve all the necessary features of the particular mode of gas explosion [2][3][4].
In this paper, we decide to focus on the problem of hydrogen self-ignition in the process of its spontaneous release from a high-pressure vessel [5].Gas release under high pressure leads to the formation of high-speed flow with shock waves and contact discontinuities [6].Herewith, the hydrogen ignites exactly at the contact surface where the diffusive and convective mixing of hydrogen with air heated by the shock wave takes place [7].Moreover, numerous experimental [8,9] and numerical [10,11] research works on hydrogen release through an orifice or a slit demonstrate a significant role of the gasdynamic structure of the contact surface.Thus, in [12], it is shown that an important role in hydrogen self-ignition belongs to the hydrogen jet evolution, its interaction with the boundary layer, and gas-dynamic instability, which gives rise to vortices at the contact surface and therefore intensifies hydrogen mixing with heated air.Much more complex patterns are observed when hydrogen is released into obstructed space, such as a channel with obstacles [13] or an open space with obstacles set up at a certain distance from the hydrogen source [14].Experimental visualization and numerical simulations show that a hydrogen jet in free space is also characterized by a complex spatial structure that can affect the self-ignition process.In view of this, there are certain demands for the utilized numerical technique."Dissipation-free" methods are considered advantageous since they possess low numerical dissipation, and high resolution of shock waves and contact discontinuities can be achieved.
The main goal of this paper is to assess the ability of the low-dissipation numerical technique CABARET [15,16] to reproduce well the features of hydrogen self-ignition when it is released into the air under high pressure.For these means, we propose several one-and two-dimensional tests related to high-pressure hydrogen release into the air and analyze the performance of three modifications of the CABARET technique.In addition, a more conventional numerical technique with first-order accuracy referred to as the "coarse particles method" (CPM) [17] is tested the same way.

Problem Setup
According to the main goals of the paper introduced above, let us consider the following problem setup.First let us propose three series of calculations describing the process of high-pressure hydrogen release into the air with different levels of detail: (1) one-dimensional non-reactive case, (2) two-dimensional non-reactive case, and (3) two-dimensional reactive case.The first series of calculations provides data on the flow structure in the most simplified case.Nevertheless, such tests give basic information on how the numerical dissipation affects the flow structure, especially in the vicinity of the contact surface between hydrogen and air.The second series provides information on the spatial structure of the jet in two dimensions.The results of the third series of calculations provide the complete pattern of the reactive flow development in the whole process of mixture preparation and its self-ignition.Here, two main results can be obtained on the ignition delay and the further development of the ignition kernel.The first one is important for the validation of numerical instruments with the use of corresponding experimental data on ignition delays and, further, for the estimation of the predictive ability of the instrument.At the same time, the data on the ignition kernel development after self-ignition would be useful for the hazardous potential evaluation of the hydrogen explosion.
In all three series of test calculations, hydrogen is stored initially under high pressure (from 50 to 700 atm) and starts to be released into ambient air (at 1 atm) after the diaphragm is ruptured instantaneously.The schematic of the two-dimensional numerical domain is presented in Figure 1.Planar geometry is used for one-dimensional calculations, while axisymmetric geometry is used for two-dimensional ones.The description of the setup parameters for the cases under consideration are mentioned in Table 1.

Mathematical Model of Reactive Gas Dynamics
The governing equations describing the reactive gas dynamics are classical Navier-Stokes equations that take into account compressibility, viscosity, thermal conductivity, multi-component diffusion, and chemical transformations [1].Their representation in the axisymmetric form is the following: where r, z-radial and axial coordinates, respectively, t-time, ρ-density, Y k = ρ k /ρ-k-th species mass fraction, p-pressure, (u r , u z )-mass velocity vector, ∂z -k-th species diffusion velocity components along the rand z-axis, respectively, E = ε + 1 2 u r 2 + u z 2 -specific total energy, ε-specific inner energy, ωk -chemical source term, κ(T)-thermal conductivity, D k (T)-diffusion coefficient of k-th species, σ ij -viscous stress tensor, h k -enthalpy of formation of k-th species, m k -k-th species molar mass.The radial spatial derivative in the left side of the equations above was converted as 1 r ∂ ∂r (r f ) = ∂ f ∂r + f r .Viscous stress tensor components are written as follows: where µ(T)-dynamic viscosity coefficient.The equation of state for a mixture of perfect gases is used with an account of the temperature dependence of the specific enthalpies and heat capacities according to the JANAF polynomials [18].The chemical kinetic mechanism of hydrogen oxidation, consisting of 19 elementary reactions between 8 elementary species (H 2 , O 2 , H 2 O, H, O, OH, HO 2 , and H 2 O 2 ) is taken from [19], where a complete kinetic mechanism for syngas oxidation is proposed.Multicomponent diffusion is calculated via zeroth-order Hirshfelder-Curtiss approximation [20].Mixture-averaged transport coefficients are obtained from the first principles of the molecular kinetic theory [21].

Description of CABARET Numerical Approach
A detailed description of the CABARET numerical algorithm can be found in [15,16,22] or in a monograph by Goloviznin et al. [23].Adaptation of the CABARET technique for combustion problems and analysis of the efficiency of the balance-characteristic form of CABARET on the problems of non-stationary combustion can be found in previous papers by the authors [2,24].Here, we provide general expressions of the CABARET procedure for the system of governing equations described above.
The solution region is split into rectangular cells with coordinates of nodes r i,j , z i,j .Following Goloviznin et al., we denote the cell centers by half-integer indices (i + 1/2, j + 1/2).Two types of variables are introduced: conservative variables and flux variables.Conservative variables ϕ = (ρ, u r , u z , E, Y k ) refer to cell centers.Flux variables ϕ f are defined at the spatial cell faces.Denoting the cell sizes as ∆r i+1/2 and ∆z j+1/2 and the temporal step between time levels t n and t n+1 as τ n+1/2 , the system (1)-( 5) can be written in implicit finite difference form: where U = U(ϕ)-conservative variables vector defined in cell centers, F = F ϕ f and geometric factor, and Q-source variables vector.The over-bar indicates the average over two time levels: η = 1 2 η n+1 + η n .Following the balance-characteristic formulation of the CABARET approach [15,25], the solution of the implicit Equation ( 7) is constructed using an explicit predictor-corrector algorithm, consisting of the following stages: 1.The predictor step, during which the conservative variables on the intermediate time level n + 1/2 are calculated.Central approximation forward in time is applied to the conservative form of the governing equations: 2. Extrapolation step, during which characteristic flux variables are extrapolated on the next time layer using a local one-dimensional characteristic form of the governing equations: where I r l (ϕ) and I z l (ϕ) are local Riemann invariants within the computational cell in the r and z directions, respectively, associated with the system of conservation Equations ( 1)-( 5), λ r l and λ z l are characteristic velocities of the l-th invariant in the r and z directions, and q k are source terms.Here, equation system (1)-( 5) is treated as Euler equations with source terms of a general form.In that case, Riemann invariants and associated characteristic velocities have the following expressions: In the r direction: In the z direction: where c is the velocity of sound c = γp/ρ.Riemann invariants for the next time step are calculated via linear extrapolation on the scale of a computational cell: After extrapolated values of the local Riemann invariants are obtained, they are corrected following the maximum principle [15].The final value of flux variables ϕ f at the cell faces is calculated from the corrected Riemann invariants.
3. At the corrector step, conservative variables on the next time step are calculated using new values of flux variables: 4. Finally, a source step is carried out, during which the source values are accounted for in conservative variables: The optimal time step between time levels t n and t n+1 is computed using the Courant stability condition: where CFL-Courant-Friedrichs-Lewy number, which regulates the stability of the scheme and should be less than 0.5 [23].

Artificial Dissipation Concept
Additional artificial dissipation can be introduced in the CABARET numerical algorithm at the extrapolation step (11) via an additional parameter σ [26]: where σ-artificial dissipation parameter variable within the range [0, 1].The value of σ = 0 corresponds to the absence of artificial dissipation, and the dissipative properties of the CABARET numerical algorithm become the same as those of the Godunov scheme at σ = 1.

Concept of the Riemann Solver
In [27], it was proposed that the Riemann solver can be used as an alternative way to calculate flux variables in the so-called "sonic points": cells where flow exhibits transition from subsonic to supersonic.The Riemann solver can be utilized in cases for which flows are characterized by high gradients in the parameter distribution.In the current paper, the Riemann solver is used to calculate the flow variables in the vicinity of contact discontinuity, where the steepest gradients of mixture composition and temperature are observed.Initial values are taken based on the values of conservative variables in adjacent cells.A schematic representation of the initial condition and the resulting solution is presented in Figure 2.
The algorithm for solving the Riemann problem implemented within the software package consists of the following stages:

•
Setting the initial conditions in the form of pressure (p L , p R ), density (ρ L , ρ R ), flow velocity (u L , u R ), and adiabatic index (γ L , γ R ) on both sides of the contact discontinuity; • Calculation of pressure and flow velocity at the contact boundary using the Newton-Raphson method [28]; • Calculation of parameters at the point corresponding to the initial position of the contact boundary based on analytical expressions for the speeds of propagation of disturbances and equations relating parameters in unperturbed regions, a compressed region, a region of steady flow, and on a rarefaction wave.
A more detailed description of the algorithm for solving the Riemann problem can be found in [29].In this study, the following condition was used as a criterion for activating the Riemann solver: where µ max and µ min -respectively, larger and smaller values among the molar masses of the mixture in neighboring cells, K-parameter greater than one.The implemented condition aims to track the contact discontinuity and activates the Riemann solver near the highest gradients of mixture composition.

One-Dimensional Non-Reactive Flow
Let us start with the analysis of a one-dimensional solution for the flow of non-reacting hydrogen released into the air under high pressure.Those calculations should provide important data on how the numerical technique reproduces the features of the flow in the vicinity of the contact surface where the self-ignition occurs.Herewith, the initial pressure in the hydrogen tank is varied in the range from 50 to 700 atm.Characteristic temperature and density profiles obtained for two different pressure values with the use of four numerical techniques are presented in Figure 3.One can see that the temperature value in the region of steady flow between the shock wave and contact discontinuity is predicted to be the same when using different numerical techniques.At the same time, the original CABARET provides much thinner contact discontinuity compared with the conventional numerical technique of first-order accuracy (CPM).That is a natural result that illustrates the role of numerical dissipation, which leads to the smoothening of the gas-dynamic discontinuities.The same effect is achieved when applying the artificial dissipation and Riemann solver concepts to CABARET (Figure 3). Figure 4 shows the evolution of the contact discontinuity thickness in time, and one can notice that the smoothening of the contact discontinuity becomes much more pronounced when using dissipative numerical approaches.So the effect of contact discontinuity smoothening can be treated as an artificial one and should be diminished.
At first sight, the low-dissipation numerical technique CABARET in its original formulation is the best choice for numerical reproduction of high-pressure hydrogen release into the air since it provides the lowest rate of contact discontinuity smoothening with time (Figure 4).However, another artifact arises in the vicinity of the contact surface-the temperature peak (Figure 3).The nature of this temperature increase seems to be related to the utilized detailed equation of state for the multicomponent mixture [30] and might be due to the disagreement between mass fluxes of mixture components and the energy fluxes, which results in erroneous temperature values calculated via table-approximated equation of state.The temperature is higher exactly at the contact surface between hydrogen and air.That can become an artificial trigger of hydrogen self-ignition, which should start first in the region where the ignition delay time achieves its local minimum (or the temperature achieves its local maximum).On the other hand, the solution is stable, and the amplitude of the temperature peak decays with time (Figure 5).Moreover, it is localized inside a narrow region of several numerical cells where the hydrogen concentration is rather low.Due to that, this feature should not affect the solution for the reactive medium crucially.Nevertheless, such fluctuations should be diminished, and both artificial dissipation and Riemann solver concepts can be fruitfully applied for this purpose (Figures 3 and 5).Based on the obtained results, the artificial viscosity concept is more efficient.It is important to note that the observed decrease in the temperature jump is accompanied by birth and rise of a density fluctuation in the same computational cell (Figure 6), which is a consequence of incorrect reproduction of the mixture composition there.Nevertheless, the effect of this artifact can also be treated as negligible since the density jump occurs in the low-temperature zone where ignition is not plausible (Figure 3).Conventionally, the artifacts can be reduced by using finer numerical grids.One can see that all the considered numerical approaches provide less increase in the contact discontinuity thickness (Figure 4) when the finer grid is used.Also, the temperature fluctuation at the contact discontinuity is reduced on the finer numerical grid (Figure 5).The results of the convergence tests are presented in Figure 7. Four numerical grids are used here: δx = 100 µm, 50 µm, 25 µm, and 12.5 µm.According to Richardson's extrapolation routine [31], all the solvers provide convergence in the considered range of numerical resolution.The convergence rates and the estimations of the limit value of contact discontinuity width for four solvers and for four values of hydrogen pressure are presented in Table 2.One can see that the CABARET solver provides a higher convergence rate, so the calculations can be carried out with the use of rougher numerical grids compared with conventional numerical techniques.Herewith, the implementation of the Riemann solver concept provides even faster convergence.

Two-Dimensional Non-Reactive Flow
Consider the axisymmetric two-dimensional non-reactive solution.That solution represents an under-expanded jet of pressurized hydrogen in the air.Herewith, a key role in the development of the jet belongs to the gas expansion in particular, leads attenuation of the shock wave.As well, one can observe the certain evolution of the spatial structure of contact discontinuity (Figure 8).The spatial structure of the jet is related to the formation and evolution of the toroidal vortex and the gas-dynamic instability of the side surface of the jet.Herewith, the large structures, such as the toroidal vortex, are reproduced well using both the conventional numerical approach (CPM) and the low-dissipation one (CABARET).At the same time, it is known that gas-dynamic instability starts its development from short-scale perturbations.So first, the short-scale perturbations rise linearly, and then the larger-scale modes arise due to non-linear effects.Due to this, the solution obtained with the use of a low-dissipation numerical technique looks much more natural since the effect of (numerical) dissipation is lower and short-scale perturbations are resolved (compare Figure 8a,b).As will be shown below, the spatial structure of the jet can play an important role in the development of hydrogen ignition, so the use of lower-dissipation solvers seems to be much more promising when reproducing two-and three-dimensional flows.However, the two-dimensional solution inherits the nonmonotonicity of the temperature profile in the region of contact discontinuity (Figure 9a).And more importantly, the high-temperature zone becomes wider with the increase in the numerical resolution.Nevertheless, as well as in the one-dimensional case, this effect can be reduced with the use of either artificial dissipation or Riemann solver concepts (Figure 9b).Herewith, artificial dissipation seems to be much more promising since it leads to less smoothening of the solution in the two-dimensional case.It is important to note that numerical dissipation (either inherent or additional) affects the two-dimensional solution stronger than the one-dimensional one.Thus, Figure 9 demonstrates that in the presence of volumetric expansion of the jet, the temperature behind the shock wave can be much lower when both the shock wave and contact discontinuity are additionally smoothened due to the dissipation.Such an underestimation of the temperature behind the shock wave can significantly affect both the ignition delay and the possibility of ignition in general.Therefore, one should accurately choose the parameters of the numerical technique in order to reduce such an artifact.Here, CABARET with artificial dissipation seems to be promising technique for two-dimensional calculations taking into consideration two abovementioned features: temperature non-monotonicity and artificial smoothening.Herewith, it should be noted that the obtained jet spatial structure, in that case, turns out to be similar to the original CABARET case (compare Figure 8a,c,d), so one can expect that the solution obtained with this solver would reproduce the original phenomenon with the highest accuracy.

Two-Dimensional Reactive Flow
Finally, let us consider the reactive two-dimensional case, which should reproduce the whole process, including the hydrogen release, mixing with air, and subsequent ignition.In the considered problem setup, the ignition takes place at the contact discontinuity, where, on the one hand, fuel and oxidizer are already mixed and, on the other hand, the temperature is high enough.Furthermore, the self-ignition starts exactly on the axis of the jet, where the effect of volumetric expansion is the weakest.The reaction wave propagates along the contact discontinuity, involving the whole volume of produced hydrogen-air mixture in the ignition process (Figure 10a,b).Here, two basic parameters of the process can be distinguished: ignition delay time and speed of reaction front propagation.The first one defines the time period during which the mixture is prepared and the ignition starts, while the latter one defines the spreading of the ignition along the whole surface of the jet.The whole temporal evolution of the reaction wave is demonstrated in Figure 10c.Here, the speed of the reaction wave is calculated in the laboratory reference frame, as shown in Figure 10a.
Figure 10c demonstrates that when using the original CABARET algorithm, the ignition process starts much earlier than in other cases under consideration.That is related to the artificial temperature jump, which is observed at the early stage of hydrogen release (see Sections 3.1 and 3.2).As soon as the ignition starts earlier, the surface of the expanding jet is smaller.Therefore, the ignition process spreads along the entire jet surface in a shorter time.When using other numerical techniques, including the CABARET versions with additional procedures such as the Riemann solver and artificial dissipation concept, the ignition starts later, and the reaction front lives longer (Figure 10c, Table 3).Herewith, the overall dynamics of the ignition process are resolved well with both the conventional low-order technique (CPM) and modified CABARET, with the only exception being that coarser numerical grids can be used when applying CABARET.When analyzing the overall evolution of the ignition process (Figure 10c), one can observe the deceleration of the reaction front.Herewith, there is no quenching of the reaction, and combustion is maintained with further expansion of the jet.That can be additionally demonstrated by the time-evolution of the volume and mass of the combustion products (Figure 11).Here, one can see that the rate of the increase in the mass of combustion products reaches a certain constant value (Figure 11b).That is related to the expansion of the combustion area with the expansion of the hydrogen jet.Those values obtained with the use of different numerical techniques are presented in Table 3.

Conclusions
When solving complex problems of gas dynamics in reacting media, the criteria for the choice of numerical techniques should be satisfied.Thus, when studying the process of pressurized hydrogen release into the air with its subsequent self-ignition, one should reproduce accurately both the gas-dynamic and kinetic features of the process especially when the hydrogen is released into obstructed spaces where the interaction between the contact surface and shock waves becomes much more important.To get recommendations on what CFD technique is more appropriate to resolve the process considered in this paper, three basic problem setups related to high-pressure hydrogen release into the atmosphere are proposed as tests for CFD techniques.All those problem setups are used to test three solvers based on the contemporary low-dissipation numerical scheme CABARET.It is demonstrated that when calculating the jet flow of pressurized hydrogen into the air, the use of the CABARET numerical technique provides a certain increase in the accuracy of contact discontinuity reproduction.At the same time, however, the non-monotonic behavior of the temperature can be observed at the early stage of the process when using CABARET.This solution feature becomes more pronounced when considering multi-dimensional cases and can lead to much earlier self-ignition of hydrogen.Corrections to the original numerical scheme can be applied to avoid this.Those corrections, in general, are related to the introduction of artificial dissipation.Nevertheless, artificial dissipation has a much smaller effect on the spatial structure of the jet compared with the numerical dissipation intrinsic to the conventional first-order numerical schemes.
The obtained results introduce the perspectives of using low-dissipation and low-order numerical schemes for risk assessments in hydrogen safety.One should expect the most appropriate solution for pressurized hydrogen self-ignition to be obtained when using the CABARET solver with the dissipator concept.At the same time, however, even low-order numerical techniques can be used sometimes when there is no need to resolve the whole spatial structure of the jet.But it is enough to get an answer to whether there is self-ignition under certain conditions and what is the ignition delay time.

Figure 1 .
Figure 1.Schematic of two-dimensional problem setup: 1-pressurized hydrogen, 2-ambient air, 3-non-permeable wall with a hole of a given diameter, 4-solid walls of the hydrogen tank, 5-open boundaries, and 6-axis of symmetry.

Figure 2 .
Figure 2. Riemann problem at the cell face between two computational cells.Dotted line is a new position of the discontinuity.

Figure 3 .Figure 4 .
Figure 3. Temperature (upper row) and density (lower row) in the vicinity of a contact discontinuity at a time instant of 20 µs.Lines with square symbols-original CABARET, lines with upper triangles-CABARET with Riemann solver, lines with lower triangles-CABARET with artificial dissipation, and lines with circles-CPM.Initial discontinuity position z = 300 mm.The calculation parameters are as follows: cell width δx = 100 µm; CFL = 0.25.

Figure 5 .
Figure 5. Temperature deviation temporal evolution.T m denotes the maximum temperature, and T s -analytical value of temperature behind the shock wave.Lines with square symbols-original CABARET, lines with upper triangles-CABARET with Riemann solver (K = 1.2), lines with lower triangles-CABARET with artificial dissipation (σ = 0.15), and lines with circles-CPM.Solid lines correspond to the data calculated with a spatial cell size of 100 µm, and dashed lines correspond to the data calculated with a spatial cell size of 25 µm.CFL = 0.25 for all cases presented.

Figure 6 .
Figure 6.Density deviation temporal evolution: ρ m denotes the maximum density, and ρ s -analytical value of density behind the shock wave.Lines with square symbols-original CABARET, lines with upper triangles-CABARET with solver (K = 1.2), lines with lower with artificial dissipation (σ = 0.15), lines with circles-CPM.correspond to the data calculated with a spatial cell size of 100 µm, and dashed lines correspond to the data calculated with a spatial cell size of 25 µm.CFL = 0.25 for all cases presented.

Figure 7 .
Figure 7. Contact discontinuity widths at time instant of 100 µs convergence plots for the case of a system with a higher pressure value of 50 atm.Lines with square symbols-original CABARET, lines with upper triangles-CABARET with Riemann solver (K = 1.2), lines with lower triangles-CABARET with artificial dissipation (σ = 0.15), and lines with circles-CPM.Dashed lines indicate values obtained by Richardson extrapolation.CFL = 0.25 for all cases presented.

5 Figure 9 .
Figure 9. Axial temperature profiles behind the leading shock wave at 10 µs obtained with original CABARET and CPM approaches (a) and with the use of artificial dissipation and Riemann solver concepts (b).CFL = 0.05 for all cases presented.

Figure 10 .
Figure 10.Temperature field at 25 µs illustrating the reaction wave propagation along the jet surface obtained with CPM (a) and original CABARET (b) methods.Black solid lines correspond to the HO 2 volume fraction of 0.01%.(c) History the reaction front speed for the case of hydrogen release under pressure of 350 atm.Black dashed lines indicate the ignition delay time.The calculation parameters are as follows: cell width δx = 100 µm; CFL = 0.05.

Figure 11 .
Figure 11.History of the the combustion products volume (a) and mass (b) for the case of hydrogen release under pressure of 350 atm.The calculation parameters are as follows: cell width δx = 100 µm; CFL = 0.05.

Table 1 .
Parameters of test calculations.

Table 2 .
Convergence rates (k) and estimated contact boundary width values (L c ).

Table 3 .
Ignition delay time and rate of combustion product mass growth.