## 2 Numerical Methods  
The Eulerian multiphase model is use for modelling the flow behaviour for each of the liquid and solid phases. This model is the most complex model, as it solves a set of momentum equations for each phase and is the most time-consuming model. For liquid-solid interaction we took into consideration only the drag force according to the Huilin-Gidaspow model and the air-liquid interaction employs the Schiller-Naumann model. We considered one dispersed phase, the solid particles, one continuous phase, the liquid, and above free surface we consider a volume of air. We did not apply surface compression and we did not include the effect of surface tension. For simulating the turbulence effects, the RNG k-\(\epsilon\) model with the dispersed option has been chosen. The dispersed turbulence model is the appropriate model when using the granular model. The geometry of the industrial stirred reactor fitted with the two new designed impellers in tandem and the mesh was generated with the help of the software Gambit 2.4. The mesh of finite elements consists of approximately 0.94 million tetrahedral elements. As regards the mesh quality, we use hybrid tetrahedral/hexahedral mesh cells as recommended for such complex geometries [2; 3]. In this study a good quality of mesh (skewness < 0.84) has been ensured throughout the computational domain. Although the mesh influence upon numerical results was not addressed and a larger mesh element number may generate more accurate results, the size of the mesh employed in this study may be assumed to be appropriate to compare the performances of the single and tandem impellers. As in our previous research papers [21; 22], we used the Sliding Mesh approach, and four fluid zones were defined: an inner rotating cylindrical volume centered on the two impellers and three other volumes containing the rest of the reactor and the baffles. The interface is located at equal distance from the blades of the two impellers and the inside edge of the baffles. An illustration of the domain and the mesh of the domain is given in Figure 2, which shows that finer meshes have been used around the two impellers where the velocity spatial gradients are expected to be large. Good care was taken in keeping the same mesh dimensions between the rotating and the stationary interface to ensure a good exchange of the hydrodynamic quantities between the two zones during calculation. A \"no-slip\" condition is applied on all walls of the geometry. Simulations of the turbulent multi-phase flow are then performed with the help of the commercial CFD code Ansys(r) Fluent, Release 16 using the RNG k-\(\epsilon\) model. There are three working materials corresponding to the three phases: liquid with the density \(\rho_{l}=690\) kg\(\cdot\)m\({}^{-3}\) and viscosity of 0.00023 kg/m\(\cdot\)s, powder with the density \(\rho_{p}=1520\) kg\(\cdot\)m\({}^{-3}\) and the particle size \(d=171\) um, and air with the density \(\rho_{a}=1.225\) kg\(\cdot\)m\({}^{-3}\) and viscosity of 1.7894 \(\times\) 10\({}^{-5}\) kg/m\(\cdot\)s. At the liquid surface, a gas zone consisting of air was added at the free surface of water, a method that has been reported to dampen the instabilities which means top surface being exposed to atmospheric pressure. Using the sliding mesh approach, one has to solve for an unsteady flow. At each time step, the position of the rotating zone relative to the stationary one is recomputed and the grid interface of the rotating zone slides along the interface of the stationary zone. For all simulations, the time step is set as a function of the two impellers rotational speed, 65 rot/min, and the number of the cells from the interface surface, so that the time step is set to a value of 5 milliseconds. At the end of each time step, after a maximum of 5 iterations, the convergence criterion reaches \(10^{-3}\) for the continuity, for momentum and turbulence quantities. Simulations are performed in double precision with the segregated implicit solver.