NUMERICAL INVESTIGATION OF LAMINAR MIXED CONVECTION IN A SQUARE CROSS-SECTIONED CYLINDRICAL ANNULAR ENCLOSURE

Steady-state laminar mixed convection of Newtonian fluids in a square cross-sectioned cylindrical annular enclosure with rotating inner wall and heated top cover has been numerically analysed based on axisymmetric incompressible flow simulations. Richardson number, Reynolds number and ri/R effects on heat and momentum transport have been investigated for the range of Richardson number 0 ≤ Ri ≤ 1, Reynolds number 500 ≤ Re ≤ 2000 and 0.25 ≤ ri/R ≤ 8 at a representative value of Prandtl number (i.e. Pr = 1.0). A scaling analysis has been also carried out in order to elucidate the possible influences of Reynolds, Richardson and Prandtl numbers and ri/R on the mean Nusselt number. It has been found that the mean Nusselt number Nu ̅̅ ̅̅ demonstrates a monotonically decreasing trend with increasing Ri whereas Nu ̅̅ ̅̅ increases with increasing ri/R and Re which is consistent with scaling estimation. It is also observed that the flow pattern in the case of purely forced convection (i.e. Ri = 0) is significantly different from those in mixed convection (i.e. Ri > 0). In the case of Ri = 0 (i.e. purely forced convection), a one-cell flow structure with two small vortexes on the top corners is observed for ri/R ≤ 1, whereas a second cell appears in the flow field for ri/R > 1 at Re = 1000. On the other hand, in the case of mixed convection (i.e. Ri > 0), two-cell and four-cell flow structures occur in the flow field depending on Ri and ri/R for the range of Ri, Re and ri/R considered here at Pr = 1.0. Based these observations, a flow regime diagram has been proposed here for mixed convection (i.e. Ri > 0) for the range of Ri, Re and ri/R analysed in this study.


INTRODUCTION
Mixed convection in a cylindrical container with rotating one of the wall has been extensively analysed due to its wide range of applications such as chemical processing, bio-chemical synthesis, polymer processing, food preparation, pharmacology. In these applications, mixed convection plays a vital role not only in heat transfer rate but also influences the mixing rate for low Reynolds number cases. Therefore, it is necessary to investigate heat transfer characteristics and flow structure in this configuration so that the rates of heat transfer and mixing can be optimised. However, heat transfer and mixing rates in this configuration depend on many parameters such as container geometry and the rotational speed of the wall. Two special cases of rotating flow problems which are the flows on top a rotating disk and inside an enclosure with a rotating one of the end covers have been commonly investigated in the current literature. Theodore von Karman pioneered the analysis of flows on top a rotating disk and such flows are commonly referred to as von-Karman flows, and an extensive review of such flows is provided in Ref. [1]. In addition to this, flow in cylindrical enclosures with a rotating cover has also been extensively analysed from various different viewpoints due to its wide range of roles in a range of different engineering applications. Vogel [2,3], Ronnenberg [4] and Bertela and Gori [5] analysed fluid flows in cylindrical enclosures with a rotating end wall, and the findings of these studies [1][2][3][4] have subsequently been extended by Escudier [6] based on an experimental analysis where the criterion for vortex breakdown has been proposed in terms of aspect ratio H/R and Reynolds number Ω 2 / . Besides these experimental studies, several numerical investigations [7][8][9][10][11][12] analysed fluid flows in cylindrical enclosures with a rotating end wall. The flow produced in a conical container by a rotating end wall has also been numerically analysed by Escudier et al. [12] who reported that vortex breakdown is suppressed beyond a certain angle of inclination of the sidewall for both convergent (increasing radius towards the rotating end wall) and divergent (decreasing radius towards the rotating end wall) geometries.
The analysis of heat transfer characteristics in cylindrical enclosures with a rotating cover received relatively limited attention [9][10][11]. Lee and Hyun [10] analysed the effects of Prandtl number on heat transfer rate in this configuration and revealed that Prandtl number has an important influence on the heat transfer characteristics and advective transport has been found to strengthen with increasing Prandtl number. Iwatsu [11] investigated the effects of Reynolds and Richardson numbers at Pr = 1, in the range of 100 ≤ Re ≤ 3000, and 0 ≤ Ri ≤ 1, on the flow pattern and heat transfer rate for swirling flows in cylindrical enclosures with an aspect ratio of unity (i.e. AR = H/R = 1), and a heated rotating top wall based on numerical simulations. The analysis by Iwatsu [11] revealed that advective (diffusive) transport strengthens (weakens) and accordingly, the mean Nusselt number increases with decreasing Richardson number. Existing analyses on flow induced by the rotation of one of the end covers in a cylindrical enclosure have been also summarised in a recent paper by Turan et al. [13] according to the boundary conditions, the governing non-dimensional parameters and the nature of the investigation (i.e. whether it is experimental or numerical).
Majority of the investigations in the existing literature related to mixed convection in enclosure spaces have been carried out for rectangular enclosure with rotating top or bottom cover. Although cylindrical annular space with rotating inner wall is more relevant to real engineering applications than rectangular enclosures with rotating top or bottom cover (especially in mixing processes), there is no study where laminar mixed convection of Newtonian fluids in cylindrical annular enclosures has been analysed in detail in existing literature. This gap in existing literature has been addressed here by numerically analysing mixed convection of Newtonian fluids in a square cross-sectioned cylindrical annular enclosure with rotating inner wall and heated top cover for the range of Richardson and Reynolds numbers 0 ≤ ≤ 1 and 500 ≤ ≤ 2000 respectively for 0.25 ≤ / ≤ 8 at a representative value of Prandtl number (i.e. = 1.0). In this respect, the main objective of the present paper is to demonstrate the influences of Richardson and Reynolds numbers and / on mixed convection induced by a rotating inner wall in a cylindrical annular enclosure with an aspect ratio of unity.
The rest of the article will be organised as follows. The necessary mathematical background and numerical implementation will be discussed in the next section, which will be followed by a detailed scaling analysis. Following these sections, results will be presented and subsequently discussed. The main findings are summarised and conclusions are drawn in the final section.

