Numerical study of hydrodynamic and thermal behaviors of agitated yield-stress ﬂuid

– Heat transfer during agitation of Bingham viscoplastic ﬂuid is studied in this paper. The ﬂuid is agitated with an anchor impeller and the heating is made by a jacketed wall of the stirred vessel. Transfers in the agitated vessel, translating hydrodynamic and thermal phenomena, are numerically predicted by means of Computational Fluid Dynamics (CFD) in transient regime. The purpose of this numerical study is to identify the rigid zones and to optimize mixing and heating performances. The Navier-Stokes and energy equations are discretized using ﬁnite volume method, and a two-dimensional analysis of the hydrodynamic and transient thermal behaviours generated in the agitated vessel are performed. Fluid rheology is modeled by the Bingham approximation and Papanastasiou’s regularization model. Results show the presence of recirculation zones and permit to explain the unpredicted Nusselt number increasing when Oldroyd number increases. This study shows also the importance of the anchor position on the size and the shape of the rigid zones and on the heating performances.


Introduction
In many industrial areas mixing is an important step such as the food, pharmaceutical or polymer industries. Mixing operation has to be controlled in order to ameliorate fina1 product properties. Non-Newtonian fluids possibly with a yield stress pose additional problems in agitation. The rheological behavior of these fluids governs the hydrodynamics found in these processes. Other difficulties for the optimization of these processes often occur with the heating of such fluids.
In the aim of controlling the mixing of viscoplastic fluids by a Rushton turbine, many numerical works were made and published: Torrez and Andre [1], Arratia et al. [2], Ein-Mozaffari and Upreti [3]. . . These works show the existence of well-mixed caverns around the impeller with near stagnant regions outside. By assuming different cavern shapes (spherical, cylindrical, and toroidal), authors have attempted to correlate cavern size to the amount of energy introduced to the system. However the stagnant regions are highly undesirable as they are poorly mixed and they also give rise to poor mass and heat transfer.
The mixing of viscoplastic fluids with anchor impeller has been studied using numerical simulations by Tanguy a Corresponding author: HammamiMoez2003@yahoo.fr et al. [4] and Bertrand et al. [5]. In their works the authors present a detailed analysis for highly viscoplastic fluids and observe a weak increase of Metzner-Otto parameter K s (proportionality coefficient between the effective shear rate and the rotation frequency). They also showed the influence of viscoplasticity on the flow pattern and on power consumption.
Later, Anne-Archard et al. [6] have investigated the hydrodynamics and power consumption in laminar stirred vessel flow using numerical computation. The mixing systems investigated are flat-bottomed vessels equipped with double helical ribbon or anchor agitators. In their work, the authors [6] have investigated the distributions of shear rates and their link to power consumption in yield stress fluids. Fully or partially sheared flow situations with plug regions are studied. Anne-Archard et al. [6] have concluded that a constant value for the Metzner-Otto parameter is not a reasonable assumption when the Bingham number varies significantly. In fact, variations in the value of Metzner-Otto coefficient induce differences in the evaluation of the generalized Reynolds number and thus in the determination of the flow regime.
Savreux et al. [7] investigated numerically the detailed flow morphology and the 2D mixing of a viscoplastic fluid in a mixer consisting of a rotating tank and a static anchor impeller. The authors [7]  and particularly the static and moving rigid zones induced by the plasticity. It is shown that the increase in rotation velocity is not generally sufficient to reduce the size of the rigid zones. The calculations are made for different types of impeller configurations: the first consisted of chamfered anchor impeller, while two aligned blades and two perpendicular blades were added respectively for the second and the third shapes. In this work power consumption is calculated and the mixing quality is determined by particle tracking.
As we can see, heat transfer to Bingham fluids in mixing process has not been investigated in the previous numerical studies; however there is some works in Newtonian cases. Baccar and Abid [8] characterized numerically the hydrodynamic and heating performances of the anchor and gate agitators utilized for the agitation of Newtonian fluid in an agitated vessel. They notice that heat transfer is specially controlled by the turbulence viscosity and is weakly affected by the clearance between the vessel wall and the blade tip of the anchor. Moreover, it appears that the anchor and gate agitators give a similar thermal performance. Also, the evolution in transient regime of the thermal state of the agitated vessel has been studied by Baccar et al. [9] in the case of radial turbine. In this work, the highest value of the local heat transfer coefficient is found at the impeller discharge region where fluid impinges on the wall at the maximum velocity. On the other hand, the change of the heat transfer coefficient is found to be very little with angular position. Nusselt numbers predicted by Baccar et al. [9] are correlated by a dimensionless equation.
Delaplace et al. [10] investigated the heat transfer for several liquids agitated in a rounded bottom vessel equipped with an atypical helical ribbon impeller. In their work, numerical simulations of heat transfer phenomena to highly viscous Newtonian at unsteady states in a jacketed vessel were made using CFD finite volume software. No temperature gradient was detected in the bulk. Also, the authors noted that the thermal boundary layer is established after 30 s from the beginning of the heating.
Recently, Tutar and Erdogdu [11,12] carried out numerical simulation for heat transfer and velocity field characteristics of two-phase flow systems in axially rotating horizontal cans. Numerical results demonstrated that a little effect of rotation on heating up time is obtained at low speeds of rotation, but a significant effect was observed at higher speeds for air-water phase. Also, the results verified the heat transfer enhancement with higher Coriolis-centrifugal buoyancy force ratio.
In this paper, we present a two-dimensional numerical modelling of the hydrodynamic and transient thermal behaviours in a stirred vessel agitated with anchor impeller. The heat transfer is applied on the jacketed wall of the stirred vessel. The aim of this work was to investigate both mixing and heating of pseudoplastic fluids possessing yield stress by studying the effect of plasticity and impeller position on the hydrodynamic and thermal structures developed in the stirred vessel. For this purpose, the resolution of dynamics and heat transfer equations, governing the transfer phenomena developed in the tank is conducted using finite-volume technique discretization.

