Herschel-Bulkley ﬂuid ﬂow within a pipe by taking into account viscous dissipation

– The present numerical study deals with the analysis of thermal characteristics during the ﬂow of an incompressible Herschel-Bulkley ﬂuid of constant physical and rheological properties. The ﬂow takes place within a pipe, of circular cross section, maintained at uniform wall temperature. Because of the viscous character of this kind of ﬂuid, viscous dissipation is taken into account. The eﬀect of this parameter is analyzed for various values of the ﬂuid index ﬂow and the Herschel-Bulkley number. The governing equations are resolved by means of the ﬁnite volume method and using the SIMPLER algorithm and the line by line method. The results reveal that heat transfer is signiﬁcantly underestimated when viscous dissipation is neglected, particularly for shear-thinning ﬂuids and for high values of the Herschel-Bulkley number. In addition, heat transfer is more important in the case of heating and the curves display a discontinuity, and negative values of Nusselt number are obtained.


Introduction
The general model of Herschel-Bulkley describes a large variety of natural and industrial fluids such as lavas, avalanche, ketchup [1], yoghurt [2], Carbopol 940 [3], waxy crude oil [4] and cement grouts [5].This kind of fluids, also called viscoplastic fluids, is characterized by a yield stress from which they start moving.If the yield stress is not reached, the fluid behaves like a rigid solid.
Many authors [6][7][8][9][10][11][12][13] studied the flow of viscoplastic fluids, especially those obeying the rheological model of Bingham while others were interested in the flow of Herschel-Bulkley fluids.Indeed, Soares et al. [14] studied, by means of finite volume method, heat transfer in the inlet region of an axial laminar flow of Herschel-Bulkley materials by taking into account axial diffusion.Their results show that, when the pipe is kept at uniform wall temperature or heat flux, heat transfer is enhanced when the Herschel-Bulkley number increases and the flow index decreases.Nouar [15] studied free and forced convection heat transfer of a Herschel-Bulkley fluid flow in a horizontal circular pipe subjected to a constant wall heat flux.He analyzed the effect of thermodependency of the fluid consistency and its density on the secondary flows and gave correlations for the Nusselt number and the wall a Corresponding author: nabilalabsi@yahoo.frshear stress.Mitsoulis and Hatzikiriakos [16] undertook an experimental study on the rolling of bread dough.The latter obeys the shear thinning Herschel-Bulkley rheological model.The authors observed that the unyielded zones are located at the outlet of the machine.
Because of the highly viscous nature of viscoplastic fluids, some researchers have introduced the dissipation function in their works in order to analyze and to understand this fluids behavior under these conditions.Vradis et al. [17] investigated numerically the simultaneous development of hydrodynamic and thermal fields of a Bingham fluid's flow in the entrance region of a circular pipe with and without viscous dissipation.Min et al. [18] studied, by means of a numerical method, hydrodynamic and thermal developing and simultaneous developing of laminar flows of a Bingham plastic in a circular pipe.The effect of viscous dissipation on the laminar forced convection flow of a Bingham fluid within a pipe has also been analyzed by Khatyr et al. [19,20].The authors considered various wall thermal conditions.Osalusi et al. [21] studied numerically, the combined effect of viscous dissipation and Joule heating on an electrically conducting Bingham fluid flow over a porous rotating disk.Boualit et al. [22]  fluid between two plane plates maintained at a constant temperature.
From the literature inspection, we note that viscous dissipation's effect for Herschel-Bulkley fluids has not been widely studied.Indeed most of studies concerned the limit case of Bingham except few of them concern the Herschel-Bulkley fluid flow, such as the interesting one released by Sayed-Ahmed [23] who used an implicit Crank-Nicolson method to obtain the temperature distributions in solving the energy equation when viscous dissipation is taken into account.The author considered the laminar heat transfer of a Herschel-Bulkley fluid in the entrance region of a square duct assuming fully developed velocity profile.He found that the Nusselt number decreases with the increase of the flow index and increases with the increase of both yield stress and Brinkman number.
The present study is purely numerical and represents an extension of previous works we carried out, adding some new results.Our contribution focuses on the effect of neglecting viscous dissipation when dealing with a Herschel-Bulkley fluid flow within a circular pipe maintained at uniform wall temperature.This investigation was motivated by the fact that only few papers deal with the Herschel-Bulkley fluid flow within a pipe, especially the shear thickening one (n > 1) while more papers investigate the limit case of a Bingham fluid such that of Min et al. [18].Furthermore, viscous dissipation has not yet been deeply investigated for shear thinning (n < 1) and shear thickening Herschel-Bulkley (n > 1) fluids, especially when considering a wall heating (Br < 0).So, we consider, for this purpose, a wide range of flow index ranging from 0.4 to 2.0 in order to draw abacuses illustrating heat transfer performances within the pipe.

