Radiation-stimulated explosive evaporation and burning of hydrogen droplets in hot aerosol mixtures

Hydrogen-oxygen (H₂/O₂) 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¹ 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₂/O₂ aerosol mixtures. These processes also take place in rocket combustion chambers and are central for understanding of the effectiveness and safety of H₂/O₂ liquid engines. Peculiarities of these processes are very low critical temperature of hydrogen, $Tc=33.145 K$⁠, and very high temperature of the combustion products, $Tg0=3500 K−3900 K$⁠.^2,3 In this letter, we explore the scenarios of evaporation and burning of H₂ droplets, induced by the IR radiation from the hot gas forming as a result of combustion of the H₂/O₂ mixture. Such extreme conditions for evaporation and burning of H₂/O₂ 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⁵ that even the simplified model cannot be solved analytically and one has to employ an involved numerical analysis. To model the H₂ 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₂ (LH₂) as incompressible and inviscid fluid and the gaseous $H2,H2O$⁠, and O₂ (⁠$GH2,GH2O$⁠, and GO₂, respectively) as ideal gases. We took into account only the diffusion of the lightest gaseous H₂ 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.⁸ An important characteristic of GH₂ is its negative thermal diffusion where the flow is directed toward the area of higher temperature.⁸ For the description of combustion of gaseous H₂/O₂ mixtures, we assumed that the burning rate is limited by the initiation reactions $H2+O2→OH+OH$ and $H2+O2→HO2+H$⁠, having the lowest rate:^3,9 $Rcomb=cH2cO2[5.5×107 exp(−19680/T)+0.74×T2.433 exp(−26926/T)] m3/mol×s$⁠. This approximation is similar to the one-step mechanism of Mitani and Williams.⁵ Stefan condition is used for the boundary

Here R is radius of the droplet; $kL(g)$ and $ρ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_ev is the evaporation mass flow rate given by well-known Hertz-Knudsen equation: $jevap=β(pH2−ps(Ts))(2πRH2Ts)−0.5$⁠, where $RH2$ is the gas constant for H₂, $β$ is the accommodation coefficient, $pH2$ is the pressure of H₂ near the droplet surface and $ps(Ts)$ is the pressure of saturated vapor; $TL(R)≡Ts$⁠. For the saturation pressure, we use the correlation^10,11$ps(T)=pc(T/Tc)λ$⁠, where $pc=1.3 MPa, Tc=33.145 K$⁠, and $λ=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₂ droplet immersed in a hot gaseous mixture of $H2O$⁠, O₂, and H₂ 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_L approaches T_c, the latent heat $qL(T)$ drops sharply to zero⁶ (see the inset to Fig. 1(a)) and, according to Eq. (1), the evaporation rate $Ṙ≡dR/dt$ diverges. Therefore, evaporation occurs explosively with a time delay $τ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₂ droplets. From analysis of experimental data,^14,15 it follows that $α≃(1.25−2.0) cm−1$ within the relevant IR range and $ε≈(0.15−0.3)$⁠.^2,3 For calculation, we used $αε=0.4 cm−1$⁠. Taking C_L and ρ_L constant, and using $Thg≫TL$ and $αR≪1$ (⁠$R≤1 mm$⁠) in Eq. (2), we obtain $Qrad≃VdrαεσThg4$ and

The value of $τdelayrad$ is found to be independent of the droplet radius. Analytical estimates for $Thg=3800 K$ and $p=1.0(1.6) atm$ give $τdelayrad=21(18) ms$⁠, which are confirmed by the numerical results, $τdelayrad≈22(18) ms$ (see Fig. 1(a)).

In stage (ii) the explosive evaporation is expected to occur at $t>τdelayrad$ as $Ṙ$ 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 ms$ and begins to increase sharply at $t>22 ms$ with $Ṙ=−0.15 m/s$ at $tst=22.2 m$ and $Ṙ=−2.2 m/s$ at $t=23.3 ms$⁠, when the droplet has evaporated almost completely (Fig. 1(a)). The duration of explosive evaporation is approximately $23 ms−22 ms=1 ms≪τdelayrad≈20 ms$⁠, which justifies the term “explosive.” The results presented in Fig. 1 were obtained using the following initial conditions: $R(0)≡R0=1 mm,$ $p(r)=1 atm,$ $TL(r)=20.4 K$⁠, and $Tgas(r≫R)≡Thg=3800 K$⁠. 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⁶ interpolation formulas for $CL(T),$ $κ(T)$ and $qL(T)$⁠, e.g., $qL(x)=(−89+478x−945x2+825x3−269x4)$ where $x=T/Tc$⁠. We note that the liquid temperature near the interface at $t=23.3 ms$ 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_L. The delay time to evaporation strongly increases as T_hg decreases: $τdelayrad≃110 ms$ at $Thg=2500 K$⁠.

The evaporation creates a cold GH₂ cloud which expands with time until its radius R_cloud stabilizes at the final stage of the explosive evaporation. The simulation showed that $Rcloud≈4 mm$ at $22.67 ms≤t≤23.30 ms$ for a droplet with initial radius $R0=1 mm$ (Fig. 1(c)). Maximal radius of the cloud is limited by the value of $Rcloud=R0ρL/ρGH23≈4R0$⁠, where $ρ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₂ droplet. When the ambient gas is pure $GH2$⁠, the local burning does not occur. The local burning occurs when the GO₂ is present in the droplet environment. The burning is seen to take place in a thin layer (the combustion front, $∼1 mm$⁠) at the distance about 20 mm from the droplet. Remarkably, the combustion front is located far from the cloud surface, $Rcloud≪20 mm$⁠. 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₂ cloud (Fig. 2). As mentioned above, burning is limited by the rate of GH₂ 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₂ 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,⁵ obtained for sufficiently small effective temperature difference $Thg−Ts$⁠. For arbitrary effective temperature difference $ΔT$⁠, we obtain the following expression:

where C_p is the isobaric specific heat of GH₂ and the approximate equality holds for $Thg≫Ts$⁠. 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) ms$ of $25μ$ droplet for $p=1(1.6) atm$ found in our simulations are close to $τevapcond=0.78(0.81) ms$ calculated with Eqs. (5) and (6). For β = 1, it was expected that $Ts=Ts(p)=const$ but for $β=0.01$ the surface temperature T_s increases to T_c, $qL(T)→0$ and the evaporation rate $Ṙ→∞$⁠. It means that the explosive evaporation for small droplets occurs only at the final stage of evaporation and only if β 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.¹⁶ The convection is relatively small when $(Nu−2)≪2$ which holds for the droplets of radius $R≤50μ, Thg≥300 K$⁠, and the velocity of ambient gas $Vg≤1 m/s$⁠.

Finally, we can establish the crossover droplet size delineating the boundary between two considered scenarios. The boundary is given by the condition $τevapcond(Rc)=τ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≃2κgCL(Tc−TL)ΔT/(αqLεσTgh4)$⁠, e.g., $Rc≃0.13(0.26) mm$ for $Tgh=3900(2500) K$⁠. The values of R_c for $Thg=2500 K$ 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₂ droplets in burning H₂/O₂ 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₂ droplets in burning H₂/O₂ aerosol mixtures. Explosive evaporation of the H₂ droplets supplies the mass flow of gaseous H₂ to the burning mixtures. As a result, the intensity of explosion of H₂/O₂ aerosol mixtures is dramatically increased. This process is crucially important for understanding the risks associated with the H₂/O₂ aerosol management, with applications to both nominal and faulty performance of H₂/O₂ liquid engines and fuel tanks.