© 2021 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
Twodimensional steady laminar natural convection of a viscoelastic fluid represented by generalized secondgrade fluid model in a square enclosure is studied. The cavity is submitted at its vertical sides to a uniform density of heat flux while the horizontal walls are insulated, without slipping conditions at all the solid boundaries. The governing conservation and constitutive equations with the corresponding boundary conditions are solved by finite volume method in a collocated grid system. The contributions of shear rate dependent and elastic characteristics of the viscoelastic fluid are investigated on momentum and heat transport. The effects of elastic number (E) in the range 0  1 on heat transfer and fluid motion are interpreted for a powerlaw index (n) in the range 1.4  0.6 and nominal values of Rayleigh number (Ra) range of 10^{3} to 10^{5}.
finite volume, generalized secondgrade model, natural convection, numerical study, square cavity, viscoelastic fluids
The study of viscoelastic flows is very important due to their wide range of applications in many areas such as geophysics, biomechanical engineering, chemical and petroleum industries. The viscoelastic fluids include both viscous and elastic effects, making their mathematical modeling complicated than the Newtonian fluids one. To describe the behavior of viscoelastic fluids different models have been used. Amongst these models, RivlinEricksen fluids or differential type fluids proposed by Rivlin and Ericksen [1], have received special attention. As a particular case of RivlinEricksen fluids we find second grade fluids, relevant issues concerning these fluids have been discussed in detail by Dunn and Fosdick [2], Fosdick and Rajagopal [3], and Fosdick and Rajagopal [4], respectively. It is clear that secondgrade fluids show only the normal stress effects. However, many realistic fluids exhibit combined effects of normal stress and shearthinning/shearthickening behavior. In such cases, generalized secondgrade fluids, proposed by Man and Sun [5], are suitable. Because of this motivation, the model considered in the present study is a generalized secondgrade fluid.
Thermal transfer of viscoelastic fluids has attained tremendous number of studies because of its presence in many areas. The previous investigations have analytically and numerically analyzed this phenomenon in order to interpret and show the effects of elasticity on the flow and heat transfer characteristics related to the viscoelastic flow in different geometries. In the precedent studies, it has been found that the elasticity affects heat transfer through the change of Nusselt number. Shenoy and Mashelkar [6] have analyzed natural convection of viscoelastic fluid by using the approximate integral method. They have noted that the effect of elasticity on Nusselt number depends on the value of Weissenberg number. Also, a set of studies [710] have been done for viscoelastic secondgrade fluids in different geometries. It has been noticed that the elasticity is considered as a resistance force, which trends to depreciate heat transfer within the flow and decelerate the fluid motion.
Natural convection flows arise from density variations with temperature or concentration within a nonisothermalfluid under the influence of gravity. This is frequently encountered in many industrial and practical engineering areas that include heat exchangers, nuclear reactors, geothermal systems, metallurgical processes, crystal growth and others. Buoyancydriven convection for nonNewtonian fluids in enclosures has been broadly studied. Such onedimensional flows often allow getting an analytical solution of the governing conservation and constitutive equations. In this context, Lamsaadi et al. [11] have analytically and numerically studied natural convection of power law fluids in a shallow cavity subjected to a uniform density of heat flux at its horizontal walls. The parallel flow approximation has been used to simplify the governing nonlinear differential equations. They have demonstrated that the flow characteristics are sensitive to the power law index but not to the Prandtl number, when this later is sufficiently large. Another model of nonNewtonian fluids has been used in the study of Allaoui and Vasseur [12], which is the CarreauYasuda nonNewtonian fluid. They have deduced a semianalytical solution for natural convection of a CarreauYasuda fluid in a vertical enclosure heated from the side walls, founding pseudorheological behavior generates a significant modification in heat transfer characteristics. Whereas 2D flows require a numerical study to solve the nonlinear governing equations. Natural convection of power law fluids inside square enclosures has been numerically studied by Turan et al. [1317] for various boundary conditions. It has been noted that Nusselt number and the flow intensity are increasing functions with rising Rayleigh number (Ra) and decreasing power law index (n), and heat transfer takes place by conduction for low values of Ra and large values of n. Buoyancydriven convection for another type of nonNewtonian fluids, which are viscoplastic fluids, confined in cavities has been numerically investigated by Hassan et al. [1820]. It has been shown that heat transfer enhances with increasing Casson parameter, for Casson viscoplastic fluids, for any value of Rayleigh number considered on their studies. While an increase in Bingham number generates a decrease in the fluid circulation and heat transfer, for Bingham fluids, and for each Rayleigh number there correspond a critical Bingham number for which the thermal transfer inside the cavity takes place by conduction mode. Furthermore, the studies of natural convection of viscoelastic fluids in 2D cavities have been numerically studied by Demir et al. [10, 21, 22]. Demir et al. [21, 22] have numerically investigated unsteady natural convection of viscoelastic CriminaleEricksonFilbey (CEF) fluid confined in a square cavity heated from below. Recently, Sheremet and Pop [10] have studied the flow of natural convection combined with thermal radiation for secondgrade fluids confined in a square enclosure, submitted to different temperatures on vertical walls. They have used a finite difference method to solve the governing nonlinear conservation and constitutive equations. Results have been obtained just in the nearly Newtonian flow (E = 0.00010.001). It has been found that the elasticity effect is more important for high values of Rayleigh number and radiation parameter as well as small values of Prandtl number.
From the above mentioned literature survey, one can find the lack of studies treating the natural convection for confined viscoelastic generalized secondgrade fluids, which exhibit both normal stress effects and variable viscosity. Consequently, the present study has considered a generalized secondgrade fluid confined in a square enclosure, submitted at its vertical sides to a uniform density of heat flux while the horizontal walls are insulated, through a numerical algorithm using a finite volume method that allows to obtain results for relatively high elasticity numbers which has been one of the main challenges in numerical simulations.
2.1 Geometry and Governing Equations
Figure 1. Sketch of the cavity and boundary conditions
The geometry of interest is a viscoelastic fluidfilled square cavity subjected at its vertical walls to a uniform density of heat flux, while the horizontal ones are insulated, without slipping conditions at all the solid boundaries (Figure 1).
The Cauchy stress tensor of an incompressible homogenous second grade fluid [1] is defined by:
$\Sigma=p I+\mu A_{1}+\alpha_{1} A_{2}+\alpha_{2} A_{1}^{2}$ (1)
where, –pI denotes the indeterminate part of the stress due to assumption of incompressibility, μ is the coefficient of viscosity, α_{1} and α_{2} are the material constants referred to normal stress modulus, the RivlinEricksen tensors A_{1} and A_{2} are given by:
$A_{1}=g r a d V^{\prime}+\left(g r a d V^{\prime}\right)^{T}$ (2)
$A_{2}=\frac{d A_{1}}{d t}+A_{1}\left(\operatorname{grad} V^{\prime}\right)+\left(\operatorname{gradV}^{\prime}\right)^{T} A_{1}$ (3)
Here V' stands for the velocity and d/dt is the material derivative defined as follow:
$\frac{d(.)}{d t}=\frac{\partial(.)}{\partial t}+[\operatorname{grad}(.)] V^{\prime}$ (4)
Then the material parameters must respect the following restrictions [2]:
μ ≥ 0; α_{1} ≥ 0 and α_{1} + α_{2} = 0.
A shortcoming of secondgrade model is that cannot predict shearthinning/shearthickening, which is the decrease/increase of viscosity with increasing share rate. The generalized secondgrade fluids exhibit both normal stress effects and variable viscosity. The stress tensor of such fluids is defined as (see [23, 24]):
$\sum=p I+\mu(\dot{\gamma}) A_{1}+\alpha_{1} A_{2}+\alpha_{2} A_{1}^{2}$ (5)
Here $\dot{\gamma}$ is the shear rate. The viscosity is represented by the powerlaw model, due to Ostwaldde Waele, which can be expressed as:
$\mu(\dot{\gamma})=k\left[2\left(\left(\frac{\partial U^{\prime}}{\partial X_{\prime}}\right)^{2}+\left(\frac{\partial V_{\prime}}{\partial Y_{\prime}}\right)^{2}\right)+\left(\frac{\partial U_{\prime}}{\partial Y^{\prime}}+\frac{\partial V^{\prime}}{\partial X_{\prime}}\right)^{2}\right]^{\frac{n1}{2}}$ (6)
where, k and n are the consistency and the power law indices, respectively.
On the basis of the assumptions commonly used in natural convection problems, the dimensionless governing equations for Boussinesq fluids, written in terms of dimensionless velocity vector components, (U, V), dimensionless pressure, P, and dimensionless temperature, T, in Cartesian coordinate system (X, Y), are:
 Continuity equation
