Numerical Simulation of Free Convection in Porous Media Modeled by Darcy-Brinkman Equation at Low Grashof Numbers

Numerical Simulation of Free Convection in Porous Media Modeled by Darcy-Brinkman Equation at Low Grashof Numbers

Betina Moundjia Zermane Samah*

Laboratoire de l’Ingénierie des Procédés de l’Environnement (LIPE), Université Constantine 3, Constantine 25000, Algeria

Corresponding Author Email:
11 September 2022
6 December 2022
16 December 2022
Available online: 
30 April 2023
| Citation

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



In this paper, we have analyzed a numerical study of free convection in porous media modeled by the Darcy-Brinkman-Forchheimer equation at low Grashof numbers; in a porous square enclosure. The right wall is maintained at a hot temperature and the left wall at a cold temperature, while the horizontal walls are adiabatic. A Fortran code based on the finite difference method and the (w-Y) algorithm was used to solve the equations that govern physical phenomena. The effects of Grashof Gr, Darcy Da and Prandlt Pr numbers were studied. The results showed that increasing the dimensionless parameters (Gr, Da, Pr) improves the heat transfer rate. Moreover, the local and mean Nusselt number is larger depending on these parameters. The horizontal and vertical velocities in the cavity in the middle of the wall are presented.


free convection, porous media, Darcy model, boussinesq equation, enclosure cavity

1. Introduction

Solid material with fine geometry, also called matrix, containing pores or small cavities and which can contain one or more fluids (liquid or gas).

This definition leads to the notion of porous medium. A porous structure is said to be closed when the pores are not connected to each other and are inaccessible to water and air since they do not emerge on the surface, and open when the pores are connected to each other, forming very narrow channels. purposes. A porous material is also characterized by its specific surface, its geometric tortuosity, and the distribution of grain and pore sizes within the solid.

The study of porous media consists in analyzing materials which have a solid structure containing pores. The fluid can flow through this medium. for example, we find sand, soil, sponges, ceramic materials ... etc. [1].

The shape and size of these pores is not uniform, it differs from one medium to another [2].

Many works have been carried out to exploit this phenomenon among them we can cite:

Garoosi et al. [3] presented a model of natural and mixed convection in porous media in a square cavity in the presence of nanofluids. Munshi and Alim [4] analyzed a digital study of mixed hydromagnetic convection in a square cavity driven by a double cover.

Hossain and Wilson [5] presented an important work of the flow of natural convection in a porous medium in cavities with non-isothermal walls in the presence of a generation of heat.

Lam and Prakash [6] simulated natural convection with entropy generation in a porous cavity with heat sources the physical model generally considered consists of a square cavity maintained on its side faces at different temperatures [7-10].

After this bibliographical study we noticed that the numerical approach method is almost always the finite volume method, we decided to work with the finite differences by eliminating the term of the pressure by (w-Y) algorithm at low Grashof numbers we will use a Fortran and Comsol code to clearly figure out the influence of different parameters on the heat transfer in a porous medium.

2. Domain Geometry

In a porous square cavity, a fluid circulates between the existing pores, whose right wall is at a hot temperature (Th) and the left wall at a cold temperature (Tc), while the horizontal is adiabatic. The flow is modeled by the Darcy –Brinkman equation in porous media and the Navier-Stokes equation in the fluid region (Figure 1).

Figure 1. Domain geometry

The heat flow is given by:

At the right wall

