*3.1. Newtonian Liquid*

The mathematical model for the Newtonian liquid regards the carrier liquid, which is pure water. Taking into account the time-averaged procedure proposed by Reynolds [60], and neglecting terms with fluctuations of density and viscosity, it is possible to obtain the conservative equations of mass, momentum, and internal energy as follows:

$$\frac{\partial}{\partial \boldsymbol{x}\_{i}} \left( \overline{\rho} \overline{\boldsymbol{U}}\_{i} \right) = 0 \tag{1}$$

$$\frac{\partial}{\partial \mathbf{x}\_{\circ}} \left( \overline{\rho} \overline{\mathbf{U}\_{i} \mathbf{U}\_{\circ}} \right) = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_{i}} + \frac{\partial}{\partial \mathbf{x}\_{\circ}} \left( \mu \frac{\partial \overline{\mathbf{U}\_{i}}}{\partial \mathbf{x}\_{\circ}} - \overline{\rho} \, \overline{u\_{i}' u\_{j}'} \right) + \overline{\rho} g\_{i} \tag{2}$$

$$\frac{\partial}{\partial \mathbf{x}\_{\circ}} \left( \overline{\rho} \overline{\mathbf{U}}\_{\circ} \overline{\mathbf{T}} \right) = \frac{\partial}{\partial \mathbf{x}\_{\circ}} \left( \frac{\mu}{Pr} \frac{\partial \overline{T}}{\partial \mathbf{x}\_{\circ}} \right) - \frac{\partial}{\partial \mathbf{x}\_{\circ}} \left( \overline{\rho} \, \overline{u'\_{\circ} t'} \right) \tag{3}$$

where the Prandtl number in Equation (3) is defined as follows:

$$\text{Pr} = \frac{\mu \ c\_p}{\lambda} \tag{4}$$

The momentum Equation (2), also named the Reynolds Averaged Navier–Stokes equation, possesses an additional term *ρ u i u j* , which appears on the right-hand side, named turbulent stress tensor or Reynold's stress tensor. It is a symmetric tensor and requires additional equations to calculate its components. To solve the set of equations, the number of dependent variables must be the same as the number of equations, which is called the closure problem. Therefore, additional equations are proposed that allow one to calculate the components of the turbulent stress tensor. Equation (3) possesses an additional term *ρ u j t* . This term is responsible for the intensification of heat exchange caused by the presence of velocity and temperature fluctuations and requires closure, too. The additional terms in Equations (2) and (3) can be modelled by a direct approach (DNS) or an indirect approach. To calculate each component of −*ρ u i u <sup>j</sup>* and −*ρ u j t* the direct or indirect methods require additional equations, either algebraic or partial differential. Although DNS and LES methods can give a high level of accuracy, they are time-consuming and computationally costly. Indirect methods can give accurate predictions with less computational power. Moreover, we deal with time-averaged quantities of velocity and temperature, instead of random quantities, which are more useful to engineers. For this study, an indirect method is employed. The indirect method employs the Boussinesque hypothesis, which introduced a new dependent variable, called turbulent viscosity, as a measure of turbulence rather than viscosity [61]. The Boussineque hypothesis assumes that the turbulent stress tensor in Equation (2) can be expressed as follows:

$$-\overline{\rho}\,\overline{u\_i'u\_j'} = \mu\_l \left(\frac{\partial\overline{\mathcal{U}\_l}}{\partial x\_j} + \frac{\partial\overline{\mathcal{U}\_l}}{\partial x\_i}\right) - \frac{2}{3}\,\overline{\rho}\,k\,\delta\_{i,j}\tag{5}$$

and the fluctuating components of velocity and temperature in the energy Equation (3), as:

$$-\overline{\rho}\,\overline{u'\_jt'} = \frac{\mu\_l}{Pr\_l}\,\frac{\partial\overline{T}}{\partial x\_i} \tag{6}$$

The turbulent Prandtl number (*Prt*) in Equation (6) was intensively examined by several researchers. Studies indicated that in the flow on a plate, *Prt* ∼= 0.5, while for the boundary layer, which exists in a pipe flow, the turbulent Prandtl number is estimated to be *Prt* = 0.9 [62].