$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$ (7)
 Momentum equations:
$\left(U \frac{\partial U}{\partial x}+V \frac{\partial U}{\partial Y}\right)=\frac{\partial P}{\partial X}+\operatorname{Pr}\left[\mu_{a}\left(\frac{\partial^{2} U}{\partial X^{2}}+\frac{\partial^{2} U}{\partial Y^{2}}\right)+\right.$
$\left.2 \frac{\partial \mu_{a} \partial U}{\partial x} \frac{\partial U}{\partial X}+\frac{\partial \mu_{a}}{\partial Y}\left(\frac{\partial U}{\partial Y}+\frac{\partial V}{\partial X}\right)\right]+\frac{\partial \tau_{x x}}{\partial x}+\frac{\partial \tau_{x y}}{\partial V}$ (8)
$\left(U \frac{\partial V}{\partial X}+V \frac{\partial V}{\partial Y}\right)=\frac{\partial P}{\partial Y}+\operatorname{Pr}\left[\mu_{a}\left(\frac{\partial^{2} V}{\partial X^{2}}+\frac{\partial^{2} V}{\partial Y^{2}}\right)+\right.$
$\left.2 \frac{\partial \mu_{a}}{\partial Y} \frac{\partial V}{\partial Y}+\frac{\partial \mu_{a}}{\partial X}\left(\frac{\partial U}{\partial Y}+\frac{\partial V}{\partial X}\right)+R a T\right]+\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \tau_{y y}}{\partial Y}$ (9)
 Energy equation:
