# Numerical simulation of combined natural convection surface radiation for large temperature gradients

Saber Hamimid, Messaoud Guellal, Madiha Bouafia

► To cite this version:

Saber Hamimid, Messaoud Guellal, Madiha Bouafia. Numerical simulation of combined natural convection surface radiation for large temperature gradients. Journal of Thermophysics and Heat Transfer, 2015, 29 (3), pp.636-645. 10.2514/1.T4437. hal-01200813

HAL Id: hal-01200813

https://hal.science/hal-01200813

Submitted on 22 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Simulation of Combined Natural Convection Surface Radiation for Large Temperature Gradients

Saber Hamimid\* and Messaoud Guellal†

Université Ferhat Abbas, Sétif 19000, Algeria

and

Madiha Bouafia‡

Université d'Évry Val-d'Essonne, 91000 Évry, France

## Nomenclature

$C_p$ = dimensionless specific heat

$F_{ij}$ = geometry view factor

$g$ = gravitational acceleration

$H$ = dimension of the enclosure

$k$ = thermal conductivity

$Nr$ = dimensionless parameter of conduction–radiation, $\sigma T_0^\frac{3}{4}H/k\Delta T$

$Nu$ = convective Nusselt number

$Pr$ = Prandtl number, $\nu/\alpha$

$p$ = dynamic pressure

$p'$ = dynamic perturbed pressure

$\bar{p}$ = mean thermodynamic pressure

$Q_r$ = dimensionless net radiative heat flux

$R_i$ = radiosity

$Ra$ = Rayleigh number, $g\beta\Delta T H^3/(\nu\alpha)$

$T$ = temperature

$t$ = dimensionless time

$U, V$ = dimensionless velocity components

$X, Y$ = dimensionless coordinates

$\alpha$ = thermal diffusivity, $k/\rho \cdot C_p$

$\beta$ = thermal expansion coefficient

$\varepsilon_b$ = temperature difference parameter (Boussinesq parameter), $\Delta T/2T_0$

$\Theta$ = dimensionless radiative temperature, $T/T_0$

$\theta$ = dimensionless temperature, $(T - T_0)/(T_h - T_c)$

$\mu$ = dimensionless dynamic viscosity

$\nu$ = kinematic viscosity

$\Pi$ = dimensionless pressure, $(p - \bar{p} + \rho_0 g y)/\rho_0(\alpha_0/H)^2$

$\rho$ = dimensionless fluid density

$\sigma$ = Stefan–Boltzmann constant

## Subscripts

$\text{avg}$ = average value

$\text{BS}$ = Boussinesq model

$\text{CP}$ = constant properties

$c$ = cold

$h$ = hot

$\text{max}$ = maximum value

$\text{mid}$ = midplane

$\text{min}$ = minimum value

$0$ = reference state

## I. Introduction

Coupled natural convection and radiation transport processes in cavities are interesting topics concerning thermal and mass transport processes involved in various types of engineering systems, such as solar energy collectors, cooling of electronic devices, multi-layered walls (facades) and double windows, thermal behavior in nuclear reactors, space technology, and processes involving high temperatures.

In most of the work reported on enclosure (coupling natural convection with surface radiation), the Boussinesq approximation has been invoked. However, for higher overheat ratios, compressibility effects need to be taken into account, since the Boussinesq approximation can become inadequate. Such an approximation leads to an accurate description of convection flows if the temperature difference in the flow domain does not exceed a few degrees. However, in many physical and industrial applications such as thermal insulation systems, chemical vapor deposition reactors [1], atmospheric flows [2,3], and combustion process [4,5], the temperature differences reach tens and hundreds of degrees; all the assumptions used to justify the Boussinesq approximation fail; and a different, more complicated modeling approach is required that accounts for realistic nonlinear fluid properties variations. The low-Mach-number (LMN) approximation suggested by Paolucci [6], and used later by Hamimid et al. [7] for pure natural convection, allows one to take into account arbitrary property variations for a general fluid.

This study illustrates a comparative investigation to analyze the physical differences between numerical simulations obtained with both the Boussinesq and LMN forms of the Navier–Stokes equations for natural convection flow coupled with surface radiation in simple geometries. The purpose of this Note is to study the combined effects of surface radiation and temperature-dependent viscosity and conductivity on the momentum and heat transfer in the enclosure, as well as to quantify how the incompressible flow assumption differs from the compressible form when solving flows driven by gravity acting upon density variations resulting from local temperature gradients.

A complete parametric study is made for different Rayleigh numbers, conductivities and viscosities of the fluid, Boussinesq parameter, and wall emissivity. From these simulations, the selected characterization parameters of the maximum convective and radiative Nusselt number, the average Nusselt number, and the normalized pressure reduction were calculated. Comparisons of these parameters were made with benchmark solutions for air with the ideal gas assumption. Additionally, we generated vertical and horizontal velocities at the midvertical plane and midhorizontal planes, respectively, streamlines, and isotherms distributions to provide a basis for further analysis.

## II. Problem Formulation

We consider a bidimensional square cavity of side $H$ (Fig. 1) filled with a variable-property Newtonian fluid of density $\rho$, molecular viscosity $\mu$, thermal conductivity $k$, and temperature $T$. The natural convection problem is described by two vertically heated walls with