If an axially symmetric flow is considered, the most suitable coordinate system is cylindrical, where the dependent variables are functions of *x*, *r*, *θ*. The assumption that the flow is axially symmetric causes the velocity component V = 0. The lack of circumferential eddies causes that component of the velocity W = 0. The last term on the right-hand side of Equation (2) is also zero because the flow is horizontal.

Taking into account the assumptions stated in the physical model, the time-averaged Equations (1)–(3) together with (5), (6) can be formulated in the final form as follows:

$$\frac{\partial}{\partial \mathcal{X}} \left( \overline{\rho} \overline{U} \right) = 0 \tag{7}$$

$$\frac{1}{r}\frac{\partial}{\partial r}\left[r\left(\mu+\mu\_l\right)\frac{\partial \overline{U}}{\partial r}\right]=\frac{\partial \overline{p}}{\partial x}\tag{8}$$

$$\overline{\rho}\overline{\mathcal{U}}\,\frac{\partial\overline{\mathcal{T}}}{\partial\mathbf{x}} = \frac{1}{r}\,\frac{\partial}{\partial r}\left[r\left(\frac{\mu}{Pr} + \frac{\mu\_l}{Pr\_l}\right)\frac{\partial\overline{T}}{\partial r}\right] \tag{9}$$

Taking into account the assumption that the flow is fully developed thermally, which means that the temperature changes linearly in the main flow direction 'ox', the convective term in Equation (9) will be delivered by an energy balance. An energy balance is developed for the slurry flowing in a tube with inner diameter *D* and unit length *L* = 1 m. Such an approach simplifies the mathematical model to one-dimensional and substantially reduces the time of computation.

Let us consider a heat flux acting between the flowing liquid and a pipe—Figure 1.

$$
\dot{Q} = \alpha \, A\_\* \, \Delta \overline{T}\_r \tag{10}
$$

where

$$A\_\* = \pi DL \tag{11}$$

and Δ*Tr* = *Tb* − *Tw* is temperature difference on radius r.

The change in the inner energy of the flowing liquid at distance *L* is the following:

$$
\dot{Q} = \dot{m} \ c\_P \,\Delta \overline{T} = \overline{\rho}\_b \,\overline{\mathcal{U}}\_b \, A \,\ c\_p \,\Delta \overline{T}\_x \tag{12}
$$

where the pipe cross section *A* = *πD*2/4, and Δ*Tx* is the slurry temperature difference at distance *x=L*.

Comparing (10) and (12), we can write:

$$
\alpha \text{ } \pi \text{ D } L \, \Delta \overline{T}\_r = \overline{\rho}\_b \, \overline{\mathcal{U}}\_b \, \text{ } \pi \, \frac{D^2}{4} \, \text{ } c\_p \, \Delta \overline{T}\_x \, \text{ } \tag{13}
$$

Dividing each side of Equation (13) by the unit length of the pipe *L*, we can write:

$$\dot{q} = \overline{\rho}\_b \,\overline{\mathcal{U}}\_b \,\,\pi \,\frac{D^2}{4} \,\, c\_p \,\, \frac{\Delta \overline{T}\_x}{L} \tag{14}$$

where the heat flux density can be written as:

$$
\dot{q} = \frac{\dot{Q}}{L} = \text{a } \pi \text{ D } \Delta \overline{T}\_r \tag{15}
$$

The final form of the convective term can be established using Equation (14), as follows:

$$\frac{\partial \overline{T}}{\partial \boldsymbol{\alpha}} = \frac{\dot{\boldsymbol{q}}}{\overline{\rho}\_b \, \overline{\mathbf{U}}\_b \, \boldsymbol{\pi} \, \boldsymbol{R}^2 \boldsymbol{c}\_p} = const \tag{16}$$