$Q=-K \cdot A \cdot \frac{d T}{d x}$[[$x=1=h \cdot A\left(T_f-T\right)$

3. Mathematical Formulation

3.1 Flow equation in porous media

The equations governing the flow are the continuity equations, momentum equation, and energy, respectively:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$             (1)

$\begin{aligned} & \frac{1}{\varepsilon} \frac{\partial u}{\partial t}+\frac{1}{\varepsilon^2}\left[u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}\right]=-\frac{1}{\rho_0} \frac{\partial p}{\partial x}+ \frac{v}{\varepsilon}\left[\frac{\partial^2 u}{\partial x^2}+\frac{\partial u^2}{\partial y^2}\right]-\frac{v}{K} u-\frac{F c}{\sqrt{K}} u \sqrt{u^2+v^2}\end{aligned}$             (2)

$\begin{aligned} & \frac{1}{\varepsilon} \frac{\partial v}{\partial t}+\frac{1}{\varepsilon^2}\left[u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}\right]= -\frac{1}{\rho_0} \frac{\partial p}{\partial y}+\frac{v}{\varepsilon}\left[\frac{\partial^2 v}{\partial x^2}+\frac{\partial v^2}{\partial y^2}\right]-\frac{v}{K} v- \\ & \frac{F c}{\sqrt{K}} v \sqrt{u^2+v^2}+g \beta\left(T-T_0\right)\end{aligned}$             (3)

$\sigma \frac{\partial T}{\partial t}+u \frac{\partial T}{\partial x}+v \frac{\partial T}{\partial y}=\alpha\left[\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right]$             (4)

In order to make the problem as general as possible we proceed to the following change of variables:

$\begin{aligned} & X=\frac{x}{L}, Y=\frac{y}{L}, U=\frac{u L}{v}, V=\frac{v L}{v}, \\ & \theta=\frac{T-T_f}{T_{c-T_f}}, P=\frac{p L^2}{\rho v^2}, \tau=\frac{t . v}{L^2}\end{aligned}$

The introduction of the preceding dimensionless variables in the differential Eqns. (1)-(4), in addition to the two numbers Pr, Gr; results in the following dimensionless mathematical model:

$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$             (5)

$\begin{aligned} & \frac{1}{\varepsilon} \frac{\partial U}{\partial \tau}+\frac{1}{\varepsilon^2}\left[U \frac{\partial U}{\partial X}+V \frac{\partial U}{\partial Y}\right]= -\frac{\partial P}{\partial X}+\frac{1}{\varepsilon}\left[\frac{\partial^2 U}{\partial X^2}+\frac{\partial U^2}{\partial Y^2}\right]-\frac{U}{D a}- \frac{F c}{\sqrt{D a}} U \sqrt{U^2+V^2}\end{aligned}$             (6)

$\begin{aligned} & \frac{1}{\varepsilon} \frac{\partial V}{\partial \tau}+\frac{1}{\varepsilon^2}\left[U \frac{\partial V}{\partial X}+V \frac{\partial V}{\partial Y}\right]=- \frac{\partial P}{\partial Y}+\frac{1}{\varepsilon}\left[\frac{\partial^2 V}{\partial X^2}+\frac{\partial V^2}{\partial Y^2}\right]-\frac{V}{D a}- \frac{F c}{\sqrt{D a}} V \sqrt{U^2+V^2}+\operatorname{Gr} \theta\end{aligned}$             (7)

$\frac{\partial \theta}{\partial \tau}+U \frac{\partial \theta}{\partial X}+V \frac{\partial \theta}{\partial Y}=\frac{1}{\operatorname{Pr}}\left[\frac{\partial^2 \theta}{\partial X^2}+\frac{\partial^2 \theta}{\partial Y^2}\right]$             (8)

By deriving Eq. (6) and Eq. (7) with respect to Y and X respectively, then by subtracting the two equations obtained, we will introduce the vorticity defined by [11, 12]:

$\omega=\frac{\partial U}{\partial Y}-\frac{\partial V}{\partial X}$             (9)

By definition the velocity vectors are given by [15]:

$\mathrm{U}=\frac{\partial \Psi}{\partial Y}, \mathrm{~V}=-\frac{\partial \Psi}{\partial X}$             (10)

By replacing U and V by their formulas in Eq. (9) we find:

$\frac{\partial^2 \Psi}{\partial Y^2}+\frac{\partial^2 \Psi}{\partial X^2}=\omega$             (11)

