Numerical method to simulate detonative combustion of hydrogen-air mixture in a containment

During an accident accompanied by core damage in a water-cooled nuclear reactor, a large amount of hydrogen is generated by fuel clad oxidation and released into the reactor containment. Hydrogen mitigation in the containment during a severe accident is important because the integrity of the structures important to safety could be threatened by an explosive combustion of the hydrogen. Even when equipped with a hydrogen mitigation system, if a highly concentrated hydrogen mixture cloud is developed in a compartment of the containment, the compartment integrity with a pressure load from a local detonation must be evaluated. In this study, a reliable numerical method to simulate the detonation of a hydrogen-air mixture was developed for the evaluation of a pressure load from a detonation. The numerical method has the capability to capture shocks robustly as well as resolve combustion phenomena stably. The Euler equations of mass, momentum, energy and gas species are employed as the governing equations of detonation. A central-upwind scheme was adopted as a shock capturing scheme and a reduced chemical kinetic model was applied for combustion. The numerical method was validated by solving experimental cases and applied for understanding detonation shock behaviour in a generic containment.

Containment building of an NPP is the last barrier retaining the radioactive materials which could be released during reactor accidents and also protects the CONTACT D. Kim dehee@kaeri.re.kr reactor and inner structures from any external impact. Events such as steam explosion, hydrogen detonation and airplane crash could bring severe impacts on the integrity of the containment building. The status report from NEA (NEA, 2018) summarized the significant progress made in the fuel coolant interaction (FCI) and the steam explosion (SE) research field. Study on the airplane crashes has been carried out to calculate the dynamic loads and its impacts by utilizing analytical tools (Saberi, Alinejad, Mahdavi, & Sepanloo, 2017;Siefert & Henkel, 2014) and experiments (Sugano et al., 1993). It is generally admitted that the containment building can withstand the impacts from the steam explosion and the airplane crashes. By the way, hydrogen release by a severe accident can develop into a detonation via hydrogen combustion. By regulation global hydrogen detonation inside the containment building should be practically eliminated on the design stage by reflecting on the containment volume and the hydrogen mass released. Nonetheless, there is still a possibility of local hydrogen combustion and it can bring impacts to compartment structures or components important to safety. Therefore, detonation wave propagation characteristics induced by the hydrogen combustion need to be investigated thoroughly for safety.
It is currently impractical to directly simulate deflagration to detonation transition (DDT) in a scale of a nuclear reactor containment because difficulties such as accurate turbulence modelling, required tremendous computational mesh to capture the small-scale physics during DDT, wide spectrum of severe accident progressions, and uncertainties of hydrogen generation in a nuclear reactor have to be resolved.
In this paper, we aim to develop and implement a numerical method to understand propagation characteristics of detonation wave such as attenuation or augmentation of strength of detonation shocks during propagating in a compartment of an NPP containment. This numerical method would be helpful for evaluation of a containment integrity from hydrogen hazard during a severe accident in an NPP containment.
During a severe accident in a water-cooled nuclear reactor, exposed fuel clads are oxidized as they meet high-temperature steam at a lowered water level inside the reactor vessel. Hydrogen gas is generated due to the oxidation of the fuel clads. The hydrogen gas released from the reactor vessel through a break or a pressurerelief valve is mixed with the air inside the containment building. The hydrogen-air mixture may lead to a supersonic explosion if the conditions are met. The explosion is called detonation and the effect of the impulsive mechanical load by the local detonation might be serious to the structures important to safety. To calculate the dynamic load induced by detonation, the peak pressure during detonation wave propagation needs to be accurately assessed.
To develop accurate and robust numerical method for simulation of detonation and understand the characteristics of detonation wave propagation, numerous research works have been conducted. Research to understand the behavior of detonation waves includes experimental works and numerical modelling (Dorofeev, Sidorov, Velmakin, & Zhernov, 1993;Pfoertner, 1991;Rao, Heidari, Wen, & Tam, 2009;SUSANA, 2016;Xiao, Breitung, Kuznetsov, Travis, & Redlinger, 2017;Yanez et al., 2011). Dorofeev et al. (1993) studied the effect of ignition location and hydrogen concentration on load for which 60 m long RUT facility was utilized and high explosive charges were given to generate detonation waves. Detonation experiments in hemispherical balloons filled with hydrogen-air mixture were conducted at Fraunhofer Institute to provide data for the CFD code verification (Pfoertner, 1991). Capabilities and limitations of CFD simulations were evaluated for the uniform hydrogenair mixture detonations in the RUT facility (Yanez et al., 2011). A numerical approach for the large scale detonation were developed and applied for the detonation tests conducted at the RUT facility (Rao et al., 2009). The COM3D code developed by KIT was utilized to analyze the RUT cases (SUSANA, 2016) and the GASFLOW-MPI code was applied to the detonation simulations in hemispherical balloons (Xiao et al., 2017).
Depending on the geometry of the experiment, wave propagation characteristics become different. An explosion in open space produces simple symmetric wave contours but an explosion inside a confined geometry makes complicated flow phenomena. Shock waves generated by the detonation can be amplified or expanded when it propagates through internal structures inside the containment building. Numerical models to resolve the detonation induced by combustion of a premixed hydrogen-air mixture should have capabilities of robust shock capturing as well as stable combustion modelling.
In this study, we utilized an OpenFOAM solver (Kroposhin, Bovtrikova, & Strijhak, 2015) to simulate the detonation phenomena. The numerical simulation model was validated through a shock tube test (Schmidt & Duffy, 1985;Toro, 1997), detonation in a hemispherical balloon (Pfoertner, 1991), and detonation in the RUT facility (Dorofeev et al., 1993). To simulate detonation accompanying chemical reactions, chemical reaction mechanism is crucial and thus several chemical reaction mechanisms were implemented and studied on the respects of numerical accuracy and stability. Through numerical validations and investigations, an accurate and robust numerical approach for combustive detonation simulation was built up. The validated numerical method was applied to a generic containment geometry so that a possible scenario of hydrogen-air detonation can be deduced from the application.
In Section 2, numerical methods including governing equations, schemes and reaction models are presented. Validation cases using the numerical method are provided in Section 3 and application to a generic containment building is given in Section 4. Concluding remarks are provided in Section 5.

