Numerical validation of XiFoam-based simulations for hydrogen explosion in RUT Facility

1. Introduction

Hydrogen presents a significant safety challenge in various industries, including nuclear power plants [1] and hydrogen fuel cell technologies [2], due to its high reactivity and wide flammability limits. Hydrogen combustion prediction tools must reliably predict the intricate pressure loads arising from diverse combustion regimes—spanning slow deflagration, turbulent combustion, and detonation—while requiring rigorous validation against experimental data at realistic scales. To address the demand for comprehensive data, the European HYCOM joint research project was initiated, focusing on establishing an experimental database for turbulent flame regimes while accounting for scale effects, multi-compartment geometries, and mixture gradients [1]. Central to this effort were large-scale tests at the Kurchatov Institute's RUT facility, which replicates reactor-relevant length scales [3]. Data from specific tests (e.g., HYC01, HYC14) provided a robust basis for the integral validation of computational fluid dynamics (CFD) and lumped-parameter codes [1]. Complementing HYCOM, the SUSANA project developed the Hydrogen Model Evaluation Protocol (HYMEP) to bolster the confidence and precision of CFD simulations in hydrogen safety assessments [2]. HYMEP encompasses all phases of safety modeling, including detonation and deflagration-to-detonation transition (DDT), with its validation database incorporating benchmark cases from the RUT facility (e.g., KI_RUT_hyd05, KI_RUT_hyd09) for thorough model evaluation [2].

Large-scale detonation data from the RUT facility have proven essential for benchmarking CFD models' predictive capabilities and identifying limitations in industrial-scale safety analyses. Yáñez et al. benchmarked multiple CFD codes against large-scale RUT detonation experiments, demonstrating strong quantitative predictions of detonation velocity, overpressure, and related parameters even on coarse grids (6–10 cm) [4]. Machniewski and Molga validated CFD simulation models by comparison with large-scale RUT tunnel detonation data, confirming higher detonation wave overpressures in confined versus unconfined conditions [5]. Huang et al. introduced the 3D DEST tool for hydrogen detonation simulations in severe accidents, showing reasonable alignment between simulated RUT facility detonation tests and experimental results [6]. Kim and Kim integrated a seven-step chemical reaction model into an OpenFOAM solver, successfully simulating RUT detonation cases (KI-RUT-HYD09, KI-RUT-HYD05) with performance matching or exceeding traditional codes [7]. Zbikowski et al. applied a ZND theory-verified LES detonation model to large-scale RUT experiments, revealing grid dependency in complex 3D geometries where stability diminishes for control volume sizes over 0.1 m [8]. Heidari et al. proposed a large-scale detonation simulation method using a single-step chemical model and coarse grids, validating through RUT tests that detonation velocity and Chapman-Jouguet (CJ) pressure predictions reasonably match experimental data [9]. Rao et al. tuned a single-step Arrhenius chemical model for coarse grids under RUT tunnel detonation conditions, achieving reasonable agreement in detonation wave propagation velocity and pressure decay [10].

This study aims to enhance confidence in hydrogen explosion simulations by proposing a XiFoam-based numerical model [11], grounded in the context of prior large-scale experiments and CFD validation activities. Specifically, it seeks to reproduce the experimental pressure data from the RUT facility. This approach will expand the accuracy and applicability of CFD analyses in the field of hydrogen safety, ultimately contributing to the safety evaluation of hydrogen energy-related facilities.

2. Numerical Simulation Model

To validate the numerical simulation model, the experimental data were obtained from the KI-RUT-HYD09 test at measurement points P7, P8, P9, P10, and P11 in Figure 1. The measured pressure data resulting from an explosion induced by igniting a 25.5% hydrogen mixture with 200g of TNT in a test facility 35.5 meters long. The experimental conditions are listed in Table 1.

Figure 1:

Configuration of RUT experiment facility [3,4]

Table 1:

Experimental condition for KI-RUT-HYD09

The simulation of the explosion phenomenon was performed using the XiFoam solver for unsteady combustion analysis of the homogeneous mixture. The general governing equations of mass, momentum, and energy conservation are applied with an ideal gas equation. The species transport equation for the turbulent premixed combustion phenomenon is modeled by applying the reaction progress variable b [12].

