Numerical Analysis of Mixed Convective Heat Transfer from a Square Cylinder Utilizing Nanoﬂuids with Multi-Phase Modelling Approach

: The present study deals with the numerical simulation of mixed convective heat transfer from an unconﬁned heated square cylinder using nanoﬂuids (Al 2 O 3- water) for Reynolds number ( Re ) 10–150, Richardson number ( Ri ) 0–1, and nanoparticles volume fractions ( ϕ ) 0–5%. Two-phase modelling approach (i.e., Eulerian-mixture model) is adopted to analyze the ﬂow and heat transfer characteristics of nanoﬂuids. A square cylinder with a constant temperature higher than that of the ambient is exposed to a uniform ﬂow. The governing equations are discretized and solved by using a ﬁnite volume method employing the SIMPLE algorithm for pressure–velocity coupling. The thermo-physical properties of nanoﬂuids are calculated from the theoretical models using a single-phase approach. The ﬂow and heat transfer characteristics of nanoﬂuids are studied for considered parameters and compared with those of the base ﬂuid. The temperature ﬁeld and ﬂow structure around the square cylinder are visualized and compared for single and multi-phase approaches. The thermal performance under thermal buoyancy conditions for both steady and unsteady ﬂow regimes is presented. Minor variations in ﬂow and thermal characteristics are observed between the two approaches for the range of nanoparticle volume fractions considered. Variation in ϕ affects C D when Reynolds number is varied from 10 to 50. Beyond Reynolds number 50, no signiﬁcant change in C D is observed with change in ϕ . The local and mean Nusselt numbers increase with Reynolds number, Richardson number, and nanoparticle volume fraction. For instance, the mean Nusselt number of nanoﬂuids at Re = 100, ϕ = 5%, and Ri = 1 is approximately 12.4% higher than that of the base ﬂuid. Overall, the thermal enhancement ratio increases with ϕ and decreases with Re regardless of Ri variation. ﬂuid. The temperature ﬁeld and ﬂow structure around the square cylinder are visualized and presented through the results obtained from numerical simulations performed based on multi-phase modelling.


Introduction
Nanofluids are the engineered colloidal suspension of nanosized particles (10-100 nm) in a base fluid [1]. They have greater thermal conductivity compared to their base fluids and are considered smart fluids because the heat transfer rates can be controlled to some extent by varying nanoparticle concentration in the base fluid [2]. Nanofluids find applications in nuclear reactors [3], electronic cooling [4], heat exchanger [5], heating buildings in cold regions [6], drying of various materials [7], and automotive applications [8,9]. Numerous experimental and numerical studies have proved that heat transfer rates augment with an increase in nanoparticle volume fractions in base fluid [10][11][12][13]. As nanofluids enhance the heat transfer rates when they flow over the surface of a hot body, it would be interesting to investigate the flow and heat transfer from a hot square cylinder placed in a uniform stream of nanofluid.
right-side walls are at different constant temperatures while the top and bottom side walls are insulated moving lids. A two-phase mixture model is adopted to analyse the thermal behaviour of nanofluids for various enclosure inclination angles ranging from −60 • to +60 • . The results show that the presence of nanoparticles and their addition in base fluid enhances the heat transfer in the cavity significantly and causes notable changes in the flow pattern. Besides, the effect of inclination angle is more pronounced at higher Ri (≥10) only. Esfe et al. [27] conducted a 3D-numerical study on free convection in a cubical cavity with porous fins using nanofluid (CuO-water) as the working fluid. The two-phase mixture model is employed, and the numerical simulations are performed. They have analyzed the result for different Rayleigh numbers (10 3 -10 5 ). The numerical results show that at lower Rayleigh number (∼ 10 3 ) the value of mean Nusselt number is relatively small since conduction dominates over convection. On the other hand, a high Rayleigh number enhances the convection term that becomes more dominant than conduction, which increases the average Nusselt number, since fluid flow irregularity intensifies consequently. Increasing the nanoparticle volume fraction shows an improvement in heat transfer as reported in the aforementioned studies. Garoosi et al. [28] performed numerical simulations considering the steady state mixed convection heat transfer (0.01 ≤ Ri ≤ 1000) of nanofluid (Cu-water, Al 2 O 3 -water, and TiO 2 -water) using a two-phase mixture model. The authors used a two-sided lid driven cavity with several pairs of heaters and coolers. It is found that at a low Ri value, the heat transfer rate increases with the increasing number of heaters and coolers (HACs). On the other hand, at high Ri, Nu M does not change significantly after a saturated number of HACs. The results reveal that the heat transfer rate increases by reducing the diameter of the nanoparticles and Ri. Although the thermal conductivity of Cu is considerably higher than that of TiO 2 , still the difference in heat transfer is small. Darzi et al. [29] studied the effect of nanofluid on combined convective heat transfer inside a finned enclosure. It is found that adding a small concentration of nanoparticles to base fluid enhances the heat transfer, especially at a low Ri value.
Recently, Ebrahimi et al. [30] investigated the heat transfer from longitudinal vortex generators placed inside microchannel heat sinks, with CuO-water and Al 2 O 3-water nanofluids as the working fluids. Improvement in heat transfer in the ranges 2.55-29.05% and 9.78-50.64% was observed for Al 2 O 3 -water and CuO-water nanofluids, respectively. The same research group (Naranjani et al. [31]) further investigated the thermal as well as hydraulic performance of heat sink with corrugated channels using nanofluids. About 22-40% enhancement in heat transfer was reported while using nanofluids compared to water as coolant. Quite recently, Saieesh and Prasad [32] investigated the laminar forced convective heat transfer from a square cylinder using the Eulerian-Eulerian mixture model considering slip velocity for Reynolds numbers 10-40 and ϕ = 0-5%. They showed the influence of slip velocity on heat transfer performance of nanofluids. Other recent works dealing with heat transfer with nanofluids include those of Arjun and Rakesh [33] and Sawicka et al. [34].
A review of the archival literature on nanofluid flow around a square cylinder reveals that no study has been made on this topic using the multi-phase modelling approach, in the presence of thermal buoyancy. In addition, no comparison has been made between the two approaches for this flow configuration. In this study, the results obtained from single and multi-phase approaches using SIMPLE algorithm are presented along with the detailed investigations. In particular, the effects of considered parameters (i.e., Re, Ri, ϕ) on the local and global flow and heat transfer characteristics have been presented and discussed. In the present study, nanofluid flow past a hot square cylinder is studied and the resulting heat transfer is compared with that of the base fluid. The temperature field and flow structure around the square cylinder are visualized and presented through the results obtained from numerical simulations performed based on multi-phase modelling.