<footnote>\*Assistant Professor, Laboratoire de Génie des Procédés Chimiques; sab_hamimid@yahoo.fr (Corresponding Author).</footnote>

<footnote>†Professor, Laboratoire de Génie des Procédés Chimiques; messaoud.guellal@gmail.com.</footnote>

<footnote>‡Professor, Laboratoire de Mécanique et d'Energétique d'Evry; madiha.bouafia@ufrst.univ-evry.fr.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 Flow configuration and coordinate system.</figcaption></figure>

prescribed temperatures $T_h$ and $T_c$ ($T_h > T_c$), which are the “hot” and “cold” wall temperatures, respectively. We assume that the horizontal surfaces are under the combined action of the convective and radiative contributions such as $-k \frac{\partial T}{\partial y} + q_r = 0$. All four walls have no-slip zero-mass flux boundary conditions ($V = 0$) and have the same emissivity value, $\epsilon$. A Cartesian $(x, y)$ coordinate system is selected with $y$ pointing vertical and $x$ horizontal; thus, the gravitational vector $g$ is applied on the $y$ coordinate. Initial conditions are composed of constant distributions of reference pressure and temperature, $p_0$ and $T_0$, and stationary flow ($V = 0$). The reference temperature $T_0$ is defined as $T_0 = (T_h + T_c)/2$. Using the mean temperature and the temperature difference between the hot and cold walls ($\Delta T = T_h - T_c$), we introduce the Boussinesq ratio $\epsilon_b = \Delta T/2T_0$. The inner surfaces, in contact with the fluid, are assumed to be gray, diffuse emitters and reflectors of radiation with identical emissivity.

In addition to assuming a Newtonian fluid, another assumption must be made in order to obtain the usual compressible Navier-Stokes equations. This assumption has to do with the relationship between the dynamic and bulk viscosity coefficients ($\mu$ and $\lambda$) that relate the strain rate and the stress: $\lambda + \frac{2}{3}\mu = 0$.

A low-Mach-number [8–14] flow model is used to describe the dynamics of the heated cavity. In this limit, the complete Navier-Stokes equations are expanded in powers of a small parameter $M^2$ (where $M$ is the Mach number), and the total pressure $P$ is decomposed into the mean thermodynamic pressure $\bar{p}(t)$ and the dynamic pressure $P'$.

### A. Governing Equations

Since we are interested in the case of large temperature differences, the low-Mach-number equations are adopted in order to describe such a flow. The equations are made nondimensional by the use of the nondimensional variables:

$$ \frac{x}{H}, \quad \frac{y}{H}, \quad \frac{uH}{\alpha_0}, \quad \frac{vH}{\alpha_0}, \quad \frac{t}{H^2/\alpha_0}, \quad \frac{pH^2}{\rho_0 \alpha_0^2}, \\ \frac{(T-T_0)}{(T_h-T_c)}, \quad \text{and} \\ \Pi = \frac{\pi}{\rho_0(\alpha_0/H)^2} = \frac{p-\bar{p}+\rho_0 gy}{\rho_0(\alpha_0/H)^2} $$

The thermophysical properties (density, dynamic viscosity, thermal conductivity, and thermal diffusivity) are scaled by $\rho_0$, $\mu_0$, $k_0$, and $\alpha_0$, where the subscript 0 denotes values at the reference temperature $T_0$. The influence of the temperature on the specific heat is assumed to be negligible so that $C_p/C_{p0} = 1$ [15].

The system Navier-Stokes equations are expressed in nondimensional conservative form [7] as

$$ \frac{\partial \rho}{\partial t} + \frac{\partial \rho U}{\partial X} + \frac{\partial \rho V}{\partial Y} = 0 \quad (1) $$

$$ \rho \left( \frac{\partial U}{\partial t} + U \frac{\partial U}{\partial X} + V \frac{\partial U}{\partial Y} \right) = - \frac{\partial \Pi}{\partial X} + Pr \left\{ \frac{\partial}{\partial X} \left( \mu \frac{\partial U}{\partial X} \right) + \frac{\partial}{\partial Y} \left( \mu \frac{\partial U}{\partial Y} \right) + \frac{1}{3} \nabla \mu \nabla \cdot V \right\} \quad (2) $$

$$ \rho \left( \frac{\partial V}{\partial t} + U \frac{\partial V}{\partial X} + V \frac{\partial V}{\partial Y} \right) = - \frac{\partial \Pi}{\partial Y} - Ra Pr \frac{\rho - 1}{2 \epsilon_b} \\ + Pr \left\{ \frac{\partial}{\partial X} \left( \mu \frac{\partial V}{\partial X} \right) + \frac{\partial}{\partial Y} \left( \mu \frac{\partial V}{\partial Y} \right) + \frac{1}{3} \nabla \mu \nabla \cdot V \right\} \quad (3) $$

$$ \rho c_p \left( \frac{\partial \theta}{\partial t} + U \frac{\partial \theta}{\partial X} + V \frac{\partial \theta}{\partial Y} \right) = \frac{\partial}{\partial X} \left( k \frac{\partial \theta}{\partial X} \right) + \frac{\partial}{\partial Y} \left( k \frac{\partial \theta}{\partial Y} \right) + \frac{\gamma - 1}{2 \epsilon_b \gamma} \frac{d\bar{P}}{dt} \quad (4) $$