MATHEMATICAL BACKGROUND Governing equations and boundary conditions
In this study, the flow is assumed to be laminar, incompressible, steady and axisymmetric (i.e. twodimensional). The conservation equations in the cylindrical coordinate system take the following form for steady-state incompressible axisymmetric swirling flows: Mass conservation equation

Momentum conservation equations
Energy conservation equation where is the dynamic viscosity, is the reference temperature for evaluating the buoyancy term ( − ) in the momentum conservation equation in the vertical direction, and here is taken to be the cold cover temperature . In addition, thermo-physical properties (thermal conductivity, specific heat, viscosity etc.) are considered to be constant and independent of temperature in this analysis for the sake of implicitly. The numerical investigation is carried out for an axisymmetric domain which is schematically shown in Figure 1. The aspect ratio (AR = H/R) of the cylindrical container is considered to be unity (i.e. AR = H/R = 1). The bottom and top covers of the cylindrical enclosures are kept at different temperatures ( < ), while the cylindrical surface is considered to be adiabatic in nature. The temperature difference between the top and bottom covers are maintained small enough to ensure that Boussinesq approximation remains valid. No-slip boundary condition and impenetrable walls indicate that all the velocity components are identically zero at the non-rotating walls. For the rotating inner wall, radial and axial velocity components are identically zero due to impenetrability and no-slip conditions, and the tangential velocity component is given by ~(1 + / ).

Non-dimensional numbers
The ratio between inertial and viscous forces represents the strength of the forced convection component in this analysis. This ratio can be quantified by the Reynolds number Re which is defined as: The natural convection component of mixed convection can be characterised by the Grashof number, which represents the ratio of buoyancy to viscous forces, in the following manner: Another important non-dimensional number is Rayleigh number, Ra, which demonstrates that the ratio of the strengths of thermal transports due to buoyancy to thermal diffusion, which is defined here in the following manner: where Pr is the Prandtl number, which is defined as: Prandtl number shows the ratio of momentum diffusion to thermal diffusion. The Prandtl number also represents the ratio of the hydrodynamic boundary layer to thermal boundary layer thicknesses. In addition to this, Richardson number is a well-known non-dimensional parameter for mixed convection, which is used to evaluate the relative importance of the natural to forced convection: The rate of convective heat transfer is generally characterised by the heat transfer coefficient h, which is expressed in a non-dimensional form in terms of the Nusselt number Nu, as: and the heat transfer coefficient h is defined as: where subscript 'wf' refers to the condition of the fluid in contact with the wall, Twall is the wall temperature and Tref is the appropriate reference temperature, which can be taken to be TC (TH) for the hot (cold) wall respectively. For this configuration the mean heat transfer coefficient ℎ ̅ and ̅̅̅̅ are given by: ℎ ̅ = ∫ ℎ2 / ( 2 − 2 ) and ̅̅̅̅ = ℎ ̅ / , respectively.

