# Numerical Simulation of Natural Gas/Hydrogen Combustion in a Novel Laboratory Combestor
## 2 Materials and Methods
### Numerical Methods and Computational Details  
The Simcenter Star CCM+ code [54], based on the finite-volume method, was used to solve the governing equations. The convective fluxes of these equations were determined using the second-order upwind scheme to calculate cell face values. The variable gradients were computed using a hybrid Gauss-least squares method and the gradient limiter proposed by [PERSON] [55]. The SIMPLE algorithm was employed for pressure-velocity coupling, and the algebraic systems of discretized equations were solved using the algebraic multigrid linear solver. Under-relaxation factors of 0.7 and 0.3 were used for the momentum and pressure, respectively, 0.8 for the turbulent kinetic energy and its dissipation rate, and 1.0 for the energy and species mass fraction equations.  
Three different meshes were tested, consisting of approximately 128, 300, and 700 thousand control volumes. This corresponds to a refinement ratio of \(r=1.33\), defined as the cubic root of the ratio of the number of control volumes between two successive grids.  
\begin{table}
\begin{tabular}{c c c c c c} \hline
**\%H\({}_{2}\)** & **Excess Air** & **Fraction of** & **Fuel Mass** & **Primary Air** & **Secondary Air** \\
**(vol.) in** & **Coefficient** & **Primary Air** & **Flow Rate** & **Mass Flow** & **Mass Flow Rate** \\
**Fuel** & **(\(\lambda\))** & **(\%), (\(\alpha\))** & **(kg/s)** & **Rate (kg/s)** & **(kg/s)** \\ \hline
0 & 1.3 & 50 & \(1.29\times 10^{-5}\) & \(1.44\times 10^{-4}\) & \(1.44\  
Figure 2 shows the temperature profiles along the centreline of the combustor and along a radial profile at the exit, in the middle axial plane of the computational domain, for pure natural gas, \(\lambda\) = 1.7 and \(\alpha\) = 50%. The estimated order of convergence was found to be \(p\) = 1.77 and \(p\) = 1.84 for two points along the centreline, at distances of 210 and 345 mm from the burner exit, respectively. These values are a little lower than the theoretical value of \(p\) = 2, which is expected for a second-order spatial discretization scheme. This deviation is not surprising, as the actual order of convergence is influenced by factors such as grid quality and the non-linearities of the governing equations.  
Additionally, the Grid Convergence Index (GCI), as defined by [PERSON] [56], was calculated for the temperature at the two points mentioned above. The following results were obtained: GCI\({}_{2,3}\) = 0.20% and GCI\({}_{1,2}\) = 0.20% for the point closer to the burner and GCI\({}_{2,3}\) = 3.43% and GCI\({}_{1,2}\) = 2.05% for the other point, where subscripts 1, 2, and 3 denote the fine, medium, and coarse grids, respectively. These values are relatively small, especially the former ones, revealing that the solution is weakly sensitive to the grid size. The ratio GCI\({}_{2,3}\) (\(t^{p}\) GCI\({}_{1,2}\)) was found to be very close to unity, namely 0.999 for one point and 0.989 for the other. This indicates that the asymptotic range of convergence was reached. Based on these results, the medium mesh, shown in Figure 3, was selected for the calculations presented in the remainder of this paper.  
Figure 3: Mesh used in the discretization of the computational domain.  
Figure 2: Influence of the grid size on the axial temperature profile.