Mathematical formulation and numerical method
In this study, the mixing system consists in a cylindrical vessel equipped with an anchor impeller (Fig. 1). Fluids with yield stress (viscoplastic fluids) are commonly described using the Bingham model: In order to simulate the flow of fluid with yield stress Papanastasiou's regularisation model [13] is used (in the dimensionless form): Since the modelling of the mixing system was done in steady state condition, computations were conducted in a reference frame rotating with the impeller. In this case, the motion of this frame requires the introduction of centrifugal and Coriolis forces to the momentum balance equation.
Then, the resulting continuity, momentum and energy coupled equations can be written in 2D dimensionless form as follow: Continuity equation: U -velocity component: V -velocity component: Energy equation: The above equations are subject to the following initial and boundary conditions. Initially, the fluid is considered at rest and has ambient uniform temperature. So the dimensionless temperature, pressure and velocities have zero initial values.
Concerning the boundary conditions, the impeller is modelled as infinitely thin surfaces and non-slip conditions are chosen on the tank and on the agitator blades. In the vessel wall, we imposed a fixed dimensionless temperature T = 1.
The vessel was always assumed to be symmetric along its axis and then it is sufficient to consider only half of the stirred vessel.
The resolution of the equations of continuity, movement and energy, is based on the finite-volume method [14]. The derivatives in the inertia terms were discretized according to the hybrid model and the pressure field was dealt with via a prediction and correction procedure, all described by Patankar [14]. Concerning temporal discretization, an implicit formulation was adopted. The temporal integration has begun using the implicit scheme of alternated directions of Douglass and Gunn [15].

Results and discussion
In this section, we will present numerical results to give a fine knowledge of the hydrodynamic characteristics and the temporal evolution of the thermal behavior induced in an agitated mixing vessel. A particular interest is provided to study the controlling parameters characterizing the fluid flow and the rigid zones sizes and forms. These parameters are Reynolds number, Oldroyd number and the impeller position.
Next, the results are turned into a fixed frame. The thermal results are given for fixed values of Prandtl number (P r = 7× 10 3 ), dimensionless impeller length (l/R = 0.2) and dimensionless rotor radius (rot/R = 0.05).