3.2 Nusselt number

Local Nusselt number, is calculated by:

$\mathrm{Nu}_{\mathrm{L}}=-\frac{\partial \theta}{\partial Y} \mid X=1$             (12)

The mean Nusselt number along the hot region is calculated by the integral of Eq. (12), which translates to:

Num $=\frac{h H}{k}=\int_0^1-\left(\frac{\partial \theta}{\partial Y}\right)_{X=1} \mathrm{dX}$             (13)

where, h represents the average heat transfer coefficient.

3.3 Boundary Conditions

Solving the system of equations obtained previously requires knowledge of the boundary conditions for each independent variable.

$\mathrm{Y}=0,0 \leq \mathrm{X} \leq 1, \theta=0$ and $\mathrm{Y}=1,0 \leq \mathrm{X} \leq 1, \theta=1$

$\mathrm{X}=0,0 \leq \mathrm{Y} \leq 1, \frac{\partial \theta}{\partial X}=0$ and $\mathrm{X}=1,0 \leq \mathrm{Y} \leq 1, \frac{\partial \theta}{\partial X}=0$

The vorticity conditions are identical on all the walls, and can be noted as follows: $\omega=\omega_p$, with [13];

$\omega_p=\frac{2\left(\Psi_p-\Psi_{p+1}\right)}{\Delta \eta^2}$             (14)

Obtained from a first order discretization where p: designates the wall and ∆η the space step in the direction normal to this one.

4. Numerical Formulation

We have used the finite difference method (explicit scheme) for the discretization of the differential equations which is conditionally stable [14], a two-dimensional numerical resolution of the vorticity and energy conservation equations has been implemented to simulate the circulation of air and temperature field in the enclosure. The resolution is carried out by the "FORTRAN" language, based on the method of Gauss Seidel [15].

The formulation (φ-ω) requires simultaneous processing of the three equations: energy equation, vorticity and current function. The main steps are summarized by the following algorithm:

  • Definition of the computation domain and data reading.
  • Mesh generation.
  • Introduction of initial and boundary conditions.
  • Beginning of the time loop.
  • Solving the equation of energy and vorticity.
  • Solving the equation of the current function using the nonlinear over relaxation (NLOR) method.
  • Calculation of speeds (U, V).
  • Calculation of wall vorticity.
  • Time increment (τ + ∆τ).
  • Repeat calculations until steady state is reached.
  • End of the time loop.
  • Printing of results.
5. Results and Discussion

A code was created in FORTRAN language, where the following simulation data are introduced (Table 1).

Table 1. Simulation data





i × j

81 × 81




Darcy number



Prandtl number



step of time



Grashof number



5.1 Grashof number effect

In reason to study the effect of the Grashof number on natural convection in porous cavity, numerical simulations were made for different values of the Grashof number (Gr =1-10-30)

The geometry considered is a square cavity containing air, the right wall of the enclosure is at (Tc) and the left at a cold temperature (Tf), while the horizontal walls (lower and upper) are adiabatic. The parameters held constant in this section are: Pr = 0.71, ε = 0.9 and Da = 0.01.

We first notice in the figure that there is a movement of the fluid in a counterclockwise direction, that is to say from the hot region to the cold region.

During its journey, the fluid particle absorbs heat from the hot wall, its temperature rises by making an upward movement to give up heat to the cold wall while descending, and after it comes back to lick the cold wall and the same phenomena occur.

There are three structures: the vertical boundary layers, the horizontal boundary layers and the cell region.

With increasing Grashof numbers the boundary layers become thin and the area of the cell widens, if one continued to increase the Grashof number of the cell would widen and break up.



Figure 2. Surfaces of (a) functions of isothermal currents (b) for Da=0.01, e=0.9, Pr=0.71