Mathematical formulation
We consider in the present paper the laminar steady flow of an incompressible Herschel-Bulkley fluid through a circular pipe of length L and diameter D maintained at constant wall temperature T w (Fig. 1).The physical and rheological properties of the fluid are constant and uniform.On the other hand, the dissipation function is taken into account in the energy equation.
Flow and forced heat transfer convection are governed by the following dimensionless equations, i.e. continuity, momentum and energy equations, which are respectively given, for a fluid of constant physical properties (ρ, C p and k), by: where: The boundary conditions consist of uniform axial velocity and temperature at the inlet (U = θ = 1, V = 0), no-slip condition and a uniform wall temperature along the wall (U = V = θ = 0) as well as fully developed velocity and temperature at the outlet ∂U ∂X = ∂V ∂X = ∂θ ∂X = 0 .Equations ( 1)-( 3) are written by considering the following dimensionless variables: axial and radial coordinates, axial and radial velocities, pressure and temperature: where V i represents the radial velocity components (i = x or r which are respectively, the axial and the radial coordinates), V 0 is the inlet velocity, p * is the pressure and finally, T 0 and T w are respectively the inlet and the wall temperatures.
The above dimensionless equations highlight the following dimensionless numbers: . The latter compares the dissipation term with the conduction term in the energy equation.It represents the viscous dissipation function.A negative value of the Brinkman number means that the fluid is heated (heating case) whereas a positive value indicates that the fluid is cooled down (cooling case).
η app , appearing in the energy Equation (3), is the dimensionless effective viscosity of the Herschel-Bulkley fluid.The general rheological model of the latter, relating the shear stress τ to the shear rate γ, is given by the following law: where K is the fluid's consistency, n is the flow index and τ 0 is the yield stress.Viscoplastic fluids in general, including those obeying the rheological model of Herschel-Bulkley, present two regions: a yielded region located near the wall and an unyielded region called "the plug flow region".This discontinuity causes numerical instabilities in the low shear rate region.In order to avoid these numerical instabilities, Equation ( 4) is replaced by the following constitutive law that's proposed by Papanastasiou [24,25]: where m is the Papanastasiou regularization parameter.
In the last decades, many researchers [26][27][28] used the Papanastasiou regularization model in their numerical and experimental investigations.Typically, a value of m = 1000 s is considered [18,29].
η is the effective viscosity of the fluid.For the Herschel-Bulkley fluid, its dimensionless form, η app , is given by: where γ * is the dimensionless shear rate, M = mV 0 /H represents the dimensionless exponential growth parameter and HB = τ 0 D n /K V n 0 is the Herschel-Bulkley number which represents the ratio of the yield stress to the nominal shear stress.
Equations ( 1)-( 3) form an elliptic system of partial differential equations.Numerical methods are needed to obtain solutions.The set of equations with the corresponding boundary conditions are, therefore, solved numerically using the finite volume method proposed by Patankar [30].This method is rapid and inexpensive.This widely used method ensures a complete conservation of mass, momentum and energy in the computational domain.Furthermore, it allows the simulation of both actual and ideal conditions.The finite volume method is based on subdividing the calculation domain in a finite number of control volumes, also called grid cells.The basic principle of this method is to integrate the differential governing equations over each cell (control volume) surrounding the calculating point, called the principal node, within a staggered grid.The scalar variables (pressure and temperature) are located at the centre of each cell while the velocity components are defined at the cell faces.Thus, the differential governing equations are transformed in algebraic equations involving the unknown values (velocities, pressure and temperature) at chosen grid points.
These equations are solved, in the present study, using a home made computer code, based on the line by line solver method.The grid adopted in the present paper is non-uniform with a mesh refinement at the inlet and near the wall of the pipe.The grid adopted in our paper consists of 250 nodes in the X direction and 50 nodes in the R direction (250 × 50).This meshing is optimal whatever the values of the Herschel-Bulkley number and the flow index and was obtained after an analysis of the sensitivity of the results to the meshing, for both hydrodynamic and thermal properties, in order to assure the accuracy of our simulations.This analysis concerned the cases where viscous dissipation was taken into account (Br = 0) and where it was neglected (Br = 0) and that, for different values of the flow index n and by considering different values of the Herschel-Bulkley number HB.The analysis of the sensitivity of the results to the meshing was supported by the comparison of the results with previous works (the validation of the computer code).However, this was possible only for the case of a Herschel-Bulkley fluid flow by neglecting viscous dissipation and the case of the limit case of a Bingham fluid (a Herschel-Bulkley fluid with n = 1) when viscous dissipation is taken into account.
The convergence criterion, which is based on the residual, is set to 10 −5 for both velocity components and temperature and to 10 −6 for pressure and 10 4 iterations were necessary to reach the solution.