$U \frac{\partial T}{\partial X}+V \frac{\partial T}{\partial Y}=\frac{\partial^{2} T}{\partial X^{2}}+\frac{\partial^{2} T}{\partial Y^{2}}$ (10)
where:
$\mu_{a}=\left[2\left(\left(\frac{\partial U}{\partial x}\right)^{2}+\left(\frac{\partial V}{\partial Y}\right)^{2}\right)+\left(\frac{\partial U}{\partial Y}+\frac{\partial V}{\partial x}\right)^{2}\right]^{\frac{n1}{2}}$ (11)
$\tau_{x x}=E\left[2 \frac{d}{d t}\left(\frac{\partial U}{\partial X}\right)+\left(\frac{\partial U}{\partial Y}\right)^{2}\left(\frac{\partial V}{\partial X}\right)^{2}\right]$ (12)
$\tau_{y y}=E\left[2 \frac{d}{d t}\left(\frac{\partial V}{\partial Y}\right)+\left(\frac{\partial V}{\partial X}\right)^{2}\left(\frac{\partial U}{\partial Y}\right)^{2}\right]$ (13)
$\tau_{x y}=E\left[\frac{d}{d t}\left(\frac{\partial U}{\partial Y}+\frac{\partial V}{\partial X}\right)+2 \frac{\partial U}{\partial X} \frac{\partial V}{\partial X}+2 \frac{\partial U}{\partial Y} \frac{\partial V}{\partial Y}\right]$ (14)
In the above equations μ_{a} is the effective viscosity, τ_{xx}, τ_{yy} and τ_{xy} are the extrastress tensor components.
The dimensionless parameters, appearing in all these equations, are the Rayleigh number $\left(R a=\frac{g \beta L^{\prime 2 n+2} q^{\prime}}{(k / \rho) \alpha^{n} \lambda}\right)$, Prandtl number $\left(P r=\frac{(k / \rho) L^{\prime 22 n}}{\alpha^{2n}}\right)$and elastic number $\left(E=\frac{\alpha_{1}}{\rho L^{2}}\right)$. The dimensionless variables used here are:
(X, Y) = (X', Y') / L', (U, V) = (U', V') / (α / L'), T = (T') / (q'L' / λ), where g represents the acceleration due to gravity, β is the thermal expansion coefficient at the reference temperature, q' stands for the constant heat flux per unit area, α is the thermal diffusivity at reference temperature, λ is the thermal conductivity at the reference temperature.
To close the above equations system, the following appropriate boundary conditions are necessary:
$U=V=\frac{\partial T}{\partial X}+1=0$, for $X=0$ and $X=1$ (15)
$U=V=\frac{\partial T}{\partial Y}=0$, for $Y=0$ and $Y=1$ (16)
The physical quantities which present the heat transfer rate by convection are the local Nusselt number, Nu, and the average Nusselt number, $\overline{N u}$ that are expressed as:
$N u=\frac{q^{\prime} L^{\prime}}{\lambda \Delta T^{\prime}}=\frac{1}{\Delta T}=\frac{1}{T(0, Y)T(1, Y)}$ (17)
$\overline{N u}=\int_{0}^{1} N u d y$ (18)
3.1 Discretization method
The above conservation and constitutive equations are solved by using a finite volume formulation in collocated grid. All of conservation equations can be written in the following general form of transport equation [25]:
$\frac{\partial}{\partial x}\left(U \phi\Gamma \frac{\partial \phi}{\partial x}\right)+\frac{\partial}{\partial y}\left(V \phi\Gamma \frac{\partial \phi}{\partial y}\right)=S_{\phi}$ (19)
Here ϕ is the working variable, г is the diffusion coefficient and S_{ϕ} is the source term, which can be linearised as:
$S_{\phi}=S_{C}+S_{P} \phi_{P}$ (20)
where, S_{c} is the constant part of S_{ϕ} that is explicitly independent on ϕ, while S_{p} is the coefficient of ϕ_{p} which is made negative to enhance the numerical stability [25].
A finite volume formulation is adopted, in collocated grid, for the special discretization. The flow domain is divided into a number of control volumes ΔV around P (Figure 2).
Figure 2. Control volume
By integrating the Eq. (19) over the control volume, the final form of discretized equations relating the variable ϕ_{p} to its neighboring grid point values can be given in every control volume as:
$A_{P} \phi_{P}=A_{W} \phi_{W}+A_{E} \phi_{E}+A_{S} \phi_{S}+A_{N} \phi_{N}+S_{\phi}$ (21)
The behavior of the numerical scheme can be improved by the discretization of source term. To compute it, we need the first gradient of τ_{xx}, τ_{yy} and τ_{xy}. For the term $\partial \tau_{\mathrm{xx}} / \partial x$, as like as the other gradients, is assuming quadratic variation of τ_{xx} along the x direction. Therefore, $\partial \tau_{\mathrm{xx}} / \partial \mathrm{x}$ is expressed as 2ax+b [26].
3.2 Solving method
The obtaining discrete system for each control volume consist of a set of linear algebraic equations, which are then easily solved by means of the line by line technique based on the tridiagonal matrix algorithm (TDMA) [25, 27].
To solve the discretized equations, an equation of pressure is clearly necessary, since it is an unknown in the momentum equations that requires the use of SemiImplicite Method for Pressure Linked Equation (SIMPLE) algorithm, in which the continuity equation is transformed to the pressure equation. To avoid checkerboard velocity and pressure distribution, cell face velocities are evaluated by Momentum Interpolation Method (MIM), firstly proposed by Rhie and Chow [28], and widely used in the literature [2931].
The solution process is reiterated until attaining the convergence criterion:
$\operatorname{MAX}\left\{\frac{\phi^{\mathrm{n}+1}\phi^{\mathrm{n}}}{\phi^{\mathrm{n}+1}}\right\} \leq 10^{7}$ (22)
where, ϕ = (U, V, P, T).
4.1 Mesh independency
To check the mesh independency of the obtaining solution, minimum stream function value, Ψ_{min}, and average Nusselt number, $\overline{N u}$, are presented for two fixed cases, in which the control parameters are given by: n=1, Ra =10^{5}, Pr = 30, E= 0.4 and for the second case n = 0.6, Ra = 10^{4}, Pr = 30 and E = 0.06. Several tests were made for four different grids. Results are tabulated in Table 1. For sufficient precision the uniform grid of 250x250 has been considered for the further study.
4.2 Validation of the numerical code
To validate the elaborated numerical code, our results, expressed in terms of minimum stream function, Ψ_{min}, and average Nusselt number, $\overline{N u}$, for both imposed temperature and flux cases, are compared with previous researches with or without elasticity effect. Hence, as can be seen from Table 2 and Figure 3 a good agreement is generally obtained.
Table 1. Minimum stream function, Ψ_{min}_{,} and average Nusselt number, $\overline{N u}$, for different meshes
n = 1, Ra = 10^{5}, Pr = 30 and E = 0.4 


$\boldsymbol{\Psi}_{\min }$ 
$\overline{\boldsymbol{N} \boldsymbol{u}}$ 
$\Delta_{\mathbf{1}}=\left\frac{\left. \left.\boldsymbol{\Psi}_{\min }\right)_{i x j}\boldsymbol{\Psi}_{\min }\right)_{250 x 250}}{ \left.\boldsymbol{\Psi}_{\min }\right)_{i x j}}\right$ 
$\Delta_{2}=\left\frac{ \left.\overline{N u}_{i x j}\overline{N u}\right)_{250 x 250}}{\overline{N u})_{i x j}}\right$ 
150x150 
 4.568350 