where T_u denotes the unburned gas temperature and T_b represents the burned gas temperature. c is defined as $$ c=\frac{T-T_{\mathrm{u}}}{T_{\mathrm{b}}-T_{\mathrm{u}}}$$. Here, c = 0 (b = 1) represents burnt products and c = 1 (b = 0) represents unburned fuel-air mixture.

where ρ is the density, u is the flow velocity, $$ {Sc}_{\mathrm{t}}=\frac{{\mu }_{\mathrm{t}}}{\rho D}$$ is the turbulent Schmidt number, μ_t is the turbulent viscosity, D is the diffusion coefficient, ρ_u is the density of unburnt mixture, and S_u is the laminar flame speed. Ξ is the flame wrinkling variable, defined by the ratio of turbulent flame speed and laminar flame speed.

where S_t represents the turbulent flame speed. In this study, the flame wrinkling factor is solved using a transport equation [12], which means a balance between the creation of wrinkles by turbulence and their destruction by diffusion.

where u_s is the flame surface velocity, GΞ is the source term representing the production of flame surface area due to turbulent eddies, R( Ξ-1) is the sink term representing the destruction of flame surface area due to molecular diffusion and chemical reactions, and (σ_s - σ_t)Ξ is the term representing the effects of strain on the flame surface, which can either increase or decrease wrinkling depending on the flame dynamics.

In this study, the laminar flame speed is also solved using a transport equation [12].

where - σ_s S_u is the negative source term that models the reduction of the laminar flame speed due to flame stretch ( σ_s ), $$ {\sigma }_{\mathrm{s}}{S}_{\mathrm{u}}^{\mathrm{\infty }}\frac{\left(S_{\mathrm{u}}^0-S_{\mathrm{u}}\right)}{\left(S_{\mathrm{u}}^0-S_{\mathrm{u}}^{\infty }\right)}$$ is the positive source term that drives the flame speed toward equilibrium, $$ S_{\mathrm{u}}^{\mathrm{\infty }}$$ is the laminar flame speed at the extinction strain rate, and $$ S_{\mathrm{u}}^0$$ ​ is the unstrained laminar flame speed, which is evaluated using Gülder’s correlation [13]:

where W(2.094), η(1.068), ξ(0.424), α(2.9), and β(-0.04) are model parameters for the hydrogen-air mixture, and Φ denotes equivalence ratio. Considering the experimental environment, the simulation conditions are selected as listed in Table 2. Details of the case setup, including fvSchemes, fvSolution, and controlDict, is available in the XiFoam tutorial by Yasari [11].

Table 2:

Initial and boundary conditions for simulation

3. Results and Discussion

To simulate ignition from a 200 g TNT detonation, equivalent temperature and pressure conditions were applied at the ignition site depicted in Figure 1. The explosion's initial pressure was set to 30.96 bar calculated using the Sadovsky correlation in atm [13]:

where Δp is overpressure (atm), m is mass of explosive (kg) and r is distance from explosion (m). And the initial temperature was set to 3,150 K based on high explosives properties [14]. The ignition region was initialized with the progress variable b = 0 and the radius of 0.5 m, considering the RUT facility geometric constraints.

Figure 2 illustrates progressively refined hexahedral meshes used to evaluate grid dependence and the influence of early ignition behavior on the subsequent explosion dynamics. Grid density was assessed by the ratio of the ignition sphere volume to the volume of a single computational cell within the ignition region, based on an ignition condition of 3,047 K and 5.84 bar. As shown in Figure 3, the pressure-time history for probes P7–P11 converged at a grid resolution of 121,980 cells, corresponding to an ignition sphere volume to single-cell volume ratio of approximately 352.47. For subsequent simulations, a numerical model with 298,688 cells was employed, achieving an ignition sphere volume to single-cell volume ratio of 890, as depicted in Figure 2(c).

Figure 2:

Mesh variation based on grid precision: (a) Volume ratio 111.30; (b) Volume ratio 352.47; (c) Volume ratio 890.00.

Figure 3:

Effect of grid resolution on pressure-time histories at sensor locations P7–P11, shown in subfigures (a)–(e), respectively

