Calculations of Flows with Deflagration Appearance in Channels with Obstacles

Abstract

Many researches devoted to deflagration appearance in hydrogen-air mixes and transition it to detonation. This is connected with the problem of detonation hydrogen engine construction. A new way of detonation engine construction connected with using different fuel mixes: CO + H₂ or kerosene-hydrogen mixes. In these cases, it seems to be preferable for mathematical modeling to use fool system of reactions for all components and radicals. Using of fool system of kinetic equation allows one to determine moment of deflagration initiation and transition it to detonation. In laboratory experiments deflagration initiation appears as a result of energy and temperature grows by successive triggering of electric discharges or near obstacles in channels. In many papers numerically investigated flows of reactive gas mixes in channels with obstacles or initiation of deflagration via hot spot as a model to electric discharges. In this chapter numerically investigated flows of reactive gas mixes in channels with obstacles with the purpose to define conditions of deflagration initiation.

1. Introduction

Several last decades idea of detonation engine was in center of attention of researchers. Detonation fuel cycle is more energetic preferable than used in aviation and rocket engine constructions Briton cycle. So, simplicity of construction and energetic effectiveness make development of detonation engine constructions very perspective.

Detonation engine construction can be used in electricity if not paid attention to size and weight of engine. An example of such energetic construction was made in Bussing [1].

During the last three decades, some experimental laboratory constructions were made [2, 3, 4] for investigation of engine constructions, and works on mixing methane-air, methane-oxygen, and kerosene-air as fuels. All these fuel mixes have low detonation possibility, which demands additional energetic income for detonation initiation. This initiation can be achieved by electric discharges [5] or by specification of region inside engine: channels with obstacles, ring or U-tube form of the channel [6], focusing of pulse jet inside sphere-shape resonator [7] or spin detonation engine, in which detonation wave propagates inside axis symmetrical region in radial direction [8].

Progress in decision of problem of detonation engine construction partly consists of numerical simulation of gas flow inside the engine. This numerical simulation demands development of sophisticated mathematical models and numerical methods of high order of accuracy.

In the first period simplified model of chemical reactions of Levin et al. [9] was used, in papers of Fujiwara [10] for example. Nowadays various full systems of chemical reactions are used, which consist of many reactions ([11] for example). In the second section of this chapter mathematical model for numerical simulation of reactive gas mixes flows is described. Also in this section numerical method of solving system of gas dynamic equations for an arbitrary number of mixed components is introduced. This numerical method is based on a variant of Shakravarthy-Osher non-linear TVD (total variation diminishing) difference scheme of third order of accuracy [12].

In the system of gas dynamics (2) exist source term of speeds of changing mix components concentrations as result of chemical reactions. Values of this term determine growth of internal energy and temperature, so defining of this term is essential part of mathematical model and numerical method. Appearance of deflagration and transition from slow deflagration to detonation in reactive gas mixes flows depend on how intensive are chemical reactions in flow.

In the third section of this chapter mathematical model of a system of ordinary differential kinetic equations for the concentration of gas mix components is described. Numerical solutions of this system of equations are discussed.

The process of detonation appearance in hydrogen-air mixes can be treated as a transition from slow chemical reactions to fast branching chain reactions in hydrogen-air (or other components) mixes. The theory of these reactions was developed by Semenov [13] (see also Denisov et al. [14]). Now this approach intensively developed in the papers of some authors: Saeid et al. [15], Liu et al. [16].

2. Numerical simulation of reactive gas mixes flows

For numerical simulation of reactive gas mixes flows the system of equations of an ideal gas for two-dimensional flows (with source terms which are velocities of changing gas mix components via chemical reaction) can be written in integral form for finite volume cell:

$$d/\mathit{dt}\underset{V}{\int }\mathit{QdV}+\underset{S}{\oint }\mathit{\text{nFdS}}+\mathrm{\Phi }=0$$E1