3.213178 
0.0048 
0.0248 
200x200 
 4.537409 
3.170324 
0.0019 
0.0116 
250x250 
 4.546154 
3.133428 
 
 
300x300 
 4.530809 
3.116613 
0.0033 
0.0053 
n = 0.6, Ra = 10^{4}, Pr = 30 and E = 0.06 

150x150 
 4.368957 
2.977868 
0.0719 
0.0386 
200x200 
 4.216745 
2.913420 
0.0384 
0.0173 
250x250 
 4.054811 
2.862851 
 
 
300x300 
 4.035047 
2.779754 
0.0048 
0.0298 
Table 2. Comparison of our simulation results with previous studies for different values of E, Pr and Ra: case of imposed constant temperature difference
Ra 
E 
Pr 
Present Work 
Sheremet and Pop [10] 
Turan et al. [13] 
De Val Davis [32] 

Ψ_{min} 
$\overline{\boldsymbol{N u}}$ 
Ψ_{min} 
$\overline{\boldsymbol{N u}}$ 
Ψ_{min} 
$\overline{\boldsymbol{N u}}$ 
Ψ_{min} 
$\overline{\boldsymbol{N u}}$ 

10^{4} 
0 
0.71 
5.073 
2.205 
 
 
 
 
5.071 
2.238 
0 
100 
5.120 
2.246 
 
 
 
2.2 
 
 

