Fluid Dynamics & Materials Processing

DOI: 10.32604/fdmp.2021.015422

ARTICLE

Mixed Convection in a Two-Sided Lid-Driven Square Cavity Filled with Different Types of Nanoparticles: A Comparative Study Assuming Nanoparticles with Different Shapes

Mostafa Zaydan1, Mehdi Riahi1,2,*, Fateh Mebarek-Oudina3 and Rachid Sehaqui1

1Laboratory of Mechanics, Faculty of Sciences Aïn-Chock, University Hassan II, Casablanca, Morocco
2Department of Mechanical Engineering, Royal Air School, Marrakech, 40002, Morocco
3Department of Physics, Faculty of Sciences, University 20 August 1955, Skikda, 21000, Algeria
*Corresponding Author: Mehdi Riahi. Email: mehdi_riahi@hotmail.fr
Received: 17 December 2020; Accepted: 01 March 2021

List of symbols:

1  Introduction

Fluid flow and mixed-convection heat transfer in nanofluid-filled enclosures represent a complicated flow phenomenon due to wall motion that involves forced convection and temperature difference that causes natural convection induced by buoyancy. This kind of flow is of great interest in many industrial sectors, such as the industrial cooling systems applied either for heating or cooling, solar collectors, chemical processing equipment, float glass production, drying technologies, etc. Conventional heat transfer fluids such as water and ethylene glycol, with their low thermal conductivity, limit the performance and compactness of many industrial and engineering electronic devices. Nanofluids are the ideal solution for these equipments since the thermal conductivity of the nanoparticles projected in suspension in the conventional fluid, improves the heat transfer. Different nature of nanoparticles are concerned by this topic including metals (Al, Cu), oxides (Al2O3, TiO2, and CuO), carbides (SiC), nitrides (AlN, SiN) or nonmetals (Graphite, carbon nanotubes) while the base fluid is usually a conductive fluid, such as water, ethylene glycol or methanol depending on application field.

Tiwari et al. [1] investigated numerically the heat transfer enhancement in mixed convection between two-sided lid-driven differentially heated square cavity utilizing nanofluids. The nanofluid used by these authors was Copper-Water nanofluid having Prandlt number Pr = 6.2 and solid volume fraction χ varied up 12%. The upper and the bottom walls are considered thermally insulated while the left and the right moving walls are maintained at different constant temperatures. It was found that the pattern formation and the heat transfer in the cavity can be greatly influenced by both the Richardson number and the direction of the moving walls.

Several studies have been interested by the same purpose in the literatures [2–5]. For instance, Rana et al. [2] have numerically studied the heat transfer enhancement in steady mixed convection flow along a vertical plate with heat source/sink utilizing nanofluids. They emphasized the effect of spherical and cylindrical shaped nanoparticles on the heat transfer using different water-based nanofluids containing Cu, Ag, CuO, Al2O3, and TiO2. It was proved that the use of nanotubes (cylindrical shaped nanoparticles) have higher heat transfer enhancement as compared to spherical shaped nanoparticles. Also, it was shown that Ag nanoparticles proved to have the highest cooling performance for this vertical plate problem where TiO2 nanoparticles have the lowest. These results have been attributed according to these authors to the high thermal conductivity of Ag and the low thermal conductivity of TiO2.

In the same framework but in the natural convection, Zaraki et al. [6] have performed a theoretical analysis of boundary layer heat and mass transfer of nanofluids with emphasis on the effects of size, shape and type of nanoparticles, type of base fluid and working temperature. The model adopted in their analysis takes into account the effects of Brownian motion and thermophoresis. It was concluded that the choice of material types of nanoparticles was crucial. Indeed, for some type of nanoparticles significant enhancement of the heat transfer rate was observed, while some others can significantly deteriorate the heat transfer from the surface. For example, the dispersion of 40 nm spherical zinc-oxide nanoparticles in water enhances the heat transfer while dispersion of 43 nm of spherical alumina nanoparticles in water leads to a deterioration of the heat transfer [6].