$$ \rho = \frac{\bar{P}}{(2\epsilon_b T + 1)} \quad (5) $$

Additional equations are needed to calculate the mean thermodynamic pressure $\bar{p}(t)$ and its derivate:

$$ \frac{d\bar{P}}{dt} = 2\epsilon_b\gamma \cdot \frac{1}{\int_{\Omega} d\Omega} \cdot \int_{S} k \frac{\partial T}{\partial n} dS \quad (6) $$

$$ \bar{P} = \frac{\int_{\Omega} d\Omega}{\int_{\Omega} [1/(2\epsilon_b\theta + 1)] d\Omega} \quad (7) $$

In this study, the fluid (air) is considered as a perfect gas at reference temperature $T_o = 300$ K (we set in this case $T_h = 480$ K and $T_c = 120$ K) and the transport coefficients $\mu(T)$ and $k(T)$ are given by Sutherland's law:

$$ \mu = (2\epsilon_b T + 1)^{3/2} \frac{(1 + S_\mu/T_0)}{2\epsilon_b T + 1 + S_\mu/T_0} \quad (8) $$

Where $S_\mu = 110.4$ [16]. The conductivity is given by

$$ k(T) = \frac{\mu(T)}{Pr} \quad (9) $$

### B. Radiative Analysis

When surface radiation transfer is considered in the analysis, the temperature boundary condition at the top and bottom adiabatic walls of the square cavity must be changed. Thus, one assumes that the solid surfaces are in thermal equilibrium under the combined action of the convective and radiative contributions, which give

$$ \left. \frac{\partial \theta}{\partial Y} \right|_{Y=0.1} - NrQ_r = 0 \quad (10) $$

where $Nr = \sigma T_0^4 H / k \Delta T$ is the dimensionless parameter of conduction-radiation, and $Q_r = q_r / \sigma T_0^4$ is the dimensionless net radiative heat flux on the corresponding insulated wall expressed as

$$ Q_{r,i} = R_i - \sum_{j=1}^{N} R_j F_{i-j} \quad (11) $$

where $R_i$ is the dimensionless radiosity of surface $A_i$, obtained by resolving the following system:

$$ \sum_{j=1}^{N} (\delta_{ij} - (1 - \epsilon_i) F_{i-j}) R_j = \epsilon_i \Theta_i^4 \quad (12) $$
where the dimensionless radiative temperature $\Theta_i$ is given by

$$ \Theta_i = \frac{T_i}{T_0} = [(T_h - T_c)\theta_i + T_0]/T_0 = \theta_i \frac{\Delta T}{T_0} + 1 = \frac{\theta_i}{\theta_0} + 1 \quad (13) $$

### C. Heat Transfer Parameters

The nondimensional heat transfer rate in terms of convective and radiative Nusselt numbers, $Nu_c$ and $Nu_r$, on the heated vertical surface are given by

$$ Nu_c = - \left. \frac{\partial \theta}{\partial X} \right|_{Y=0.1} \quad (14) $$

$$ Nu_r = N_r Q_r |_{Y=0.1} \quad (15) $$

The average convective Nusselt number was calculated by integrating the temperature gradient over the vertical wall as

$$ Nu_{c\_avg} = \frac{1}{A} \int_{0}^{A} -\frac{\partial \theta}{\partial X} dX \quad (16) $$

The average radiative Nusselt number is obtained by integrating the net radiative flux

$$ Nu_{r\_avg} = N_r \frac{1}{A} \int_{0}^{A} Q_r dX \quad (17) $$

The total average Nusselt number is calculated by summing the average values of convective and radiative Nusselt numbers [17]:

$$ Nu_{avg} = \frac{1}{A} \int_{0}^{A} \left( -\left. \frac{\partial \theta}{\partial X} \right|_{0,Y} + N_r Q_r(0, Y) \right) dY \quad (18) $$

### D. Numerical Modeling

The numerical solution of the governing differential equations for the velocity, pressure, and temperature fields is obtained by using a finite volume technique. A power scheme was also used in approximating advection-diffusion terms. The SIMPLER algorithm, for which the details can be found in work by Patankar [18], with a staggered grid is employed to solve the coupling between pressure and velocity. The radiosities of the elemental wall surfaces are expressed as a function of elemental wall surface temperature, emissivity, and the view factors. The radiosity $R_i$ and temperature $\theta_i$ are connected by a matrix of the type

$$ [A_{i,j}]\{R_i\} = \{\sigma\theta_i^4\} \quad (19) $$

The inverse of the matrix $[A_{i,j}]$ is determined (only once) by the Gauss elimination method. The coefficients of $[A]$ are constants and depend only on the emissivity and view factors.

In a two-dimensional formulation, the view factors are analytic [19]:

$$ F_{ij} = \frac{-1}{2(x_2 - x_1)} \left[ \sqrt{x_2^2 + y^2} \Big|_{y_1}^{y_2} - \sqrt{x_2^2 + y^2} \Big|_{y_1}^{y_2} \right] \quad (20) $$