Numerical implementation, grid-independency, and bench-marking
In this study, the conservation equations of mass, momentum and energy have been solved in the framework of a finite-volume method using a computational fluid dynamics (CFD) software ANSYS-FLUENT, which was previously successfully used for simulations of the flows induced by rotating one of the end covers of a cylindrical container [11,13]. A second-order central difference scheme is used for the discretisation of the diffusive terms and a second-order up-wind scheme is used for the convective terms. Coupling of pressure and velocity is achieved using the well-known SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm [14]. The convergence criteria have been taken to be 10 -7 for all the relative (scaled) residuals.
Three different non-uniform meshes M1 (100 × 100), M2 (200 × 200) and M3 (300 × 300) have been investigated for each configuration, and the details of these meshes have been provided in Table 1 where the normalised minimum grid spacing ∆ , / and grid expansion ratio are provided. The numerical uncertainties for the mean Nusselt number ̅̅̅̅ for = 1.0 and / = 1 at = 1000 and = 1 are shown in Table 1. It indicates that the maximum relative error levels (ea) among the meshes are under 1.0%. Based on this analysis, the simulations have been conducted using mesh M2 (200 × 200) for this configurations, which is found to be sufficient for providing high accuracy and computational efficiency. In addition to a grid-independency study, the simulation results have also been compared to the study by Iwatsu [11], in which mixed convection of Newtonian fluids in a cylindrical enclosure with heated rotating top cover, for different Richardson and Reynolds number values at = 1.0. As shown in Figure 2, the present simulations results remain in excellent agreement with the corresponding benchmark data reported by Iwatsu [11].

SCALING ANALYSIS
A scaling analysis has been carried out in order to elucidate the possible influences of Reynolds, Richardson and Prandtl numbers and / on the mean Nusselt number. The wall heat flux can be scaled as: Using Eq. (11), the Nusselt number can be scaled as: where fı is a function of Re, Ri, and Pr, which accounts for the ratio of hydrodynamic to thermal boundary layer thicknesses (i.e. ℎ ⁄~1( , , )). In order to determine the hydrodynamic boundary thickness , the order of magnitudes of inertial and viscous forces in the radial direction can be equated: The shear stress can be scaled as ~ ( / ) and thus Eq. (13) can be rewritten as: Using Eq. (14), hydrodynamic boundary thickness can be estimated as: Here, and are the characteristic velocity scales in radial and tangential directions respectively. Eq.
Equation (18) offers important physical insights into the influences of Re, Ri, and / on the mean Nusselt number ̅̅̅̅ . In the following section, this scaling predictions will be used for discussing Ri and / effects on ̅̅̅̅ .