Governing equations
Since propagation speeds of detonation waves are very high, flow physics such as molecular diffusion, thermal conduction and viscous effects can be ignored for detonation simulation (Berkenbosch, 1995;Fickett & Davis, 1979). Reactive Euler equations can be written as follows: where ρ, U, p, S, Y i ,ω i , E,Q are the mass density of the mixture, velocity vector, pressure, momentum sources, mass fraction of species i, reaction rate species i, specific total energy, heat sources, respectively.

Numerical schemes for convection terms and reaction models
Numerical schemes for the convection terms in the governing equations need to have shock capturing capability to resolve detonation waves. Historically, numerical schemes for solving Riemann problems have been mainly developed in two categories: upwind schemes and central difference schemes (Hirsch, 1988;Toro, 1997). Upwind schemes utilize characteristic information for robust and accurate calculation, and among them, Roe's flux difference splitting family and van Leer's flux vector splitting family are well known. Upwind schemes originated from Godunov's scheme (Toro, 1997) in which variables are approximated as constant values in a computational cell and accurate calculation is carried out at the cell interface. Godunov's scheme is accurate but time-consuming due to calculation at the cell interface. Therefore, approximate numerical schemes such as Roe's and van Leer's schemes were derived.
On the other hand, the Lax-Friedrichs scheme is a representative central difference scheme where characteristic information is not required. Jameson's scheme with an artificial dissipation (Hirsch, 1988) is a well-known central difference scheme. Recently, Kurganov-Tadmor developed a second-order central-upwind semi-discrete KT scheme (Kurganov & Tadmor, 2000) which can be considered a descendant of the Lax-Friedrichs scheme.
Generally speaking, upwind schemes are more robust and accurate but complicated to implement because the Riemann solver and characteristic decomposition are required. On the other hand, central difference schemes are simple but more diffusive than upwind schemes, especially in the discontinuity region.
The KT scheme is a central difference scheme but it can give comparable results to upwind schemes. Its simple structure makes it a good candidate for multi-fluid, multi-phase simulations where characteristic decomposition of the system becomes complicated or unavailable. We want to extend this research to multi-fluid, multiphase flow simulations. To this end, the KT scheme is employed and validated through benchmark problems and a generic NPP case.
For general hyperbolic conservation laws which leads, where x is spatial variable, w is conserved quantity and F is convection flux, discretization by the semi-discrete KT scheme is expressed as, where numerical flux H is written as The speed of wave propagation c i+1/2 and reconstructed values w + i+1/2 , w − i+1/2 of left and right faces at cell boundary i + 1/2 are defined as wherew n i , (w x ) n i are cell averaged values and spatial derivatives at time t n and space x i , respectively. Spatial derivatives are obtained by using slope limiters such as 'minmod.' A chemical reaction is solved by the finite-rate chemical kinetics models for which simplified kinetic models or elementary kinetic models can be employed. Results of the reacting flows are sensitive to the finite-rate mechanism. Therefore, 1-step, 2-step, and 7-step mechanisms were implemented and studied. While finite-rate mechanisms using 1-step and 2-step gave numerically unstable or diffusive results, it was found out that the Shang's 7-step chemical reaction mechanism (Shang, Chen, Liaw, Chen, & Wang, 1995) can produce accurate and stable physics.
An OpenFOAM solver which implemented the KT scheme and is capable of solving chemical reactions was developed by Kroposhin et al. (Kroposhin et al., 2015). For the OpenFOAM solver, we applied a finite rate Arrhenius reaction model with a reduced hydrogen-air 7 step chemistry (Shang et al., 1995) for detonative combustion simulations. Model constants are given in Table 1 and were used for the simulations in Sections 3 and 4.  (Shang et al., 1995).

