## Governing Equation and Simulation Method
### Governing Equation
Cell-averaged quantities of the finite volume method are defined as follows:  
\[\mathbf{\bar{Q}} =\frac{1}{\alpha V}\int_{z_{1}}^{z_{2}}\int_{y_{1}}^{y_{2}}\int_ {x_{1}}^{x_{2}}H\mathbf{Q}\mathrm{d}x\mathrm{d}y\mathrm{d}z, \alpha =\frac{1}{V}\int_{z_{1}}^{z_{2}}\int_{y_{1}}^{y_{2}}\int_{x_{1}}^ {x_{2}}H\mathrm{d}x\mathrm{d}y\mathrm{d}z,\] \[V =(x_{2}-x_{1})(y_{2}-y_{1})(z_{2}-z_{1}),\] \[\mathbf{\mathbf{\mathsf{E}}}_{\mathrm{a}}-\mathbf{\mathsf{E}}_{\mathrm{d}} =\frac{1}{\beta_{x}S_{x}}\int_{z_{1}}^{z_{2}}\int_{y_{1}}^{y_{2}}H( \mathbf{\mathsf{E}}_{\mathrm{a}}-\mathbf{\mathsf{E}}_{\mathrm{d}})\mathrm{d}y \mathrm{d}z, \beta_{x} =\frac{1}{S_{x}}\int_{z_{1}}^{z_{2}}\int_{y_{1}}^{y_{2}}H\mathrm{ d}y\mathrm{d}z,\] \[S_{x} =(y_{2}-y_{1})(z_{2}-z_{1}),\] \[\bar{\mathbf{\mathsf{F}}}_{\mathrm{a}}-\bar{\mathbf{\mathsf{F}}}_{ \mathrm{d}} =\frac{1}{\beta_{y}S_{y}}\int_{z_{1}}^{z_{2}}\int_{x_{1}}^{x_{2}}H( \mathbf{\mathsf{F}}_{\mathrm{a}}-\mathbf{\mathsf{F}}_{\mathrm{d}})\mathrm{d}x \mathrm{d}z, \beta_{y} =\frac{1}{S_{y}}\int_{z_{1}}^{z_{2}}\int_{x_{1}}^{x_{2}}H\mathrm{ d}x\mathrm{d}z, \tag{8}\] \[S_{y} =(x_{2}-x_{1})(z_{2}-z_{1}),\] \[\bar{\mathbf{\mathsf{G}}}_{\mathrm{a}}-\bar{\mathbf{\mathsf{G}}}_{ \mathrm{d}} =\frac{1}{\beta_{z}S_{z}}\int_{y_{1}}^{y_{2}}\int_{x_{1}}^{x_{2}}H( \mathbf{\mathsf{G}}_{\mathrm{a}}-\mathbf{\mathsf{G}}_{\mathrm{d}})\mathrm{d}x \mathrm{d}y, \beta_{z} =\frac{1}{S_{z}}\int_{y_{1}}^{y_{2}}\int_{x_{1}}^{x_{2}}H\mathrm{ d}x\mathrm{d}y,\] \[S_{z} =(x_{2}-x_{1})(y_{2}-y_{1}),\] \[\bar{\mathbf{\sigma}} =\frac{1}{S_{\text{cut-plane}}}=\iint_{\text{cut-plane}}\mathbf{ \sigma}\mathrm{d}S, S_{\text{cut-plane}} =\iint_{\text{cut-plane}}\mathrm{d}S,\]  
where \(\alpha\) is the volume fraction of the air in the cell; \(\beta_{x}\), \(\beta_{y}\), \(\beta_{z}\) are the area fractions of the air in each cell interface; \(V\) is the volume of the cell; and \(S_{x}\), \(S_{y}\), and \(S_{z}\) are the area of each cell interfaces. Using these cell-averaged values, governing Equation (6) is reduced to the following form:  
\[\frac{\mathrm{d}\alpha\bar{\mathbf{\mathsf{Q}}}}{\mathrm{d}t}+\frac{[\beta_{x}(\bm {\mathsf{E}}_{\mathrm{a}}-\mathbf{\mathsf{E}}_{\mathrm{d}})]_{x_{1}}^{x_{2}}}{x_{ 2}-x_{1}}+\frac{[\beta_{y}(\bar{\mathbf{\mathsf{F}}}_{\mathrm{a}}-\bar{\mathbf{\mathsf{ F}}}_{\mathrm{d}})]_{y_{1}}^{y_{2}}}{y_{2}-y_{1}}+\frac{[\beta_{z}(\bar{\mathbf{\mathsf{G}}}_{ \mathrm{a}}-\bar{\mathbf{\mathsf{G}}}_{\mathrm{d}})]_{z_{1}}^{z_{2}}}{z_{2}-z_{1} }=\frac{S_{\text{cut-plane}}}{V}\bar{\mathbf{\sigma}}. \tag{9}\]  
The projection of the cut plane to each cell interface is approximated in the following form:  
\[V =(x_{2}-x_{1})(y_{2}-y_{1})(z_{2}-z_{1}),\]0  
These approximation tend to be exact equations at the limit of small curvature of the cut plane. Thus in the case of the sufficiently small control volume, \(S_{\text{cut-plane}}\) and \(\mathbf{n}\) are evaluated as follows:  
\[V =(x_{2}-x_{1})(y_{2}-y_{1})(z_{2}-z_{1}),\]1  
\[V =(x_{2}-x_{1})(y_{2}-y_{1})(z_{2}-z_{1}),\]2The advection flux at the cell interface is calculated using the Simple Low-dissipation AUSM (SLAU) scheme [17]. The primitive variables on the cell interface are reconstructed using Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) of the fifth-order accuracy [18]. The velocity components are obtained using the modification of the MUSCL scheme proposed by [PERSON] [19]. The time integration to second-order accuracy is calculated using the Total Variation Diminishing (TVD) Runge-Kutta method [20].