Results and discussion
In the present paper, the Reynolds and the Prandtl numbers are taken equal to 20 and 50, respectively.These values give a moderate value 1000 for the Peclet number since most fluids encountered in industries have a high Prandtl number [31] and also since P e > 100 allows to neglect axial diffusion [32].We note that all the curves, except the ones corresponding to HB = 0 (shear thinning power-law fluid, Newtonian fluid and shear thickening power-law fluid for n = 0.5, n = 1 and n = 1.5, respectively), present two regions: a region close to the wall characterized by a parabolic velocity profile called the "yielded region" and a zone around the centreline which represents a uniform velocity distribution and is called the "core region" (unyielded region).In this region, the shear stress is less than the yield stress.Consequently, the fluid resists to deformations and moves like a rigid solid.

Validation of the computer code
Figures 2a-2c compare also our numerical results with those obtained by the following analytical solution for the fully developed velocity profile [33], after adaptation to our dimensionless parameters: r 0 is the dimensionless plug size which corresponds to our dimensionless writing: r 0 = r * 0 /2.r * 0 is the dimensionless plug size in the paper of Peixinho et al. [33].
Here, the dimensionless plug size r * 0 is calculated from the following equation proposed by Peixinho et al. [33] but adapted to our dimensionless parameters: The figures show that the numerical and analytical solutions are in good agreement.Indeed, for the case of a  shear thinning Herschel-Bulkley fluid (n = 0.5), the difference between the two results is about 1.75%, 2.60% and 3.89% for HB = 2, 5 and 10, respectively.For the case of a Bingham fluid (n = 1), it does not exceed 1%, 2.06% and 1.05% for HB = 2, 5 and 10, respectively and finally, for the case of a shear thickening Herschel-Bulkley fluid (n = 1.5), it is less than 0.19%, 0.12% and 0.23% for HB = 2, 5 and 10, respectively.
For HB = 0, the fluid can be described by the powerlaw model to which corresponds the fully developed velocity profile given in Equation ( 9) [34]: The comparison between our numerical results and those calculated from Equation ( 9) for a shear thinning fluid (Fig. 2a), a Newtonian fluid (Fig. 2b) and a shear thickening fluid (Fig. 2c) shows a good agreement.The difference between both results is about 2.9%, 1.14% and 2.43% for n = 0.5, n = 1 and 1.5, respectively.It is to notify that the plug size (i.e. the extend of the plug flow region in the half of the pipe) was calculated using Equation ( 8) by an iterative procedure in Excel.The extend of the plug flow region is assembled in Table 1.
We can notice from Table 1 that the increase of the Herschel-Bulkley number leads to the increase of the extend of the plug flow region (unyielded region) for which the values are more significant for a shear thinning Herschel-Bulkley fluid.This can be explained by the fact that the wall velocity gradient (Figs.2a-2c) is more important for weak values of the flow index.However, large wall velocity gradients are characterized by a decrease of the extend of the yield region, which offers to the plug flow the opportunity to stretch out in a large distance of the pipe.
The Table 1 shows also that the great values of the extend of the unyielded region are obtained in the case of a shear thickening Herschel-Bulkley fluid.