Spatial operators and linear algebra solver
Governing equations are discretized by the finite volume method in which computational domain is divided by finite volumes. Over each finite volume, the governing equations are integrated and volume integrals are converted into surface integrals by Gauss theorem. Discretized forms of convection, diffusion terms of the governing equations require flux terms at faces between two neighboring cells. Variables, its gradients and divergence terms at the faces have to be interpolated from the information of neighboring cellcentered values. Discretized equations are solved iteratively by a segregated manner for which PIMPLE method was employed (Moukalled, Mangani, & Darwish, 2016). Matrices resulted from the discretized equations over the whole domain were solved by a conjugate gradient-based method with a matrix preconditioner. For time advancement, the Euler implicit method was applied. Information about numerical schemes and methods for the OpenFOAM solver are summarized in Table 2.

1-D shock tube problem
The shock tube problem is a fundamental test case for which the exact solution is known (Toro, 1997). Initially, a 1D long tube is divided into two regions by a thin diaphragm. Each region is filled with the same gas, but the thermodynamic parameters such as pressure, density and temperature are different from each other. As the diaphragm breaks abruptly, the gas in the higher pressure region moves fast into the low-pressure gas region. An expansion wave is generated when the high-pressure gas expands and a shock wave propagates as the low-pressure gas is compressed. A contact discontinuity moves behind the shock wave between the expansion wave and the shock wave. The 1D shock tube problem was chosen to see if the OpenFOAM solver can resolve key characteristics of a shock wave such as expansion fan, contact discontinuity and shock. Initial conditions were set to the left side and the right side separated by a center position inside the shock tube. Pressure, density, velocity, and temperature conditions were given as follows: (10000, 0.1, 0, 298), where the units of the properties are Pa, kg/m 3 , m/s, K, respectively. Computational domain ranges over −5 m ≤ x ≤ 5 m which is discretized with regular-sized cells.    Grid convergence test was carried out with different grid systems having 100, 200, 400, 800 and 1600 cells. As shown in Figure 1, case with more number of grids produces closer result to the exact solution. It is found that the resolution of the contact discontinuity is more sensitive to the number of grids than that of the shock. Figure 1 also presents time step sensitivity of the solver with 400 cells in which the shock location is more sensitive to the CFL number than the contact discontinuity location. CFL number lower than 0.2 gives good resolution but results with CFL number greater than 0.4 shows dispersive solutions.
Based on the time step sensitivity result, CFL number was set to 0.02 and zero gradient boundary conditions for flow variables were applied to both tube ends.
As the diaphragm is broken at time t = 0 ms, shock, contact discontinuity and expansion fan are formed. The shock and the contact discontinuity move right and the expansion fan moves left. At t = 7 ms, the pressure, density, and velocity profiles are displayed on the top of Figure 2 for which 100 grid points were used. For comparison with the OpenFOAM solver, calculations utilizing the 1st order accurate Roe scheme (Hirsch, 1988) and the second order accurate upwind TVD scheme (Hirsch, 1988) were also carried out. Three solvers can resolve the propagation of shock, contact discontinuity, and expansion fan, but Roe scheme's numerical dissipation smeared the shock more severely than the other two solvers. A grid refinement test with 400 grid points was performed and its result is shown on the bottom of Figure 2. The Roe scheme is diffusive, especially at the contact discontinuity and expansion fan regions. It was found that the OpenFOAM solver can give comparable solution to the second order upwind TVD scheme.
For the shock tube problem, the OpenFOAM solver resolved the key flow patterns and its accuracy corresponded with a second order scheme.