To evaluate the influence of wall heat transfer on simulation outcomes, based on an ignition condition of 3150 K and 30.57 bar, simulations were conducted with two wall boundary conditions: a fixed temperature of 300 K and an adiabatic condition using a zero-gradient boundary. These simulations enabled a quantitative assessment of wall heat transfer effects on the pressure dynamics. As shown in Figure 4, pressure measurements at all sensor locations exhibited negligible differences between the two configurations. Consequently, wall heat transfer was determined to have minimal impact on the results, and the adiabatic condition, implemented with a zero-gradient boundary for wall temperature, was adopted for subsequent simulations.

Figure 4:

Comparison of wall boundary conditions: fixed wall temperature (300 K) and adiabatic wall condition (zero-gradient), with subfigures (a)-(e) showing pressure-time histories at sensor locations P7–P11, respectively

To identify the turbulence model most consistent with experimental overpressure data, two approaches were compared: a large-eddy simulation (LES) k-equation model and a shear-stress transport detached-eddy simulation (DES) k-ω SST-DES model. As a result, Figure 5 illustrates the temporal evolution of a representative flame surface. The k-ω SST-DES exhibits greater flame surface wrinkling compared to the k-equation LES as the flame propagates along a curved path, resulting in an expanded combustion region where the progress variable ranges from 0.1 < b < 0.9. However, Figure 6, the k-ω SST-DES overpredicts pressure-time histories at sensor locations P7 and P11 and underpredicts them at P8, P9, and P10, compared to the LES k-equation model, which demonstrates closer agreement with experimental data. This suggests that non-uniformity in turbulent flame propagation may influence wall pressure measurements.

Figure 5:

Comparison of turbulence models: (a) k-equation (LES) and (b) k-ω SST-DES (SSTDES). Flame surface contours on the x-y plane at z = 5 m are shown, with blue representing b = 0.9 and red representing b = 0.1, respectively

Figure 6:

Comparison of turbulence models: k-equation (LES) and k-ω SST-DES (SSTDES), with subfigures (a)-(e) showing pressure-time histories at sensor locations P7–P11, respectively

In addition to pressure-time histories, Figure 7 presents local pressure field distributions visualized using two iso-surface types: shock wave and flame fronts. The shock front was detected by evaluating the first and second derivatives of the density gradient projected onto the normalized flow direction, as described by:

Figure 7:

Shock wave and flame front iso-surfaces with normal vectors and pressure field distribution on the x-y plane (z = 5 m)

where ε is a filtering threshold [15]. The flame front was extracted using the contour of b = 0.1. Figure 7 illustrates the representative temporal evolution of these iso-surfaces along a curved path, with the shock wave propagating ahead of the flame front and elevated pressure distributed non-uniformly between them. This also suggests that the curvature and propagation behavior of the shock wave and flame front influence the pressure distribution.

Figure 8 presents the propagation speeds of the flame and shock fronts, calculated based on the time-dependent displacement of their mean surface coordinates. The mean position of each front was determined by averaging the x, y, and z coordinates of all surface points at each time step. The instantaneous propagation speed was computed by dividing the displacement of the mean coordinates by the time interval, as given by:

Figure 8:

Propagation speeds of the flame and shock fronts.

where X_mean​ is the mean coordinates of all points on the shock or flame front at each time step, and Δt is the sampling time step, set to 5×10^-6 s.

The shock wave initially propagates at approximately 3,000 m/s at ignition, but decelerates as it advances into the unburned hydrogen mixture, indicating not a self-sustained detonation.

This deceleration is attributed to the transition from a high-energy initiator, such as TNT with a detonation velocity ranging from 2,000-7,000 m/s depending on mixture properties [16,17], to hydrogen mixture, which exhibits a lower detonation velocity of 1,000-2,000 m/s [18,19]. In contrast, the flame front velocity is initially lower than that of the shock wave but gradually approaches it over time. This behavior is likely due to turbulence-enhanced combustion driven by interactions between the flame front and surrounding flow structures, particularly along curved paths and in canyon regions.