Rigid zone characteristics
When the applied stress is beneath a threshold value called the yield stress, rigid zone is obtained and moves like a solid. We note that three types of rigid zones can be detected (Fig. 2): -Zones Z1: cling to the impeller on both sides.  To obtain a precise determination of the rigid zones, the influence of the mesh on the results and in particular on these zones was studied. The grid-independent test is used to determine the minimum grid resolution that achieves a solution independent of the mesh quantities. In this test three different mesh sizes were evaluated: 6480 (coarse mesh), 14 688 (fine mesh), and 25 920 (very fine mesh). Figure 3 shows that, a computational domain consisting of 108 × 136 grid points with uniform grid spacing respectively in radial and angular directions was found to be sufficient for producing accurate results at reasonable computed time. However, to have more accuracy we must use the very fine mesh (25 920) or the mesh must be highly refined near the impeller [7] to obtain the rigid zones Z1 stuck to the impeller. Nevertheless, quantitatively, the overall area of unyielded zones Z1 can be neglected (less than 3% of the total rigid zone area [7]). On the other hand, it can also be seen in Figure 4 that the effect of these zones on the hydrodynamic behavior is not significant.

Yield stress effect on the hydrodynamic behaviour
For R 1 = 0.8 and Re = 1, Figure 5 shows the effect of the Oldroyd number on the flow pattern inside the tank. A similar hydrodynamic behaviour to the Newtonian case is observed for Od = 5, a vortice is formed between the impeller and the wall of the stirred vessel. In these cases the flow is strongly dominated by the tangential component. When Od = 50, the shape of the hydrodynamic structure changes and a second recirculation zone takes birth nearby the other side of the agitator. Also far from the impeller, when the yield stress increases, the fluid velocity drastically decreases.
To further investigate the effects of viscoplasticity, Figure 6 presents the non-dimensional radial profiles of tangential velocities taken at the median and the impeller planes for R 1 = 0.8 and Re = 1. As soon as plasticity effects are present, a linearity of the tangential velocities profiles is observed in the median plane. This indicates that the fluid exhibits no deformation and solid structures Z2 are formed. Also for Od = 50, it is noteworthy that two zones Z3 are formed in the peripheral regions where there is a poor transfer. On the other hand, it can be seen at the impeller plane (Fig. 6a) that recirculation flow becomes. Figure 7 shows the rigid zones for R 1 = 0.8 and different Re and Od numbers. It can be noted that two solid cores Z2 appear on either side of the tank center even for small yield stress. As can be seen in Figure 7, the solid cores Z2 size increases considerably with Od. However, the two parts of the peripheral rigid zone Z3 take birth only for high values of Od. This can be explained  by the fundamental variation in the structure of the flow for Od = 50: the fluid skirts around the impeller and the flow is clearly braked in the peripheral zones (Fig. 5).

Yield stress and inertia effects on the rigid zones
If Od number has a strongly effect on the conformation of the rigid zones, however, as it can be seen in Figure 7, Re number has a negligible influence on the unyielded zones. Also for Od = 5, the increase in the inertia is not sufficient to reduce the size of the rigid zones. We note that in some specific configurations, especially for a low yield stress value, when Re = 20 the rigid zones become asymmetrical.
To verify the accuracy of our simulation, a comparison is conducted in Figure 7, where we have reproduced the rigid zones obtained by Savreux et al. [7] for Re = 0.05 and R 1 = 0.8. We can notice a good agreement between our results and those of the previous study.

