## 3 Numerical Method and Validation  
### Governing Equation  
For the calculation, which was based on a solid-liquid two-phase flow, the governing equations of the liquid phase and the solid phase were, respectively taken into account. For the fluid phase, which was treated as incompressible, the Reynolds-Averaged Navier-Stokes (RANS) was employed in the simulation. The mass equation and momentum Equation were as follows [25]:  
\[\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x_{i}}(\rho\overline{ u_{i}})=0, \tag{9}\]  
\[\frac{\partial}{\partial t}(\rho\overline{u_{i}})+\frac{\partial}{\partial x_ {i}}(\rho\overline{u_{i}u_{i}})=\frac{\partial\rho}{\partial x_{i}}+\frac{ \partial\rho}{\partial x_{j}}\left(\mu\frac{\partial\overline{u_{i}}}{\partial x _{j}}-\rho\overline{u^{\prime}_{i}u^{\prime}_{j}}\right)+F_{i}. \tag{10}\]  
where \(\rho\) and \(\mu\) represent the density and dynamic viscosity of the fluid, respectively; \(\overline{u_{i}}\) and \(\overline{u^{\prime}_{i}}\) are the mean components of velocity and pulsation component of the velocity in the \(i\)-directions, respectively; \(p\) is the total pressure; \(F_{i}\) denotes the external body force components in the \(i\)-direction; and \(-\rho\overline{u^{\prime}_{i}u^{\prime}}\) denotes the Reynolds stress.  
To simulate the turbulent flow, the realizable \(k-\epsilon\) model was employed thanks to the stability and precision it demonstrates when dealing with pressure gradients [26]. The transport Equations were expressed as follows [25]:  
\[\frac{\partial(\rho k)}{\partial t}+\frac{\partial(\rho ku_{i})}{\partial x_ {i}}=\frac{\partial}{\partial x_{i}}\Bigg{[}\bigg{(}\mu+\frac{\mu_{t}}{ \sigma_{k}}\bigg{)}\frac{\partial k}{\partial x_{j}}\Bigg{]}+G_{k}+G_{b}-\rho \varepsilon-Y_{M}, \tag{11}\]  
\[\frac{\partial(\rho\varepsilon)}{\partial t}+\frac{\partial(\rho\varepsilon u _{i})}{\partial x_{i}}=\frac{\partial}{\partial x_{i}}\Bigg{[}\bigg{(}\mu+ \frac{\mu_{t}}{\sigma_{\varepsilon}}\bigg{)}\frac{\partial\varepsilon}{ \partial x_{j}}\Bigg{]}+\rho C_{1}S\varepsilon-\rho C_{2}\frac{\varepsilon^{ 2}}{k+\sqrt{\
u\varepsilon}}+C_{1\varepsilon}\frac{\varepsilon}{K}C_{3 \varepsilon}G_{b}. \tag{12}\]  
where \(k\) is the turbulent kinetic energy; \(\varepsilon\) is the dissipation rate; \(G_{k}\) and \(G_{b}\) represent the turbulent kinetic energy generated from the average speed gradients and buoyancy effect, respectively; \(S\) is the average strain rate; \(Y_{M}\) is the overall dissipation rate due to the fluctuating dilatation; \(C_{1}=\max[0.43,\eta/(\eta+5)]\) and \(\eta=S\cdot k/\varepsilon\); \(C_{3\varepsilon}\) is a function representing the effect of buoyancy on the dissipation rate; \(C_{2}\) and \(C_{1\varepsilon}\) are constants; \(\sigma_{k}\) and \(\sigma_{\varepsilon}\) denote the turbulent Prandtl number for \(k\) and \(\varepsilon\), respectively; and, generally, \(C_{1\varepsilon}=1.44\), \(C_{2}=1.9\), \(\sigma_{k}=1.0\) and \(\sigma_{\varepsilon}=1.2\) were the default values in the software.  
For the governing equation of the discrete phase particle, the motion trajectory was predicted by integrating its force balance in a Lagrangian reference system. The force balance equation can be written as [27]:  
\[m_{p}\frac{d\vec{u}_{p}}{dt}=m_{p}\frac{\vec{u}-\vec{u}_{p}}{\tau_{r}}+m_{p} \frac{\vec{g}\big{(}\rho_{p}-\rho\big{)}}{\rho_{p}}+\vec{F}. \tag{13}\]  
where \(m_{p}\) is the mass of the discrete solid particle, \(\vec{u}\) is the velocity of the fluid phase, \(\vec{u}_{p}\) is the velocity of the discrete solid particle, \(\rho_{p}\) and \(\rho\) are the densities of discrete particles and continuous fluid, respectively, \(m_{p}\frac{\vec{u}-\vec{u}_{p}}{\tau_{r}}\) is the drag force, and \(\vec{F}\) is the additional force. \(\vec{F}\) includes the virtual mass force, thermophoretic force, Brownian force, Saffman's lift force, etc.