## 2 Methodology
### Multiphase Flow Model and Turbulent Model  
The Volume of Fluid (VOF) method is a widely used multiphase flow method. This method provides a sharp representation of fluid interfaces, which is particularly valuable for capturing details involving complex flows. The specific governing equations are as follows [17]:  
The mass conservation equation is  
\[\frac{\partial\rho}{\partial t}+\
abla\cdot\left(\rho\overline{U}\right)=0 \tag{2}\]  
The momentum equation can be expressed as  
\[\frac{\partial(\rho\overline{U})}{\partial t}+\
abla\cdot\left(\rho\overline{ U}\overline{U}\right)=-\
abla P+\
abla\Big{[}\mu\Big{(}\
abla\overline{U}+ \
abla\overline{U}^{T}\Big{)}\Big{]}+\rho\overline{g} \tag{3}\]  
where \(\rho\) is the density of the fluid, \(\overline{U}\) is the velocity vector of the fluid, \(\mu\) represents the kinematic viscosity coefficient, \(P\) represents pressure, and \(\overline{g}\) is the gravitational acceleration.  
The VOF method exploits a phase fraction \(F_{P}\) to represent the volume fraction. It ranges from zero to one, where \(F_{P}=0\) represents air, \(F_{P}=1\) represents slurry, and \(0<F_{P}<1\)  
Figure 1: Schematic of the numerical model and pressure type.  
represents a two-phase interface. The phase transformation equation and the properties of the fluid can be written as  
\[\frac{\partial F_{P}}{\partial t}+\mathbf{U}\cdot\
abla F_{P}=0 \tag{4}\]  
\[\rho=(1-F_{P})\rho_{A}+F_{P}\rho_{P}\text{ and }\mu=(1-F_{P})\mu_{A}+F_{P}\mu_{P} \tag{5}\]  
Due to the fluctuating pressure, turbulent behavior occurs within the computational domain.  
The realizable k-\(\varepsilon\) model enhances the physical realism of the turbulence predictions by incorporating more accurate terms and coefficients [18]. It also takes into account the anisotropic nature of turbulence, which is especially important for flows with significant directional variations. Therefore, this model is adopted in this study, and the expressions are as follows [18]:  
\[\frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_{j}}(\rho ku_{ j})=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{t}}{\sigma_{k}} \right)\frac{\partial k}{\partial x_{j}}\right]+G_{k}+G_{b}-\rho\varepsilon-Y _{M}+S_{k} \tag{6}\]  
\[\frac{\partial}{\partial t}(\rho\varepsilon)+\frac{\partial}{\partial x_{j}} \big{(}\rho\varepsilon u_{j}\big{)}\delta\varepsilon=\frac{\partial}{ \partial x_{j}}\Bigg{[}\left(\mu+\frac{\mu_{t}}{\sigma_{x}}\right)\frac{ \partial\varepsilon}{\partial x_{j}}\Bigg{]}+\rho C_{i}S\varepsilon-\rho C _{2}\frac{\varepsilon^{2}}{k+\sqrt{\upsilon\varepsilon}}+C_{1}\varepsilon \frac{\varepsilon}{k}C_{3}G_{b}+S_{\varepsilon} \tag{7}\]  
In the equation, \(G_{k}\) represents the generation term of turbulent kinetic energy caused by Reynolds-averaged velocity gradients. \(\sigma_{k}\) and \(\sigma_{\varepsilon}\) are the Prandtl numbers for turbulent kinetic energy and dissipation rate, respectively. \(G_{b}\) represents the turbulent kinetic energy (J) generated due to buoyancy, and \(Y_{M}\) represents the contribution of fluctuation expansion in compressible turbulence to the overall dissipation rate. \(S_{k}\) and \(S_{\varepsilon}\) are user-defined source terms. \(C_{1\varepsilon}\), \(C_{2\varepsilon}\), and \(C_{3\varepsilon}\) represent the numerical constants.