The convective term (16) delivered from the energy balance replaces the convective term on the left-hand side of Equation (9). Taking into account Equation (16), one can say that the convective term *∂*T/*∂*x = 0 if q = 0. The heat flux acting on the pipeline is positive if Tw > Tb or negative if Tw < Tb. The thermal properties of the slurry are treated as known quantities and will be presented. The influence of temperature on a carrier liquid density and a thermal conductivity was approximated using the Taylor expansion as follows:

$$
\Phi = a + b\overline{T} + c\overline{T}^2 \tag{17}
$$

while slurry properties, such as density, coefficient of heat conduction, and specific heat, were calculated as follows:

$$
\Phi\_{SL} = \mathbb{C}\Phi\_S + (1 - \mathbb{C})\Phi\_L \tag{18}
$$

Finally, the mathematical model of heat transfer in liquid flow constitutes two partial differential equations, namely, (8) and (9), since the continuity Equation (7) is included in the momentum Equation (8), and the complementary Equation (16). Equation sets possess the following dependent variables: U(r), μt(r), p(x), T(r). To calculate the velocity and temperature profiles across a pipe, a set of equations requires closure. Turbulent viscosity μt(r) will be calculated using an indirect method together with a suitable turbulence model, while the pressure gradient *∂*p/*∂*x will be treated as a known value. For the preset value of *∂*p/*∂*x, the velocity and temperature profiles will be calculated. In other words, for known values of *∂*p/*∂*x the time-averaged profiles of U and T will be calculated.

To calculate the turbulent viscosity μt(r), the Launder and Sharma turbulence model was chosen [63]. The chosen model fulfills the requirements formulated by Spalding [64,65], an honorable founder of Computational Fluid Dynamics (CFD). The Launder and Sharma model was successfully examined for comprehensive types of flow, including homogeneous slurries. The researchers emphasized that it is one of the first and most widely used models and has been shown to agree well with the experimental and DNS data for a wide range of turbulent flow problems, performing better than many other k–ε models [66–69].

$$k = \frac{\overline{u\_i' u\_i'}}{2} \tag{19}$$

$$
\varepsilon = \nu \frac{\overline{\partial u\_i'} \, \overline{\partial u\_i'}}{\partial \boldsymbol{\alpha}\_k} \tag{20}
$$

The Launder and Sharma turbulence model assumes that the turbulent viscosity depends on the kinetic energy of the turbulence, defined by Equation (19), and the dissipation rate of the kinetic energy of the turbulence, defined for homogeneous turbulence by Equation (20) [70]. On the basis of dimensional analyzes, Launder and Sharma assumed that the turbulent viscosity depends on *k*, *ε* and *ρ* as follows:

$$
\mu\_t = f\_\mu \frac{\overline{\rho}}{\varepsilon} k^2 \tag{21}
$$

where *fμ* is called the damping of the turbulence function or the wall function. The function causes a reduction in turbulent viscosity if y → 0. The wall function was developed empirically by matching the predictions with measurements and is the following:

$$f\_{\mu} = 0.09 \exp\left[\frac{-3.4}{\left(1 + \frac{Rc\_l}{50}\right)^2}\right] \tag{22}$$

where the turbulent Reynolds number (*Ret*) was developed from dimensionless analysis, as follows:

$$Re\_l = \frac{\rho \, k^2}{\mu \, \varepsilon} \tag{23}$$

Equations for the kinetic energy of turbulence and its dissipation rate are obtained from the Navier–Stokes equations, using a time-averaged procedure [63].

Taking into account the assumptions made in the physical model, the final form of the equations for the kinetic energy of turbulence and its dissipation rate, in cylindrical coordinates, are as follows:

$$\frac{1}{r}\frac{\partial}{\partial r}\left[r\left(\mu + \frac{\mu\_t}{\sigma\_k}\right)\frac{\partial k}{\partial r}\right] + \mu\_l \left(\frac{\partial \overline{U}}{\partial r}\right)^2 = \rho \varepsilon + 2\mu \left(\frac{\partial k^{1/2}}{\partial r}\right)^2\tag{24}$$

$$\frac{1}{r}\frac{\partial}{\partial r}\left[r\left(\mu+\frac{\mu\_l}{\sigma\_\varepsilon}\right)\frac{\partial \varepsilon}{\partial r}\right]+\mathbb{C}\_1\frac{\varepsilon}{k}\mu\_l\left(\frac{\partial \overline{U}}{\partial r}\right)^2=\mathbb{C}\_2\left[1-0.3\exp\left(-\mathrm{Re}\_l^2\right)\right]\frac{\rho\varepsilon^2}{k}-2\frac{\mu}{\rho}\mu\_l\left(\frac{\partial^2 \overline{U}}{\partial r^2}\right)^2\tag{25}$$

The set of partial differential Equations (8), (9), (24) and (25) together with the complementary Equations (16), (21)–(23) is dedicated to a Newtonian liquid.