Open-ended shock tube
Some experiments and simulations of blast flow fields aim to investigate flow physics and predict the dynamic pressure load from a blast wave which are normally developed from firing military weapons (Bai et al., 2009) or rocket engines. Since the blast loading can have serious effects on neighboring structures or personnel nearby, refined experiments or simulations are required.
A cylindrical shock tube with one open end was tested to investigate the blast flow fields (Schmidt & Duffy, 1985) and a numerical study was carried out to reproduce the detailed experiments (Wang & Widhopf, 1990). Figure 3 illustrates the cylindrical shock tube configuration and pressure measurement locations. When planar shock waves propagate from inside the shock tube and reach its open end, oval-shaped shock waves are generated and complicated flow patterns such as the following are formed: vortex, blast wave, reflected wave, jet shear layer, jet shock, Mach disc, recompression shocks, triple point, contact surface, and slipstream (Schmidt & Duffy, 1985;Wang & Widhopf, 1990).  Since flow features are axisymmetric, a 2D axisymmetric domain was created for simulation. Tube inner diameter, D, is 0.152 m and streamwise and radial direction length scales of the computational domain are 6 × D ( = 0.912 m), 3 × D ( = 0.456 m), respectively. The computational domain was discretized with a Cartesian grid system. The centerline of the tube to the tube inner surface was discretized with 20 uniform cells, and the same cell size was applied in the entire flow domain. The initial pressure ratio of the pressure inside the tube and the ambient air is 3.42. Initial conditions were set so that the high-pressure field inside the tube produced a planar shock of Mach strength 1.76.
Grid convergence test was performed with mesh size, dx = 0.05 × D, 0.025 × D, 0.0125 × D and 0.00625 × D. Figure 4 shows the result at the measuring point 9 in which finer meshes give more accurate solutions. For   further detailed analysis, grid system with mesh size dx = 0.025 × D was applied to compare the solution with a reference CFD solution (Wang & Widhopf, 1990) with the same mesh size.
CFL number was set to 0.2 and total pressure, total temperature, zero gradient were given for the pressure, temperature and velocity inlet condition, respectively. Zero gradient conditions of pressure and temperature were applied for wall and far-field boundaries. Slip and pressure inlet/outlet velocity condition were applied for wall and far-field boundaries, respectively.
The developed blast flow field is displayed in Figure  5 which shows pressure contours at t = 0.2, 0.5, 1.0 and 1.5 ms following the time stamps recorded in Refs (Schmidt & Duffy, 1985;Wang & Widhopf, 1990). At t = 0.2 ms, a blast wave was formed and a vortex was generated at the tube exit due to a pressure difference between the jet stream and the ambient air. It was observed that the Mach disc is apparently generated at t = 1.5 ms and shocks, jet shock and vortex are developed as a consequence of complicated flow interactions.
To measure the pressure fluctuations, pressure transducers were mounted on the experiment. At the same locations, numerical results using the OpenFOAM solver were obtained and compared with the experimental results (Schmidt & Duffy, 1985) and the numerical simulation by Wang et al. (Wang & Widhopf, 1990). Figure  6 shows overpressure transient histories at several measurement locations. The results using the OpenFOAM solver are in good agreement with the experimental results (Schmidt & Duffy, 1985) and numerical results (Wang & Widhopf, 1990). At point 4, the numerical results deviate from the experimental results which might be influenced by the measuring device as mentioned by Wang et al. (Wang & Widhopf, 1990). In some results, time delay is shown but the pressure peak and wave shape is similar to the experimental data and compares well with the numerical simulation (Harten, 1983;Wang & Widhopf, 1990).

