## 3 Numerical Model and Method
### Numerical Model
\[f_{n,ij}=-\kappa_{n}\delta_{n,ij}-\sigma_{n}v_{n,ij} \tag{7}\]  
and  
\[f_{t,ij}=\kappa_{t}\delta_{t,ij}-\sigma_{t}v_{t,ij} \tag{8}\]  
where \(\delta_{n,ij}\) and \(\delta_{t,ij}\) express the normal and tangential interpenetration of the two particles \(i\) and \(j\), \(\kappa_{n}\),\(\kappa_{t}\) are respectively the stiffness at normal and tangential contact, \(\sigma_{n}\), \(\sigma_{t}\) are respectively the normal and tangential viscous damping coefficient and where normal and tangential relative velocities \(v_{n,ij}\) and \(v_{t,ij}\) are given by the relations:  
\[m_{p}\frac{du_{p}}{dt}=F_{P}+F_{D}+F_{C} \tag{1}\]0  
\[m_{p}\frac{du_{p}}{dt}=F_{P}+F_{D}+F_{C} \tag{1}\]1  
The friction slider limits the tangential force governed by Coulomb's friction law [24]. \[m_{p}\frac{du_{p}}{dt}=F_{P}+F_{D}+F_{C} \tag{1}\]2  
Water phase is modeled in CFD-FLUENT as a continuous fluid phase. The equations of continuity and momentum conservation are solved for each time step related to the fluid phase:  
The continuum equation for the \(q\) phase is:  
\[m_{p}\frac{du_{p}}{dt}=F_{P}+F_{D}+F_{C} \tag{1}\]3  
The momentum balance for the \(q\) phase is:  
\[m_{p}\frac{du_{p}}{dt}=F_{P}+F_{D}+F_{C} \tag{1}\]4  
In which expressions, \(\rho_{q}\) is the \(q\) phase density, \(g\) is the acceleration of gravity, \(\overset{\rightarrow}{u}_{p}\) and \(\overset{\rightarrow}{u}_{q}\) are the phase velocities, \(\overline{\tau}_{q}\) is the \(q\) phase stress-strain tensor given by:  
\[m_{p}\frac{du_{p}}{dt}=F_{P}+F_{D}+F_{C} \tag{1}\]5  
Here, \(\mu_{q}\) and \(\lambda_{q}\) are the \(q\) phase's shear and bulk viscosity. \(k_{sl}\) is the interphase momentum exchange coefficient. For granular flows, momentum exchange between phases is based on the fluid-solid exchange coefficients value. In CFD-ANSYS Fluent, the interphase exchange coefficient models are empirically based. The Gidaspow model, which is a combination of the Wen and Yu model [25] and the Ergun equation [26], is used in this study. This model is recommended for dense fluidized beds. \[k_{sl}=\left\{ \begin{array}{cc}\frac{3}{4}C_{D}\frac{\varepsilon_{s} \varepsilon_{l}\rho_{l}\left|\stackrel{{\rightarrow}}{{u}}_{s} \stackrel{{\rightarrow}}{{u}}_{l}\right|}{d_{s}}\varepsilon_{l}^{- 265},&\varepsilon_{l}>0.8\\ 150\frac{\varepsilon_{s}(1-\varepsilon_{l})\pi_{l}}{\varepsilon_{l}d_{s}^{2} }+1.75\frac{\rho_{l}\varepsilon_{l}\left|\stackrel{{ \rightarrow}}{{u}}_{s}\stackrel{{\rightarrow}}{{u}}_{l}\right|} {d_{s}}\end{array} ; \begin{array}{cc}\varepsilon_{l}>0.8\\ \varepsilon_{l}\leq 0.8\end{array} \right\} \tag{15}\]  
where \(\mu_{l}\) is the water dynamic viscosity, \(d_{s}\) is the particle diameter and \(C_{D}\) is drag force coefficient of a single particle. \(C_{D}\) is expressed by:  
\[C_{D}=\frac{24}{\varepsilon_{l}R_{\varepsilon_{s}}}\Big{[}1+0.15(\varepsilon_ {l}R_{\varepsilon_{s}})^{0.687}\Big{]} \tag{16}\]  
With the Reynolds number as:  
\[R_{\varepsilon_{s}}=\frac{\varepsilon_{l}\rho_{l}d_{s}\Big{|}\stackrel{{ \rightarrow}}{{u}}_{s}-\stackrel{{\rightarrow}}{{u}}_{l} \Big{|}}{\mu_{l}} \tag{17}\]  
More detailed information on the Gidaspow model can be found in Reference [27].