The CYCAS code solves transient 3D Navier–Stokes (N-S) equations. The conservation equations are listed below: V o l u m e   e q u a t i o n : ∂ V ∂ t = ∇ V ⋅ ( − u ⇀ ) , (1) M a s s   e q u a t i o n : ∂ ρ ∂ t = ∇ ⋅ [ − ρ u ⇀ ] + S ρ , (2) M o m e n t u m   e q u a t i o n : ∂ ( ρ u ⇀ ) ∂ t = ∇ ⋅ [ − ρ u ⇀ ( u ⇀ ) ] − ∇ p − ∇ ⋅ σ ⇀ + ρ g ⇀ + S m , (3) I n t e r n a l   e n e r g y : ∂ ( ρ I ) ∂ t = ∇ ⋅ [ − ρ I ( u ⇀ ) ] − p ∇ ⋅ u − ∇ ⋅ q ⇀ + S I , (4) where u ⇀ is the fluid velocity vector, V is the discretized fluid control volume, I is the internal energy, and p , σ ⇀ , and g ⇀ are the fluid pressure, viscous stress tensor, and gravity acceleration, respectively. S ρ , S m , and S I are the external source terms for each equation caused by phase change, chemical reaction, and so on.

The secondary order accuracy discrete scheme, i.e., van Leer MUSCL algorithm is adopted. The all-speed implicit method, i.e., Implicit Continuous Eulerian–Arbitrary Lagrangian Eulerian (Hirt et al., 1997), is adopted to solve the transport equations. The detailed information on the transport equations of CYCAS can be found in Chen et al. (2016).

The models for the phenomenon that are relevant to hydrogen safety in the containment are validated with both domestic and international projects, namely the turbulence, convection, and condensation model, as listed in Table 1. Two validation cases, the HyJet test on the BMC facility and the ISP-23 test on the HDR facility are presented in Chen et al. (2016). The other tests are included in the internal reports. Figure 2 gives the geometry models and comparison between the experimental data of the three cases, i.e., OECD SETH Test 9, PANDA HYMERES, and H2PAR, the experimental details of which can be found in Auban et al. (2007), Mimouni et al. (2011), and Paladino et al. (2014), respectively.

The BVM describes the development of the flame with the combustion process variable ξ by solving the transportation equation for ξ , as shown in Eq. 7. The second-order accuracy scheme is adopted in the spatial domain, while the first order in the time domain. The combustion source term ∫ ​ ρ u n S ξ d V is described as a function of the turbulence flame speed, as shown in Eq. 8. The Zimont model is selected for the turbulence flame speed model, as given in Eq. 9 (Zimont and Lipatnikov, 1995; Zimont, 2000b). The turbulence flame speed S t takes the higher value between the laminar flame speed S l and the Zimont correlation, which means that the flame is laminar for low turbulence intensity. d d t ∫ ​ ρ ξ d V = ∮ ​ [ − ρ u → ξ + ( μ S c + μ t S c t ) ∇ ξ ] d S + ∫ ​ ρ u n S ξ d V , (7) ∫ ​ ρ u n S ξ d V = ρ u n | ∇ ξ | S t , (8) S t = 0.52   u ' D a 0.25 , (9) where ξ is the combustion process variable: ξ =   0 , indicating the unburnt gas mixture and ξ =   1 , indicating the completely burnt gas. μ is the viscosity and S c is the Schmidt number; the subscript t indicates the turbulence value. ρ u n is the density for the unburnt gas mixture and S t is the turbulence flame speed. Da is the Damköhler number, of which the definition equation is shown in Eq. 10. D a = ( l / u ' ) ( l f / S l ) , (10) where l is the integral turbulence length, u ' is the fluctuation velocity, l f is the flame thickness, and S l is the laminar flame speed, with the formation as shown in Eq. 11. S l ( T , x H 2 O ) = S l ( 300 ) F ( T ) F ( x h 2 o ) F ( P ) . (11)

The terms of S l ( 300 ) is taken from Szabó et al. (2012), the terms F ( T ) and F ( x h 2 o ) are from Koroll et al. (1993), and the term of F ( P ) is from Hu et al. (2009). By doing so, the validation range of Eq. 11 can envelop the atmospheric condition during the accident conditions, as listed below:

Stoichiometric ratio [ 0.1252 ,   6.5 ] ;

The spark ignition is modeled as an external source term to the transport equation of ξ , as shown below (Zimont and Lipatnikov, 1995). The source term indicates that the gas mixture in the ignition cell is ignited by a pulse of spark. Moreover, the quasi-laminar flame is assumed during the ignition stage according to Zimont and Lipatnikov (1995). d d t ∫ ρ ξ d V = ρ u n W i g n = ρ u n ⁡ exp { − [ ( t − t 0 σ t ) 2 ] } ( 1 − ξ ) , (12) where t 0 is the time duration of one spark ignition, σ t = t 0 / 5 .

The geometry model of the THAI facility is built in cylindrical coordination, with the computational region radius of 1.6 m and height of 9.2 m, as shown in Figure 3. The structure is modeled with four layers with different materials for each layer according to NEA (2011). The models of convective, condensation, and radiation heat transfers are calculated for the heat loss on the structure surface. Two cases from the THAI facility (NEA, 2011) are calculated, one is the HD-12 test without steam, and the other is the HD-22 test with steam.

The initial conditions are listed as follows (NEA, 2011):

HD-12: atmospheric pressure at 1.485 bar, gas temperature at 291 K, hydrogen at 7.98 vol%.

The ignition is located at 0.7 m from the bottom of the facility. The initial turbulence condition is set as k o = 1.5 e − 4 m 2 / s 2 and ε o = 4.8 e − 4 m 2 / s 3 , based on Sathiah et al. (2016). Both the convective and radioactive heat transfers to the structures are considered in this calculation.

Firstly, a mesh sensitivity analysis is conducted with three cases with meshes of 35 × 24 × 292, 35 × 24 × 146, and 16 × 24 × 73 grids. The result shows that the first two cases are in accord with each other, while the third case underestimates the peak pressure when compared with the first two cases as shown in Figure 5. Therefore, the second mesh, which is 35 × 24 × 146 grids in radial, peripheral, and axial directions, respectively, with a total mesh number of 122,640 is used for the next HD-22 analysis.

The r-z view of the temperature at 2 s after the ignition is shown in Figure 4 The pressure and flame position are compared with the experimental data as shown in Figures 5, 6. Figure 5 shows that the peak pressure is overestimated in the CYCAS calculation. The pressure decreases at a smaller slope after the peak, indicating that the heat loss of the structure is underestimated. This could be the reason of the overestimation of the peak pressure. The combustion model in CYCAS overestimates the flame speed in both cases: it takes about 5.8 s to reach the ceiling in the experiment, while the time duration is 4.73 s in the CYCAS calculation for the HD-12 calculation. It takes about 3.9 s to reach the ceiling in the experiment, while the time duration is 4.56 s in the CYCAS calculation for the HD-22 calculation. The same overestimates also occur in the FLUENT calculation for HD-12 in the study by Sathiah et al. (2016), as well as the calculation for HD-22 in NEA (2011). This is mainly due to the model limitation on the flame development from the laminar flame to turbulent flame (NEA, 2011). Despite this limitation, the combustion model in CYCAS manages to have reasonable prediction on the slow deflagration for gas mixture, both with and without steam.