This phenomenon is explained by the increase in inertia, and therefore in convection speeds. The thermal fields presented by the temperature contours for the different Grashof numbers are shown in Figure 2, we can see that at Gr=1, the isothermal lines are almost vertically aligned. and this can be explained by the dominance of conduction over convection (Figure 2).

5.2 Effect of Darcy number

To study the influence of the Darcy number on free convection, calculations will be made for different values of Darcy number (Da = 10-3, 10-2, 10-1).



Figure 3. Surfaces of (a) functions of isothermal currents (b) for Gr=10, e=0.9, Pr=0.71

Figure 3 shows the current lines and isotherms, in the porous square cavity.

At low Darcy number, Da = 10-3, fluid movement in the cavity is low, when we increase the Darcy number from 10-3 to 10-1, the viscous effects become more important and therefore the heat transfer by convection will be more significant.

5.3 Effect of Prandtl number

The study of the influence of this parameter is made by several authors, we quote [15, 16] whose showed that the number of Pr is the most important parameter in the study of heat transfer

In this section, we analyze the influence of Prandtl number in a porous square cavity.

A series of calculations were performed by scanning the entire range of the Prandtl number,). The values of the number of Prandtl chosen are: Pr = 0.1- 0.71 - 10. The parameters held constant in this section are: Gr=10, e=0.9, Da=0.01.

As Prandtl's number is increased from 0.1 to 10, the heat transfer by free convection effects become more prominent.

5.4 Profiles of velocity



Figure 4. Profiles de Vitesse (a) horizontals (b) verticals for e=0.9, Da=0.01, Pr= 0.71

The profiles of the velocity components in the middle of the enclosure for different Grashof numbers are shown in the Figure 4.

The results obtained are in good agreement with the previous results, the horizontal and vertical speed profiles show very remarkable peaks for the large values of the parameters studied.

5.5 Number of Nusselt

Figure 5 presents the variation of the local Nusselt number Nu, for different parameters.

The Nusselt value clearly indicates that, the heat transfer conduction mode is dominant for low values of Pr, Da, Gr numbers (Figure 5).

Figure 5. Local Nusselt number

5.6 Mean Nusselt number

Table 2. Mean Nusselt number


Effect of Gr

Effect of Da

Effect of Pr


1   10   30

10-3 10-2 10-1

0.1 0.71 10

Nu m

4.19 7.35   9.5

2.29 7.35 8.65

3.49 7.35 19.15

Table 2 shows the variation of average Nusselt number, for different values of porosity, Grashof number Darcy and Prandtl numbers.

It can be noted that the Nusselt number increases with increasing of these parameters.

5.7 The sensitivity of the grid

The value of the mesh steps is an important factor on which it directly depends on the precision of the numerical results. In order to minimize the influence of this mesh on the solution we carried out several simulations by comparing the values of the average Nusselt number. The results obtained are given by Table 3.

Table 3. Comparison of Nusselt number of different grids


Nu m

Error (%)










The number of Nusselt for a mesh of (101×101) is taken as reference for the calculation of the error expressed in percentage.

According to the results, we decided to work with a mesh of (81×81). Where the Error is calculated by:

Error $=\left|\frac{N u_m(101 \times 101)-N u_m(i \times m)}{N u_m(i \times m)}\right| \times 100$

6. Results of Comsol Code

Comsol is used by several authors such as [17, 18] who have simulated finite element heat transfer

6.1 Validation of Comsol code

In the first place we carried out the work of article [5] for the validation of the work. We entered their data their parameters in detail in the code of Comsol without any changes.

Figure 6. Domain geometry and boundary conditions [5]

Figure 7. Result of temperature gradients

The geometry of the study [5] is square and contains the following boundary conditions: (Figure 6).

The results obtained for the same geometries by Comsol are shown in the Figure 7:

The Conjugate Heat Transfer in porous media, Laminar Flow Multiphysics interface is used to simulate the coupling between heat transfer and water flow.