where $Q=\left(\rho ,\overrightarrow{m},\mathit{\rho e},\rho c_i\right),i=1,\dots n$, is vector of m conservative unknown, n is number of mix components, $\overrightarrow{m}=\left(\mathit{\rho u},\mathit{\rho v}\right)$ is vector of impulse, and $c_i={\rho }_i/\rho$ are mass concentrations of mix components (where i = 1, …, n is number of mix component). The source term has the form:

$$\mathrm{\Phi }={\left(0,0,0,0,\mathit{d\rho }{c}_i/{\mathit{dt}}_{\mathit{\text{chem}}}\right)}_{.}$$E2

Tensor of flows is equal to: $\tilde{F}=\left(\overrightarrow{m},\left(\overrightarrow{m}\cdot \overrightarrow{m}\right)/\rho +P\tilde{I},\overrightarrow{m}\left(e+P\right)/\rho ,\overrightarrow{m}{c}_i\right)$.

For implementation non-linear TVD difference scheme, it is necessary to define flows on the border of finite volume cell on the basis of some approximate Riemann solvers. In this chapter those approximate Riemann solvers where defined (initially find in Martyushov [17]) for arbitrary number of gas mix components on the basis of Roe-Pike method, explanation of which can be found, for example, in Toro [18]. Flows on the border of finite volume cell for approximate Riemann solvers are:

$$F_L+\sum _{{\tilde{\lambda }}_k\le 0}{\mathrm{\Delta }}_{i+1/2}{W}_k{\tilde{\lambda }}_k{r}_k=F_R+\sum _{{\tilde{\lambda }}_k\ge 0}{\mathrm{\Delta }}_{i+1/2}{W}_k{\tilde{\lambda }}_k{r}_k\mathrm{\Delta }F=\sum _{k=1}^m{\mathrm{\Delta }}_{i+1/2}{W}_k{\tilde{\lambda }}_k{r}_k,$$E3

where $\mathrm{\Delta }{W}_k$ intensity of characteristic wave. For calculation of flow on the bound of control volume we have the next relations:

$$F_{i+1/2}=\frac{1}{2}\left(F_L+F_R\right)-\frac{1}{2}\sum _{k=1}^m{\mathrm{\Delta }}_{i+1/2}{W}_k\mid {\tilde{\lambda }}_k\mid r_k$$E4

where $r_k{I}_k$ are right and left eigen vectors of Jacobian matrix $A\left(Q\right)=\partial F/\partial Q$. $\mathrm{\Delta }W=L\mathrm{\Delta }Q$. For calculation of $\mathrm{\Delta }W$ it is sufficient to use this equality, where $\mathrm{\Delta }Q={\left(\mathrm{\Delta }\rho ,\mathrm{\Delta }\left(\mathit{\rho U}\right),\mathrm{\Delta }\left(\mathit{\rho V}\right),\mathrm{\Delta }E,\mathrm{\Delta }\left(\rho c_1\right),\dots \mathrm{\Delta }\left(\rho c_n\right)\right)}^T,i=1,n.$ where right eigen vectors are:

$$\begin{array}{lll}{r}_1& =& {\left(0,0,1,V,0,\dots ,0\right)}^T,r_2={\left(1,U,V,\omega ,0,\dots ,0\right)}^T,r_3={\left(1,U+a,V,H+\mathit{Ua},c_1,\dots ,c_n\right)}^T,\\ r_4& =& {\left(1,U-a,V,H-\mathit{Ua},c_1,\dots ,c_n\right)}^T,r_5={\left(0,0,0,h_1,1,0,\dots ,0\right)}^T,\\ r_6& =& {\left(0,0,0,h_2,0,1,0,\dots ,0\right)}^T,r_n={\left(0,0,0,h_n,0,0,\dots ,1\right)}^T.\end{array}$$E5

The next eigen values where chosen: $\lambda =\left(U,U,U+a,U-a,U,\dots ,U\right)$.