Problem Description and Geometrical Configuration
A two-dimensional square cylinder (height, D) placed in an unbounded domain is considered as shown in Figure 1a. The cylinder is hot with its surface maintained at a constant temperature θ w greater than that of the ambient fluid that flows in the positive x-direction. The ambient fluid is Al 2 O 3-water nanofluid flowing with a uniform velocity (U 0 ) and has a temperature θ 0 . The thermal buoyancy is considered with the gravity acting opposite to the flow direction and parallel to the flow inlet. The heated cylinder exchanges heat with the nanofluid flowing past its surface, resulting in heat transfer. The boundaries of the computational domain are placed far from the heated cylinder and appropriate boundary conditions are used while performing the calculations. The simulations are carried out for both steady and unsteady flow regimes. The heat transfer taking place from the cylinder can be modelled with appropriate governing equations to predict the flow and thermal field. Since nanofluids contain solid particles dispersed in liquid medium, two kinds of approaches can be followed. One is the single-phase modelling approach and the other is the two-phase modelling approach. In the next section, we describe both the approaches.

Problem Description and Geometrical Configuration
A two-dimensional square cylinder (height, D) placed in an unbounded domain is considered as shown in Figure 1a. The cylinder is hot with its surface maintained at a constant temperature greater than that of the ambient fluid that flows in the positive x-direction. The ambient fluid is Al2O3-water nanofluid flowing with a uniform velocity ( 0 ) and has a temperature 0 . The thermal buoyancy is considered with the gravity acting opposite to the flow direction and parallel to the flow inlet. The heated cylinder exchanges heat with the nanofluid flowing past its surface, resulting in heat transfer. The boundaries of the computational domain are placed far from the heated cylinder and appropriate boundary conditions are used while performing the calculations. The simulations are carried out for both steady and unsteady flow regimes. The heat transfer taking place from the cylinder can be modelled with appropriate governing equations to predict the flow and thermal field. Since nanofluids contain solid particles dispersed in liquid medium, two kinds of approaches can be followed. One is the single-phase modelling approach and the other is the two-phase modelling approach. In the next section, we describe both the approaches.

