**4. Numerical Procedure**

The mathematical model assumes the following boundary conditions:

$$\text{Let the pipe wall for } \mathbf{r} = \mathbf{R} \text{: } \mathbf{T} = \mathbf{T}\_{\text{W}} \text{: } \mathbf{q} = \text{const}; \text{ U = 0; } \mathbf{k} = \mathbf{0}; \text{ } \varepsilon = \mathbf{0}; \tag{35}$$

$$\text{symmetry axis for } \mathbf{r} = 0; \,\partial\Gamma/\partial\mathbf{r} = 0; \,\partial\mathbf{U}/\partial\mathbf{r} = 0, \,\partial\mathbf{k}/\partial\mathbf{r} = 0 \text{ and } \partial\varepsilon/\partial\mathbf{r} = 0. \tag{36}$$

The boundary condition (35) indicates that there is no sleep velocity, k and ε on the pipe wall (r = R). Very near the wall, the dissipation rate is equal to 2 μ(*∂*k0.5/*∂*r)2. This term was added to Equation (24) to allow to be set to zero on the pipe wall (y = 0). The boundary condition (36) indicates that axially symmetric conditions were applied to all dependent variables.

The set of partial differential Equations (8), (9), (24) and (25) was computed using a finite difference scheme and its own computer code. The set of equations was solved taking into account the TDMA approach, with an iteration procedure, and using the control volume method [74]. The control volume was obtained for a pipe of length equal to L=1m and rotating the radius of the pipe around the symmetry axis at an angle of 1 radian.

The calculations were carried out in a proper order. For the predetermined value of *∂p*/*∂x*, the following inlet conditions were used:


Attention was paid to the effect of the number of grid points localized across a radius of a pipe. It is well known that the number of nodal points strongly affects the accuracy of the predictions [74,75]. Therefore, the radius of the pipe was divided into 80 nodal points that were not uniformly distributed across the pipe. Most of the nodal points were located near the wall of the pipe. The number of nodal points was experimentally set to provide nodally independent predictions.

The iteration cycles were repeated until the convergence criterion, defined by Equation (37), was achieved.

$$\sum\_{j} \left| \frac{\mathcal{Q}\_j^n - \mathcal{Q}\_j^{n-1}}{\mathcal{Q}\_j^n} \right| \le 0.001\tag{37}$$

where ∅*<sup>n</sup> <sup>j</sup>* is a general dependent variable <sup>∅</sup> = U, T, k, <sup>ε</sup>; the jth is the nodal point after the nth iteration cycle and the ∅*n*−<sup>1</sup> *<sup>j</sup>* is the (n − 1)th iteration cycle.

Finally, the hydrodynamically and thermally developed turbulent flow of the fine dispersive slurry was calculated and the following profiles were obtained: T(r), U(r), k(r), ε(r) and the following dependent variables were calculated: Ub, Tb, Re, Pr, Nu.