Figure 9 compares numerical simulation pressure data with experimental results, incorporating a +9 ms time offset to account for uncertainties in ignition delay between the experiment and simulation. The maximum pressure distribution from the cross-section shows closer agreement with experimental data than point-based measurements from discrete sensors, with the largest discrepancies observed at sensor P7. Additional simulations conducted at initial ignition temperatures of 3,000 K and 3,300 K (±150 K from baseline) revealed minimal deviations in peak overpressure, demonstrating XiFoam's capability to accurately reproduce the time-dependent overpressure evolution of hydrogen explosions. Combining point-based measurements and the maximum pressure distribution from the cross-section improves the evaluation of experimental conditions.

Figure 9:

Comparison of point-based and surface-based pressure measurements, with subfigures (a)-(e) showing pressure-time histories at sensor locations P7–P11, respectively. The first peak time includes interpolated maximum pressure distributions from the cross-section and wall-mounted sensor values

4. Conclusion

This study employed the XiFoam solver to simulate large-scale hydrogen explosions, validated against the KI-RUT-HYD09 hydrogen–TNT experiments at the RUT facility. By systematically optimizing simulation parameters, including mesh resolution, wall boundary, turbulence model, and pressure evaluation methods, the numerical model achieved robust agreement with experimental data. Key findings are: (1) maximum overpressure convergence at ignition sphere-to-cell volume ratios above 352.47; (2) negligible effects of wall temperature variations on pressure and flame propagation; (3) superior performance of the k-equation LES turbulence model over the k-ω SST-DES model in reproducing pressure-time histories; (4) closer alignment of maximum pressure distributions from cross-sections with experimental data compared to point-based sensor measurements, particularly at sensors P7–P11; and (5) shock front propagation at an initial 3,000 m/s, decelerating to a hydrogen detonation velocity of 1,400 m/s, with the flame front trailing due to turbulence-enhanced combustion along curved paths and canyon regions.

Applying a +9 ms time offset to account for ignition delay uncertainties between the experiment and simulation, XiFoam accurately captured the time-dependent overpressure evolution and shock–flame separation dynamics. The optimized parameter set provides a reliable framework for future CFD-based hydrogen explosion modeling, enhancing safety assessments for hydrogen energy facilities.

Acknowledgments

This research was supported by Korea Institute of Marine Science & Technology Promotion (KIMST) funded by the Ministry of Ocean and Fisheries (20220603), Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2022R1I1A3063360), and the University innovation support project Research Fund of Hoseo University in 2025 (2024-0282-01).

References