It combines the Heat Transfer in Solids and Laminar Flow interfaces. The Non-Isothermal Flow Multiphysics coupling is automatically added.

We notice that the heat transfer propagates over time in a remarkable way. View the scale of temperature we can clearly see that the transfer by conduction convection is more important in the case of systems with Cartesian coordinates.

After having validated our work in relation to the study [5] we have changed the geometry, the boundary conditions the parameters, and after carrying out all the steps, we modified the mesh, and finally we had the following results.

6.2 Velocity

The field of the velocity vector of the flow is shown in Figure 8.

Figure 8. Domain of velocity of present geometry

The flow is symmetrical by compared to the half height of the channel. Away from the entrance, the flow is developed following the horizontal direction. The developed flow is unidirectional (with a totally horizontal speed) and unidimensional because this speed depends only on the vertical coordinate. In the vertical direction, the horizontal speed is zero at the walls. Moving away from these, the speed horizontal increases rapidly and becomes uniform from a distance equal to 0.2 of each wall. The uniform speed value is 1.

We have found it useful to show this velocity vector field in order to visualize the development of the flow. Indeed, we note that, in the case where the natural convection predominates the flow develops mainly along the hot wall and propagates along the upper wall before the exit. We also notice that in the rest of the cavity the water forms low intensity recirculation zones.

6.3 Pressure

Figure 9. Domain of pressure of present geometry

The flow obtained is uniform. The horizontal speed is constant and maintains its input value. The vertical speed is everywhere zero. The pressure drop is axial but the pressure does not undergo any vertical variation. It has been found and verified that Darcy's regime is dominant.

The uniform flow of the Darcy regime is due to the low value of the Darcy number.

Thus, depending on the direction of the flow, the convective terms, the Forchheimer term and the Brinkman terms are very weak compared to the Darcy terms and the pressure gradient. When the stationary regime is established, the temporal variation is zero and the Darcy and pressure gradient terms are balanced (Figure 9).

6.4 Temperature

The hot air at the inlet undergoes rapid axial and transverse cooling. This is due to the axial convective effect of the uniform flow combined with the diffusion of heat transverse increased by the thermal conductivity of the porous medium.

The inlet temperature drops axially and approaches asymptotically from that of the walls, from the axial position of the middle of the geometric configuration (Figure 10).

Figure 10. Domain of temperature of present geometry

6.5 Isothermal contours

The isotherms are shown in Figure 11. The distribution of heat in the cavity is in accordance with the circulation of the water revealed by the is currents (Figure 12). Indeed, we note a heating of the water starting from the entry, all along the right wall to the left.

We observe that the isotherms represent an almost linear evolution starting from the cold wall inside the cavity. So the temperature distribution is linear as we observe a strong temperature gradient in the vicinity of the left isothermal faces.

Figure 11. Isotherm contours of present geometry



Figure 12. Surfaces of (a) functions of isothermal currents (b) for Gr=10, e=0.9, Da=0.01

7. Conclusion

The simulation of natural convection in a porous square cavity containing a fluid was presented the right and left walls of the cavity are at different temperatures, while the top and bottom walls are kept adiabatic.

The influence of Grashof, Darcy, and Prandtl numbers, were studied.

We have found that Nusselt number increases by increasing porosity and also Grashof, Darcy and Prandtl numbers.

The results obtained are in good agreement with the previous results, the horizontal and vertical velocity profiles show very remarkable peaks for the large values of the parameters studied.

By comparing our results with those of the works already mentioned in the introduction, we find that the influence of the dimensional parameters is more noticeable in the case of low Grashof numbers.



Average Nusselt number



local Nusselt number






Dimensionless Pressure, $p L^2 / \rho v^2$



Prandtl number, ν/α



Grashof number, g β(Th-Tc)L32






Hot wall Temperature



Cold wall Temperature



Horizontal velocity Component



Dimensionless Horizontal velocity Component



vertical velocity Component



Dimensionless vertical velocity Component