10^{5} 
0.001 
23 
11.028 
4.630 
10.8 
4.6 
 
 
 
 
Figure 3. Comparison of our simulation results (left) with those of Turan et al. [14] (right) for Pr = 100 and different values of Ra and n: case of imposed constant heat flux
4.3 Flow and heat transfer in viscoelastic fluids
The main objective of this study is to evaluate the effects of Ra, n and E on the flow and heat transfer generated by natural convection for a generalized secondgrade fluid. To do this, numerical trials are conducted for various values of the governing parameters E (0  1), n (0.6, 1, 1.4), Ra (10^{3 } 10^{5}) and Prandtl number was kept at a constant value Pr = 30.
4.3.1 Combined effects of Rayleigh and elastic numbers
In this part, the power law index was kept at a constant value, n = 1, while the Rayleigh number was varied between 10^{3} and 10^{5} and the elasticity parameter was ranged between 0 and 1. It was found that the attainable maximum elastic number, at which a converged solution was obtained, strongly depends on Rayleigh number. At Ra = 10^{3} the maximum elastic number reached was 1, while at Ra = 10^{5} it was 0.4.
Typical streamlines and isotherms illustrating elasticity effects on the flow structure and temperature patterns are presented in Figures 4a4b for two different values of Ra (10^{3} and 10^{5}) and various values of E. At Ra = 10^{3} (Figure 4a), a circular clockwise flow is formed inside the enclosure for all values of elastic number. In fact, the buoyancy forces generated due to the temperature difference between the vertical walls cause the fluid to rise along the hot wall and to descend along the cold wall. At the same time, the corresponding isotherms are nearly parallel to the isothermal walls indicating that the pseudoconductive regime dominates the heat transfer mechanism. A growth of elastic number shows a nonsignificant modification of streamlines and isotherms. An increase in Rayleigh number to Ra = 10^{5} (relatively high convection), which indicates a strengthening of buoyancy forces comparing to viscous forces, generates, in the Newtonian case (E = 0), an elongation of isotherms and cell core of streamlines along the horizontal axis, which indicates that convection dominates heat transfer. However, for a viscoelastic fluid, a more important flow resistance is generated as the value of elastic number increases. Therefore, the cell core of streamlines becomes less elongated, and the corresponding isotherms seem to be less nonlinear with increasing E. Hence, the phenomenon proves that the elasticity reduces the convection process.
Additionally, Figure 5 shows the distribution of nondimensional temperature T along the centerline at Y = 0.5 for two values of Ra (10^{3} and 10^{5}) and different values of E. It is clear to notice that, in general, the distribution of T exhibits an increase in boundary layer thickness with rising E. This increasing trend of thermal boundary layer thickness with a growth of E indicates that the effects of convection weaken with elasticity. An augmentation of the thermal boundary layer thickness gives rise to a decrease in the magnitude of heat flux at the vertical walls, which acts to decrease the average Nusselt number $\overline{N u}$ as can be shown in Figure 6, where the variation of the average Nusselt number $\overline{N u}$ with E is presented. It is clear that the effect of elasticity is more pronounced for Ra=10^{5} than Ra = 10^{3} (This behavior is consistent with results of Sheremet and Pop [10]), which means that the parameter E plays an inhibiting role against the mixing role of Ra.
Figure 4. Stream lines (right) and isotherms (left) for different values of E and various values of Ra: (a) Ra = 10^{3} and (b) Ra =10^{5}
Figure 5. Variations of nondimensional temperature along the midlength of the cavity for different values of E and various values of Ra: Ra = 10^{3} (left) and Ra = 10^{5} (right)
Figure 6. Average Nusselt number $\overline{N u}$ for different values of E and various values of Ra: Ra = 10^{3} (left) and Ra = 10^{5} (right)
Figure 7. Stream lines (top) and isotherms (bottom) for n = 1.4, Ra = 10^{5} and (a) E = 0, (b) E = 0.4 and (c) E = 0.6
Figure 8. Nondimensional temperature (left) and vertical velocity (right) along the line Y = 0.5 for n = 1.4, Ra = 10^{5} and different values of E
4.3.2 Combined effects of power law index and elastic number
In this part of the study, mutual effects of the rheological parameters have been discussed. During the variation of power law index, the governing equations become heavily nonlinear, mainly for shearthinning fluids (n<1). Therefore, the attainable maximum elastic number, at which a converged solution was obtained, strongly depends on Rayleigh number and power law index. For n = 0.6, the value of Ra was taken equal to 10^{4 }and E was ranged from 0 to 0.06. Despite this restriction to the small values of elastic number, it was noticed that interesting structural changes in both streamlines and isotherms. For shearthickening fluids (n>1), it was considered n = 1.4, Ra = 10^{5} and E was varied from 0 to 0.6.
Results presented in the previous section are given for n =1. Now let us increase the value of n. consequently, the effects of viscous forces become increasingly strong in comparison to buoyancy forces because of shearthickening. Figure 7 presents typical streamlines and isotherms for Ra = 10^{5}, n= 1.4 and different values of E. A very limited qualitative effect can be observed in the flow structure and temperature patterns by changing E. Furthermore, the distributions of nondimensional temperature T and the vertical velocity component V along the midlength of the cavity are presented in Figure 8; it is observed a nonsignificant modification in the temperature distribution, and a slight decrease in the magnitude of vertical velocity. Moreover, Table 3 shows the dependence of the flow intensity (absolute value of minimum stream function), Ψ_{min}, and the mean heat transfer (average Nusselt number), $\overline{N u}$, on E for Ra = 10^{5} and two different values of n (n = 1.4 and n = 1). It is clear from Table 3 that the increase in E causes the flow intensity, Ψ_{min}, and average Nusselt number, $\overline{N u}$, to decrease because of the flow resistance offered by elasticity. Table 3 further shows that, at n = 1.4, elasticity influences slightly the mean heat transfer and flow intensity. Here, an increase in E from 0 to 0.6 generates a diminution of Ψ_{min} and $\overline{N u}$ at 2.08% and 0.77%, respectively. Whereas in the case where n = 1 and Ra = 10^{5}, a growth of elastic number from 0 up to 0.4 leads to a decrease in Ψ_{min} and $\overline{N u}$ at 10.76 % and 12.57%, respectively. In other words, the influence of elasticity is less important for shearthickening fluids. It is probably due to the fact that this parameter does not affect the flow when the viscous forces are more important.
Table 3. Absolute value of minimum stream function, average Nusselt number for Pr = 30, Ra = 10^{5}, n = 1.4 and n=1 versus E
n = 1.4 