$$ F_{i-k} = -\frac{1}{2(x_2 - x_1)} \left[ \sqrt{(x_2-x)^2 + H^2} \Big|_{x=x_1}^{x=x_2} - \sqrt{(x_1-x)^2 + H^2} \Big|_{x=x_1}^{x=x_2} \right] \quad (21) $$

For the most computations reported in this study, 242 × 242 grid points were chosen to optimize the relation between the accuracy required and the computing time; nevertheless, in the case of variable properties, a 322 × 322 grid was used. The grid is nonuniform and very fine in the *near-wall* region.

### III. Verification

The validity of the solution procedure described previously is demonstrated in this section by presenting a solution in the classical case of natural convection and surface thermal radiation in a Boussinesq approximation. A good agreement was found between the present calculations and those reported in the work of [20] (Table 1).

### IV. Results and Discussion

In this section, the low-Mach-number code is applied to analyze the effect of the Boussinesq parameter $\epsilon_b$, temperature-dependent fluid properties, and surface radiation on heat transfer and fluid flow in steady-state conditions.

The fixed parameters of simulation are $Pr = 0.71$, $T_0 = 300$ K, $A = 1$, and $\epsilon$. However, to examine the effect of the Boussinesq parameter and emissivity, we consider different values of $\epsilon_b$ (0...0, 6) and $\epsilon$ (0...1).

When the radiation exchange is taken into account, the characteristic dimension of the cavity may be calculated for a specific Rayleigh number ($H = 0.097$ in the case $Ra = 10^6$). For the Boussinesq solution, the difference temperature $\Delta T$ is fixed to 10 K. Substituting values of $H$ and $\Delta T$ into the expression of $Nr$ gives

$$ Nr = \frac{\sigma T_0^4}{k \Delta T} \left( \frac{\nu \alpha}{g \beta \Delta T} \right)^{1/3} Ra^{1/3} $$

#### A. Effect of Boussinesq Parameter $\epsilon_b$

Figures 2–10 allow comparison of the results obtained from Boussinesq, weak non-Boussinesq, and strong non-Boussinesq cases

<table><caption>Table 1 Nusselt number at active walls obtained with $T_0 = 293.5$ K and $\Delta t = 10$ K<sup>a,b</sup></caption><thead><tr><th rowspan="3">Ra</th><th rowspan="3">H</th><th rowspan="3">ε</th><th colspan="6">Wang et al. [20]</th><th colspan="6">Present work</th></tr><tr><th colspan="3">Hot wall</th><th colspan="3">Cold wall</th><th colspan="3">Hot wall</th><th colspan="3">Cold wall</th></tr><tr><th>Nu<sub>c</sub></th><th>Nu<sub>r</sub></th><th>Nu<sub>t</sub></th><th>Nu<sub>c</sub></th><th>Nu<sub>r</sub></th><th>Nu<sub>t</sub></th><th>Nu<sub>c</sub></th><th>Nu<sub>r</sub></th><th>Nu<sub>t</sub></th><th>Nu<sub>c</sub></th><th>Nu<sub>r</sub></th><th>Nu<sub>t</sub></th></tr></thead><tbody><tr><td>10<sup>4</sup></td><td>0.021</td><td>0</td><td>2.246</td><td>0</td><td>2.246</td><td>2.246</td><td>0</td><td>2.246</td><td>2.246</td><td>0</td><td>2.246</td><td>2.246</td><td>0</td><td>2.246</td></tr><tr><td>10<sup>4</sup></td><td>0.021</td><td>0.2</td><td>2.260</td><td>0.507</td><td>2.767</td><td>2.268</td><td>0.499</td><td>2.767</td><td>2.262</td><td>0.507</td><td>2.769</td><td>2.271</td><td>0.498</td><td>2.769</td></tr><tr><td>10<sup>4</sup></td><td>0.021</td><td>0.8</td><td>2.249</td><td>2.401</td><td>4.650</td><td>2.278</td><td>2.372</td><td>4.650</td><td>2.255</td><td>2.401</td><td>4.656</td><td>2.284</td><td>2.371</td><td>4.656</td></tr><tr><td>10<sup>5</sup></td><td>0.045</td><td>0</td><td>4.540</td><td>0</td><td>4.540</td><td>4.540</td><td>0</td><td>4.540</td><td>4.532</td><td>0</td><td>4.532</td><td>4.532</td><td>0</td><td>4.532</td></tr><tr><td>10<sup>5</sup></td><td>0.045</td><td>0.2</td><td>4.394</td><td>1.090</td><td>5.484</td><td>4.411</td><td>1.073</td><td>5.484</td><td>4.398</td><td>1.090</td><td>5.489</td><td>4.417</td><td>1.072</td><td>5.489</td></tr><tr><td>10<sup>5</sup></td><td>0.045</td><td>0.8</td><td>4.189</td><td>5.196</td><td>9.385</td><td>4.247</td><td>5.137</td><td>9.384</td><td>4.200</td><td>5.196</td><td>9.397</td><td>4.261</td><td>5.136</td><td>9.397</td></tr><tr><td>10<sup>6</sup></td><td>0.097</td><td>0</td><td>8.852</td><td>0</td><td>8.852</td><td>8.852</td><td>0</td><td>8.852</td><td>8.863</td><td>0</td><td>8.863</td><td>8.863</td><td>0</td><td>8.863</td></tr><tr><td>10<sup>6</sup></td><td>0.097</td><td>0.2</td><td>8.381</td><td>2.355</td><td>10.736</td><td>8.417</td><td>2.319</td><td>10.736</td><td>8.379</td><td>2.355</td><td>10.734</td><td>8.416</td><td>2.318</td><td>10.734</td></tr><tr><td>10<sup>6</sup></td><td>0.097</td><td>0.8</td><td>7.815</td><td>11.265</td><td>19.080</td><td>7.930</td><td>11.150</td><td>19.078</td><td>7.861</td><td>11.265</td><td>19.126</td><td>7.971</td><td>11.151</td><td>19.126</td></tr></tbody></table>