Influences of
The variations of the mean Nusselt number ̅̅̅̅ with Richardson number are presented in Figure 4 for different / values at = 1000 and = 1. Figure 4 shows that ̅̅̅̅ monotonically decreases with increasing which is consistent with the scaling estimation given by Eq. (18). This can be explained in the following manner. For mixed convection (i.e. > 0), the relative strengths of inertial, buoyancy and viscous forces determine the flow behaviour. For = 0, which corresponds to purely forced convection, the flow is governed by the inertial and viscous forces. The influence of buoyancy force starts to strengthen with increasing , and therefore the competition between buoyancy and viscous forces becomes increasingly important with increasing . However, this configuration, where the top cover is hotter than the bottom one, represents a stable condition where the lighter hot fluid sits on top of heavier cold fluid so the effects of natural convection remain significant only close to the heated top cover and thus the net advective transport weakens and the mean Nusselt number decreases with increasing Ri. In addition, it is also observed from Figure 4 that the influence of Richardson number on the mean Nusselt number weakens with increasing / and a decreasing in the mean Nusselt number with increasing Ri remains relatively limited levels in the case of high / values. This is also expected from the scaling estimation given by Eq. (18). Eq. (18) explicitly indicates that the ratio of /(1 + / ) decreases with increasing / . As a results of this, the influence of Richardson number on the mean Nusselt number decreases with increasing / . The observations made from Figure 4 can further be explained from the non-dimensional temperature countors , which are presented in Figure 3 for different values at = 1000 and = 1. It can be observed from Figure 3 that isotherms become increasingly parallel to the horizontal walls with increasing Ri in the case of small / (i.e. / = 0.25). This indicates that convection weakens with increasing Ri and thermal transport takes place predominantly due to conduction for large values of Ri. This leads to a decrease in the mean Nusselt number with increasing Ri as observed in Figure 4. On the other hand, it is also observed from Figure 3 that the influence of Richardson number on the isotherm weakens with increasing / (i.e. / = 4) which is also consistent with the observation in Figure 4.  Figure 5. It can be seen from Figure 5 that ̅̅̅̅ increases with increasing / for both low and high Ri values. This is consistent with the scaling estimate of the mean Nusselt number given by Eq. (18) which also suggests that ̅̅̅̅ is expected to increase with increasing / . An increase in Nusselt number with increasing / bears the signature of strengthening of advective transport. This can be substantiated from Figure 6 Figure 6 demonstrates that the magnitude of increases with increasing / values, which is an evident of strengthening of advective transport in the flow field. This also leads to an increase in the mean Nusselt number observed in Figure 5. The strengthening of advective transport with increasing / , which is observed in Figs. 5 and 6, could alternatively be explained by integrating convective heat transport through the boundary layer thickness on the bottom cover: where is the hydro-dynamic boundary layer thickneess on the horizontal walls. Substituting

Influences of
The variation of the mean Nusselt number ̅̅̅̅ with Reynolds number for different / values at = 0.1 (i.e. a representative mixed convection case) and Pr = 1 are shown in Figure 7. It can be seen from Figure  7 that ̅̅̅̅ increases with increasing Re for all / case. This is consistent with the scaling estimate of the mean Nusselt number given by Eq. (18) which also suggests that ̅̅̅̅ is expected to increase with increasing . In addition to this, the distributions of non-dimensional swirl velocity component (= / ) along the vertical mid-plane ( / = 0.5) are shown in Figure 8 for different values / at = 0.1 (i.e. a representative mixed convection case) at Pr = 1. It can be seen from Figure 8 that the magnitude of increases with increasing for all / cases. This indicate a strengthening of advective transport with increasing Reynolds number, which leads to an increase in heat transfer rate in the enclosure as observed in Figure 7. This can also be explained by Eq. (21a) which indicates that the convective heat transport strengthens with increasing as observed in Figs. 7 and 8. For this reason, the mean Nusselt number ̅̅̅̅ increaes with increasing for all / cases (Figure 7).

Flow regimes
The variation of non-dimensional stream functions Ψ = / are shown in Figure 9 for different and / values at = 1000 and = 1. It is apparent from Figure 9 that the flow pattern in the case of purely forced convection (i.e. = 0) is significantly different from those in mixed convection. In the case of = 0 (i.e. purely forced convection), a one-cell flow structure with two small vortexes on the top corners is observed for / ≤ 1, whereas a second cell appears in the flow field for / > 1 at = 1000. On the other hand, in the case of mixed convection (i.e. > 0), two-cell and four-cell flow structures occur in the flow field depending on , and / for the range of and / considered here at = 1000. It is also worth noting that similar behaviour has also been observed for different Reynolds number values (i.e. Re = 500 and 2000). Based these observations, a flow regime diagram has been proposed here for mixed convection (i.e. > 0) for the range of and / analysed in this study. According to this diagram, the flow patterns in the mixed convection for this configuration are classified into two different zones in terms of Richardson number and / as shown in Figure 10. The region, where + 1 ≤ 1.75 ( / ) 0.22 ( + 1 ≤ 1.50 ( / ) 0.22 ) for ≤ 1000 ( > 1000) is termed as the Regime 1, which exhibits two-cell flow pattern. The Regime 2 is characterised by + 1 > 1.75 ( / ) 0.22 ( + 1 > 1.50 ( / ) 0.22 ) for ≤ 1000 ( > 1000) where a four-cell flow pattern occurs. It needs to be highlighted that the boundaries which distinguish one regime from another in Figure 10 are based on the observations made from simulation results. As such these boundaries should not be treated rigidly but need to be considered only in an order of magnitude sense.