Raza et al. [7] have developed a house code to study the effect of different base fluids on the hydrodynamic and thermal characteristics of Titania nanofluids in cylindrical annulus with discrete heat source. He found that the impact of TiO2 nanofluid heat transfer is related to the basic fluid types.

In this paper, we extend the analysis performed by Tiwari and Das to investigating the heat transfer in this system when different types and shapes of nanoparticles are incorporated in different base fluids. We consider the case where the right wall moves down while the left wall moves to the top which represents the case where the shear and the buoyancy forces act in opposite directions. Using highly validated Fourth-order compact numerical method for the resolution of the nonlinear equations governing the complex physics of the system, we attempt to shed light on the effect of the nature and the shape of the nanoparticles dispersed in different types of base fluids on the heat transfer enhancement in the cavity.

This paper is organized as follows: The studied configuration and its related mathematical formulation are defined in Section 2. Section 3 is devoted to presenting details of the numerical formulation used to solve the heat transfer problem and comparison with the existing results in specific configurations. In Section 4, pertinent results are discussed while the conclusion is addressed in Section 5.

2  Mathematical Formulation

The geometry of the present problem is shown in Fig. 1. It consists of a two-dimensional square cavity with the height of L. The temperature of the right wall is considered to be maintained at high temperature of Th as the left sidewall is kept at low temperature of Tc. The horizontal walls are adiabatic and impermeable. The vertical walls are moving with the same velocities in different directions such as the right wall moves downwards with the constant velocity Vp while the left one moves upwards with –Vp. The cavity is filled with a nanofluid with different nanoparticles immersed separately in different base fluids where their pertinent thermophysical properties are given in Tabs. 1 and 2.

Figure 1: Configuration of the problem in the case of the mixed convection

Table 1: Thermophysical properties of nanoparticles

Table 2: Thermophysical properties of base fluids

The flow is considered incompressible, steady and laminar and governed by the laws of conservation of mass, momentum and energy. Taking into account these assumptions, these equations take the following form:

∂u~∂x~+∂v~∂y~=0 (1)

u~∂u~∂x~+v~∂v~∂y~=−1ρnf,0∂P~∂x~+μnfρnf,0Δ2u~ (2)

u~∂v~∂x~+v~∂v~∂y~=−1ρnf,0∂P~∂y~+[χβsρs+(1−χ)ρfβf]ρnf,0g(T~−T~0)+μnfρnf,0Δ2v~ (3)

u~∂T~∂x~+v~∂T~∂y~=αnfΔ2T~ (4)

αnf is the thermal diffusivity of the nanofluid, given by

αnf=κeff(ρCp)nf.

The density ρnf , the heat capacity (ρCp)nf , and the thermal expansion coefficient (ρβ)nf , of the nanofluid are obtained from the following equations, respectively:

ρnf=χρs+(1−χ)ρf (5)

(ρCp)nf=χ(ρCp)s+(1−χ)(ρCp)f (6)

(ρβ)nf=χ(ρβ)s+(1−χ)(ρβ)f (7)

where χ being the volume fraction of the solid particles and subscripts f, nf and s designate base fluid, nanofluid and solid nanoparticles, respectively. The effective thermal conductivity of the nanofluid is approximated by the Maxwell–Garnetts model [8] and μnf is the effective dynamic viscosity of the nanofluid given respectively by:

κnfκf=(κs+2κf)−2χ(κf−κs)(κs+2κf)+χ(κf−κs) (8)

μnf=μf(1−χ)2.5 (9)

In order to proceed to the numerical solution of the system Eqs. (1)–(4) the following dimensionless variables are introduced:

x=x~L,y=y~L,u=u~Vp,v=v~Vp,P=P~ρnfVp2,t=t~VpL,T=T~−TCTH−TC

By substitution of Eqs. (5) into Eqs. (1)–(4), the following system of equation is derived:

∂u∂x+∂v∂y=0 (10)

u∂u∂x+v∂u∂y=−∂P∂x+μnfνfρnf,01ReΔ2u (11)

u∂v∂x+v∂v∂y=−∂P∂y+[χβsρs+(1−χ)ρfβf]ρnf,0βfRi.T+μnfνfρnf,01ReΔ2v (12)