Effect of viscous dissipation on heat transfer
Viscous dissipation conveys the deterioration of a fluid mechanical energy in heat.This phenomenon is introduced in the heat balance by adding a source term to the energy equation (the dissipation function).Its intensity depends on both the local velocity gradient and the fluid's effective viscosity [35].In addition, it converses, in general, the contribution of all viscous stresses: normal and shear stresses.However, in most of cases (for incompressible fluids in particular), the contribution of shear stresses during the fluid's deformation is dominant.
Viscous dissipation is quantified by the Brinkman dimensionless number, which represents the ratio between the heat generated by viscosity and that dissipated by conduction.Hence, viscous dissipation effects are more significant when the Brinkman number exceeds the unity.Therefore, taking into account this function in the energy equation leads to significant modifications in the thermal properties of the fluid flow.
The dissipation function is often neglected in studies.It is actually valid for non viscous Newtonian fluids such as air, and for some non-Newtonian fluids but not for viscoplastic fluids which are very viscous.This has been confirmed by many authors [17,18,36] who considered the flow of a viscoplastic fluid in a circular pipe, by taking into account viscous dissipation.Their studies reveal also that viscous dissipation affects only heat transfer characteristics.This is predictable since the momentum and energy equations are decoupled.For this reason, the analysis in the present paper concerns only the thermal structure of the flow, characterized mainly by the local Nusselt number variations given by the following expression: Figures 3a-3 c show the effect of the Brinkman number on the axial evolution of Nusselt number for a shear thinning Herschel-Bulkley fluid (n = 0.5), a Bingham fluid (n = 1) and a shear thickening Herschel-Bulkley fluid (n = 1.5), respectively.The results concern both heating (Br < 0) and cooling (Br > 0) cases.For the limit case of a Bingham fluid (Fig. 4b), the results are compared to those obtained from the numerical study of Min et al. [18].The agreement is good whether or not viscous dissipation is neglected or taken into account, since the deviation between the two studies does not exceed 3%.The curves show a decrease of the Nusselt number near the inlet, before reaching an asymptotic value which corresponds to a thermal fully developed flow.It is interesting to note that for both heating and cooling (Br = 0), this asymptotic value is independent of Brinkman number.Indeed, it is equal to 16.44, 10.66 and 9.40 for n = 0.5, n = 1 and n = 1.5, respectively.These values are more important than the ones obtained when viscous dissipation is neglected (Br = 0).This means that heat transfer is underestimated when viscous dissipation is neglected.Actually, for viscous fluids, the fluid's internal heating, caused by viscous dissipation, leads to the increase of the temperature within the fluid, which exceeds in this case the unity unlike to the case of negligible viscous dissipation (Figs.4a-4c).
It is to be noted that heat transfer is more significant for the heating case (Br < 0).The Nusselt number increases consequently when the Brinkman number increases.As a result, neglecting viscous dissipation leads to underestimate heat transfer by about 256%, 180% and 164% for n = 0.5, n = 1 and n = 1.5, respectively.It is interesting to note, that in the case of heating (Br < 0), the curves present a discontinuity and negative values of the Nusselt number are obtained.This represents a change in heat direction such as reported by Jambal et al. [34] for the case of a power-law fluid.
In order to explain, more in detail, the discontinuities as well as negative values of the Nusselt number, we consider the case of a shear thinning fluid (n = 0.5).This explanation is based on the expression of the Nusselt number (Eq. ( 10)).
Thus, negative values of the Nusselt number are due, as mentioned before, to the change in heat direction.Indeed, for Br = −10, for instance, the Nusselt number decreases from the inlet until the value 0 at a position X = 0.05, since the dimensionless temperature gradient becomes equal to 0 (Fig. 5) which means that the fluid's dimensional temperature becomes equal to the wall temperature (T = T w ).
The curve continues its decrease between the axial positions X = 0.05 and X = 8.42 (Fig. 3a) with negative values of the Nusselt number, because of negative values of the dimensionless bulk temperature, i.e. the fluid's dimensional bulk temperature increases asymptotically along the pipe but remains less than the wall temperature (T m < T w ).At the axial position X = 8.42, the dimensionless bulk temperature is equal to 0 (θ m = 0) according to Figure 4a, it means that the bulk temperature is equal to the wall temperature.This corresponds, according to Equation (10), to an infinite value of the Nusselt number (Fig. 4a) which will be followed by a brusque change in the sign since from this position, the dimensionless bulk temperature becomes negative (i.e.T m > T w ) and the temperature gradient is still positive (Fig. 5).Therefore, heat transfer will occur in the opposite direction, i.e. the fluid will heat the wall and not the opposite.
It is to be noted that minimum and maximum values of the Nusselt number are observed at the axial positions X = 3.80 and X = 320, respectively.It is worth noting that for a shear thickening fluid of a flow index equal to 1.6, the discontinuity of the Nusselt number curves still exists whatever the value of the Herschel-Bulkley number.However, it still exists for n = 1.8only for Br = −2 whatever the value of the Herschel-Bulkley number but it disappears completely for n = 2 and the behavior becomes similar to that obtained for the cooling case (Br > 0).This can be explained by the fact that the fluid tends to have a solid behavior.