Yield stress effect on the thermal behaviour
To understand better the thermal treatment of mixed Bingham viscoplastic fluid, we have reproduced the thermal behaviour evolution with time for fixed Prandtl number (P r = 7 × 10 3 ). The fluid is heated by a jacketed wall of the stirred vessel. For Od = 50, Re = 1 and R 1 equal to 0.8, Figures 8 and 9 show the evolution with time of the temperature distribution in the agitated mixing vessel. The results show the unexpected trend of rapid temperature increasing in the peripheral region. In fact, in the  Z3 zones we have a nearly one-dimensional conduction regime, but the vortices, formed between the impeller and the wall of the stirred vessel, enhance the heat transfer. The same phenomenon is observed in the middle of tank where the second recirculation zone improves the transfer.
In order to give more accurate information about the Od number effect on the temperature distribution, Figure 10 represents temperature radial profiles in the median plane at t = 1500 for Re = 1. It is found that there is negligible difference between the Newtonian case (Od = 0) and Od = 5. Also for Od = 50, the effect of the rigid zones Z2 is not clear. This is can be explained by the fact that the angular rotation speed of the confined fluid is different of the impeller one. So these zones are unlikely to have any impacts for mixing efficiency: there is not a stagnant zone, but the fluid is only sometimes confined in the rigid zones and then the particles can migrate outside these zones and be mixed.
In Figure 11, showing the temporal evolution of the averaged temperature in the mixed tank for R 1 = 0.8, we can see that the temperature increases more rapidly when Od = 50. This unpredicted tendency can be explained by the generation of important vortexes on both sides of the impeller. In fact, conjugated cells improve the fluid homogenisation, and subsequently, convective transfer will be performed.
In order to characterize the instantaneous mean convective heat transfer coefficient through the wall of the agitated vessel, obtained by integrating of the local coefficient over its correspondent exchange surface, we calculate the instantaneous average Nusselt number, defined as follows: To study the Oldroyd number effect on the transient heat transfer, the averaged Nusselt number evolution with time is plotted in Figure 12 for various Re and Od numbers and a radius R 1 = 0.8. In general, it is found that N u(t) decreases at first and then varies slowly to attend a constant value. The results show the trend of increasing N u with Re number. On the other hand, when Od = 50 heat transfer is promoted. This tendency is in adequacy with the velocity profile ( Fig. 6) that shows that the recirculation flow is intense when the yield stress increases. Figure 13 resumes the effect of the impeller position R 1 on the location and shape of the rigid zones. Three configurations are studied in this section: R 1 = 0.7, 0.8 and 0.9. It can be noticed that the decrease of the radius R 1 (equal to 0.7) can reduce the percentage of rigid zone for the case Od = 5. However, for Od = 50 the Z3 zones tend to join together to form a continuous annulus and this configuration prevents the fluid homogenisation. For high Od number, it appears that the configuration with R 1 = 0.9 is effective for reducing the rigid zones, in particular Z3.

Effect of the impeller position on the rigid zones
To verify this behavior, we have reproduced in Figure 14, the average temperature evolution for different radii R 1 (0.7, 0.8 and 0.9) and Od numbers (5 and 50) with Re = 1. This figure clearly demonstrates the importance of the impeller position effect on the heat transfer performance. Moreover, it was noted that for Od = 5 the reduction in rigid zone size when R 1 = 0.7 is not a criterion to indicate the heating quality improving. Also, although the same zone rigid percentage is almost obtained for R 1 = 0.7 and R 1 = 0.8 when Od = 50, there is a large difference in thermal behaviours. However, it seems that the configuration R 1 = 0.9 was the best for agitating and heating viscoplastic fluids. Also, for this configuration there is not a considerable difference between Od = 5 and 50 cases.

Conclusion
In this study, the mixing and heating of Bingham viscoplastic fluid in a vessel agitated by an anchor impeller  were simulated to analyze the complex flow structures velocities and temperature distributions in transient regime.
Indeed, resolution of coupled momentum and heat transfer equations give interesting local information concerning the hydrodynamic characteristics and the evolution with time of the thermal behaviours. The variation of the Nusselt number with the rheological parameters is also presented. Firstly, we have put in evidence the importance of the yield stress effect on the flow structure. We have demonstrated that the area of the unyielded zones increases drastically when the Oldroyd number increases. However, the recirculation flow is intense near the impeller when the yield stress becomes important.
Secondly, it is found that the convective heat transfer is improved when the yield stress is important as a consequence of the apparition of recirculation zones. In fact, in spite of the increases of the rigid zone area, these zones are unlikely to have any repercussion for heating efficiency: the fluid is only sometimes confined in the rigid zones and then the particles can migrate outside these zones and be heated.
Thirdly, we have demonstrated that the increasing of the ratio R 1 improves heat transfer between the jacketed wall of the stirred vessel and the mixed fluid. In the particular case of R 1 = 0.7 and Od = 5, we have noted that the decreasing in the rigid zone percentage masks the decreasing of the heating speed. Thus, it was noted that the reduction in rigid zone size is not obligatory the criterion of improving heating performance.