Detonation in a hemispherical balloon
In the preceding two test cases, fluid flows without chemical reactions were solved to validate the shock capturing capability of the solver in one and multi-dimensional domains. Flows coupled with chemical reactions give rise to a new stiffness problem in the equation solving procedure. Therefore, when solving the chemical reaction mechanism, numerical stability must be guaranteed. Figure 11. Simulation of the spherical balloon test.
The first case accompanying chemical reactions is a simple test without complex geometrical effects. Experiments on detonation by a premixed hydrogen-air mixture inside a hemispherical balloon was performed at Fraunhofer Institute to provide experimental data for verification of the CFD models developed at KIT (Xiao et al., 2017). Initially, the hemispherical balloon was filled with 29.05% hydrogen and 70.95% air at 1 bar and 304 K. The atmosphere pressure and temperature were the same for both inside and outside the balloon. The configuration is shown in Figure 7. Along the path of the moving shock, pressure sensors were mounted on the ground to measure the pressure transients in the experiment (Xiao et al., 2017).
The 3D computational domain was set to be a quarter of a hemisphere with 8 m radius due to symmetrical conditions. Sensitivity of number of grids was studied using three sets of grid systems, i.e. a spherical coarse mesh with 731,190 cells, a spherical fine mesh with 5,771,190 cells, and a Cartesian fine mesh with 8,000,000 cells. Figure 8 displays that the Cartesian fine mesh gives more diffusive result than the spherical coarse mesh even though about 10 times more cells were used and the spherical coarse mesh produces comparable accuracy to the spherical fine mesh. Therefore, reasonable grid system with 731,190 cells were used in which tetrahedral meshes for the ignition region near the center and hexahedral meshes for the remaining region were generated. Xiao et al. (2017) used GASFLOW-MPI code with 15,625,000 cells for the same computational domain.
CFL number was set to 0.1. Symmetry boundary conditions for all the flow variables were given for all the symmetric planes and wave transmissive pressure, zero gradient temperature, pressure inlet/outlet velocity boundary condition were applied for the far-field boundary.
The mixture was ignited from the center of the balloon with 50 g of high explosives. The onset of detonation is generally modeled by a pressure pulse instead of simulating detailed chemical reactions (SUSANA, 2016;Xiao et al., 2017). The initial pressure pulse by the high explosive was modeled by the following equation in which the pressure is a function of the explosive mass m and the distance of explosion r (Bajic, Bogdanov, & Jeremic, 2009 (11) Effects of chemical reaction mechanism were investigated implementing 1-, 2-, and 7-step models. Three chemical reaction models were coupled with the Open-FOAM solver and the simulation results are shown in Figure 9. The 7-step model gives much accurate and stable performance compared with 1-, and 2-step models.
The 1-step model shows a little unstable behavior after pressure peak and the 2-step model is even more diffusive than the 1-step model.
Pressure contours at several elapsed times are shown in Figure 10. In Figure 10, the white half circle displays the hemispherical balloon position and shock front passes through the balloon boundary between 10.8 and 11.2 ms. In Figure 11, we compared results using the current method with the experimental data as well as the results from GASFLOW-MPI (Xiao et al., 2017). Pressure profiles at the locations of 1.5, 2.75, 3.25 and 4.0 m off the balloon center are obtained from the OpenFOAM solver with the 7-step chemistry. The current numerical method predicts the detonation quite well and gives similar resolution to GASFLOW-MPI (Xiao et al., 2017) even though current solver used 1/20 less number of grids than the GASFLOW-MPI.
A weak shock appears at 12.6 ms in Figure 11. It was reported that a partial detonation shock was reflected Figure 15. Detonation shock propagation in the RUT test (KI-RUT-HYD05), pressure contours from 101.325 kPa to 3,049.14 kPa with 45 levels.  from the hydrogen-air mixture/air interface at the balloon boundary as the balloon is torn by the shock (Xiao et al., 2017). The weak shock moving inwardly is apparently seen in Figure 10 at time = 11.6 ms.

Detonation in RUT
The OpenFOAM solver was also applied for the 3D RUT benchmark problem (see Figure 12) to simulate complicated shock interactions by the complex geometry of the RUT facility (SUSANA, 2016). Among several experiments, KI-RUT-HYD09 and KI-RUT-HYD05 were chosen for detonation simulation (SUSANA, 2016).
In the case of KI-RUT-HYD09, 25.5% volume hydrogen in air was included. The ignition point is positioned at the end of the curved channel and is located 80 cm above the floor and 50 cm from the wall. Pressure was measured by pressure transducers mounted on the rear wall in the tunnel and canyon as shown in Figure 12 (point 7 ∼ 12). In the experiment, a 200 g high explosive charge was used for ignition. Similar to the KI-RUT-HYD09 case, the KI-RUT-HYD05 case was set to have a volume of 20% volume hydrogen in air. The ignition point is installed at the corner of the canyon whose location is 80 cm above the floor and 50 cm from the wall. Pressure was measured by instruments mounted on the front wall as shown in Figure 12 (point 1 ∼ 6). The high explosive charge was the same as KI-RUT-HYD09. For numerical simulation, the computational domain was discretized with hexahedral cells of 2,211,900 which is much less than 7.7 million cells used (SUSANA, 2016).
CFL number was set to 0.1 and zero gradient conditions for the flow variables were applied for all the wall boundaries except slip boundary condition for the velocity. Figure 13 shows the wave propagation for the KI-RUT-HYD09 case. Detonation starts in the end of the curved channel. Complicated shock interactions appear as time increases. Shock waves reflect from the side walls and form shock stems near the walls and expand when the shock front enters the canyon area. In Figure 14, data from the experiment (Dorofeev et al., 1993) and COM3D (SUSANA, 2016) are compared with those obtained from the OpenFOAM solver at measuring points 11 and 8. The pressure peak is lower than the experiment but the Open-FOAM solver gives similar results to COM3D at point 11 and better results than COM3D at point 8 (SUSANA, 2016). Figure 15 displays the wave propagation of the KI-RUT-HYD09 case. Figure 16 shows comparison results between the current method, experiment and COM3D. It also shows discrepancy between the experiment and the numerical simulations is large than the KI-RUT-HYD09 case.
Discrepancy between the experiment and the numerical simulations might be caused by ignoring support structures. Through the channel and canyon areas, several support structures were installed to protect the facility from the detonation, but for the simulations, the detailed geometry couldn't be reflected due to the absence of the geometrical data. When the shock wave passes through obstacles, the strength of the shock wave can increase due to the change of flow area but the OpenFOAM and COM3D simulations were carried out without reflecting such geometries.

Application to detonation in a generic NPP containment
If there is a possibility of DDT (deflagration to detonation transition) in a compartment filled with a hydrogen mixture due to its geometric characteristic and hydrogen concentration, how severely a detonation impulse can affect the structural integrity must be evaluated. In order to assess the applicability of the developed numerical method to a reactor containment it was applied to a generic containment which has a large free volume of 95,400 m 3 with a simplified internal structure to understand the characteristics of a detonation shock behavior in a dome region. The generic containment is composed of a cylindrical part and a hemi-spherical dome with a diameter of 46 m and a height of 79 m. Figure 17 shows the internal structure and free volume surrounded by its containment wall.
In this simulation, it is assumed that a hydrogen-air mixture cloud with a concentration of 20 volume% is developed only in the dome region and the other region is filled with air as shown in Figure 17. A detonation is initiated near the corner of the dome region. In order to develop an initial detonation in the simulation, a pressure load from 200 g of TNT calculated by Equation (11) is applied at the assumed point.
For the simulation of a detonation shock behavior in a dome region of the generic plant, a cluster of 16 Intel I7 desktop PCs were used. The CPU time spent for the physical time of 0.03 s on a mesh with 4,750, 299 cells (98.5% hexahedral cells) was 14,038 s (about 3.9 h) using 64 CPU cores.  CFL number was set to 0.1 and zero gradient conditions for the flow variables were applied for all the wall boundaries except slip boundary condition for the velocity.
In Figure 18, the behavior of the detonation shock from an initiation is displayed. Initially, a spherical shock builds up, but it becomes an arc-shaped shock by interacting with the dome wall. Hence, a reflected shock forms from the detonation shock. In the case when a hydrogenair mixture cloud is generated in the containment dome region, a detonation shock develops at the corner of the dome, runs along the dome wall, and finishes at the opposite corner with a full combustion of the hydrogen cloud. This propagation characteristic of the detonation shock along the concave wall of the dome creates a reflected shock. Generally, it is known that the strength of a detonation shock becomes bigger when it is reflected on a surface. In order to evaluate the strength of the pressure load along the dome wall, pressures at some locations on the dome wall, depicted in Figure 19, was sampled during the simulation. Figure 20 is the pressure-time histories at the locations. The pressure load on the dome wall at point 3 reaches 11 bars but the pressure load at point 12 reaches 22 bars. It is thought that the main cause of the doubled pressure load at point 12 is strongly related to the development of the reflected shock.

Conclusion
An accurate and robust numerical method for detonation simulation was presented which utilizes the OpenFOAM solver and a 7-step chemical reaction model. Comparing 1-, 2-, 7-step chemical reaction mechanisms, it was figured out that the OpenFOAM solver coupled with the Shang's 7-step model is very accurate and stable for detonative combustion simulation. To validate the capability of the presented numerical method, various cases with detonation have been calculated with the solver and compared with experiment and other simulation results. From the validations, it was figured out that our numerical method can resolve complicated shock structures and flow physics very well with moderate grid points. Therefore, the solver can give accurate pressure peak values for dynamic loading to structures inside the containment building under hydrogen-air mixture detonation. In order to understand the characteristics of detonation shock behavior in a dome region of a containment, our numerical method was applied to a generic containment. It was found that a pressure load from a detonation is magnified by a reflected shock which is made by the concave shape of the dome geometry. We discovered that the simulation of a detonation shock should be carried out considering real geometry and increased strength of the shock due to shock interaction can cause severe impact on the structural integrity of safety-related structures.
In the future, other effects such as steam and turbulence will be implemented to our numerical method for further realistic analysis inside the containment.