u∂T∂x+v∂T∂y=αnfαf1Pr.ReΔ2T (13)

The dimensionless parameters appeared in the problem under consideration are:

Gr=βfgH3(Th−Tc)νf2,Grashofnumber

Pr=νfαf,Prandtlnumber

Re=VpHνf,Reynoldsnumber

Ri=GrRe2,Richardsonnumber

and the boundary conditions are:

x=1,0≤y≤1,u=0,v=−1,T=1Rightwallx=0,0≤y≤1,u=0,v=1,T=0Leftwally=0,0≤x≤1,u=0,v=0,∂T∂y=0Bottomwally=1,0≤x≤1,u=0,v=0,∂T∂y=0Topwall (14)

To go from the pressure-velocity formulation (u, v, p) to the stream function-vorticity formulation (ψ, ω) where

u=∂ψ∂y,v=−∂ψ∂xandω=∂v∂x−∂u∂y (15)

the system of equations governing the dynamics of the physical problem is written in the following general form

∂2ψ∂x2+∂2ψ∂y2=−ω (16)

(∂ψ∂y)(∂Γ∂x)−(∂ψ∂x)(∂Γ∂y)=ε(∂2Γ∂x2+∂2Γ∂y2)+η(∂T∂x) (17)

where

Γ=ω{ε=μnfνfρnf,01Reη=[χβsρs+(1−χ)ρfβf]βfρnf,0RiΓ=T{ε=αnfαf1Pr.Reη=0

3  Numerical Method

For the first and second x-derivatives, the following discretizations with fourth order precision are considered:

∂Φ∂x=Φx−Δx26∂3Φ∂x3+O(Δx4) (18)

∂2Φ∂x2=Φxx−Δx212∂4Φ∂x4+O(Δx4) (19)

where Φx and Φxx are the standard second order centered discretizations in the x -direction.

Substituting the expressions (18) and (19) and their corresponding y -derivatives into the Eqs. (16) and (17) yield the following form of the kinematic equation:

ψxx+ψyy−Δx212∂4ψ∂x4−Δy212∂4ψ∂y4+O(Δx4,Δy4)=−ω (20)

while the general conservation equation becomes:

εΓxx+εΓyy−εΔx212(∂4Γ∂x4)−εΔy212(∂4Γ∂y4)+O(Δx4,Δy4)=ψyΓx−ψxΓy−Δy26Γx(∂3ψ∂y3)−Δx26ψy(∂3Γ∂x3)+Δx26Γy(∂3ψ∂x3)+Δy26ψx(∂3Γ∂y3)−ηTx+ηΔx26(∂3T∂x3)+O(Δx4,Δx2,Δy2,Δy4). (21)

Considering the obtained form of the Eq. (20), one can remark that the fourth x and y-derivatives of the stream function should be replaced by their expressions resulting from Eq. (16) where

∂4ψ∂x4=−∂2ω∂x2−∂4ψ∂x2∂y2 (22)

∂4ψ∂y4=−∂2ω∂y2−∂4ψ∂y2∂x2 (23)

And by using the standard form of second-order centered discretization, the Eqs. (22) and (23) can be written as:

∂4ψ∂x4=−ωxx−ψxxyy+O(Δx2,Δy2) (24)

∂4ψ∂y4=−ωyy−ψxxyy+O(Δx2,Δy2) (25)

Taking into account these expressions, we obtain finally the following finite difference form of the kinematic equation:

Ψxx+Ψyy=−ω−Δx212ωxx−Δy212ωyy−(Δx212+Δy212)ψxxyy+O(Δx4,Δx2,Δy2Δy4) (26)

whose solution is a solution of the stream function Eq. (16) with a fourth-order spatial precision.

The same approach is applied to obtain the fourth-order approximation of the general conservation equations along the ( Ox) and (Oy) axes. For instance, the third and the fourth x -derivative are written as:

∂3Γ∂x3=1ε∂2ψ∂x∂y∂Γ∂x+1ε∂ψ∂y∂2Γ∂x2−1ε∂2ψ∂x2∂Γ∂y−1ε∂ψ∂x∂2Γ∂x∂y−ηε∂2T∂x2−∂3Γ∂x∂y2 (27)

∂4Γ∂x4=1ε∂3ψ∂x2∂y∂Γ∂x+1ε∂2ψ∂x∂y∂2Γ∂x2+1ε∂2ψ∂x∂y∂2Γ∂x2+1ε∂ψ∂y∂3Γ∂x3−1ε∂3ψ∂x3∂Γ∂y−1ε∂2ψ∂x2∂2Γ∂x∂y−1ε∂2ψ∂x2∂2Γ∂x∂y−1ε∂ψ∂x∂3Γ∂x2∂y−ηε∂3T∂x3−∂4Γ∂x2∂y2 (28)

The final form of our fourth order compact formulation (FOCF) scheme of the equations governing the dynamics of the system is as follows:

ψxx+ψyy=−ω+A (29)

ε(1+B)Γxx+ε(1+C)Γyy=(ψy+D)Γx−(ψx+E)Γy−ηΥ1Tx+F (30)

with:

A=−Δx212ωxx−Δy212ωyy−(Δx212+Δy212)Ψxxyy,B=−1εΔx26Ψxy+1ε2Δx212ΨyΨy

C=1εΔy26Ψxy+1ε2Δy212ΨxΨx

D=(Δx212+Δy212)Ψxxy−1εΔx212ΨyΨxy+1εΔy212ΨxΨyy

E=(Δx212+Δy212)Ψxyy−1εΔx212ΨyΨxx+1εΔy212ΨxΨxy (31)

F=(Δx212+Δy212)ψyΓxyy−(Δx212+Δy212)ψxΓxxy−1εΔx26ψxxΓxy+Δy26ψyyΓxy+1ε(Δx212+Δy212)ψxψyΓxy−1ε(Δx212−Δy212)ΓxΓy−(Δx212+Δy212)Γxxyy+η.G

where

G=Υ2Υ1Δx212ψxyTx+(Υ2Υ1+1Υ12)Δx212ψyTxx−Υ2Υ1Δx212ψxxTy−(Υ2Υ1Δx212+1Υ12Δy212)ψxTxy−1Υ1(Δx212+Δy212)Txyy

And

Υ1=μnfνfρnf,01.Re;Υ2=αnfαf1Pr.Re

Note that if A, B, C, D, E and F vanishes, the Eqs. (29) and (30) are reduced to the second order formulation. These equations associated to the boundary conditions are solved numerically by the implicit scheme with alternating directions (ADI: Alternating directions Implicit) method established by Peaceman et al. [9] in 1955 which consists in advancing half a step of time following x and then the second half step following y .

Such iterative numerical algorithm requires pseudo time derivatives to be assigned to the equations presenting the physical system such as:

∂2ψ∂x2+∂2ψ∂y2=−ω+∂ψ∂t (32)

∂Γ∂t+(∂ψ∂y)(∂Γ∂x)−(∂ψ∂x)(∂Γ∂y)=ε(∂2Γ∂x2+∂2Γ∂y2)+η(∂T∂x) (33)

The solution of this system will not depend on the time derivative introduced by the numerical scheme since we are dealing with a steady flow configuration where the source of unsteadiness does not exist. Consequently, a first order (Δt) accurate which is an implicit Euler step is added to these pseudo time derivatives as well as to the non-linear terms in the general conservation equation. In this framework, the system of Eqs. (32) and (33) is rewritten as follows:

ψn+Δt∂2ψn+1∂x2+Δt∂2ψn+1∂y2=−Δtωn+ψn+1 (34)

Γn+1+Δt(∂ψn∂y)(∂Γn+1∂x)−Δt(∂ψn∂x)(∂Γn+1∂y)=εΔt(∂2Γn+1∂x2+∂2Γn+1∂y2)+ηΔt(∂Tn+1∂x)+Γn (35)

This approach has been used recently by Erturk et al. [10] where the solution of the Eqs. (34) and (35) converge always toward that of the steady Eqs. (16) and (17). For avoiding any kind of repetition in this paper, we will not present the different steps used by these authors and the reader should be referred to [10] for more details.

In order to check the convergence of the obtained results during the iterations, several residuals factors are calculated. We define the first residual expression related to the stream function and the general conservation function respectively as follows:

RES1ψ=max(|ψi−1,jn+1−2ψi,jn+1+ψi+1,jn+1Δx2+ψi,j−1n+1−2ψi,jn+1+ψi,j+1n+1Δy2+ωi,jn+1−Ai,jn+1|) (36)

RES1Γ=max(|ε(1+Bi,jn+1)(Γi−1,jn+1−2Γi,jn+1+Γi+1,jn+1Δx2)+ε(1+Ci,jn+1)(Γi,j−1n+1−2Γi,jn+1+Γi,j+1n+1Δy2)−(ψi,j+1n+1−ψi,j−1n+12Δy+Di,jn+1)(Γi+1,jn+1−Γi−1,jn+12Δx)+(ψi+1,jn+1−ψi−1,jn+12Δx+Ei,jn+1)(Γi,j+1n+1−Γi,j−1n+12Δy)+(ηΥ1Ti+1,jn+1−Ti−1,jn+12Δx−Fi,jn+1)|) (37)

The value of RES1 is an indication related to the time derivatives introduced by the numerical procedure. It would generally vanish since the solution converges to the steady state.

A second residual parameter, RES2, is defined by the following expressions presenting the difference in absolute value between two successive iteration steps for the variables stream function, vorticity and energy:

RES2ψ=max(|ψi,jn+1−ψi,jn|) (38)

RES2Γ=max(|Γi,jn+1−Γi,jn|) (39)

RES2 gives an indication of the significant digit on which the code is iterated.

The third residual parameter, RES3, is similar to RES2 except that it is normalized by the representative value of the previous time step. This makes possible to get an indication later on the maximum percentage variation of change in ψ and Γ at each iteration step. RES3 is defined as follows:

RES3ψ=max(|ψi,jn+1−ψi,jnψi,jn+1|) (40)

RES3Γ=max(|Γi,jn+1−Γi,jnΓi,jn+1|) (41)

In our calculations and for all Rayleigh numbers, we consider that convergence is obtained when both RES1ψ≤10−6 and RES1Γ≤10−6 are reached. This value is chosen to ensure an accurate precision of the solution. At these convergence levels, the second residual parameters are of the order of RES2ψ≤10−10 and RES2Γ≤10−11 which means that the stream function, vorticity and the energy equations are accurate to the precision of 9 and 10 digits respectively after the decimal point, to a grid point and even more accurate to the rest of the meshes. Similar convergence criteria is adopted for the other residuals where RES3ψ≤10−10 and RES3Γ≤10−9 . These small residues guarantee that our solutions are indeed very precise.

3.1 Boundary Condition for the Vorticity

Dimensionless boundary conditions for ψ are:

ψ=0forX=0,X=1and0≤Y≤1 (42)

ψ=0forY=0,Y=1and0≤X≤1 (43)

For vorticity, Störtkuh et al. [11] have presented an analytical asymptotic solution near the corners of cavity and using finite element bilinear shape functions and they also have presented a singularity-removed boundary condition for vorticity at the corner points as well as at the wall points.

In this study, we follow Störtkuh et al. [11] for the boundary conditions and use the following expressions for calculating vorticity values at the walls:

13Δh2[...12−412111]ψ+19[...1221214114]ω=−VΔh (44)

13Δh2[....−212.121]ψ+19[....112.1214]ω=−V2Δh (45)

and at the corner points:

where V is the velocity of the walls. In our case, this velocity equals 1 for the left vertical wall, –1 for the right vertical wall whereas 0 for the other horizontal walls.

In explicit notation, for the wall points shown in Fig. 2a, the vorticity is calculated as following:

ωb=−32Δh2(ψd+ψe+ψf)−18(2ωa+2ωc+ωd+4ωe+ωf) (46)

Figure 2: Grid points at the walls and at the corners: (a) wall points and (b) corner points

Similarly, for the corner points shown in Fig. 2b, the vorticity is calculated as the following:

ωb=−32Δh2ψf−14(2ωc+2ωe+ωf) (47)

The reader should be referred to Störtkuh et al. [11] for more details on the boundary conditions.

Nulocal=hnfLkf=−(knfkf)∂θ∂Y|Y=0 (48)

Nu¯=∫01(−(knfkf)∂θ∂Y|Y=0)dx (49)

The averaged and the local heat transfer rates at the left cold wall of the cavity are presented in terms of the averaged Nusselt numbers, Nu¯ and the local Nusselt number Nulocal defined respectively as follows.

3.2 Code Validation

A good concordance is found by comparing the results obtained by our numerical code with those of de Vahl Davis [12], Markatos et al. [13], Hadjisophcleous et al. [14], Tiwari et al. [1], El Bouihi et al. [15], Khanafer et al. [16], Barakos et al. [17], Fusegi et al. [18] and Goodarzi et al. [19]. This comparison is performed on the problem of natural convection in a differentially heated square enclosure filled with air having the Prandlt number Pr = 0.71. In Tab. 3, we report the numerical values of the average Nusselt number obtained for different values of the Rayleigh number.

Table 3: Average Nusselt number comparison with previous works

For giving more reliability to our code, we attempt to compare our numerical results with previous works concerned by mixed convection in differentially heated cavity containing air where the top wall is moving with a constant velocity and kept hot while the bottom wall is cooled. The two vertical walls are insulated. At Gr =⟬ 10 ⟭^2 and for different values of the Reynolds number the obtained results are summarized in Tab. 4. As one can notice, excellent agreement is found by comparing the results obtained by Tiwari et al. [1], Sharif [20], Waheed [21], Khanafer et al. [22], Iwatsu et al. [23] and Cheng [24].

Table 4: Comparison of the Nusselt number with previous works (C.V.M)

4  Results and Discussion

4.1 Grid Test

Computations with reasonable accuracy require an optimum choose of the mesh. For that, we start by shedding light on the influence of the size and nodes distribution on the average Nusselt number Nu¯ and the stream function |Ψmin| . We report in Tab. 5 results related to these physical concepts for an assigned value of the Richardson number Ri = 10² and the volume fraction χ = 10% with a uniform distribution of nodes. As one can notice, Nu¯ and |Ψmin| become independent to the number of nodes when the grid 104 × 104 is reached. Consequently, this latter mesh size will be adopted in the rest of our calculations.

Table 5: Values of the physical quantities and for different grid values

This grid test is performed on the framework of a mixed convection configuration using square cavity filled with air (Pr = 0.71) with a vertical thermal gradient and moving upper horizontal wall.

4.2 Heat Transfer and Pattern Formation in Ag-Water Nanofluid

Our attention is focused on the mixed convective regimes where both forced and natural convections coexist corresponding to a Richardson number range 0.1≤Ri≤10 . Typically, the natural convection is negligible, compared to the forced one, when Ri < 0.1 whereas it becomes predominant when Ri > 10. However, when Ri = 0.1 we are in presence of mixed convection regime since both the Rayleigh and the Reynolds numbers are involved in the system. In this case, the forced convection is slightly more important in comparison to the natural one, which is also present, since Ra<Re2(Ra=0.1Re2) . The mixed convective regime is also considered when Ri = 10 but in this case the natural convection is more important in comparison to the forced one since Ra>Re2(Ra=10Re2) . The pure forced convection is reached when Ri = 0 (Ra=0) while the pure natural one is reached when the Richardson number tends to infinity (Re=0) .

The configuration we considered in this study is when the left vertical wall moves constantly upward while the right one moves downwards with the same velocity. Here, the buoyancy and shear forces are opposing each other. The Grashof number is fixed at Gr=104 and variation of the Richadson number is produced by a variation of the Reynolds number.

Fig. 3 illustrates the streamlines (left part) and the isotherms (right part) for Water and Water-Ag nanofluid when gradually increasing the volume fraction (χ = 0.04, 0.2) when Ri = 0.1. As one can see from the Figs. 3a–3c and 3e–3g, the flow pattern is not influenced by the presence of the nanoparticles where only one convective roll is formed in the whole cavity. The direction of rotation of this cell is clockwise as it can be deduced from the vertical velocity profile related to this case plotted at the mid-plane (x=12) in Fig. 6a. This latter exhibits a symmetry point located at the center of the cavity, where

u(x=12,y)=−u(x=12,1−y)∀y∈[0,1]. (50)

Figure 3: Streamlines (A, C, E, G) and isotherms (B, D, F, H) for Ri=0.1 : (A) χ=0% ; (B) χ=0% ; (C) χ=8% ; (D) χ=8% ; (E) χ=16% ; (F) χ=16% ; (G) χ=20% ; (H) χ=20%

Figure 4: Streamlines (A, C, E, G) and isotherms (B, D, F, H) for Ri=1 : (A) χ=0% ; (B) χ=0% ; (C) χ=8% ; (D) χ=8% ; (E) χ=16% ; (F) χ=16% ; (G) χ=20% ; (H) χ=20%

Figure 5: Streamlines (A, C, E, G) and isotherms (B, D, F, H) for Ri=10 : (A) χ=0% ; (B) χ=0% ; (C) χ=8% ; (D) χ=8% ; (E) χ=16% ; (F) χ=16% ; (G) χ=20% ; (H) χ=20%

Figure 6: Vertical mid-plane u-velocity: (a) Ri=0.1 ; (b) Ri=1 ; (c) Ri=10

This symmetry property means that the flow, at this mid-plan, changes its direction from the upper half of the cavity 12≤y≤0 to the lower one 0≤y≤12 . It is to be noted, in addition, that the rotation’s direction of the cell affects considerably the fluid flow in the cavity especially in the corners. Here, indeed, only the left-upper and right-lower corners are occupied by the fluid in contrast to the other corners.

The isotherms related to this configuration are depicted in Figs. 3b–3d and 3f–3h where a strong convection is noticed for all the values of the nanoparticle fractions. Similarly to the streamlines, the isotherms are present in the left-upper and right-lower corners ensuring heat transfer in these regions. The temperature at the horizontal mid-plane (y=12) is constant inside the cavity except near the vertical boundaries due to the presence of two thermal boundary layers near the vertical walls as shown in the Fig. 7a. The temperature function also exhibits the centre of the cavity (x=12,y=12) as a symmetry point and it verifies the following law:

T(x,y=12)=1−T(1−x,y=12)∀x∈[0,1]/[x1,x2]. (51)

Figure 7: Horizontal mid-plane temperature profile: (a) Ri=0.1 ; (b) Ri=1 ; (c) Ri=10

The interval [x1,x2] is a part of the cavity where the temperature is constant having the value T=0.5 , an average value between the temperatures of the vertical walls. In this Richardson number, x1≈0.15 and x2=0.85 .

Decreasing the Reynolds number until Ri =1 is reached, the flow pattern and the heat transfer characteristics are changed as shown in Fig. 4. Indeed, secondary flow composed by three cells is observed for χ≤8% : two cells near the left and the right walls rotating in clockwise direction and a main cell placed at the middle of the cavity rotating in the anti-clockwise direction, see Figs. 4a–4c and 4e–4g. In addition, an increase in the nanoparticles concentration beyond χ=8% leads to a stabilization of the flow in the sense that it reduces the three cells into one cell rotating in clockwise direction. This feature is responsible of the appearance of phenomenon commonly known as flow reversal observed in several flow configurations as the Couette-Taylor system with modulated rotation of the cylinders [25–28]. Indeed, it concerns a change in the direction of the vortices and results in discontinuity in the first derivative of the critical Taylor number as well as a jump in critical wave number. This phenomenon can be clearly deduced from the velocity profile plotted at the vertical mid-plan (x=12) in Fig. 7b. Meanwhile, we recover a dynamics similar to that where Ri < 1, a configuration in which the forced convection is predominant than the natural one, and therefore, one can conclude that an increase in the nanoparticle amount leads to a change in the convection mechanism in the sense that the secondary flow becomes mainly due to forced convection. Moreover, the symmetry property given by the expression (50) for Ri = 0.1 remains also valid for Ri = 1.