E 
Ψ_{min} 
$\overline{N u}$ 
0 
3.747493 
2.224913 
0.2 
3.739863 
2.221718 
0.4 
3.726862 
2.219437 
0.6 
3.669501 
2.207776 
n = 1 

E 
Ψ_{min} 
$\overline{N u}$ 
0 
5.094553 
3.584049 
0.2 
4.841439 
3.414926 
0.4 
4.546154 
3.133428 
Figure 9. Stream lines (top) and isotherms (bottom) for n = 0.6, Ra = 10^{4} and (a) E = 0, (b) E = 0.04 and (c) E = 0.06
Figure 10. Nondimensional temperature (left) and vertical velocity (right) along the line Y = 0.5 for n = 0.6, Ra = 10^{4} and different values of E
For n = 0.6, the effects of convection become stronger. As expected, the flow structure is more sensitive to the change of the elastic parameter although the variation interval of E is smaller. According to Figure 9, with the increase in E, streamlines become more curved in the vicinity of the walls, and the central cell starts to be less elongated. Also, the isotherms are less nonlinear due to the fact that the elasticity incapacitates the convection process. On the other hand, Figure 10 shows the variation of nondimensional temperature T and the vertical velocity component V along the midlength of the cavity. It is clear to notice that, an increase in elastic parameter E leads to an increase in thermal boundary layer thickness and a decrease in magnitude of the vertical velocity, which in turn indicates that the effects of convection weaken with elasticity. Moreover, Table 4 shows the dependence of Ψ_{min} and $\overline{N u}$ on E for Ra = 10^{4} and two different values of n (n = 0.6 and n = 1). It is evident from Table 4 that the flow intensity and average Nusselt number decrease with the enhancement of elastic number as a result of the flow resistance generated by elasticity. Table 4 also shows that, at n = 0.6, the mean heat transfer and flow intensity are significantly affected by the rise of elastic number. An increase in E from 0 to 0.06 leads to a reduction of Ψ_{min} and $\overline{N u}$ at 19.01% and 19.06%, respectively. This change is more important comparing to the case when n = 1 and Ra = 10^{4}, where an increase in E from 0 to 0.6 leads to a reduction of Ψ_{min} and $\overline{N u}$ at 5.81 % and 3.39%, respectively.
Table 4. Absolute value of minimum stream function, average Nusselt number for Pr = 30, Ra = 10^{4}, n = 0.6 and n=1 versus E
n = 0.6 

E 
Ψ_{min} 
$\overline{N u}$ 
0 
5.006324 
3.536906 
0.02 
4.646111 
3.319593 
0.04 
4.208706 
2.995955 
0.06 
4.054811 
2.862851 
n = 1 