Abacuses
In order to make our results useful, we propose abacuses which give an idea on how heat transfer is undervalued when viscous dissipation is neglected.These abacuses, shown in Figures 9a-9c for HB = 2, HB = 5 and HB = 10, respectively, display the variation of the asymptotic values of the Nusselt number according to both the Herschel-Bulkley number and the flow index.The figures display a decrease revealing that heat transfer is more significant for low values of the flow index whether viscous dissipation is neglected or taken into account.
Figure 10 shows, as for it, the variations of the ratio between the asymptotic values of the Nusselt number when viscous dissipation is taken into account and those when viscous dissipation is neglected.We notice that the increase of the Herschel-Bulkley number leads to the increase of heat transfer for low values of the flow index.This increase becomes less and less important when the flow index increases.For n = 0.6, for instance, heat transfer when taking into account viscous dissipation is 3.40, 4.10 and 5.00 times more important than that obtained when viscous dissipation is neglected for HB = 2, 5 and 10, respectively.However, it does not exceed 2.50 for n = 2 which corresponds to a shear thickening Herschel-Bulkley fluid.
It is to be noticed that only asymptotic values of the Nusselt number were represented instead of mean values, since the representation of the latter may cause errors because of the curves discontinuities and negative values of the Nusselt number.

Conclusion
The present numerical study, based on finite volume method, enables us to analyse the impact of neglecting  The calculation and the fully developed velocity profiles show that the extend of the core region increases when both the Herschel-Bulkley number and the flow index increase.Furthermore, the results reveal that for this kind of fluids, it is required to take into account viscous dissipation in calculation since when it is neglected, significant undervaluation of heat transfer occurs.Indeed, the asymptotic value of the Nusselt number is notably greater than the one corresponding to the case of a neglected viscous dissipation, particularly for a Herschel-Bulkley fluid of a weak index flow (shear thinning) and high Herschel-Bulkley number.It has been noted also that in the case of heating (Br < 0), heat transfer is more important and the curves display a discontinuity, and negative values of Nusselt number are obtained.

Figures
Figures 2a-2c show the fully developed velocity profiles for various values of the Herschel-Bulkley number when viscous dissipation is neglected (Br = 0), for a shear thinning Herschel-Bulkley fluid (n = 0.5), a Bingham fluid (n = 1) and a shear thickening Herschel-Bulkley fluid (n = 1.5), respectively.We note that all the curves, except the ones corresponding to HB = 0 (shear thinning power-law fluid, Newtonian fluid and shear thickening power-law fluid for n = 0.5, n = 1 and n = 1.5, respectively), present two regions: a region close to the wall characterized by a parabolic velocity profile called the "yielded region" and a zone around the centreline which represents a uniform velocity distribution and is called the "core region" (unyielded region).In this region, the shear stress is less than the yield stress.Consequently, the fluid resists to deformations and moves like a rigid solid.Figures 2a-2c compare also our numerical results with those obtained by the following analytical solution for the fully developed velocity profile[33], after adaptation to our dimensionless parameters:

Table 1 .
Extend of the plug flow region to both Herschel-Bulkley number and index flow.P r = 50; Re = 20; Br = 0. Extend of the plug flow region r 0 HB n = 0.5 n = 1

Fig. 5 .
Fig. 5. Axial evolution of the dimensionless wall temperature gradient for different values of the Brinkman number.Re = 20, P r = 50, HB = 2.

Figures
Figures 3a-3c reveal that when the flow index increases, the discontinuity observed for the heating case (Br < 0) begins to disappear.This result led us to deal with our analysis in deep and extend the study to other values of both flow index and Herschel-Bulkley number.Thus, we consider the case of a shear thickening Herschel-Bulkley fluid (n > 1) and three values of the Herschel-Bulkley number: HB = 2, HB = 5 and HB = 10.The results are shown in Figures 6, 7 and 8 for n = 1.6, n = 1.8 and n = 2, respectively.It is worth noting that for a shear thickening fluid of a flow index equal to 1.6, the discontinuity of the Nusselt number curves still exists whatever the value of the Herschel-Bulkley number.However, it still exists for n = 1.8only for Br = −2 whatever the value of the Herschel-Bulkley number but it disappears completely for n = 2 and the behavior becomes similar to that obtained for the cooling case (Br > 0).This can be explained by the fact that the fluid tends to have a solid behavior.