• W. Breitung, S. Dorofeev, A. Kotchourko, R. Redlinger, W. Scholtyssek, A. Bentaib, J.-P. L’Heriteau, P. Pailhories, J. Eyink, M. Movahed, K.-G. Petzold, M. Heitsch, V. Alekseev, A. Denkevitse, M. Kuznetsov, A. Efimenko, M.V. Okune, T. Huld, and D. Baraldi, “Integral large scale experiments on hydrogen combustion for severe accident code validation-HYCOM,” Nuclear Engineering and Design, vol. 235, pp. 253–270, 2005. [https://doi.org/10.1016/j.nucengdes.2004.08.063]
• D. Baraldi, D. Melideo, A. Kotchourko, K. Ren, J. Yanez, O. Jedicke, S.G. Giannissi, I.C. Tolias, A.G. Venetsanos, J. Keenan, D. Makarov, V. Molkov, S. Slater, F. Verbecke, and A. Duclos, “Development of a model evaluation protocol for CFD analysis of hydrogen safety issues the SUSANA project,” International Journal of Hydrogen Energy, vol. 42, pp. 7633–7643, 2017. [https://doi.org/10.1016/j.ijhydene.2016.05.212]
• S. B. Dorofeev, V. P. Sidorov, W. Breitung, A. Kotchourko, “Large scale combustion test in the RUT facility: Experimental study, numerical simulations and analysis on turbulent deflagration and DDT,” Proceedings of SMiRT-14, Lyon, France, pp. 2463–2471, 1997.
• J. Yáñez, A. Kotchourko, A. Lelyakin, A. Gavrikov, A. Efimenko, M. Zbikowski, D. Makarov, and V. Molkov, “A comparison exercise on the CFD detonation simulation in large-scale confined volumes,” International Journal of Hydrogen Energy, vol. 36, pp. 2613–2619, 2011. [https://doi.org/10.1016/j.ijhydene.2010.04.133]
• P. Machniewski and E. Molga, “CFD analysis of large-scale hydrogen detonation and blast wave overpressure in partially confined spaces,” Process Safety and Environmental Protection, vol. 158, pp. 537–546, 2022. [https://doi.org/10.1016/j.psep.2021.12.032]
• T. Huang, Y.P. Zhang, W.X. Tian, G.H. Su, and S.Z. Qiu, “A 3-D simulation tool for hydrogen detonation during severe accident and its application,” Annals of Nuclear Energy, vol. 104, pp. 113–123, 2017. [https://doi.org/10.1016/j.anucene.2017.02.011]
• D. Kim and J. Kim, “Numerical method to simulate detonative combustion of hydrogen-air mixture in a containment,” Engineering Applications of Computational Fluid Mechanics, vol. 13, no. 1, pp. 938–953, 2019. [https://doi.org/10.1080/19942060.2019.1660219]
• M. Zbikowski, D. Makarov, and V. Molkov, “Numerical simulations of large-scale detonation tests in the RUT facility by the LES model,” Journal of Hazardous Materials, Vol. 181, pp. 949–956, 2010. [https://doi.org/10.1016/j.jhazmat.2010.05.105]
• A. Heidari, S. Ferraris, J.X. Wen, and V.H.Y. Tam, “Numerical simulation of large scale hydrogen detonation,” International Journal of Hydrogen Energy, vol. 36, pp. 2538–2544, 2011. [https://doi.org/10.1016/j.ijhydene.2010.05.093]
• V.C. Madhav Rao, A. Heidari, J.X. Wen, and V.H.Y. Tam, “Numerical study of large scale hydrogen explosions and detonation,” SYMPOSIUM SERIES NO. 155 Hazards XXI (IChemE), pp. 630–637, 2009.
• E. Yasari, Tutorial XiFoam, 2010, Available from http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2010/ehsanYasari/ehsanYasariReport.pdf, (Accessed 2025-10-07).
• H. G. Weller, G. Tabor, A. D. Gosmanand C. Fureby, “Application of a Flame-Wrinkling LES Combustion Model to a Turbulent Mixing Layer,” Twenty-Seventh Symposium (International) on Combustion/The Combustion Institute, pp. 899–907, 1998. [https://doi.org/10.1016/S0082-0784(98)80487-6]
• Z. Bajic, J. Bogdanov, R. Jeremic, “Blast effects evaluation using TNT equivalent, Scientific Technical Review,” Vol. 59, No. 3–4, pp. 50–52, 2009.
• B. Gelfand, M. Silnikov, “Blast Effects Caused by Explosions,” European Research Office of the U.S. Army, London, England, Final Technical Report, N62558-04-M-0004, 2004.
• Z. Wu, Y. Xu, W. Wang, R. Hu, “Review of shock wave detection method in CFD post-processing,” Chinese Journal of Aeronautics, vol. 26, no. 3, pp. 501–513, 2013. [https://doi.org/10.1016/j.cja.2013.05.001]
• W. B. Cybulski, W. Payman, D. W Woodhead, “Explosion waves and shock waves. VII.The velocity of detonation in cast T.N.T,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 197, no. 1048, 1949. [https://doi.org/10.1098/rspa.1949.0050]
• G. D. Kozak, B. N. Kondrikov, A. I. Sumin, “Dependence of detonation velocity on charge density for foamed alumotol (Al/TNT) and TNT mixtures,” Combustion, Explosion and Shock Waves, vol. 34, pp. 448–452, 1998. [https://doi.org/10.1007/BF02675614]
• Y. Gao, B. Zhang, H. D. Ng, J. H. S. Lee, “An experimental investigation of detonation limits in hydrogen-oxygen-argon mixtures,” International Journal of Hydrogen Energy, vol. 41, no. 14, pp. 6076–6083, 2016. [https://doi.org/10.1016/j.ijhydene.2016.02.130]
• W. Rudy, R. Porowski, A. Teodorczyk, “Propagation of hydrogen-air detonation in tube with obstacles,” Journal of Power Technologies, vol. 91, no. 3, pp. 122–129, 2011.