E 
Ψ_{min} 
$\overline{N u}$ 
0 
3.124333 
1.977504 
0.2 
2.991091 
1.946583 
0.4 
2.979582 
1.980614 
0.6 
2.942718 
1.910477 
In the present work, the problem of natural convection of Generalized SecondGrade fluids, which exhibit both normal stress and shearthinning/shearthickening effects, confined in a square cavity has been numerically solved using a finite volume method and SIMPLE algorithm in a nonstaggered grid system. Numerical simulations have been conducted for different values of Rayleigh number, power law index and elastic number. It has been found that the flow and thermal structures as well as the thermohydraulic characteristics of the fluid are sensitive to its rheological behavior, since the flow intensity and heat transfer rate decrease with the elasticity effect, which represents a resistance force. Besides, the effect of elasticity is more pronounced for high values of Rayleigh number and low values of power law index.
The authors like to express their thankfulness to the National Center of Scientific and Technical Research (CNRSTMorocco) for providing a computing infrastructure during this work.
A_{1,2} 
RivlinEricksen tensors 
E 
elastic number 
g 
gravitational acceleration, m.s^{2} 
k 
consistency index for a powerlaw fluid 
L' 
height or width of the enclosure, m 
n 
power law index 
Nu 
local Nusselt number 
$\overline{N u}$ 
mean Nusselt number 
P 
dimensionless pressure 
Pr 
generalized Prandtl number 
q' 
constant heat flux, W. m^{2} 
Ra 
generalized Rayleigh number 
T 
dimensionless temperature 
(U,V) 
dimensionless horizontal and vertical velocities 
(U',V') 
horizontal and vertical velocities, m.s^{1} 
(X,Y) 
dimensionless horizontal and vertical coordinates 
(X',Y') 
horizontal and vertical coordinates, m 
Greek symbols 

$\alpha_{l, 2}$ 
normal stress modulus 
$\beta$ 
thermal expansion coefficient, K^{1} 
λ 
thermal conductivity, W.m^{1}. K^{1} 
τ_{ij} 
elastic stress tensor 
$\mu$ 
dynamic viscosity, kg. m^{1}.s^{1} 
$\mu_{a}$ 
dimensionless effective viscosity of fluid 
$\rho$ 
density of fluid, Kg. m^{3} 
Ψ 
dimensionless stream function 
$\Gamma$ 
diffusion coefficient 
$\phi$ 
working variable 
Subscripts 

$a$ 
effective variable 
min 
minimum value 
Superscripts 

