## 2 Methodology
### Governing Equations
\tag{9}\]  
The stress tensor of each phase is calculated as follows:  
\[\overline{\overline{\tau}_{\mathrm{I}}}=\mu_{\mathrm{I}}\bigg{\{}\bigg{[} \
abla\overset{\rightarrow}{u_{\mathrm{I}}}+(\
abla\overset{\rightarrow}{u _{\mathrm{I}}})^{T}\bigg{]}-\frac{2}{3}(\
abla\cdot\overset{\rightarrow}{u_ {\mathrm{I}}})^{\overset{=}{}{I}}\bigg{\}} \tag{10}\]  
\[\overline{\overline{\tau}_{\mathrm{s}}}=\mu_{\mathrm{s}}\bigg{[}\
abla \overset{\rightarrow}{u_{\mathrm{s}}}+(\
abla\overset{\rightarrow}{u_{ \mathrm{s}}})^{T}\bigg{]}+(\xi_{\mathrm{s}}-\frac{2}{3}\mu_{\mathrm{s}})( \
abla\cdot\overset{\rightarrow}{u_{\mathrm{s}}})^{\overset{=}{}{I}} \tag{11}\]  
where \(\overset{=}{I}\) is the unit tensor, \(\mu\) is the shear viscosity, \(\xi_{\mathrm{s}}\) is the bulk viscosity of solids, and \(p_{\mathrm{s}}\) is the solids pressure. The solids pressure is defined as follows:  
\[p_{\mathrm{s}}=\alpha_{\mathrm{s}}\rho_{\mathrm{s}}\rho_{\mathrm{s}}+2\rho_{ \mathrm{s}}(1+e)\alpha_{\mathrm{s}}^{2}\rho_{\mathrm{0}}\rho_{\mathrm{s}} \tag{12}\]  
where \(e\) is the restitution coefficient, \(\rho_{\mathrm{s}}\) is the granular temperature, and \(g_{\mathrm{0}}\) is the radial distribution function and can be written as follows:  
\[g_{\mathrm{0}}=\left[1-\left(\frac{\alpha_{\mathrm{s}}}{\alpha_{\mathrm{smax}} }\right)^{\frac{1}{3}}\right]^{-1} \tag{13}\]  
The bulk viscosity of solids is defined as follows:  
\[\zeta_{\mathrm{s}}=\frac{4}{3}\alpha_{\mathrm{s}}^{2}\rho_{\mathrm{s}}d_{ \mathrm{s}}g_{\mathrm{0}}(1+e)\sqrt{\frac{\Theta_{\mathrm{s}}}{\pi}} \tag{14}\]  
The solids shear viscosity is defined as follows:  
\[\mu_{\mathrm{s}}=\frac{4}{5}\alpha_{\mathrm{s}}\rho_{\mathrm{s}}d_{\mathrm{s}}g_ {\mathrm{0}}(1+e)\sqrt{\frac{\Theta_{\mathrm{s}}}{\pi}}+\frac{10\rho_{\mathrm{s}}d _{\mathrm{s}}\sqrt{\pi\Theta_{\mathrm{s}}}}{96(1+e)\alpha_{\mathrm{s}}g_{ \mathrm{0}}}\left[1+\frac{4}{5}g_{\mathrm{0}}\alpha_{\mathrm{s}}(1+e)\right]^{2} \tag{15}\]The granular temperature transport equation is calculated as follows:  
\[0=(-p_{\mathrm{s}}\overline{\overline{1}}+\overline{\overline{\tau_{\mathrm{s}}}}): \
abla\overset{\rightarrow}{u_{\mathrm{s}}}-\gamma_{\Theta_{\mathrm{s}}}+\phi_{ \mathrm{ls}} \tag{16}\]  
where \((-p_{\mathrm{s}}\overline{\overline{1}}+\overline{\overline{\tau_{\mathrm{s}}}}) :\
abla\overset{\rightarrow}{u_{\mathrm{s}}}\) is the generation of energy by the solid stress tensor, \(\gamma_{\Theta_{\mathrm{s}}}\) is the collisional dissipation of energy, and \(\phi_{\mathrm{ls}}\) is the energy exchange between the fluid and solid phases. In this study, the realizable \(k\)-\(\epsilon\) turbulence model was used as a turbulence closure model. The turbulence kinetic energy \(k\) and the dissipation rate \(\epsilon\) were obtained from the following transport equations [21; 22; 23]. \[\frac{\partial(\rho_{\mathrm{l}}k)}{\partial t}+\frac{\partial(\rho_{\mathrm{ l}}\overset{\rightarrow}{u_{j}})}{\partial x_{j}}=\frac{\partial}{\partial x_{j}} \Bigg{[}(\mu_{\mathrm{l}}+\frac{\mu_{\mathrm{t}}}{\sigma_{k}})\frac{\partial k }{\partial x_{j}}\Bigg{]}+G_{\mathrm{k}}-\rho_{\mathrm{l}}\epsilon \tag{17}\]  
\[\frac{\partial(\rho_{\mathrm{l}}\epsilon)}{\partial t}+\frac{\partial(\rho_{ \mathrm{l}}\overset{\rightarrow}{\epsilon}\overset{\rightarrow}{u_{j}})}{ \partial x_{j}}=\frac{\partial}{\partial x_{j}}\Bigg{[}(\mu_{\mathrm{l}}+\frac {\mu_{\mathrm{t}}}{\sigma_{\epsilon}})\frac{\partial\epsilon}{\partial x_{j}} \Bigg{]}+\rho_{\mathrm{l}}C_{1}S\epsilon-\rho_{\mathrm{l}}C_{2}\frac{\epsilon ^{2}}{k+\sqrt{\frac{\mu_{\mathrm{l}}}{\rho_{\mathrm{l}}}}\epsilon} \tag{18}\]  
where \(G_{\mathrm{k}}\) represents the generation of turbulence kinetic energy due to the mean velocity gradients and \(\mu_{\mathrm{t}}=\rho_{\mathrm{l}}C_{\mathrm{u}}k^{2}/\epsilon\) is the turbulent viscosity of the fluid. The related empirical constants are \(C_{1}=\max[0.43,\frac{\eta}{\eta+5}]\), \(\eta=S\frac{k}{\epsilon}\), \(\sigma_{\mathrm{k}}=1.0\), \(\sigma_{\mathrm{\varepsilon}}=1.2\) and \(C_{2}=1.9\).