Darcy number, $k / L^2$



Forchheimer coefficient, (=1.75/$\left.\sqrt{150} \varepsilon^{3 / 2}\right)$






Acceleration of gravity



Height of the cavity



Local heat transfer coefficient


Greek symbols


Dimensionless time



thermal diffusion



dimensionless temperature



Coefficient of thermal expansion



kinematic Viscosity of fluid



fluid density



current function









basis function



[1] Su, Y., Davidson, J.H. (2015). Modeling approaches to natural convection in porous media. Springer Briefs in Applied Sciences and Technology, USA, 2015. 

[2] Nield, D.A., Bejan, A. (2006). Convection in porous media, 3: New York: Springer.

[3] Garoosi, F., Bagheri, G., Rashidi, M.M. (2015). Two phase simulation of natural convection and mixed convection of the nanofluid in a square cavity. Powder Technology, 275: 239-256.

[4] Munshi, M.J.H., Alim, M.A. (2017). Effect of hydromagnetic mixed convection double lid-driven square cavity with inside elliptic heated block. Journal of Scientific Research, 9(1): 1-11.

[5] Hossain, M.A., Wilson, M. (2002). Natural convection flow in a fluid-saturated porous medium enclosed by non-isothermal walls with heat generation. International Journal of Thermal Sciences, 41(5): 447-454.

[6] Lam, P.A.K., Prakash, K.A. (2014). A numerical study on natural convection and entropy generation in a porous enclosure with heat sources. International Journal of Heat and Mass Transfer, 69: 390-407.

[7] Chan, B.K.C., Ivey, C.M., Barry, J.M. (1970). Natural convection in enclosed porous media with rectangular boundaries. Journal of Heat Transfer, 92(1): 21-27.

[8] Walker, K.L., Homsy, G.M. (1978). Convection in a porous cavity. Journal of Fluid Mechanics, 87(3): 449-474. 

[9] Bejan, A., Khair, K.R. (1985). Heat and mass transfer by natural convection in a porous medium. International Journal of Heat and Mass Transfer, 28(5): 909-918.

[10] Trevisan, O.V., Bejan, A. (1986). Mass and heat transfer by natural convection in a vertical slot filled with porous medium. International Journal of Heat and Mass Transfer, 29(3): 403-415.

[11] Wiley, J. (1969). Applied Numerical Methods. WileyNew York.

[12] Husain, S., Mustafa, J. (2014). Effects of Prandtl number on heat transfer and fluid flow in an internally heated vertical annulus. Conference: “International Conference on Advancements in Mechanical Engineering” (ICAME-2014) At: Al-Falah University, Dhauj Faridabad (India). 

[13] Bazargan, M., Mohseni, M. (2009). Effect of turbulent Prandtl number on convective heat transfer to turbulent upflow of supercritical carbon dioxide. In Heat Transfer Summer Conference, 43574: 295-302.

[14] Steiner, T.R. (2022). High temperature steady-state experiment for computational radiative heat transfer validation using COMSOL and ANSYS. Results in Engineering, 13: 100354.

[15] Ramachandran, A.V., Zorzano, M.P., Martín-Torres, J. (2022). Numerical heat transfer study of a space environmental testing facility using COMSOL Multiphysics. Thermal Science and Engineering Progress, 29: 101205.

[16] Shi, Y., Rui, S., Xu, S., Wang, N., Wang, Y. (2022). COMSOL Modeling of Heat Transfer in SVE Process. Environments, 9(5): 58. 

[17] Labed, A., Zermane, S. (2020). Simulation of heat transfer in different geometries immersed in a solar pond using Fortran and COMSOL codes. International Journal of Heat and Technology, 38(1): 231-239.

[18] Ramachandran, A.V., Zorzano, M.P., Martín-Torres, J. (2022). Numerical heat transfer study of a space environmental testing facility using COMSOL Multiphysics. Thermal Science and Engineering Progress, 29: 101205.