' 
Dimensional variables 
[1] Rivlin, R.S., Ericksen, J.L. (1955). Stress deformation relations for isotropic materials. Collected Papers of RS Rivlin, 9111013. https://doi.org/10.1007/9781461224167_61
[2] Dunn, J.E., Fosdick, R.L. (1974). Thermodynamics, stability, and boundedness of fluids of complexity 2 and fluids of second grade. Archive for Rational Mechanics and Analysis, 56(3): 191252. https://doi.org/10.1007/BF00280970
[3] Fosdick, R.L., Rajagopal, K.R. (1979). Anomalous features in the model of second order fluids. Archive for Rational Mechanics and Analysis, 70(2): 145152. https://doi.org/10.1007/BF00250351
[4] Fosdick, R.L., Rajagopal, K.R. (1980). Thermodynamics and stability of fluids of third grade. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 339(1738): 351377. https://doi.org/10.1098/rspa.1980.0005
[5] Man, C.S., Sun, Q.X. (1987). On the significance of normal stress effects in the flow of glaciers. Journal of Glaciology, 33: 268273. https://doi.org/10.3189/S0022143000008832
[6] Shenoy, A.V., Mashelkar, R.A. (1978). Laminar natural convection heat transfer to a viscoelastic fluid. Chemical Engineering Science, 33: 769776. https://doi.org/10.1016/00092509(78)800566
[7] Mustafa, N., Asghar, S., Hossain, M.A. (2010). Natural convection flow of secondgrade fluid along a vertical heated surface with variable heat flux. Int. J. Heat Mass. Transfer, 53(2526): 58565862. https://doi.org/10.1016/j.ijheatmasstransfer.2010.07.060
[8] Prasad, V.R., Bhuvanavijaya, R., Bandaru, M. (2016). Natural convection on heat transfer flow of nonnewtonian second grade fluid over horizontal circular cylinder with thermal radiation. Journal of Naval Architecture and Marine Engineering, 13(1): 6378. https://doi.org/10.3329/jname.v13i1.20703
[9] Ewis, K.M. (2019). Natural convection flow of Rivlin–Ericksen fluid and heat transfer through nonDarcy medium with radiation. Advances in Mechanical Engineering, 11(8). https://doi.org/10.1177/1687814019866033
[10] Sheremet, M.A., Pop, I. (2018). Natural convection combined with thermal radiation in a square cavity filled with a viscoelastic fluid. International Journal of Numerical Methods for Heat and Fluid Flow, 28(3): 624640. https://doi.org/10.1108/HFF0220170059
[11] Lamsaadi, M., Naimi, M., Hasnaoui, M. (2005). Natural convection of nonNewtonian powerlaw fluids in a shallow horizontal rectangular cavity uniformly heated from below. Heat Mass Transfer, 41(3): 239249. https://doi.org/10.1007/s0023100405308
[12] Alloui, Z., Vasseur, P. (2015). Natural convection of CarreauYasuda nonNewtonian fluids in a vertical cavity heated from the sides. International Journal of Heat and Mass Transfer, 84: 912924. https://doi.org/10.1016/j.ijheatmasstransfer.2015.01.092
[13] Turan, O., Sachdeva, A., Chakraborty, N., Poole, R.J. (2011). Laminar natural convection of powerlaw fluids in a square enclosure with differentially heated side walls subjected to constant temperatures. Journal of NonNewtonian Fluid Mechanics, 166(1718): 10491063. https://doi.org/10.1016/j.jnnfm.2011.06.003
[14] Turan, O., Sachdeva, A., Poole, R.J., Chakraborty, N. (2012). Laminar natural convection of powerlaw fluids in a square enclosure with differentially heated sidewalls subjected to constant wall heat flux. Journal of Heat Transfer, 134(12). https://doi.org/10.1115/1.4007123
[15] Kaddiri, M., Naïmi, M., Raji, A., Hasnaoui, M. (2012). RayleighBénard convection of nonNewtonian Powerlaw fluids with temperaturedependent viscosity. ISRN Thermodynamics. https://doi.org/10.5402/2012/614712
[16] Raisi, A. (2016). Natural convection of nonNewtonian fluids in a square cavity with a localized heat source. Strojniski Vestnik/Journal of Mechanical Engineering 62(10). https://doi.org/10.5545/svjme.2015.3218
[17] Horimek, A., Noureddine, B., Benkhchiba, A., AitMessaoudene, N. (2017). Laminar natural convection of powerlaw fluid in a differentially heated inclined square cavity. Annales de ChimieScience des Matériaux, 41(34): 261281. https://doi.org/10.3166/acsm.41.261280
[18] Hassan, M.A., Pathak, M., Khan, M. (2013). Natural convection of viscoplastic fluids in a square enclosure. Journal of Heat Transfer, 135(12). https://doi.org/10.1115/1.4024896
[19] Pop, I., Mikhail, S. (2017). Free convection in a square cavity filled with a Casson fluid under the effects of thermal radiation and viscous dissipation. International Journal of Numerical Methods for Heat & Fluid Flow, 27(10): 23182332. https://doi.org/10.1108/HFF0920160352
[20] Devi, T.S., Lakshmi, C.V., Venkatadri, K., Prasad, V.R., Bég, O.A., Reddy, M.S. (2020). Simulation of unsteady natural convection flow of a Casson viscoplastic fluid in a square enclosure utilizing a MAC algorithm. Heat Transfer, 49(4): 17691787. https://doi.org/10.1002/htj.21690
[21] Demir, H., Akyoldoz, F.T. (2000). Unsteady thermal convection of a nonNewtonian fluid. International Journal of Engineering Science, 38(17): 19231938. https://doi.org/10.1016/S00207225(00)000112
[22] Demir, H. (2001). Thermal convection of viscoelastic fluid with Biot boundary conduction. Mathematics and Computers in Simulation (MATCOM), 56(3): 277296.
[23] Walicki, E., Walicka, A. (2002). Convergent flows of molten polymers modeled by generalized secondgrade fluids of powerlaw type. Mechanics of Composite Materials, 38(1): 8994. https://doi.org/10.1023/A:1014017125466
[24] Carapau, F. (2010). Onedimensional viscoelastic fluid model where viscosity and normal stress coefficients depend on the shear rate. Nonlinear Analysis: Real World Applications, 11(5): 43424354. https://doi.org/10.1016/j.nonrwa.2010.05.020
[25] Patankar, S. (1980). Heat Transfer and Fluid Flow. Hemisphere, New York.
[26] Raghay, S., Hakim, A. (2001). Numerical simulation of White–Metzner fluid in a 4: 1 contraction. International Journal for Numerical Methods in Fluids, 35(5): 559573. https://doi.org/10.1002/10970363(20010315)35:5<559::AIDFLD102>3.0.CO;2P
[27] Van Doormaal, J.P., Raithby, G.D. (1984). Enhancement of the SIMPLE method for preciding incompressible fluid flows. Numerical Heat Transfer, 7: 147163. https://doi.org/10.1080/01495728408961817
[28] Rhie, C.M., Chow, W.L. (1983). Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal, 21(11): 15251532. https://doi.org/10.2514/3.8284
[29] Rahman, M.M., Siikonen, T. (2000). An improved simple method on a collocated grid. Numerical Heat Transfer: Part B: Fundamentals, 38(2): 177201. https://doi.org/10.1080/104077900750034661
[30] Darwish, M., Sraj, I., Moukalled, F. (2007). A coupled incompressible flow solver on structured grids. Numerical Heat Transfer, Part B, 52: 353371. https://doi.org/10.1080/10407790701372785
[31] Kolmogorov, D.K., Shen, W.Z., Sørensen, N.N., Sørensen, J.N. (2015). Fully consistent SIMPLElike algorithms on collocated grids. Numerical Heat Transfer, Part B: Fundamentals, 67(2): 101123. https://doi.org/10.1080/10407790.2014.949583
[32] de Vahl Davis, G. (1983). Natural convection of air in a square cavity: A bench mark numerical solution. International Journal for Numerical Methods in Fluids, 3(3): 249264. https://doi.org/10.1002/FLD.1650030305