## 2 Methodology
### Simulation Domain and Mesh  
Explosion simulations were conducted using the FLACS-CFD, which employs a structured Cartesian grid for computational fluid dynamics (CFD) analysis. The computational domain was divided into two regions: a core domain surrounding the fuel preparation room and a stretch domain extending outward to minimize boundary reflections and ensure numerical stability.  
A uniform grid with an aspect ratio of 1:1 was applied within the core region to achieve high spatial resolution. The outer stretch domain was given a gradually increasing mesh size to reduce the computational load. The overall domain measured 30 m \(\times\) 45 m \(\times\) 15 m in the X-, Y-, and Z-directions, respectively, extending more than three times the dimensions of the FPR, to allow the complete dissipation of pressure waves before they reached the boundaries. The mesh strategy and domain layout are shown in Figure 4.  
Figure 4: Generation of grids for explosion scenario.  
A mesh convergence study was performed using four grid sizes: 0.1 m, 0.2 m, 0.5 m, and 1.0 m. As shown in Figure 5, the peak explosion pressures varied depending on the mesh size, with a maximum difference of approximately 3.0 bar between the coarsest and finest grids. When the grid size was refined to 0.2 m or smaller, this variation became negligible. Therefore, 0.2 m was selected as the final resolution, resulting in 292,110 control volumes, a maximum aspect ratio of 12.45, and a maximum relative difference of 19.24%, all of which are within acceptable limits for engineering applications.  
The governing equations were solved using the finite volume method on a staggered grid for velocity and a cell-centred grid for scalar quantities, such as pressure, temperature, and species mass fractions. The model is based on the compressible Favre-averaged Navier-Stokes equations, which include the conservation laws for mass, momentum, energy, and species. The turbulence was modelled using a standard two-equation k-\(\epsilon\) approach, and combustion was represented with a premixed turbulent flame model that incorporated laminar flamelet correlations and flame surface area enhancement mechanisms.  
The numerical scheme offers second-order spatial accuracy and first-order (optionally, second-order) temporal accuracy. The SIMPLE algorithm was employed for pressure-velocity coupling and extended to account for compressibility. The time stepping followed the Courant-Friedrichs-Lewy (CFL) criteria with CFLC = 5 and CFLV = 0.5, and simulations proceeded until the mass residual decreased below \(10^{-4}\).  
Euler-type open boundaries were applied in all directions to allow free pressure wave propagation and to suppress artificial reflections, whereas solid wall boundaries were defined for the enclosure. Complex internal geometries and unresolved obstructions were modelled using the distributed porosity approach, while accounting for localized flow resistance, turbulence generation, and flame-folding effects. Standard wall functions were used for the near-wall treatment.  
To evaluate the explosion dynamics, 11 monitoring points were installed within the FPR, as illustrated in Figure 6 and listed in Table 4. These sensors recorded the pressure, temperature, and species mass fractions over time, enabling the quantitative analysis of flame propagation, overpressure characteristics, and combustion behaviour.  
Additionally, radiative heat loss sub-models were activated to account for the quenching effects near the wall-adjacent ignition points, in accordance with the FLACS-CFD recommendations. Although direct flame extinction was not modelled, the use of appropriate sub-models and a fine mesh resolution (0.2 m) indirectly represented thermal losses  
Figure 5: Maximum pressure according to the grid size.  
near solid boundaries while satisfying the minimum requirement for sub-grid modelling (\(\geq\)1 cm) [33].