<sup>a</sup>$Nu_t = Nu_c + Nu_r$.

<sup>b</sup>Comparison with values published in [20].
<figure><img src="image_3.png" /><figcaption>Fig. 2 Streamlines (top) and temperature contours (bottom) at Ra = 10⁶ and ε = 0.2. From left to right, Boussinesq model (BS), εb = 0.017, and εb = 0.6.</figcaption></figure>

at Ra = 10⁶ in the presence of radiation with a wall emissivity of ε = 0.2. The properties are assumed to be variable.

From a general point of view, it can be seen that, for a low temperature difference (ε = 0.017), the solutions for compressible and incompressible models did not differ significantly. Isotherms and streamlines (Fig. 2) indicate a symmetric structure with respect to the cavity center with the presence of two secondary flows symmetrically arranged at the cavity center. At the horizontal walls, the isotherms show the existence of temperature gradients due to radiative exchange.

For a large temperature difference (ε<sub>b</sub> = 0.6), the symmetric nature of the flow is broken. A rapid look in Fig. 2 allows us to note that the isotherms' slopes near the lower wall are different from those of the upper wall; moreover, the secondary rolls inside the primary flow are no longer symmetric.

Considering the convective heat transfer along the hot and cold walls, there is a remarkable change in their profiles according to the

parameter ε<sub>b</sub> (Fig. 3). This behavior is especially observed at the bottom, where the Nusselt number for ε<sub>b</sub> = 0.6 is larger than for ε<sub>b</sub> = 0.017 on the hot wall. An opposed comportment is noticed on the cold wall.

Profiles of the local radiative Nusselt number on the vertical walls are depicted on Fig. 4. The quasi-linear character of these profiles reflects a uniform temperature on each side. Note, however, that the difference between Nu<sub>r</sub>(h) and Nu<sub>c</sub>(c) is even more important, as the parameter ε<sub>b</sub> is small. This remark is also observed on the horizontal walls (Fig. 5). The shape of these profiles reflects the temperature evolution on the top and bottom walls, which results from a balance between convective and radiative exchanges. Temperature profiles on the horizontal walls are plotted in Fig. 6. With increasing ε<sub>b</sub>, the bottom wall is cooling while the top wall is heating.

Horizontal profiles of air temperature at the midheight of the cavity are displayed in Fig. 7. One can see that, for ε<sub>b</sub> = 0.6, the thermal boundary layer gets thicker near the hot wall and thinner near the cold

<figure><img src="image_4.png" /><figcaption>Fig. 3 Distributions of local Nusselt numbers on the hot and cold walls at Ra = 10⁶ and ε = 0.2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Distribution of the local radiative Nusselt number on vertical walls at Ra = 10⁶ and ε = 0.2.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5 Radiative Nusselt number distribution on the horizontal walls at Ra = 10⁶ and ε = 0.2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 Distributions of upper and lower wall temperatures for Ra = 10⁶ and ε = 0.2.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Cross-section of the temperature at midheight cavity for Ra = 10⁶ and ε = 0.2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8 Cross section of the temperature at x = 0.5 for Ra = 10⁶ and ε = 0.2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 Horizontal velocity at the vertical cross section at x = 0.5 for Ra = 10⁶ and ε = 0.2.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10 Vertical velocity at the horizontal cross section at y = 0.5 for Ra = 10⁶ and ε = 0.2.</figcaption></figure>

wall, showing a deflection of peaks toward the cold wall. This behavior is due entirely to non-Boussinesq effects. In the center region, the fluid heats up with increasing ε<sub>b</sub> and the stratification decreases, as shown in Fig. 8.

The effect of the parameter ε<sub>b</sub> on the flow is presented in Figs. 9 and 10, where horizontal and vertical velocity components are plotted at x = 0.5 and y = 0.5, respectively. Compared to the case with a low value (ε<sub>b</sub> = 0.017), profiles of the U component show that the boundary layer contracts in the vicinity of the bottom wall and

