We present results of analytical and numerical investigation of explosive evaporation and burning scenarios of hydrogen droplets in hydrogen/oxygen aerosols. The following two scenarios have been elucidated. The first scenario, corresponding to sufficiently large droplets, is characterized by three stages: (i) an essentially homogeneous heating of a droplet to a near-critical temperature by IR radiation from the hot gas; (ii) explosive evaporation; and (iii) burning of hydrogen cloud formed by evaporation. The second scenario, corresponding to small droplets, differs in that a droplet is heated mainly by thermal conduction from the hot gas. The heating is accompanied by evaporation which can become explosive at the final stage of evaporation. The crossover droplet size separating the two scenarios is calculated. Conservative finite-difference numerical analysis is used to explore the predicted scenarios and verify analytical estimates.

Hydrogen-oxygen (*H*_{2}/*O*_{2}) mixtures are the most effective propellant for modern vehicles and rockets, such as the Space Launch System. These mixtures are explosive and the Challenger disaster of 1986^{1} has been associated with their self-ignition.^{2,3} The risk of explosion is controlled by complex processes of evaporation and burning of droplets in the *H*_{2}/*O*_{2} aerosol mixtures. These processes also take place in rocket combustion chambers and are central for understanding of the effectiveness and safety of *H*_{2}/*O*_{2} liquid engines. Peculiarities of these processes are very low critical temperature of hydrogen, $Tc=33.145\u2009K$, and very high temperature of the combustion products, $Tg0=3500\u2009K\u22123900\u2009K$.^{2,3} In this letter, we explore the scenarios of evaporation and burning of *H*_{2} droplets, induced by the IR radiation from the hot gas forming as a result of combustion of the *H*_{2}/*O*_{2} mixture. Such extreme conditions for evaporation and burning of *H*_{2}/*O*_{2} aerosol mixtures have not been explored in the literature to date.^{4,5}

Theoretical arguments allow to develop a simplified model which captures the most essential characteristics of the processes and distinguishes between different scenarios of burning. Still, the phenomena are so complex^{5} that even the simplified model cannot be solved analytically and one has to employ an involved numerical analysis. To model the *H*_{2} droplet evaporation and burning, we used transient equations of gas dynamics under the assumption of spherical symmetry. The equations follow from the conservation of gas species masses (with sources and sinks due to the burning), energy, and momenta.^{7,8} The boundary conditions at the liquid-gas interface are determined by the mass balance and heat release due to evaporation (Stefan condition) and Hertz-Knudsen equation (see below). We treated the liquid *H*_{2} (*LH*_{2}) as incompressible and inviscid fluid and the gaseous $H2,H2O$, and *O*_{2} ($GH2,GH2O$, and *GO*_{2}, respectively) as ideal gases. We took into account only the diffusion of the lightest gaseous *H*_{2} that can be described by the coefficient $DH2=(T/T0)32(p0/p)DH2(T0,p0)$, where *T* and *p* are temperature and pressure of the mixture.^{8} An important characteristic of *GH*_{2} is its negative thermal diffusion where the flow is directed toward the area of higher temperature.^{8} For the description of combustion of gaseous *H*_{2}/*O*_{2} mixtures, we assumed that the burning rate is limited by the initiation reactions $H2+O2\u2192OH+OH$ and $H2+O2\u2192HO2+H$, having the lowest rate:^{3,9} $Rcomb=cH2cO2[5.5\xd7107\u2009exp(\u221219680/T)+0.74\xd7T2.433\u2009exp(\u221226926/T)]\u2009m3/mol\xd7s$. This approximation is similar to the one-step mechanism of Mitani and Williams.^{5} Stefan condition is used for the boundary

Here *R* is radius of the droplet; $kL(g)$ and $\rho L(g)$ are thermal-diffusion coefficient and density of liquid (gas), respectively; *q _{L}* is the evaporation heat; $TL(g)$ is the liquid (gas) temperature;

*j*is the evaporation mass flow rate given by well-known Hertz-Knudsen equation: $jevap=\beta (pH2\u2212ps(Ts))(2\pi RH2Ts)\u22120.5$, where $RH2$ is the gas constant for

_{ev}*H*

_{2}, $\beta $ is the accommodation coefficient, $pH2$ is the pressure of

*H*

_{2}near the droplet surface and $ps(Ts)$ is the pressure of saturated vapor; $TL(R)\u2261Ts$. For the saturation pressure, we use the correlation

^{10,11}$ps(T)=pc(T/Tc)\lambda $, where $pc=1.3\u2009MPa,\u2009Tc=33.145\u2009K$, and $\lambda =5.3$.

The numerical solution is based on conservative implicit finite-difference schemes using balance method for transformed original dynamical equations. The transformation of the equation was attained as follows. At each time step, the mapping of liquid and gas regions into two unit length computational segments was performed by an appropriate coordinate transformation. A non-uniform mesh clustered near the liquid/vapor interface was used in the computational region. A conservative implicit finite-difference schemes were developed for approximating the equations. Coefficients of the equations were calculated using the predicted value of instantaneous velocity of the phase change interface. The resulting equations are further linearized. At each time step, the system of algebraic equations, including the interface velocity, was solved by Gauss method after the LU-decomposition of the linear system. We note that similar approach was used for a problem of crystal growth.^{12,13} Theoretical and numerical analyses reveal the following scenarios and stages of burning.

Scenario I describes burning of a large droplet. We consider an *LH*_{2} droplet immersed in a hot gaseous mixture of $H2O$, *O*_{2}, and *H*_{2} at $T=Thg$. We show that under the action of IR radiation from the hot gas, the droplet is heated almost homogeneously to the temperature close to *T _{c}* (stage (i) – the “bulk” effect). As

*T*approaches

_{L}*T*, the latent heat $qL(T)$ drops sharply to zero

_{c}^{6}(see the inset to Fig. 1(a)) and, according to Eq. (1), the evaporation rate $R\u0307\u2261dR/dt$ diverges. Therefore, evaporation occurs explosively with a time delay $\tau delayrad$ which can be estimated using the energy balance equation

where *C _{L}* is the specific heat of liquid hydrogen,

*σ*is the Stefan-Boltzmann constant,

*α*is the absorption coefficient, and

*ε*is the overlap coefficient of the radiation spectrum of the gas and the absorption spectrum of

*H*

_{2}droplets. From analysis of experimental data,

^{14,15}it follows that $\alpha \u2243(1.25\u22122.0)\u2009cm\u22121$ within the relevant IR range and $\epsilon \u2248(0.15\u22120.3)$.

^{2,3}For calculation, we used $\alpha \epsilon =0.4\u2009cm\u22121$. Taking

*C*and

_{L}*ρ*constant, and using $Thg\u226bTL$ and $\alpha R\u226a1$ ($R\u22641\u2009mm$) in Eq. (2), we obtain $Qrad\u2243Vdr\alpha \epsilon \sigma Thg4$ and

_{L}The value of $\tau delayrad$ is found to be independent of the droplet radius. Analytical estimates for $Thg=3800\u2009K$ and $p=1.0(1.6)\u2009atm$ give $\tau delayrad=21(18)\u2009ms$, which are confirmed by the numerical results, $\tau delayrad\u224822(18)\u2009ms$ (see Fig. 1(a)).

In stage (ii) the explosive evaporation is expected to occur at $t>\tau delayrad$ as $R\u0307$ drastically decreases due to the properties of the evaporation heat (see inset in Fig. 1(a)) and Eq. (1). Indeed, the evaporation rate is found to be negligible at $t<22\u2009ms$ and begins to increase sharply at $t>22\u2009ms$ with $R\u0307=\u22120.15\u2009m/s$ at $tst=22.2\u2009m$ and $R\u0307=\u22122.2\u2009m/s$ at $t=23.3\u2009ms$, when the droplet has evaporated almost completely (Fig. 1(a)). The duration of explosive evaporation is approximately $23\u2009ms\u221222\u2009ms=1\u2009ms\u226a\tau delayrad\u224820\u2009ms$, which justifies the term “explosive.” The results presented in Fig. 1 were obtained using the following initial conditions: $R(0)\u2261R0=1\u2009mm,$ $p(r)=1\u2009atm,$ $TL(r)=20.4\u2009K$, and $Tgas(r\u226bR)\u2261Thg=3800\u2009K$. The initial profiles of temperature and $GH2$ concentration are shown by curves 1 in Fig. 1(c). We assumed that the initial droplet temperature is equal to the saturation temperature at ambient pressure and used NIST^{6} interpolation formulas for $CL(T),$ $\kappa (T)$ and $qL(T)$, e.g., $qL(x)=(\u221289+478x\u2212945x2+825x3\u2212269x4)$ where $x=T/Tc$. We note that the liquid temperature near the interface at $t=23.3\u2009ms$ is very close to *T _{c}* and $qL(T)$ has decreased by several orders of magnitude (see the inset to Fig. 1(a)) from its value at

*T*. The delay time to evaporation strongly increases as

_{L}*T*decreases: $\tau delayrad\u2243110\u2009ms$ at $Thg=2500\u2009K$.

_{hg}The evaporation creates a cold *GH*_{2} cloud which expands with time until its radius *R _{cloud}* stabilizes at the final stage of the explosive evaporation. The simulation showed that $Rcloud\u22484\u2009mm$ at $22.67\u2009ms\u2264t\u226423.30\u2009ms$ for a droplet with initial radius $R0=1\u2009mm$ (Fig. 1(c)). Maximal radius of the cloud is limited by the value of $Rcloud=R0\rho L/\rho GH23\u22484R0$, where $\rho GH2$ corresponds to saturation temperature.

We note that the preceding discussion and the results presented in Figs. 1(a) and 1(c) are insensitive to the initial composition of the gas surrounding the *H*_{2} droplet. When the ambient gas is pure $GH2$, the local burning does not occur. The local burning occurs when the *GO*_{2} is present in the droplet environment. The burning is seen to take place in a thin layer (the combustion front, $\u223c1\u2009mm$) at the distance about 20 mm from the droplet. Remarkably, the combustion front is located far from the cloud surface, $Rcloud\u226a20\u2009mm$. This effect is due the fact that the burning is fueled by the negative thermal diffusion of the light $H2$ into the heavy $H2O$ vapor containing a small fraction of $GO2$ increasing away from the cloud (Fig. 1(b)).

Stage (iii) is the main phase of burning which occurs after the droplet's complete evaporation and formation of the *GH*_{2} cloud (Fig. 2). As mentioned above, burning is limited by the rate of *GH*_{2} diffusion into the heavy $H2O$ vapor containing a small fraction of $GO2$. As a consequence, the combustion front propagates away from the cloud surface as the burning continues. The cloud shrinks due to the *GH*_{2} diffusion, leading to the shift of $H2O$ boundary toward the cloud center (Fig. 2(b)).

Scenario II describes burning of a small droplet. We show that evaporation of small droplets occurs mainly due to heat conduction from the hot gas at the droplet surface (“surface” effect). For heat-conduction-driven evaporation, the dynamics of *R*(*t*) is described by the equation

where Eq. (5) defines the effective temperature difference between the liquid surface and the ambient gas. We note that Eq. (4) formally corresponds to earlier results,^{5} obtained for sufficiently small effective temperature difference $Thg\u2212Ts$. For arbitrary effective temperature difference $\Delta T$, we obtain the following expression:

where *C _{p}* is the isobaric specific heat of

*GH*

_{2}and the approximate equality holds for $Thg\u226bTs$. The analytical prediction of Eqs. (4)–(6) is supported well by the numerical simulations for sufficiently small droplets. The evaporation times $0.77(0.8)\u2009ms$ of $25\mu $ droplet for $p=1(1.6)\u2009atm$ found in our simulations are close to $\tau evapcond=0.78(0.81)\u2009ms$ calculated with Eqs. (5) and (6). For

*β*= 1, it was expected that $Ts=Ts(p)=const$ but for $\beta =0.01$ the surface temperature

*T*increases to

_{s}*T*, $qL(T)\u21920$ and the evaporation rate $R\u0307\u2192\u221e$. It means that the explosive evaporation for small droplets occurs only at the final stage of evaporation and only if

_{c}*β*is small enough. Similar to Scenario I, the burning occurs mainly at the vapor cloud surface after the droplet has evaporated completely. We note in passing that Eq. (5) can be used to describe evaporation due to the convective heat transfer if the factor of 2 is replaced by the appropriate Nusselt number.

^{16}The convection is relatively small when $(Nu\u22122)\u226a2$ which holds for the droplets of radius $R\u226450\mu ,\u2009Thg\u2265300\u2009K$, and the velocity of ambient gas $Vg\u22641\u2009m/s$.

Finally, we can establish the crossover droplet size delineating the boundary between two considered scenarios. The boundary is given by the condition $\tau evapcond(Rc)=\tau delayrad(Rc)$. Using Eqs. (3) and (5), we find that the radiation-stimulated explosive evaporation is dominant for the droplets of (initial) radius $R0>Rc\u22432\kappa gCL(Tc\u2212TL)\Delta T/(\alpha qL\epsilon \sigma Tgh4)$, e.g., $Rc\u22430.13(0.26)\u2009mm$ for $Tgh=3900(2500)\u2009K$. The values of *R _{c}* for $Thg=2500\u2009K$ and 3800 K are presented in Fig. 3, where the evaporation time-scales are shown as functions of the droplet initial radius, in accordance with Eqs. (3) and (5).

Numerical and analytical analyses of evaporation scenarios of cryogenic *H*_{2} droplets in burning *H*_{2}/*O*_{2} lead to the following results. Evaporation of sufficiently small droplets is driven by thermal conductance from the hot ambient gas. However, for sufficiently large droplets the main mechanism of heat transfer is radiation from the combustion products. It is found that radiation leads to explosive evaporation of cryogenic *H*_{2} droplets in burning *H*_{2}/*O*_{2} aerosol mixtures. Explosive evaporation of the *H*_{2} droplets supplies the mass flow of gaseous *H*_{2} to the burning mixtures. As a result, the intensity of explosion of *H*_{2}/*O*_{2} aerosol mixtures is dramatically increased. This process is crucially important for understanding the risks associated with the *H*_{2}/*O*_{2} aerosol management, with applications to both nominal and faulty performance of *H*_{2}/*O*_{2} liquid engines and fuel tanks.

## References

_{2}/O

_{2}: Rate constants for $H2+O2\u2192H+HO2$ at high temperature