## 3 Numerical Method  
For two dimensional cases, the CFD model for gas phase has the following component form:  
\[\frac{\partial\varepsilon}{\partial t}+\frac{\partial(\varepsilon u_{x})}{ \partial x}+\frac{\partial\Big{(}\varepsilon u_{y}\Big{)}}{\partial y}=0 \tag{47a}\] \[ \begin{array}{c}\rho_{g}\frac{\partial(\varepsilon u_{x})}{ \partial t}+\rho_{g}\frac{\partial(\varepsilon u_{x}u_{y})}{\partial x}+\rho_ {g}\frac{\partial(\varepsilon u_{x}u_{y})}{\partial y}\\ \qquad\qquad=-\varepsilon\frac{\partial p}{\partial x}-S_{px}+\frac{\mu_{g}}{ 3}\bigg{[}\frac{\partial}{\partial x}\big{(}\varepsilon\frac{\partial u_{x}} {\partial x}\big{)}+\frac{\partial}{\partial x}\bigg{(}\varepsilon\frac{ \partial u_{y}}{\partial y}\bigg{)}\bigg{]}+\mu_{g}\bigg{[}\frac{\partial}{ \partial x}\bigg{(}\varepsilon\frac{\partial u_{x}}{\partial x}\bigg{)}+ \frac{\partial}{\partial y}\bigg{(}\varepsilon\frac{\partial u_{x}}{\partial y }\bigg{)}\bigg{]}\\ \rho_{g}\frac{\partial(\varepsilon u_{y})}{\partial t}+\rho_{g}\frac{\partial( \varepsilon u_{y}u_{x})}{\partial x}+\rho_{g}\frac{\partial(\varepsilon u_{y }u_{y})}{\partial y}\\ \qquad\qquad=-\varepsilon\frac{\partial p}{\partial x}-S_{py}+\frac{\mu_{g}}{ 3}\bigg{[}\frac{\partial}{\partial y}\bigg{(}\varepsilon\frac{\partial u_{x} }{\partial x}\bigg{)}+\frac{\partial}{\partial y}\bigg{(}\varepsilon\frac{ \partial u_{y}}{\partial y}\bigg{)}\bigg{]}+\mu_{g}\bigg{[}\frac{\partial}{ \partial x}\bigg{(}\varepsilon\frac{\partial u_{y}}{\partial x}\bigg{)}+ \frac{\partial}{\partial y}\bigg{(}\varepsilon\frac{\partial u_{y}}{\partial y }\bigg{)}\bigg{]}-\varepsilon\rho_{g}g\end{array}  \tag{47c}\]  
where the subscripts \(x\) and \(y\) refers to the components in \(x\) and \(y\) directions, respectively. The finite volume method (FVM) with non-staggered grids is used to discretize Equation (47). Figure 3 illustrates the control volume for a grid node \(P\). By using this control volume, the discrete forms of Equations (47b) and (47c) at the grid node \(P\) can be uniformly written as:  
\[a_{P}\mathbf{u}_{P}=a_{E}\mathbf{u}_{E}+a_{W}\mathbf{u}_{W}+a_{N}\mathbf{u}_{ N}+a_{S}\mathbf{u}_{S}+\mathbf{B} \tag{48}\]where the components of the source term \(\mathbf{B}\) are given by  
\[B_{x}=\frac{\varepsilon_{P}\Delta x\Delta y}{\Delta t}(u_{x}^{0})_{P} -\frac{\varepsilon_{P}\Delta y}{\rho_{g}}(p_{e}-p_{w})-\frac{(S_{P _{P}})_{p}\Delta x\Delta y}{\rho_{g}} \tag{49a}\] \[+\frac{\mu_{g}}{\beta_{g}}\frac{\Delta y}{\Delta x}\left[ \varepsilon_{e}((u_{x})_{P}-(u_{x})_{P})-\varepsilon_{w}((u_{x})_{P}-(u_{x})_{ W})\right]\] \[+\frac{\mu_{g}}{\beta_{g}}\left[\varepsilon_{e}\big{(}(u_{y})_{P_{4}}-(u_{ y})_{P_{3}}\big{)}-\varepsilon_{w}\big{(}(u_{y})_{P_{1}}-(u_{y})_{P_{2}}\big{)}\right]\] \[B_{y}=\frac{\varepsilon_{P}\Delta x\Delta y}{\Delta t}(u_{y}^{0})_{P} -\frac{\varepsilon_{P}\Delta x}{\rho_{g}}(p_{n}-p_{s})-\frac{(S_{P _{P}})_{p}\Delta x\Delta y}{\rho_{g}}-\varepsilon_{P}\delta\Delta x\Delta y +\frac{\mu_{g}}{\beta_{g}}\frac{\Delta y}{\Delta y}\left[\varepsilon_{n} \big{(}(u_{y})_{N}-(u_{y})_{P}\big{)}-\varepsilon_{s}\big{(}(u_{y})_{P}-(u_{y })_{S}\big{)}\right]\] (49b) \[+\frac{\mu_{g}}{\beta_{g}}\left[\varepsilon_{n}\big{(}(u_{x})_{P_{4}}-(u_{x})_{ P_{1}}\big{)}-\varepsilon_{s}\big{(}(u_{x})_{P_{3}}-(u_{x})_{P_{2}}\big{)}\right]\]  
Here \(u_{x}^{0}\) and \(u_{y}^{0}\) are the gas velocities at the previous time step.