<table><caption>Table 2 Variation of thermophysical properties as a function of Boussinesq parameter</caption><thead><tr><th>ε<sub>b</sub></th><th>ΔT</th><th>Δρ, %</th><th>Δk, %</th><th>Δμ, %</th></tr></thead><tbody><tr><td>0.0083</td><td>5</td><td>1.59</td><td>1.19</td><td>1.26</td></tr><tr><td>0.017</td><td>10</td><td>3.24</td><td>2.57</td><td>2.50</td></tr><tr><td>0.025</td><td>15</td><td>4.87</td><td>3.73</td><td>3.87</td></tr><tr><td>0.033</td><td>20</td><td>6.47</td><td>4.98</td><td>4.95</td></tr><tr><td>0.042</td><td>25</td><td>8.05</td><td>6.20</td><td>6.28</td></tr><tr><td>0.05</td><td>30</td><td>9.59</td><td>7.42</td><td>7.46</td></tr><tr><td>0.058</td><td>35</td><td>11.02</td><td>8.62</td><td>8.50</td></tr><tr><td>0.083</td><td>50</td><td>15.48</td><td>12.04</td><td>12.05</td></tr><tr><td>0.167</td><td>100</td><td>28.61</td><td>22.89</td><td>22.84</td></tr><tr><td>0.25</td><td>150</td><td>40.00</td><td>32.74</td><td>32.78</td></tr><tr><td>0.33</td><td>200</td><td>49.97</td><td>41.84</td><td>41.87</td></tr><tr><td>0.417</td><td>250</td><td>58.81</td><td>50.39</td><td>50.44</td></tr><tr><td>0.5</td><td>300</td><td>66.67</td><td>58.59</td><td>58.53</td></tr><tr><td>0.6</td><td>360</td><td>75.01</td><td>67.95</td><td>67.97</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Fig. 11 Variation of the average convection Nusselt number on the hot wall as a function of Boussinesq parameter.</figcaption></figure>

expands near the opposite wall. The same trend is observed for the V component near the vertical walls but not in the same proportions. It is interesting to note that, in the stratified region, the flow remains invariant for the component V, and the peak value near the cold vertical wall is slightly lower for $ε_b = 0.6$ than for $ε_b = 0.017$.

In Fig. 11, the low variations of Nusselt number for the weak values of the Boussinesq parameter $ε_b$ show that, under the condition $ε_b ≤ 0.05$, the incompressible model based on the Boussinesq approximation can be used to simulate both pure natural convection ($ε = 0$) and combined natural convection surface radiation ($ε > 0$). Furthermore, one may note that this interval of Boussinesq parameter corresponds to a temperature difference $ΔT ≤ 30°C$ and maximal variations of thermophysical properties less than 9.5% for the density and 7.5% for the conductivity and viscosity (see Table 2). It is also

<figure><img src="image_13.png" /><figcaption>Fig. 13 Distributions of local convective Nusselt numbers on the hot and cold walls for Ra = 10⁶ and ε = 0.5.</figcaption></figure>

interesting to note that Gray and Giorgini [21] showed in their study that, for air, the Boussinesq approximation is valid for $ΔT ≤ 28.6 °C$ corresponding to the density variation of 10%.

## B. Effect of Conductivity and Viscosity Variations

Another aim of this study is to examine the effect of fluid properties on the flow and heat transfer in the case of large temperature gradient. Simulation parameters for the compressible model are $ε_b = 0.6$, Ra = 10⁶, $ε = 0.5$, and a reference temperature $T_0 = 300$ K.

A comparison between results obtained with constant and variable properties is first illustrated by isotherms and streamlines plotted in Fig. 12. The most visible effect on isotherms occurs at the top wall.

<figure><img src="image_14.png" /><figcaption>Fig. 12 Streamlines (top) and temperature contours (bottom) for Ra = 10⁶ and ε = 0.5: a) constant properties, and b) Sutherland's law.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14 Distributions of local radiative Nusselt numbers on the vertical walls for Ra = 10⁶ and ε = 0.5.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15 Radiative Nusselt number distribution on the horizontal walls for Ra = 10⁶ and ε = 0.5.</figcaption></figure>

The center region of the cavity is marked by a sensible change in the fluid flow. Note that the symmetrical nature of the flow is strongly broken in the case of the constant.

Distributions of convective and radiative Nusselt numbers along the vertical walls are displayed in Figs. 13 and 14. By taking into account variable properties, the most noted change concerns the decrease of the convective Nusselt number on the cold wall at the upper part of the cavity. Unlike convection, $Nu_r$ profiles undergo a weak translation, showing a slight decrease of radiative heat transfer for variable properties. On the horizontal walls, the fluid properties weakly affect the radiative exchanges for which the effects are more visible on the upper wall (Fig. 15). Consideration of variable

<figure><img src="image_17.png" /><figcaption>Fig. 17 Cross section of the temperature at midheight cavity for Ra = 10⁶ and ε = 0.5.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 18 Horizontal velocity at the vertical cross section at x = 0.5 for Ra = 10⁶ and ε = 0.5.</figcaption></figure>

properties slightly enhances the radiative exchange at the top wall. The temperature profiles on the horizontal walls are higher for variable properties than for constant properties (Fig. 16).

The effect of fluid properties on the air temperature is evidenced in Fig. 17. By assuming variable properties, the temperature considerably increases in both boundary layers and at the core cavity. For the velocity components, the influence of fluid properties essentially appears near the bottom wall, at the cavity center for $U$, and near the hot wall for $V$. Variable properties lead to a decrease of the magnitude of both components, as seen in Figs. 18 and 19. It is interesting to note that the boundary layers near the hot wall remain thicker for both cases (constant or variable properties), reflecting that this behavior is mainly due to the nonlinear density variation.

