Numerical Study of Natural Convection for Generalized Second-Grade Fluids Confined in a Square Cavity Subjected to Horizontal Heat Flux

Numerical Study of Natural Convection for Generalized Second-Grade Fluids Confined in a Square Cavity Subjected to Horizontal Heat Flux

Hamza Daghab Mourad KaddiriSaid Raghay Mohamed Lamsaadi Hassan El Harfi 

Industrial Engineering Laboratory, Sultan Moulay Slimane University, B.P. 523, Béni-Mellal 23000, Morocco

Laboratory of Applied Mathematics and Computing, Cadi Ayyad University, B.P. 549, Marrakech 40000, Morocco

Polydisciplinary Faculty, Sultan Moulay Slimane University, B.P. 592, Béni-Mellal 23000, Morocco

Corresponding Author Email:
24 August 2020
8 December 2020
17 December 2020
Available online: 
30 April 2021
| Citation

© 2021 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (



Two-dimensional steady laminar natural convection of a viscoelastic fluid represented by generalized second-grade 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 power-law index (n) in the range 1.4 - 0.6 and nominal values of Rayleigh number (Ra) range of 103 to 105.


finite volume, generalized second-grade model, natural convection, numerical study, square cavity, viscoelastic fluids

1. Introduction

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, Rivlin-Ericksen fluids or differential type fluids proposed by Rivlin and Ericksen [1], have received special attention. As a particular case of Rivlin-Ericksen 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 second-grade fluids show only the normal stress effects. However, many realistic fluids exhibit combined effects of normal stress and shear-thinning/shear-thickening behavior. In such cases, generalized second-grade fluids, proposed by Man and Sun [5], are suitable. Because of this motivation, the model considered in the present study is a generalized second-grade 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 [7-10] have been done for viscoelastic second-grade 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 non-isothermal-fluid 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. Buoyancy-driven convection for non-Newtonian fluids in enclosures has been broadly studied. Such one-dimensional 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 non-linear 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 non-Newtonian fluids has been used in the study of Allaoui and Vasseur [12], which is the Carreau-Yasuda non-Newtonian fluid. They have deduced a semi-analytical solution for natural convection of a Carreau-Yasuda fluid in a vertical enclosure heated from the side walls, founding pseudo-rheological behavior generates a significant modification in heat transfer characteristics. Whereas 2-D flows require a numerical study to solve the non-linear governing equations. Natural convection of power law fluids inside square enclosures has been numerically studied by Turan et al. [13-17] 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. Buoyancy-driven convection for another type of non-Newtonian fluids, which are viscoplastic fluids, confined in cavities has been numerically investigated by Hassan et al. [18-20]. 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 2-D 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 Criminale-Erickson-Filbey (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 second-grade fluids confined in a square enclosure, submitted to different temperatures on vertical walls. They have used a finite difference method to solve the governing non-linear conservation and constitutive equations. Results have been obtained just in the nearly Newtonian flow (E = 0.0001-0.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 second-grade fluids, which exhibit both normal stress effects and variable viscosity. Consequently, the present study has considered a generalized second-grade 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. Geometry and Governing Equations

2.1 Geometry and Governing Equations

Figure 1. Sketch of the cavity and boundary conditions

The geometry of interest is a viscoelastic fluid-filled 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 Rivlin-Ericksen tensors A1 and A2 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 second-grade model is that cannot predict shear-thinning/shear-thickening, which is the decrease/increase of viscosity with increasing share rate. The generalized second-grade 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 power-law model, due to Ostwald-de 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{n-1}{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)


$\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{n-1}{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 extra-stress 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 2-2 n}}{\alpha^{2-n}}\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. Numerics

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, Sc is the constant part of Sϕ that is explicitly independent on ϕ, while Sp 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 Semi-Implicite 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 [29-31].

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. Results and Discussion

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 =105, Pr = 30, E= 0.4 and for the second case n = 0.6, Ra = 104, 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 = 105, 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|$ 


- 4.568350





- 4.537409





- 4.546154





- 4.530809




n = 0.6, Ra = 104, Pr = 30 and E = 0.06


- 4.368957





- 4.216745





- 4.054811





- 4.035047




Table 2. Comparison of our simulation results with previous studies for different values of E, Pr and Ra: case of imposed constant temperature difference




Present Work

Sheremet and Pop [10]

Turan et al. [13]

De Val Davis [32]


$\overline{\boldsymbol{N u}}$


$\overline{\boldsymbol{N u}}$


$\overline{\boldsymbol{N u}}$


$\overline{\boldsymbol{N u}}$

































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 second-grade 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 (103 - 105) 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 103 and 105 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 = 103 the maximum elastic number reached was 1, while at Ra = 105 it was 0.4.

Typical streamlines and isotherms illustrating elasticity effects on the flow structure and temperature patterns are presented in Figures 4a-4b for two different values of Ra (103 and 105) and various values of E. At Ra = 103 (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 pseudo-conductive regime dominates the heat transfer mechanism. A growth of elastic number shows a non-significant modification of streamlines and isotherms. An increase in Rayleigh number to Ra = 105 (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 non-linear with increasing E. Hence, the phenomenon proves that the elasticity reduces the convection process.

Additionally, Figure 5 shows the distribution of non-dimensional temperature T along the centerline at Y = 0.5 for two values of Ra (103 and 105) 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=105 than Ra = 103 (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 = 103 and (b) Ra =105

Figure 5. Variations of non-dimensional temperature along the mid-length of the cavity for different values of E and various values of Ra: Ra = 103 (left) and Ra = 105 (right)

Figure 6. Average Nusselt number $\overline{N u}$  for different values of E and various values of Ra: Ra = 103 (left) and Ra = 105 (right)

Figure 7. Stream lines (top) and isotherms (bottom) for n = 1.4, Ra = 105 and (a) E = 0, (b) E = 0.4 and (c) E = 0.6

Figure 8. Non-dimensional temperature (left) and vertical velocity (right) along the line Y = 0.5 for n = 1.4, Ra = 105 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 non-linear, mainly for shear-thinning 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 104 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 shear-thickening fluids (n>1), it was considered n = 1.4, Ra = 105 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 shear-thickening. Figure 7 presents typical streamlines and isotherms for Ra = 105, 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 non-dimensional temperature T and the vertical velocity component V along the mid-length of the cavity are presented in Figure 8; it is observed a non-significant 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 = 105 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 = 105, 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 shear-thickening 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 = 105, n = 1.4 and n=1 versus E

n = 1.4



$\overline{N u}$













n = 1



$\overline{N u}$










Figure 9. Stream lines (top) and isotherms (bottom) for n = 0.6, Ra = 104 and (a) E = 0, (b) E = 0.04 and (c) E = 0.06

Figure 10. Non-dimensional temperature (left) and vertical velocity (right) along the line Y = 0.5 for n = 0.6, Ra = 104 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 non-linear due to the fact that the elasticity incapacitates the convection process. On the other hand, Figure 10 shows the variation of non-dimensional temperature T and the vertical velocity component V along the mid-length 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 = 104 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 = 104, 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 = 104, n = 0.6 and n=1 versus E

n = 0.6



$\overline{N u}$













n = 1



$\overline{N u}$













5. Conclusion

In the present work, the problem of natural convection of Generalized Second-Grade fluids, which exhibit both normal stress and shear-thinning/shear-thickening effects, confined in a square cavity has been numerically solved using a finite volume method and SIMPLE algorithm in a non-staggered 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 thermo-hydraulic 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 (CNRST-Morocco) for providing a computing infrastructure during this work.



Rivlin-Ericksen tensors


elastic number


gravitational acceleration, m.s-2


consistency index for a power-law fluid


height or width of the enclosure, m


power law index


local Nusselt number

$\overline{N u}$

mean Nusselt number


dimensionless pressure


generalized Prandtl number


constant heat flux, W. m-2


generalized Rayleigh number


dimensionless temperature


dimensionless horizontal and vertical velocities


horizontal and vertical velocities, m.s-1


dimensionless horizontal and vertical coordinates


horizontal and vertical coordinates, m

Greek symbols

$\alpha_{l, 2}$

normal stress modulus


thermal expansion coefficient, K-1


thermal conductivity, W.m-1. K-1


elastic stress tensor


dynamic viscosity, kg. m-1.s-1


dimensionless effective viscosity of fluid


density of fluid, Kg. m-3


dimensionless stream function


diffusion coefficient


working variable



effective variable


minimum value



Dimensional variables


[1] Rivlin, R.S., Ericksen, J.L. (1955). Stress deformation relations for isotropic materials. Collected Papers of RS Rivlin, 911-1013. 

[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): 191-252.

[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): 145-152.

[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): 351-377. 

[5] Man, C.S., Sun, Q.X. (1987). On the significance of normal stress effects in the flow of glaciers. Journal of Glaciology, 33: 268-273.

[6] Shenoy, A.V., Mashelkar, R.A. (1978). Laminar natural convection heat transfer to a viscoelastic fluid. Chemical Engineering Science, 33: 769-776.

[7] Mustafa, N., Asghar, S., Hossain, M.A. (2010). Natural convection flow of second-grade fluid along a vertical heated surface with variable heat flux. Int. J. Heat Mass. Transfer, 53(25-26): 5856-5862.

[8] Prasad, V.R., Bhuvanavijaya, R., Bandaru, M. (2016). Natural convection on heat transfer flow of non-newtonian second grade fluid over horizontal circular cylinder with thermal radiation. Journal of Naval Architecture and Marine Engineering, 13(1): 63-78.

[9] Ewis, K.M. (2019). Natural convection flow of Rivlin–Ericksen fluid and heat transfer through non-Darcy medium with radiation. Advances in Mechanical Engineering, 11(8).

[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): 624-640.

[11] Lamsaadi, M., Naimi, M., Hasnaoui, M. (2005). Natural convection of non-Newtonian power-law fluids in a shallow horizontal rectangular cavity uniformly heated from below. Heat Mass Transfer, 41(3): 239-249.

[12] Alloui, Z., Vasseur, P. (2015). Natural convection of Carreau-Yasuda non-Newtonian fluids in a vertical cavity heated from the sides. International Journal of Heat and Mass Transfer, 84: 912-924. 

[13] Turan, O., Sachdeva, A., Chakraborty, N., Poole, R.J. (2011). Laminar natural convection of power-law fluids in a square enclosure with differentially heated side walls subjected to constant temperatures. Journal of Non-Newtonian Fluid Mechanics, 166(17-18): 1049-1063. 

[14] Turan, O., Sachdeva, A., Poole, R.J., Chakraborty, N. (2012). Laminar natural convection of power-law fluids in a square enclosure with differentially heated sidewalls subjected to constant wall heat flux. Journal of Heat Transfer, 134(12).

[15] Kaddiri, M., Naïmi, M., Raji, A., Hasnaoui, M. (2012). Rayleigh-Bénard convection of non-Newtonian Power-law fluids with temperature-dependent viscosity. ISRN Thermodynamics.

[16] Raisi, A. (2016). Natural convection of non-Newtonian fluids in a square cavity with a localized heat source. Strojniski Vestnik/Journal of Mechanical Engineering 62(10).

[17] Horimek, A., Noureddine, B., Benkhchiba, A., Ait-Messaoudene, N. (2017). Laminar natural convection of power-law fluid in a differentially heated inclined square cavity. Annales de Chimie-Science des Matériaux, 41(3-4): 261-281.

[18] Hassan, M.A., Pathak, M., Khan, M. (2013). Natural convection of viscoplastic fluids in a square enclosure. Journal of Heat Transfer, 135(12).

[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): 2318-2332.

[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): 1769-1787.

[21] Demir, H., Akyoldoz, F.T. (2000). Unsteady thermal convection of a non-Newtonian fluid. International Journal of Engineering Science, 38(17): 1923-1938.

[22] Demir, H. (2001). Thermal convection of viscoelastic fluid with Biot boundary conduction. Mathematics and Computers in Simulation (MATCOM), 56(3): 277-296.

[23] Walicki, E., Walicka, A. (2002). Convergent flows of molten polymers modeled by generalized second-grade fluids of power-law type. Mechanics of Composite Materials, 38(1): 89-94.

[24] Carapau, F. (2010). One-dimensional viscoelastic fluid model where viscosity and normal stress coefficients depend on the shear rate. Nonlinear Analysis: Real World Applications, 11(5): 4342-4354. 

[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): 559-573.<559::AID-FLD102>3.0.CO;2-P

[27] Van Doormaal, J.P., Raithby, G.D. (1984). Enhancement of the SIMPLE method for preciding incompressible fluid flows. Numerical Heat Transfer, 7: 147-163.

[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): 1525-1532.

[29] Rahman, M.M., Siikonen, T. (2000). An improved simple method on a collocated grid. Numerical Heat Transfer: Part B: Fundamentals, 38(2): 177-201. 

[30] Darwish, M., Sraj, I., Moukalled, F. (2007). A coupled incompressible flow solver on structured grids. Numerical Heat Transfer, Part B, 52: 353-371. 

[31] Kolmogorov, D.K., Shen, W.Z., Sørensen, N.N., Sørensen, J.N. (2015). Fully consistent SIMPLE-like algorithms on collocated grids. Numerical Heat Transfer, Part B: Fundamentals, 67(2): 101-123.

[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): 249-264.