## 3 Simulation Setup
The filtered transport equation for instantaneous species is described as follows [29]:  
\[\frac{\partial\overline{\rho}Y_{k}}{\partial t}+\
abla\cdot(\overline{\rho} \overline{u}\gamma_{k})=\
abla\cdot(\overline{\rho}\Xi_{\Delta}FD_{k}\
abla Y _{k})+\frac{\Xi_{\Delta}\dot{\omega}_{k}(Q)}{F}, \tag{1}\]  
where \(\rho\) and \(u\) are the density and velocity vector, respectively; \(\dot{\omega_{k}}\) is the chemical source term of \(k\); \(F\) is the thickening factor. The flame can be resolved in the LES grid with a proper thickening factor. \(D_{k}\) is the mixture-averaged molecular diffusion coefficient, calculated from the following [31]:  
\[D_{k}=\frac{1-x_{k}}{\sum_{j\
eq k}^{N}\left(\frac{x_{j}}{D_{kj}}\right)}, \tag{2}\]  
where \(x_{k}\) is the mole fraction of species \(k\) and \(D_{jk}\) is the binary diffusion coefficient of species \(j\) and \(k\). \(\Xi_{\Delta}\) is the sub-grid flame wrinkling factor. The wrinkling factor is used to compensate for the lost flame surface due to the thickening process. This study employed the power-law flame wrinkling model shown in [32; 33]:  
\[\Xi_{\Delta}=\bigg{\{}1+min\bigg{[}max\bigg{(}\frac{\Delta}{\delta_{l}}-1,0 \bigg{)},\Gamma\bigg{(}\frac{\Delta}{\delta_{l}},\frac{u_{\Delta}^{\prime}}{S _{l}},Re_{\Delta}\bigg{)}\frac{u_{\Delta}^{\prime}}{S_{l}}\bigg{]}\bigg{\}}^{ \beta}, \tag{3}\]  
where \(u_{\Delta}^{\prime}\) is the sub-grid-scale turbulent velocity and \(Re_{\Delta}\) is the sub-grid-scale Reynolds number; \(\beta\) is the model coefficient. The thickening factor \(F\) is often associated with the mesh size \(\Delta_{x}\), given by [34]:  
\[F(n)=\frac{n\Delta_{x}}{\delta_{l}^{0}}, \tag{4}\]  
where \(\delta_{l}^{0}\) is the laminar flame thickness, which can be calculated from a 0D simulation using CHEMKIN-PRO [35]. The number of grids \(n\) is typically 5-10 [36]. The heat release rate method was used to select the flame region to determine the thickening factor dynamically. The dynamical control method of \(F\) is shown by the following [37]:  
\[F_{dyn}=1+(F(n)-1)tanh\bigg{\{}const\text{-}min\bigg{[}1,max\bigg{(}0,(\frac{ HRR-minH}{maxH-minH})\bigg{)}\bigg{]}\bigg{\}}, \tag{5}\]  
where \(minH\) and \(maxH\) are the minimum and maximum values of the heat release rate (calculated by one-dimension laminar simulation), respectively. The \(F_{dyn}\) depends on time and space, which are recalculated in each time step. Table 1 shows the parameters used in the thickened flame model. TF modeling was performed using OpenFOAM 2.3.0 open-source software, which has been discussed in detail and validated in one of our previous works [38]. The governing equations were solved using a pressure-based finite volume method code in OpenFOAM 2.3.0. The Pressure Implicit Split Operator (PISO) algorithm was used to handle the pressure-velocity coupling term [39]. The transient term was discretized by the implicit Euler scheme. Second-order linear schemes were used to discretize the convection and diffusion terms. The fixed time step \(\Delta t\) was set at \(1\times 10^{-6}\) s, and the maximum Courant number was always less than 0.3 during the simulations. This study employed the 30-species skeletal reaction mechanism proposed by [PERSON] et al. [40]. The validation of this mechanism is described in detail in Appendix B.