Single-Phase Model (SPM)
The single-phase modelling approach assumes that the nanoparticles in the base fluid can be easily fluidized and reach the base fluid velocity. Therefore, the nanofluid is considered a homogeneous fluid. Further, it is assumed that both the liquid and particle The single-phase modelling approach assumes that the nanoparticles in the base fluid can be easily fluidized and reach the base fluid velocity. Therefore, the nanofluid is considered a homogeneous fluid. Further, it is assumed that both the liquid and particle phases are in thermal equilibrium and move with the same velocity [35]. The governing equations of continuity, momentum, and energy [36] are given as follows: • Continuity equation: • Momentum equation: • Energy equation: where, and K n f are the velocity vector, density, viscosity, thermal expansion coefficient, specific heat, and thermal conductivity, respectively. The difference between the solid surface temperature and free-stream average temperature is used as the characteristic temperature difference.

Nanofluids Modelling
The volume fraction of nanofluid (ϕ) is the volumetric concentration of the nanoparticles in the nanofluid. The effective properties of nanofluids such as the effective density, viscosity, thermal expansion coefficient, and thermal conductivity are given by • Effective viscosity (Brinkman [37]): • Effective specific heat (Xuan and Roetzel [38]): • Effective thermal conductivity (Xie et al. [39]): where, g is ratio of the nanolayer thickness (δ) to the original particle radius (r p ). The nanoparticle diameter (2r p ) and nanolayer thickness (δ) are taken as 30 nm and 2 nm, respectively.
The thermal expansion coefficient of nanofluids can be estimated utilizing the volume fraction of nanoparticles on a weight basis and is given by [40] (ρβ) n f = 1 − ϕ np (ρβ) b f + ϕ np (ρβ) np . Unlike the other models (i.e., Eulerian model, VOF model), the mixture model is based on a single fluid but two-phase approach. The phases have their velocity vector (i.e., interpenetrating) closely following the flow, and coupling between them is substantial [26]. In the control volume, primary and secondary phases have their separate volume fractions. The primary phase influences the secondary phase via drag and turbulence, while the secondary phase, in turn, influences the primary phase via a reduction in mean momentum and turbulence. The mixture model is based on the following assumptions [20]: (1) all phases are allocated a single pressure; (2) the secondary phase is assumed to be spherical in shape with uniform particle size and their interactions between different dispersed phases are neglected; (3) the concentrations of the secondary dispersed phases are solved from scalar equations taking into account the correction due to phase slip. It is to be noted that turbulence generation in the secondary phases is not accounted for, nor is the turbulence of the primary phase directly affected by the presence of the secondary phase [41]. Instead of utilizing the governing equations of each phase separately, the continuity, momentum, and energy equations for the mixture are employed and written in the dimensional form as

•
Continuity equation: • Momentum equation: • Energy equation: • Volume fraction equation: In the MPM model, each phase has its own velocity vector field and within a given control volume there exists a certain fraction for each phase. In addition to Equations (9)-(11), the mixture model solves the volume fraction equation for the secondary phase. It then uses an algebraic expression to calculate the relative velocity between the phases where, V m , ρ m , µ m , K m , and P m are velocity, density, viscosity, thermal conductivity, and pressure of mixture; V dr,np and V dr,b f are the drift velocity of nanoparticle and drift velocity of base fluid; ρ b f and c (p,b f ) are the density and specific heat of the base fluid; ρ b f and c (p,n f ) are the density and specific heat of nanofluids; φ is the solid volume fraction of the nanoparticles, respectively. The mixture velocity → V m is determined as follows: . In this equation, → V dr,np is the drift velocity for the secondary phase and is expressed as, The velocity of the secondary phase in relation to the primary phase is known as the relative or slip velocity and it is defined as The drift velocity is related to the slip velocity as Energies 2021, 14, 5485 7 of 26 The following equations are proposed by Manninen et al. [42] and Schiller [43] to calculate the slip velocity → V (np,b f ) and drag function (f drag ), respectively: In the above equation, and the acceleration is given as For mixture model calculations, the thermo-physical properties of nanofluids have been taken from their respective models [37][38][39]. The physical properties of nanoparticle and base fluid are shown in Table 1.