<figure><img src="image_19.png" /><figcaption>Fig. 16 Distributions of upper and lower wall temperatures for Ra = 10⁶ and ε = 0.5.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19 Vertical velocity at the horizontal cross section at y = 0.5 for Ra = 10⁶ and ε = 0.5.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 20 Streamlines (top) and temperature contours (bottom) for Ra = 10⁶ and ε_b = 0.6.</figcaption></figure>

### C. Effect of Surface Radiation

In this section, computations are carried out for different values of emissivity ($0.1 \le \varepsilon \le 1.0$), Ra = 10⁶, ε_b = 0.6, and T₀ = 300 K. The properties are assumed to be variable according to Sutherland's law.

Typical isotherms and streamlines are plotted in Fig. 20 for different values of emissivity. Compared to the case without radiation, the temperature contours show the existence of temperature gradients near the top and bottom walls due to the radiative exchange. At the center region, the stratification is changed and the flow is modified but remains globally with the same structure: a primary flow along the walls and inside two asymmetrical secondary flows.

Figure 21 shows the variations of local convective Nusselt numbers on the vertical walls as a function of emissivity ε. Convective

exchanges on the hot wall are greatly influenced by the wall radiation, especially on the lower part of the hot wall. Increasing the wall emissivity ε produces a decreasing of Nu_h, whereas it is less significant for Nu_c. For the local radiative exchange, an opposite trend is observed. Increasing ε leads to an increase of Nu_r on the hot wall and a decrease on the cold wall (see Fig. 22). On the horizontal walls, it is clearly seen that an increase of the wall emissivity leads to the cooling of the top wall and the heating of the bottom wall (Fig. 23). This behavior is explained by the fact that the top wall loses heat while the bottom wall receives heat. The distribution of the radiative Nusselt number on the top and bottom walls, presented in Fig. 24, shows that the net radiative heat flux is essentially positive at the top and negative at the bottom. Vertical profiles of air temperature

<figure><img src="image_22.png" /><figcaption>Fig. 21 Variations of convective Nusselt numbers at the hot and cold walls.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 23 Distribution of upper and lower wall temperatures for Ra = 10⁶ and ε_b = 0.6.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22 Variations of radiative Nusselt number at the hot and cold walls.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 24 Distribution of radiative Nusselt number on upper and lower walls for Ra = 10⁶ and ε_b = 0.6.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 25 Cross section of the temperature at $x = 0.5$ for $Ra = 10^6$ and $\epsilon_b = 0.6$.</figcaption></figure>

show that the emissivity reduces the stratification in the core (Fig. 25).

Figure 26 shows the variations of the average convective and radiative Nusselt numbers ($Nuc_{avg}$, $Nur_{avg}$) according to the emissivity for $Ra = 10^5$, $10^6$, and $10^7$. With increasing of the wall emissivity $\epsilon$, the radiative Nusselt number on the hot wall increases and the convective Nusselt number $Nuc_{avg}$ decreases. On the cold wall, an opposite trend is observed for the radiative Nusselt number, whereas the convective Nusselt number keeps the same shape of variation.

Velocity components are plotted in Figs. 27 and 28 across the midplanes. Unlike the vertical velocity where $\epsilon$ has almost no effect, the influence of the wall emissivity on the horizontal component indicates that the flow is intensified near the horizontal walls and the stratified region is almost at rest.

<figure><img src="image_27.png" /><figcaption>Fig. 27 Vertical velocity at $y = 0.5$ for $Ra = 10^6$ and $\epsilon_b = 0.6$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 28 Horizontal velocity at $x = 0.5$ for $Ra = 10^6$ and $\epsilon_b = 0.6$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 26 Variation of the average convection and radiation Nusselt numbers on hot (top) and cold (bottom) walls as a function of emissivity (PC: constant properties, PV: variable properties).</figcaption></figure>
## V. Conclusions

In this Note, computations have been made in the context of a low-Mach-number formulation to study the combined natural convection and surface radiation in a square enclosure using a control volume approach. The present study shows the effects of the temperature difference parameter, surface emissivity, and temperature-dependent fluid properties on the flow and heat transfer. In view of the results presented, the main points can be summarized as follows:

1) Under non-Boussinesq conditions, and when the temperature difference increases, a significant variation between incompressible and compressible solutions is observed. In this case, the compressible model based on a low-Mach-number approximation can be used with a sufficient accuracy to simulate the combined natural convective surface radiation.

2) The variation of thermophysical properties could put a significant influence on the convective and radiative heat transfer and on the fluid flow.

3) The radiation exchange reduces the stratification at the cavity center and intensifies the flow at the horizontal walls.

4) The presence of radiation increases the temperatures of the fluid, and it considerably modifies the fluid flow and temperature distribution.

## References

[1] Ibrahim, J., and Paolucci, S., "The Modeling of Realistic Chemical Vapor Infiltration/Deposition Reactors," *International Journal for Numerical Methods in Fluids*, Vol. 64, No. 5, 2010, pp. 473-516.
doi:10.1002/fld.v64:5

[2] Klein, R., "Multiple Spatial Scales in Engineering and Atmospheric Low Mach Number Flows," *ESAIM: Mathematical Modeling and Numerical Analysis (M2AN)*, Vol. 39, No. 3, 2005, pp. 537-559.
doi:10.1051/m2an:2005022