Denoting $\mathrm{\Delta }W=\left({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,\dots ,{\alpha }_n\right)$ for finding $\mathrm{\Delta }W=L\mathrm{\Delta }Q$ we can solve system of equation for finding $\mathrm{\Delta }\mathit{WR}=\mathrm{\Delta }Q$:

$$\begin{array}{lll}{\alpha }_2& +& {\alpha }_3+{\alpha }_4=\mathrm{\Delta }\rho \\ U{\alpha }_2& +& \left(U+a\right){\alpha }_3+\left(U-a\right){\alpha }_4=\mathrm{\Delta }\left(\mathit{\rho U}\right)\\ {\alpha }_1& +& V{\alpha }_2+V{\alpha }_3+V{\alpha }_4=\mathrm{\Delta }\left(\mathit{\rho V}\right)\\ V{\alpha }_1& +& {\mathit{\omega \alpha }}_2+\left(H+\mathit{Ua}\right){\alpha }_3+\left(H-\mathit{Ua}\right){\alpha }_4+\sum _{i=1}^n{\alpha }_6{h}_i=\mathrm{\Delta }\left(E\right)\\ c_i{\alpha }_3& +& c_i{\alpha }_4+{\alpha }_i=\mathrm{\Delta }\left(\rho c_i\right),i=1,n.\end{array}$$E6

System (6) can be solved by consequent excluding of unknowns:

$$\begin{array}{lll}{\alpha }_2& =& \left(\mathrm{\Delta }\rho \left(H-U^2-\sum _{i=1}^n{c}_i{h}_i\right)+\mathrm{\Delta }\left(\mathit{\rho U}\right)U+\sum _{i=1}^n{h}_i\mathrm{\Delta }\left(\rho c_i\right)-\overline{E}\right)/\omega \\ {\alpha }_4& =& \left(\left(U+a\right)\mathrm{\Delta }\rho -\mathrm{\Delta }\left(\mathit{\rho U}\right)-a{\alpha }_2\right)/2a\\ {\alpha }_3& =& \left(-\left(U-a\right)\mathrm{\Delta }\rho +\mathrm{\Delta }\left(\mathit{\rho U}\right)-a{\alpha }_2\right)/2a\\ {\alpha }_1& =& \mathrm{\Delta }\left(\mathit{\rho V}\right)-V\mathrm{\Delta }\rho \\ {\alpha }_{i+4}& =& \mathrm{\Delta }\left(\rho c_i\right)-c_i\left(\mathrm{\Delta }\rho -{\alpha }_2\right),i=1,n\\ \overline{E}& =& \mathrm{\Delta }E-V\left(\mathrm{\Delta }\left(\mathit{\rho V}\right)-V\mathrm{\Delta }\rho \right).\end{array}$$E7

The TVD difference scheme of Chacravarthy-Osher [12] was used for solving system (1). This scheme is one of third order of accuracy and can be written in such form (for flow on the border of finite volume in (4)):

Preliminary step:

$$\begin{array}{lll}{\alpha }_{1,m+1/2}^i& =& l_{m+1/2}^i\left(U_m-U_{m-1}\right);{\alpha }_{2,m+1/2}^i=l_{m+1/2}^i\left(U_{m+1}-U_m\right);{\alpha }_{3,m+1/2}^i=l_{m+1/2}^i\left(U_{m+2}-U_{m+1}\right);\\ {\lambda }_{m+1/2}^{i+}& =& \left\{\begin{array}{c}{\lambda }_{m+1/2}^i,.{\lambda }_{m+1/2}^i>0\\ 0,.{\lambda }_{m+1/2}^i<0\end{array}\right.,{\lambda }_{m+1/2}^{i-}=\left\{\begin{array}{c}{\lambda }_{m+1/2}^i,.{\lambda }_{m+1/2}^i<0\\ 0,.{\lambda }_{m+1/2}^i>0\end{array}\right.;\end{array}$$E8

$$\begin{array}{l}{\overline{\alpha }}_{1,m+1/2}^i=\text{min mod}\left\{{\alpha }_{1,m+1/2}^i,b{\alpha }_{2,m+1/2}^i\right\};{\tilde{\alpha }}_{2,m+1/2}^i=\text{min mod}\left\{{\alpha }_{2,m+1/2}^i,b{\alpha }_{1,m+1/2}^i\right\};\\ {\overline{\alpha }}_{2,m+1/2}^i=\text{min mod}\left\{{\alpha }_{2,m+1/2}^i,b{\alpha }_{3,m+1/2}^i\right\};{\tilde{\alpha }}_{3,m+1/2}^i=\text{min mod}\left\{{\alpha }_{3,m+1/2}^i,b{\alpha }_{2,m+1/2}^i\right\}\\ \text{where min mod}\left(x,y\right)=sgn\left(x\right)max[0,min(\left|x\right|,ysgn\left(x\right)].\end{array}$$

Main step:

$$\begin{array}{lll}{F}_{m+1/2}& =& H_{m+1/2}+\sum _i\left(\frac{1+\delta }{4}{\tilde{\alpha }}_{1,m+1/2}^i+\frac{1-\delta }{4}{\overline{\alpha }}_{2,m+1/2}^i\right){\lambda }_{m+1/2}^{i+}{r}_{m+1/2}^i\\ & & -\sum _i\left(\frac{1+\delta }{4}{\overline{\alpha }}_{2,m+1/2}^i+\frac{1-\delta }{4}{\tilde{\alpha }}_{3,m+1/2}^i\right){\lambda }_{m+1/2}^{i-}{r}_{m+1/2}^i\end{array}$$E9

The first term in (9) is:

$$H_{m+1/2}=\frac{1}{2}\left[F\left(Q_{m+1}\right)+F\left(Q_m\right)\right]-\frac{1}{2}\sum _i\left({\lambda }_{m+1/2}^{i+}-{\lambda }_{m+1/2}^{i-}\right){\alpha }_{2,m+1/2}^i{r}_{m+1/2}^i$$E10

Parameter $b=\left(3-\delta \right)/\left(1-\delta \right)$, and parameter $\delta$ determine order of approximation of the scheme. Meaning $\delta$ = 1/3 gives to scheme third order of approximation.

Time approximation was made on the basis of Runge-Kutta scheme of fourth order of accuracy.

Implementation of the Roe-Pike method for calculation of flows of reactive gas mixes flows was used by author in Martyushov [17]. Modification of numerical method [17] in this chapter consists of using meaning of parameter $\delta$ = 1. This meaning of $\delta$ preserve pure upwind variant of this scheme. This variant of scheme made it possible to avoid oscillations of values of radicals components concentration in regions of deflagration appearance.

The velocity of gas mix in gas dynamic flows is of order of sound speed 300 m/s in air. Pressure in the CI dimension system is of 10⁵ pascal. For the aim of minimizing truncation errors procedure of transition to non-dimensional unknowns was provided. Non- dimensional gas dynamic values were used: $P^{\prime }=P/\left({\rho }_0{c}_0^2\right),{\rho }^{\prime }=\rho ,{\overrightarrow{V}}^{\prime }=\overrightarrow{V}/c_0^2,T^{\prime }=T/c_0^2,t^{\prime }=t/c_0$—time of process is of $sec\times {10}^{-5}$ order.

For this reason, in Section 5 of this chapter physical parameters in results have dimension only temperature is in Kelvin degrees (pressure can be found from equation of state for ideal gas mix), time of process is of order $sec\times {10}^{-5}$ and concentration of components are measured in percent of density or molar units.

3. Kinetic model

For process of deflagration in reactive gas mixes the most important term in system (1) is $\mathit{d\rho }{c}_i/{\mathit{dt}}_{\mathit{\text{chem}}}$ (in source tern $\mathrm{\Phi }$) which determine, appearance of deflagration and transitions it to detonation, and as a result intensive growth of energy and temperature. The algorithm for calculation of this term will be described in this section of the chapter.

Equations of chemical reactions can be presented as follows:

$$\sum _{i=1}^n{\alpha }_{\mathit{ij}}{A}_i=\sum _{i=1}^n{\beta }_{\mathit{ij}}{B}_i,j=1,\dots ,M$$E11

where n is number of components of the mix and M is number of reactions, coefficients ${\alpha }_{\mathit{ij}}$, ${\beta }_{\mathit{ij}}$ are numbers of molecules (or moll) of mix components in direct and inverse reactions, $A_i$ is reagents and $B_i$ is products of direct reactions. Arrhenius low was used for calculating speeds of changing of mix components concentration $c_i$:

$$\frac{{\mathit{dc}}_i}{\mathit{dt}}=\sum _{j=1}^M\left({\beta }_{\mathit{ij}}-{\alpha }_{\mathit{ij}}\right)w_j\left(c,T\right),$$E12

$$w_j\left(c,T\right)=k_f\left(T\right)\prod _{i=1}^n{c}_i^{{\alpha }_{\mathit{ij}}}-k_b\left(T\right)\prod _{i=1}^n{c}_i^{{\beta }_{\mathit{ij}}},$$E13

$$k_{f,b}=A_{f,b}{T}^{l_{f,b}}exp\left(-E_{f,b}/\mathit{RT}\right)),$$E14

where meanings of ${\alpha }_{\mathit{ij}}$, ${\beta }_{\mathit{ij}}$ are taken from (11). In the present chapter gas mix of nine component: Н₂, О₂, Н, О, Н₂О, ОН, НО₂, Н₂О₂, N₂ was treated. The next 11 reactions were chosen:

Two algorithms were used for numerical solving of system (12)–(14). The first one was solving of fool system of stiff ODE (12)–(14) on the basis of implicit numerical Gear method. Constants of reactions in (12)–(14) were taken from Ibragimova et al. [19]. The other meanings of this constant can be found in Refs. [11, 14].

The model problem of initiation of deflagration as a result of sudden temperature increase in hot spot in closed volume was numerically simulated. In Figure 1(a) and (b) (initiation of deflagration as result of temperature meaning increasing to 2.4 and 2.7 times from initial atmospheric one (293 K degree) in part of closed volume), graphics of molar concentration of mix components are drawn, consequently: line 1 is graphic of Н₂О molar concentration, line 2 is graphic of Н₂ molar concentration, line 3 is graphic of О₂, line 4 (line up to the mark “4”) is graphic of 1/Т (inverse to temperature), line 5 is graphic of H₂О₂ molar concentration, line 6 is graphic of H molar concentration, and line 7 is graphic of HO₂ molar concentration (graphics 5–7 multiplied by 10²). On x-axes is time in seconds, and on y-axes is concentration of components in unity of part of mixed mass. In Figure 1(b) graphics 1–5 have the same meanings, line 6 is graphic of НО₂, and line 7 is graphic of OH molar concentration (graphics 5–7 multiplied by 10²).

Figure 1.

(a)–(c) Graphics of molar concentrations of mixed components from time (time is measured in seconds, meanings varies from 1 to 30, and concentrations measured in part of mass units).

Intensive changes of О₂ and Н₂ (that means deflagration initiation) begin after grows of H₂О₂ value, so that means that the influence of radical H₂О₂ appearance is important just in the same degree as influence of radical H appearance for initiation of deflagration.

Idea of the second calculation algorithm for the same problem is the next [13, 14]. In process of deflagration appearance in hydrogen-air gas mix sudden explosion appears after long period of induction. In this induction period concentrations of radicals Н, О, HO_2, H₂О_2, and ОН grows. This explosion mechanism is branching chain reactions. Theory of branching chain reactions was introduced by Semenov [13]. In this theory system of two differential equations for concentration of H and H₂O₂ (in Semenov [13] only one equation for H) is solving separately, concentration of radicals H, O, OH, and HO₂ are found from algebraic relation from assumption of “quasi-stationary” concentration, and for “large unknowns” system of three differential equations is solved with bigger time step. System of two differential equations for concentration of H and H₂O₂ take s form:

$$\begin{array}{l}d\left[\mathrm{H}\right]/\mathit{dt}=2k_{\mathit{0}}\left[{\mathrm{H}}_2\right]\left[{\mathrm{O}}_2\right]+\mathrm{Ф}\left[\mathrm{H}\right]–k_7{\left[\mathrm{H}\right]}^2\\ \mathrm{Ф}=2k_2\left[{\mathrm{O}}_2\right]-k_5\left[{\mathrm{O}}_2\right]–k_6\left[{\mathrm{H}}_2\mathrm{O}\right]–k_{11}\left[{\mathrm{H}}_2{\mathrm{O}}_2\right].\\ d\left[{\mathrm{H}}_2{\mathrm{O}}_2\right]/\mathit{dt}=k_8{\left[{\mathrm{HO}}_2\right]}^2-k_9\left[{\mathrm{H}}_2{\mathrm{O}}_2\right]\left[\mathrm{M}\right]+k_{10}{\left[\mathrm{OH}\right]}^2\left[\mathrm{M}\right]-k_9\left[{\mathrm{H}}_2{\mathrm{O}}_2\right]\left[\mathrm{H}\right].\end{array}$$E15

Molar concentrations of “quasi -stationary” components [O], [HO₂], and [OH] can be found in algebraic equations:

$$\begin{array}{l}\left[\mathrm{O}\right]=\left[\mathrm{H}\right]\left(k_2\left[{\mathrm{O}}_2\right]\right)/\left(k_3\left[{\mathrm{H}}_2\right]\right);\\ \left[\mathrm{OH}\right]=(2\left[{\mathrm{O}}_2\right]\left[{\mathrm{H}}_2\right]k_0+k_6\left[{\mathrm{H}}_2\mathrm{O}\right]\left[\mathrm{H}\right]+2k_9\left[{\mathrm{H}}_2{\mathrm{O}}_2\right]\left[\mathrm{M}\right]\\ +k_{11}\left[{\mathrm{H}}_2{\mathrm{O}}_2\right]\left[\mathrm{H}\right])(k_1\left[{\mathrm{H}}_2\right]+2k_{10}\left[\mathrm{M}\right];\\ {\left[{\mathrm{HO}}_2\right]}^2=k_5\left[\mathrm{H}\right]\left[{\mathrm{O}}_2\right]\left[\mathrm{M}\right]/\left(2k_8\right).\end{array}$$E16

The other components of the mix are found from equations for the rest unknowns of system (12)–(14).

The same model problem was numerically simulated using branching chain reaction theory (15) and (16).

Results of calculation 0n the basis of this algorithm are shown in Figure 1c, namely graphics of mass concentration of components (initial temperature meaning is increasing to 2.5 times from initial atmospheric one (293 K degree) in part of closed volume), line 1 is graphic of Н₂О mass concentration, line 2 is graphic of Н₂ mass concentration, line 3 is graphic of О₂ and line 4 is graphic of H₂О₂ multiplied by 10².

Calculations of systems of ODE (12)–(14) and (15)–(16) are provided by Goer difference scheme with time scale of order sec. Calculation of system of gas dynamic (4) by algorithm (9) provided with time scale of order $sec\times {10}^{-5}$. For reason of different time steps in whole algorithm (gas dynamic plus kinetic), it demands to provide 20 gas dynamic steps for one chemical step.

4. Construction of calculation grids

Curvilinear structured grids, used in numerical simulation of this chapter, were constructed by algorithm of author [20] and consist of solution of system of PDE equation of parabolic type. Construction of such grids is not a simple problem, sometimes calculation grid, constructed absolutely writes theoretically, produces oscillations, and wrong results in gas dynamics calculation. Examples of curvilinear grids for axis symmetrical channels with obstacles are shown in Figure 2(a)–(с).

Figure 2.

(a)–(c) Calculation grids for axis symmetrical channels with obstacles, near which initiation of deflagration appears (every tenth coordinate line is drawn).

5. Numerical simulation of gas mixes flows with deflagration

Numerical simulations of two-dimensional axis symmetrical flows of hydrogen-air mixes with deflagration were provided for some configurations of axis symmetrical channels with obstacles.

Due to high speed of gas dynamics flows and chemical reactions for mathematical model (1) normalizing performed similar to ones for Reynolds number definition for viscous liquid flows performed. The next characteristic values were used for normalizing. So non-dimension meanings for time, linear scale, and physical parameters were used. On graphics lower only temperature was in Kelvin degree and concentration of components of mixes in molar units. For other physical parameters (include time scale and linear sizes of calculation region), non-dimension parameters were taken in the next form:

$$t^{\prime }=t\cdot l/c_0,P^{\prime }=P/c_0^2{\rho }_0,E^{\prime }=E/c_0^2,$$E17

where ${\rho }_0$ is density of atmosphere gas, c₀² is square of atmosphere sound speed, and l is character size of calculation region. In non-dimension time scale calculation were provided for 3000–5000 time step (in non-dimension time units till 200–300 units).

Formulation of the mathematical problem for all geometry regions (Figure 2(a)–(c)) is the next. Initially calculation region was filled by hydrogen-air mix with a shock wave and the region after it near the left boundary of the region. On the left boundary of calculation regions boundary conditions with parameters of gas mixes behind the shock wave were established. On the top boundary of calculation region non-penetration conditions were used, on the bottom boundary conditions of axes of symmetry or non-penetration conditions were established. On the right boundary non-reflection conditions were used. In right part of the region before the shock wave meanings of gas dynamics parameters were atmospheric ones.

For all flows calculated in this section stoichiometric hydrogen-air composition of mix (20% of H₂, 10% of О₂, and 70% of N₂₎ was used in whole region with small (less than 2%) include of water as catalyst.

5.1 Flow in non-reflecting nozzle

Axis symmetrical channel with non-reflecting nozzle (calculation grid drawn in Figure 2(b)) initially filled with hydrogen-air mix with little amount of water (less than 2%, it involved as collision partner in reactions 5 and 7 of Table 1) with atmosphere meanings of gas dynamics values. At the initial moment shock wave of intensity M_sh.w. = 1.9 (Mach number of sh0ck wave is the relation of shock wave speed to sound speed before shock) began to move from the left side of the region.

H2 + O2 – k0 – 2OH	H2 + OH – k1 – H + H2O	H + O2 – k2 –O + OH
H2 + O – k3 – H + OH	HO2 + M – k5 –H + O2 + M	H + H2O – k6 – H2 + OH
2H + M – k7 – H2 + M	2HO2 – k8 –H2O2 + O2	H2O2 + M– k9 –2OH + M
2OH + M – k10 –H2O2 + M	H + H2O2 – k11 –OH + H2O

Table 1.

Chemical reactions (component M in reactions means some collision partner of reaction, k_i—speed of appropriate direct reaction $k_f$).

During the spreading of shock wave on the top bound of the region (conical surface) deflagration of gas mix appears. On Figure 3(a)–(c) level lines of temperature and meanings of water concentration in percentage (in number column) in consequent time moments are drawn. Change of this concentration inside the flow (Figure 3(c), from 5% to 15% of whole mass) demonstrates deflagration appearance. On Figure 3(d) graphics of water concentration along top boundary of the channel in consequent time moments are drawn (along y-axis concentration in percents of mass, along x-axes number of grid point in x-direction). Temperature of deflagration appearance approximately 700–800 K degree.

Figure 3.

(a)–(c) Vectors of velocity, level lines of temperature and in right number column corresponding to level lines meanings of concentration in percent of H₂O in consecutive time moments; (d) graphic of concentration of water along top boundary of the channel in consequent moments (along y-axis graphics of concentration in moll units, along x-axes is numbers of grid point in the x-direction).

5.2 Flow in cone-form channel

Axis symmetrical calculation region is a channel between cone surface and cylinder. The shock wave of intensity M_sh.w. = 3.1 spreading from the left boundary of the region. There is no one point of compression of flow in this problem. Compression takes place along all top boundaries of the region and strengthens with time. In Figure 4(a) and (b) level lines of temperature and temperature meanings (in right number column) in two consequent time moments are drawn. In Figure 4(c), graphics of concentration of water along top boundary of the channel in consequent time moments are drawn (along y-axis concentration in percents of whole mass, along x-axes number of grid point in x-direction).

Figure 4.

(a) and (b) Vectors of velocity, level lines of temperature and in right number column corresponding to level lines meanings of temperature in Kelvin in two consecutive time moments; (c) graphics of molar concentration of water along top boundary of the channel in consequent time moments (along y-axis graphics of concentration in moll units, along x-axes is numbers of grid point in x-direction).

5.3 Flow in complex form channel with obstacles

Numerical simulation of flows with deflagration appearance in axis-symmetric channels of complex form with obstacles was provided. Curvilinear structural grids for region of channel were constructed on the basis of algorithm [20] (Figure 2(c)). In Figure 5(a) and (b) level lines of temperature and meanings of molecular hydrogen in percents of mass concentration multiply to 10³ (in right number column) in two consequent moments are drawn. Areas of concentration of level lines correspond to fronts of deflagration waves.

Figure 5.

(a) and (b) Vectors of velocity, level lines of temperature and in write column corresponding to level lines meanings of molar concentration of H₂ in percentage multiplied by 10³ inside the axis-symmetric channel in consecutive time moments are drawn.

5.4 Flow in channel form which consist of combination of cylinders and cones

The axis symmetrical calculation region consists of combination of cylinders and cones. The shock wave of intensity M_sh.w. = 2.7 spreading from the left boundary of the region. Deflagration appears at the boundary point between second cone and second cylinder. In Figure 6(a)–(c) level lines of temperature and meanings of temperature in Kelvin degree (in right column) in three consequent moments are drawn. In Figure 6(d), graphics of concentration of water along top boundary of the channel in consequent moments are drawn (along y-axis concentration in molar units, along x-axes is numbers of grid point in x-direction).

Figure 6.

(a)–(c) Vectors of velocity, level lines of temperature and in right column corresponding to level lines meanings of temperature in Kelvin in consecutive three moments; (d) graphic of molar concentration of water along top boundary of the channel in consequent time moments (along y-axis concentration in moll units, along x-axes is numbers of grid point in x-direction).

5.5 Flow in cylinder with obstacles

The calculation region (calculation grid is drawn in Figure 1(a)) consists of cylinder with three obstacles. The shock wave of initial intensity M_sh.w. = 2.2 moves from the left boundary of the region. Deflagration appears at the obstacles and, later near axes of symmetry. In Figure 7(a)–(d) level lines of temperature and water concentration (at the same picture) in four consequent moments are drawn. In Figure 7(e), graphics of molar concentration of water along ax of symmetry of the channel in consequent time moments are drawn (along y-axis concentration of water in moll units, along x-axes number of grid point in x-direction).

Figure 7.

(a)–(d) Vectors of velocity, level lines of temperature and water concentration (on the same picture); on (e) graphics of concentration of water along the ax of symmetry of the channel in consequent time moments are drawn (along y-axis concentration in moll units, along x-axes is numbers of grid point in x-direction).

6. Conclusions

In the chapter algorithm of calculation flows of reactive gas mixes with arbitrary number of component was provided. This algorithm consists of modification of one, developed in Martyushov [17] and based on pure upwind variant of Chakravarthy-Osher difference scheme. This modification allows to avoid oscillations of values of radicals concentrations in regions of deflagration appearance during calculations of hydrogen-air mixes flows.

Calculations of flows with the initiation of deflagration in hydrogen-air mixes in some geometries of channels with obstacles by this algorithm modification were provided.

For getting solution of system of kinetic equations as part of gas dynamic calculation two algorithms were used: Algorithm of solution of fool system of kinetic equation for multistage hydrogen-air reaction, and algorithm of separate solution of two parts of system of kinetic equation with different size of time step based on Semenov’s theory of branching chain reaction.

Structured curvilinear calculation grids for calculation performed were constructed by author’s algorithm.