Boundary Conditions
At the inlet boundary, a uniform flow profile (i.e., U 0 = 1, and V = 0) is assumed. A zero-shear boundary condition is specified along the top and bottom boundaries, (i.e., ∂U/∂Y = 0, and V = 0) of the domain. The right-side boundary is designated as the outlet. This boundary is located sufficiently far downstream from the cylinder, and it is considered as the pressure outlet (i.e., default option in FLUENT, known as 'PRESSURE OUTLET'), which assumes a zero-gauge (static) pressure (P = 0) for the operating pressure. For the velocities, the following boundary conditions are used: ∂U/∂Y = 0, and ∂V/∂X = 0 [46]. The no-slip condition (i.e., U = 0, and V = 0) is applied on the cylinder surface. For thermal boundary conditions, the top and bottom boundaries of the domain are assumed to be adiabatic, except on the heated cylinder where the non-dimensional temperature (i.e., Θ = (Θ − θ 0 )/ (θ w − θ 0 )) is unity, and the inlet boundary is kept as Θ = 0. The physical properties of the fluid are assumed to be constant except for the density in body force, which varies linearly with temperature (i.e., Boussinesq's hypothesis).

Grid Sensitivity Analysis and Code Verification
Sensitivity of the grid to the obtained results is extensively tested. A rectangular domain is employed as shown in Figure 1a. In order to minimise the influence of boundary effects, the top, bottom, inlet, and outlet boundaries are placed sufficiently far away from the square cylinder. As shown in Figure 1b, a structured and non-uniform grid system is used in the entire computational domain. A fine mesh is placed close to the surface of the cylinder. Tests are carried out with four different grid sizes close to the cylinder viz., δ = 0.001D, 0.003D, and 0.01D. Among these, δ = 0.003D is found to be the optimum size of the grid close to the cylinder surface. The grids are evenly distributed around the cylinder surface. A non-uniform structured grid, with ∆ = 0.25D, is applied elsewhere. The grids are stretched by smooth transition using different bias factors i.e., (growth rate) × (number of divisions −1 ). Table 2  point are found to be the best choice as they predict the flow features best while incurring comparatively less computational time. Table 2. Grid sensitivity and downstream length (L D ) dependence test on drag coefficient (C D ) and Nusselt number (Nu M ) of cylinder at Re = 50 and volume fraction, ϕ = 5%.

Numerical Method
The governing equations of fluid flow and heat transfer, namely Equations (1)-(3) (for single-phase simulations), and Equations (10)-(12) (for two-phase simulations) are solved using the commercial CFD software-ANSYS FLUENT (service pack; 15.0.7) [47]. The QUICK scheme is utilized for discretizing the convection terms, while the second-order central difference scheme is used for the diffusion terms. The SIMPLE algorithm is used for pressure-velocity coupling. It is found that the SIMPLE algorithm shows good agreement between experimental and numerical results [23]. The node-based method is adopted to find gradients on the mesh surface. Then, 2D numerical simulations are carried out for an isothermally heated cylinder immersed in nanofluid when both the imposed flow and the buoyancy induced motion are in the same direction, i.e., the so-called buoyancy aiding configuration. In the Boussinesq approximation context, the suitable forms of the momentum and thermal energy equations for the mixture model are solved numerically. In all the simulations, solutions are assumed to be converged when the residual in each cell dropped to 10 −6 .

Results and Discussions
Mixed convective heat transfer from a square cylinder in a uniform flow, with nanofluids as working fluid, is investigated for the following operating parameters: number (Re cr ) for wake instability is observed to be equal to 44.7, 45, 46, 47 ± 2, and 49.5 in the work of Park and Yang [48], Saha et al. [49], Jiang et al. [50], Sohankar et al. [51], and Abdelhamid et al. 2021 [52], respectively. It can be seen in Figure 2b that the flow separates at the leading edge of the cylinder at Re = 10. Consequently, a separation bubble (steady recirculating region) consisting of twin symmetric vortices forms at the leeward side of the cylinder. This separation bubble increases in size with Re, and the flow remains steady at Re = 30 [53]. At a higher Reynolds number, i.e., Re ∼ 80, the flow separates at the leading edge of the cylinder and reattaches at a short distance downstream, thus forming a small recirculation region on the side faces of the cylinder. Overall, distinct flow patterns such as steady flow separation at the trailing edge with a separation bubble at the leeward side, flow separation at the trailing edge with vortex shedding, separation at the leading edge and reattachment on the sides of the cylinder, and separation at the leading edge with no reattachment, can be identified, as seen in subplots (i-v) of Figure 2b.  The vorticity contours are presented in Figure 2c to gain further insights into fluid flow, especially near the cylinder. At Re = 50, the vortices in the separation bubble start to separate alternately from the trailing edge of the square cylinder. As a result, the positive (i.e., anticlockwise rotation of the fluid remarked by solid lines) and negative vortices (i.e., clockwise rotation of the fluid remarked by dashed lines) grow periodically and start to shed from the cylinder and move downstream due to the Bènard-von Kàrmàn instability phenomena. This phenomenon is referred to as vortex shedding, which can be clearly seen in subplot (iii) of Figure 2c.
Time-averaged streamlines obtained by averaging the stream function during a shedding cycle are presented in Figures 3 and 4. The recirculation bubble at the leeward side of the cylinder increases with Re when Ri is maintained constant (see, subplots (i-iv) of Figure 3a,b). The opposite trend is observed with increasing Ri when Re is maintained constant (see, subplots (i-iii) of Figure 3a,b). This is most likely because buoyancy increases the velocity gradient at the cylinder surface and reduces the pressure over the surface of the cylinder that affects the size of the recirculation bubble. The density of the velocity vectors is high at Ri = 1 compared to Ri = 0 due to aiding buoyancy as seen in Figure 5b. Therefore, at any fixed value of Re, the wake length at Ri = 0 would be higher than that at Ri = 0.5 and 1. Similar findings are also reported by Sharma et al. [54] in their study on the mixed convection heat transfer from a square cylinder under thermal buoyancy at low Re values.
Further, the downstream stretching of vorticity contours increases with Re at a fixed value of Ri. The magnitude of the vorticity near the surface of the cylinder increases with an increase in Re and/or Ri, as seen in Figure 6a,b. Further, the effects of thermal buoyancy on heat transfer are discussed in Section 5.2.1. Comparison of the results between single-phase and multi-phase approaches is shown through streamline plots and vorticity contours in Figure 7I,II at different values of Re and Ri. It is observed that the wake and vorticity patterns show qualitatively similar trends for both the approaches for the considered parameters.

Time-Averaged Pressure Coefficient
The variation of time-averaged pressure coefficient (C P ) on the face of the cylinder is presented for Re = 10, 50, 80, and 100 at Ri = 0, 0.5, and 1.0 for nanoparticle volume fraction ϕ = 0% and 5% in Figure 8a-f. The maximum value of C P is observed (near the front stagnation point) at the windward side (i.e., CD) of the cylinder compared to other surfaces (i.e., AB, BC, and DA). This is consistent with the results of Gupta et al. [55] who investigated the flow and heat transfer from a semi-circular cylinder in a confined domain in the presence of buoyancy. The difference in magnitude of C P on the windward (CD) and leeward (AB) surfaces decrease with Re. As seen in subplots (a-c) of Figure 8, the pressure coefficient increases with Ri at a fixed value of Re. Still, C P has no appreciable change with an increase of ϕ value in both buoyancy and non-buoyancy cases as observed in subplots (d−f) of Figure 8. Overall, C P varies with Re and Ri regardless of the ϕ variation.     Further, the downstream stretching of vorticity contours increases with Re at a fixed value of Ri. The magnitude of the vorticity near the surface of the cylinder increases with an increase in Re and/or Ri, as seen in Figure 6a,b. Further, the effects of thermal buoyancy on heat transfer are discussed in Section 5.2.1. Comparison of the results between singlephase and multi-phase approaches is shown through streamline plots and vorticity contours in Figure 7I,II at different values of Re and Ri. It is observed that the wake and vorticity patterns show qualitatively similar trends for both the approaches for the considered parameters.

Time-Averaged Pressure Coefficient
The variation of time-averaged pressure coefficient ( ) on the face of the cylinder is presented for Re = 10, 50, 80, and 100 at Ri = 0, 0.5, and 1.0 for nanoparticle volume fraction φ = 0% and 5% in Figure 8a  (i.e., AB, BC, and DA). This is consistent with the results of Gupta et al. [55] who investigated the flow and heat transfer from a semi-circular cylinder in a confined domain in the presence of buoyancy. The difference in magnitude of on the windward (CD) and leeward (AB) surfaces decrease with Re. As seen in subplots (a-c) of Figure 8, the pressure coefficient increases with Ri at a fixed value of Re. Still, has no appreciable change with an increase of φ value in both buoyancy and non-buoyancy cases as observed in subplots (d−f) of Figure 8. Overall, varies with Re and Ri regardless of the φ variation.

Time-Averaged Drag Coefficients
The variation of time-averaged drag coefficient for the flow past a square cylinder is presented as a function of Re for different values of φ at Ri = 0, 0.5, and 1, in Figure 9a-c. The drag coefficient ( ) decreases with an increase of Re from 10 to 150, as can be seen in Figure 9a. It is well documented that at Ri = 0, in the steady flow regime, is mostly due to viscous drag and it decreases with an increase in Re. In the unsteady vortex shedding regime, as eddies are continuously shed from the cylinder, is mainly due to pressure drag. A slight increment in is observed with an increase in Re from 100 to 150, as seen in the zoomed view in subplot (a) of Figure 9. This decrement of with Re in steady flow regime is significantly high compared to unsteady flow regime. For lower values (i.e., 10 ≤ Re ≤ 50) of Re at Ri = 0, increases with an increase in the nanoparticles volume fraction (φ) in base fluid. This is noticeable at Re = 10 in the zoomed view of subplot (a) of Figure 9. It is noted that the effective viscosity ( ) increases with φ, and this leads to an increase in viscous drag force that exerts a more retarding force to the shear layers and consequently the magnitude of is increased. For higher Re values (>50), effects of an increment of φ on is found to be negligible. Similarly, at higher values of Ri (i.e., 0.5 and 1), effects of φ on are not significant (see subplots (b,c) of Figure 9). It can also be concluded from Figure 9a-c that increases with Ri in the presence of buoyancy as

Time-Averaged Drag Coefficients
The variation of time-averaged drag coefficient for the flow past a square cylinder is presented as a function of Re for different values of ϕ at Ri = 0, 0.5, and 1, in Figure 9a-c. The drag coefficient (C D ) decreases with an increase of Re from 10 to 150, as can be seen in Figure 9a. It is well documented that at Ri = 0, in the steady flow regime, C D is mostly due to viscous drag and it decreases with an increase in Re. In the unsteady vortex shedding regime, as eddies are continuously shed from the cylinder, C D is mainly due to pressure drag. A slight increment in C D is observed with an increase in Re from 100 to 150, as seen in the zoomed view in subplot (a) of Figure 9. This decrement of C D with Re in steady flow regime is significantly high compared to unsteady flow regime. For lower values (i.e., 10 ≤ Re ≤ 50) of Re at Ri = 0, C D increases with an increase in the nanoparticles volume fraction (ϕ) in base fluid. This is noticeable at Re = 10 in the zoomed view of subplot (a) of Figure 9. It is noted that the effective viscosity (µ n f ) increases with ϕ, and this leads to an increase in viscous drag force that exerts a more retarding force to the shear layers and consequently the magnitude of C D is increased. For higher Re values (>50), effects of an increment of ϕ on C D is found to be negligible. Similarly, at higher values of Ri (i.e., 0.5 and 1), effects of ϕ on C D are not significant (see subplots (b,c) of Figure 9). It can also be concluded from Figure 9a-c that C D increases with Ri in the presence of buoyancy as more forces are exerted on the cylinder. The findings are also consistent with the literature reported for the mixed convection heat transfer from a circular cylinder [56].

Isotherms
Contours of isotherms around the heated square cylinder are shown in Figures 10-12 for different values of Re, Ri, φ where a comparison between single-phase and multiphase models are also shown. Figure 10a-c shows the isotherms at Ri = 0, 0.5, and 1 when Re = 10, 50, and 100. At Ri = 0, when Re is increased from 10 to 100, the clustering of the isotherms increases around the cylinder. One can understand that at lower values of Re heat transfer occurs mainly due to diffusion process and at higher values of Re the thermal boundary layer spreads along the flow direction due to the dominance of convective transport (see, Figure 10a). Under aiding buoyancy conditions (Ri = 0.5 and 1) within the same Re range, clustering and lateral thinning of the thermal boundary layer can be seen in the wake region towards the downstream direction as Ri increases from 0.5 to 1. This clustering of the thermal boundary layer on the windward side is highest, followed by the top/bottom face and then the rear face of the cylinder and accordingly heat transfer rate changes along the surface (see, Section 5.2.2). The edging of the thermal boundary layer along the centre line and towards the downstream direction increases with increasing value of Ri and/or Re as observed in Figure 10b,c.

Isotherms
Contours of isotherms around the heated square cylinder are shown in Figures 10-12 for different values of Re, Ri, ϕ where a comparison between single-phase and multi-phase models are also shown. Figure 10a-c shows the isotherms at Ri = 0, 0.5, and 1 when Re = 10, 50, and 100. At Ri = 0, when Re is increased from 10 to 100, the clustering of the isotherms increases around the cylinder. One can understand that at lower values of Re heat transfer occurs mainly due to diffusion process and at higher values of Re the thermal boundary layer spreads along the flow direction due to the dominance of convective transport (see, Figure 10a). Under aiding buoyancy conditions (Ri = 0.5 and 1) within the same Re range, clustering and lateral thinning of the thermal boundary layer can be seen in the wake region towards the downstream direction as Ri increases from 0.5 to 1. This clustering of the thermal boundary layer on the windward side is highest, followed by the top/bottom face and then the rear face of the cylinder and accordingly heat transfer rate changes along the surface (see, Section 5.2.2). The edging of the thermal boundary layer along the centre line and towards the downstream direction increases with increasing value of Ri and/or Re as observed in Figure 10b,c.
In Figure 11a,b, a comparison between the isotherms for base fluid (ϕ = 0%) and the nanofluid (ϕ = 5%) is presented at Ri = 0 and 1. In the subplot, the upper half of the isotherms is for nanofluid while the lower half shows the results for base fluid. From the figure, it is observed that the effect of ϕ on the thickness of a thermal boundary layer is more noticeable only at a lower value of Re (see, subplot (i) of Figure 11a,b). It is known that adding nanoparticles to base fluid increases the effective viscosity (µ n f ) and the effective thermal conductivity (K n f ) of fluid. Consequently, an increase in µ n f reduces the convection effect while an increase in K n f enhances the heat transfer. Still, fluid momentum produced due to the buoyancy and inertial forces are high enough to overcome the decrement of convection induced by viscosity. Furthermore, minor changes are found on the thermal boundary layer with increment in ϕ at higher values of Ri and Re, as seen in subplots (ii,iii) of Figure 11b. A comparison of the isotherms for single and multi-phase models are presented in Figure 12a,b. A significant change in isotherm patterns can be observed at a lower value of Re only with the change in Ri (see, subplot (i) of Figure 12a,b). In Figure 11a,b, a comparison between the isotherms for base fluid (φ = 0%) and the nanofluid (φ = 5%) is presented at Ri = 0 and 1. In the subplot, the upper half of the isotherms is for nanofluid while the lower half shows the results for base fluid. From the figure, it is observed that the effect of φ on the thickness of a thermal boundary layer is more noticeable only at a lower value of Re (see, subplot (i) of Figure 11a,b). It is known that adding nanoparticles to base fluid increases the effective viscosity ( ) and the effective thermal conductivity ( ) of fluid. Consequently, an increase in reduces the convection effect while an increase in enhances the heat transfer. Still, fluid momentum produced due to the buoyancy and inertial forces are high enough to overcome the decrement of convection induced by viscosity. Furthermore, minor changes are found on the thermal boundary layer with increment in φ at higher values of Ri and Re, as seen in subplots (ii,iii) of Figure 11b. A comparison of the isotherms for single and multi-phase models are presented in Figure 12a,b. A significant change in isotherm patterns can be observed at a lower value of Re only with the change in Ri (see, subplot (i) of Figure 12a,b). In Figure 11a,b, a comparison between the isotherms for base fluid (φ = 0%) and the nanofluid (φ = 5%) is presented at Ri = 0 and 1. In the subplot, the upper half of the isotherms is for nanofluid while the lower half shows the results for base fluid. From the figure, it is observed that the effect of φ on the thickness of a thermal boundary layer is more noticeable only at a lower value of Re (see, subplot (i) of Figure 11a,b). It is known that adding nanoparticles to base fluid increases the effective viscosity ( ) and the effective thermal conductivity ( ) of fluid. Consequently, an increase in reduces the convection effect while an increase in enhances the heat transfer. Still, fluid momentum produced due to the buoyancy and inertial forces are high enough to overcome the decrement of convection induced by viscosity. Furthermore, minor changes are found on the thermal boundary layer with increment in φ at higher values of Ri and Re, as seen in subplots (ii,iii) of Figure 11b. A comparison of the isotherms for single and multi-phase models are presented in Figure 12a

Local and Mean Nusselt Number of the Cylinder
The local Nusselt number ( ) variation over the square cylinder is presented in Figure 13a-i. The effects of Re, Ri, and φ values on for the considered parameters are shown. Peak values of occur at the corners of the square cylinder due to large temperature gradients. The maximum value of is noticed at the windward side (i.e., CD) of the cylinder, which increases with Re, for given values of Ri and φ (see, subplot (a-c) of Figure 13). A significant enhancement is noticed for magnitude when Ri is increased for a given Re and φ (see, subplot (d-f) of Figure 13). Similarly, effects of the addition of nanoparticles are measurable at the windward surface (CD) of the cylinder at the given value of Ri and Re.
The mean Nusselt number ( ) variation is shown in Figure 14a-c. Since φ changes the thermophysical properties of the nanofluids, accordingly the mean Nusselt number increases with an increase in φ. A significant increment in is observed at higher values of φ and Re, which justifies the use of nanofluids. An increment of Re and Ri rises the convective heat transfer rate and fluid momentum that creates a temperature gradient in the vicinity of the cylinder resulting in enhanced . The quantitative comparison of value obtained from the single-phase and multi-phase models is presented in Table  3. It is found that the use of the MPM approach indicates a higher value of than the SPM approach for the same operating parameters. A similar finding is reported in the literature for Reynolds number in the range of 10 ≤ Re ≤ 40 [24].

Local and Mean Nusselt Number of the Cylinder
The local Nusselt number (Nu l ) variation over the square cylinder is presented in Figure 13a-i. The effects of Re, Ri, and ϕ values on Nu l for the considered parameters are shown. Peak values of Nu l occur at the corners of the square cylinder due to large temperature gradients. The maximum value of Nu l is noticed at the windward side (i.e., CD) of the cylinder, which increases with Re, for given values of Ri and ϕ (see, subplot (a-c) of Figure 13). A significant enhancement is noticed for Nu l magnitude when Ri is increased for a given Re and ϕ (see, subplot (d-f) of Figure 13). Similarly, effects of the addition of nanoparticles are measurable at the windward surface (CD) of the cylinder at the given value of Ri and Re.
The mean Nusselt number (Nu M ) variation is shown in Figure 14a-c. Since ϕ changes the thermophysical properties of the nanofluids, accordingly the mean Nusselt number increases with an increase in ϕ. A significant increment in Nu M is observed at higher values of ϕ and Re, which justifies the use of nanofluids. An increment of Re and Ri rises the convective heat transfer rate and fluid momentum that creates a temperature gradient in the vicinity of the cylinder resulting in enhanced Nu M . The quantitative comparison of Nu M value obtained from the single-phase and multi-phase models is presented in Table 3. It is found that the use of the MPM approach indicates a higher value of Nu M than the SPM approach for the same operating parameters. A similar finding is reported in the literature for Reynolds number in the range of 10 ≤ Re ≤ 40 [24].        Further, there are no significant effects on E under thermal buoyancy conditions. Moreover, higher heat transfer enhancement occurs at low Reynolds number and at high volume fraction. At the same time, the drag coefficient of the cylinder and viscosity of fluid is also greater in magnitude. Table 4 presents the Nusselt number (Nu M ) deviation between the single-phase modelling and multi-phase modelling approaches at Ri = 1.0. A percentage difference less than 1% in Nu M is observed between the two models.

Conclusions
Mixed convective heat transfer from a heated square cylinder placed in a uniform flow is studied numerically with Al2O3-water nanofluids as the working fluid using the two-phase mixture model. This model gives better consistency due to the inclusion of the multi-phase approach while considering slip velocity between nanoparticles and base

Conclusions
Mixed convective heat transfer from a heated square cylinder placed in a uniform flow is studied numerically with Al 2 O 3 -water nanofluids as the working fluid using the two-phase mixture model. This model gives better consistency due to the inclusion of the multi-phase approach while considering slip velocity between nanoparticles and base fluid. The fluid flow and heat transfer behaviour are presented through vorticity, streamlines, and thermal contours. By increasing the solid volume fraction, minor variations in the flow and thermal patterns are observed for the base fluid under thermal buoyancy condition. For a fixed Ri, the recirculation bubble increases with Re. Magnitude of vorticity increases with an increase in Re and Ri. Streamlines and vorticity contours qualitatively show similar trends for both the single-phase and multi-phase approaches. The magnitude of the pressure coefficient decreases with an increase in Re from 10 to 100. Minor variation in C p is observed with variation in nanoparticle volume fraction from 0 to 5%. A slight rise in C D is observed for nanofluids when Re varies from 10 to 50. Beyond this value, no significant change in C D is observed for any value of ϕ. A remarkable change in isotherm patterns can be observed at a lower value of Re only with the change in Ri. A significant enhancement in heat transfer is noticed for Nu l magnitude when Ri is increased for a given Re and ϕ. The local and mean Nusselt numbers increase with Reynolds number, Richardson number, and nanoparticle volume fraction. For instance, the mean Nusselt number of nanofluids at Re = 100, ϕ = 5%, and Ri = 1 is approximately 12.4% higher than that of the base fluid. A comparison of the calculated value of Nu M is made between the single-phase and multi-phase models. In the mixture model, effective conductivity and viscosity of nanofluids are found to be sensitive parameters for heat transfer calculation. The overall thermal enhancement ratio increases with ϕ and decreases with Re, and almost remains constant at a lower value of ϕ regardless of the variation in Ri.
To expand this study as future work, the effect of nanolayer thickness and nanoparticles diameter on the overall heat transfer rate can be studied for different shapes of bluff bodies at higher thermal buoyancy.

Data Availability Statement:
The data that support the findings of this study are available within the article.

Conflicts of Interest:
The authors declare no conflict of interest. Nusselt number, −∂Θ/∂n p dimensional pressure, N/m 2 P non-dimensional pressure, p/ρU 2