[3] Klein, R., "Asymptotic Analyses for Atmospheric Flows and the Construction of Asymptotically Adaptive Numerical Methods," *Zeitschrift für Angewandte Mathematik und Mechanik*, Vol. 80, Nos. 11-12, 2000, pp. 765-777.
doi:10.1002/(ISSN)1521-4001

[4] Madja, A., and Sethian, J., "The Derivation and Numerical Solution of the Equations for Zero Mach Number Combustion," *Combustion Science and Technology*, Vol. 42, No. 3-4, 1985, pp. 185-205.
doi:10.1080/00102208508960376

[5] Chen, S., Liu, Z., Zhang, C., He, Z., Tian, Z., Shi, B., and Zheng, C., "A Novel Coupled Lattice Boltzmann Model for Low Mach Number Combustion Simulation," *Applied Mathematics and Computation*, Vol. 193, No. 1, 2007, pp. 266-284.
doi:10.1016/j.amc.2007.03.087

[6] Paolucci, S., "On the Filtering of Sound from the Navier-Stokes Equations," Sandia National Labs. TR-82-8257, Albuquerque, NM, 1982.

[7] Hamimid, S., Guellal, M., and Bouafia, M., "Numerical Study of Natural Convection in a Square Cavity Under Non-Boussinesq Conditions," *Thermal Science* [online], 2014, Paper 84.
doi:10.2298/TSCI130810084H

[8] Suslov, S. A., and Paolucci, S., "Nonlinear Analysis of Convection Flow in a Tall Vertical Enclosure Under Non-Boussinesq Conditions," *Journal of Fluid Mechanics*, Vol. 344, Aug. 1997, pp. 1-41.
doi:10.1017/S0022112097005971

[9] Suslov, S. A., "Analysis of Coupling Between Hydrodynamic and Thermal Instabilities in Non-Boussinesq Convection," *Nonlinear Analysis*, Vol. 63, Nos. 5-7, 2005, pp. e1435-e1443.
doi:10.1016/j.na.2005.01.057

[10] Codina, R., "A Stabilised Finite Element Method for Generalised Stationary Incompressible Flows," *Computer Methods in Applied Mechanics and Engineering*, Vol. 190, Nos. 20-21, Feb. 2001, pp. 2681-2706.
doi:10.1016/S0045-7825(00)00260-7

[11] Vierendeels, J., Merci, B., and Dick, E., "A Multigrid Method for Natural Convective Heat Transfer with Large Temperature Differences," *Journal of Computational and Applied Mathematics*, Vol. 168, Nos. 1-2, July 2004, pp. 509-517.
doi:10.1016/j.cam.2003.08.081

[12] Liu, W., and Makhviladze, G., "An Implicit Finite Element Solution of Thermal Flows at Low Mach Number," *Journal of Computational Physics*, Vol. 227, No. 5, Feb. 2008, pp. 2743-2757.
doi:10.1016/j.jcp.2007.10.025

[13] Bouafia, M., and Daube, O., "Natural Convection for Large Temperature Gradients Around a Square Solid Body Within a Rectangular Cavity," *International Journal of Heat and Mass Transfer*, Vol. 50, Nos. 17-18, Aug. 2007, pp. 3599-3615.
doi:10.1016/j.ijheatmasstransfer.2006.05.013

[14] Chenoweth, D. R., and Paolucci, S., "Natural Convection in an Enclosed Vertical Air Layer with Large Horizontal Temperature Differences," *Journal of Fluid Mechanics*, Vol. 169, Aug. 1986, pp. 173-210.
doi:10.1017/S0022112086000587

[15] Suslov, S. A., and Paolucci, S., "Stability of Natural Convection Flow in a Tall Vertical Enclosure Under Non-Boussinesq Conditions," *International Journal of Heat and Mass Transfer*, Vol. 38, No. 12, Aug. 1995, pp. 2143-2157.
doi:10.1016/0017-9310(94)00348-Y

[16] White, F. M., *Viscous Fluid Flow*, McGraw-Hill, New York, 1974, pp. 28-29.

[17] Bouali, H., Mezrhaba, A., Amaouia, H., and Bouzidib, M., "Radiation-Natural Convection Heat Transfer in an Inclined Rectangular Enclosure," *International Journal of Thermal Sciences*, Vol. 45, No. 6, June 2006, pp. 553-566.
doi:10.1016/j.ijthermalsci.2005.10.001

[18] Patankar, S. V., *Numerical Heat Transfer and Fluid Flow*, McGraw-Hill, New York, 1980, pp. 131-132.

[19] Howell, J. R., *A Catalog of Radiation Configuration Factors*, McGraw-Hill, New York, 1982, p. 182.

[20] Wang, H., Xin, S., and Quéré, P. L., "Étude Numérique du Couplage de la Convection Naturelle avec le Rayonnement de Surfaces en Cavité Carrée Remplie d'Air," *Comptes Rendus Mécanique*, Vol. 334, No. 1, 2006, pp. 48-57.
doi:10.1016/j.crme.2005.10.011

[21] Gray, D. D., and Giorcini, A., "The Validity of the Boussinesq Approximation for Liquids and Gases," *International Journal of Heat and Mass Transfer*, Vol. 19, No. 5, May 1976, pp. 545-551.
doi:10.1016/0017-9310(76)90168-X