# Numerical simulation of a double skin with secondary ventilation flow on adiabatic wall

Mohamed Bouraoui, M.S. Rouabah, Aissa Abidi-Saad, A. Korichi, Catalin Viorel Popa, Guillaume Polidori

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

## To cite this version:

Mohamed Bouraoui, M.S. Rouabah, Aissa Abidi-Saad, A. Korichi, Catalin Viorel Popa, et al.. Numerical simulation of a double skin with secondary ventilation flow on adiabatic wall. Frontiers in Heat and Mass Transfer, 2017, 8, 10.5098/hmt.8.18. . hal-03162618

HAL Id: hal-03162618

https://hal.science/hal-03162618v1

Submitted on 21 Sep 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

# NUMERICAL SIMULATION OF A DOUBLE SKIN WITH SECONDARY VENTILATION FLOW ON ADIABATIC WALL

M. Bouraoui<sup>a</sup>, M. S. Rouabah<sup>a</sup>, A. Abidi-Saad<sup>b,c,d,*</sup>, A. Korichi<sup>e</sup>, C. Popa<sup>b</sup>, G. Polidori<sup>b</sup>

## ABSTRACT

Study concerning a double skin flow with secondary ventilation was conducted numerically, in order to understand the basic mechanisms of the free convection in an open channel asymmetrically heated with uniform heat flux density (510 W/m²). The vertical channel corresponds to a double skin façade, which was immersed in a tank filled with water. The tank corresponds to the environment which allows us to overcome pressure conditions at the inlet and the outlet of the channel. The use of water allows neglecting radiation effect. The mass conservation equations of momentum and energy are solved using the finite volume method and numerical simulations are performed using Ansys Fluent CFD software. Our interest, in this work, is the study of the influence of the position and the size of an opening in the adiabatic wall, which represents secondary ventilation, on the flow dynamics and the heat transfer.

**Keywords:** Natural convection; asymmetric heating; vertical channel; CFD;

## 1. INTRODUCTION

In this work, we are interested in the study of laminar flow of natural convection in a vertical open channel heated asymmetrically. The benefit of this study is that the plane channel is representative of many problems such as chimney, Trombe wall, solar collector or double-skin façade. The flow in an asymmetrically heated vertical plane channel at constant flux can be considered as a prototype of semi-confined natural convection flows.

Elenbaas (1942) experimentally observed two flow dynamic regimes within a vertical channel. The appearance of these regimes depends on the modified Rayleigh number Ra<sup>+</sup>. Indeed, for a relatively low Rayleigh number the author founds a fully developed regime in which all the fluid entering from below exits through the top of the channel with a profile similar to that of Poiseuille. However, for high Rayleigh number, a boundary layer is observed near the heated wall; the fluid simultaneously enters from below and exits through the top of the channel via a reversed flow along the unheated wall. Aung et al. (1972) highlighted also two regimes in a numerical study concerning a natural convection induced flow in a channel heated symmetrically or asymmetrically: the first is developed at low Rayleigh numbers and the second is of boundary layer type for high Rayleigh numbers. Azevedo and Sparrow (1986) interested in the heat transfer intensification in natural convection flows in a vertical plane channel by studying the influence of the opening on the unheated wall. In this configuration, it has been shown that the average Nusselt number on the heated wall is almost independent of the location and size of the opening.

Webb and Hill (1989) have observed an upward flow of boundary layer type near the heated wall, coming from the bottom of the channel. This latter is accompanied by a downward flow developing in the vicinity of the opposite unheated wall coming from the top of the channel. The dynamic structure of the flow has been evidenced by Ospir et al. (2012) using a visualization technique by laser tomography with discrete tracers. These authors focused specifically on the influence of both the channel aspect ratio (A/b) and Rayleigh number on the flow features. Their results show that the reversed flow occupies different proportions of the channel width for different modified Rayleigh numbers. This finding has been verified in a comparative study between simulations and experimental results performed by Popa et al. (2012).

Concerning the study of the radiation contribution in a vertical channel, Li et al. (2013) have shown that the heat radiation causes the disappearance of the recirculation in the upper part of the channel. Finally, a recent experimental study based on laser tomography by Polidori et al. (2015) in a vertical plane channel asymmetrically heated, allowed to observe simultaneously the appearance of the dynamic boundary layer and a reversed flow coming from the top of the channel at the very first times of the heating process.

The present numerical work deals with the study of the influence of both the position and the size of an opening at the adiabatic wall, considered as a secondary ventilation, on the flow dynamics and heat transfer. The vertical plane channel was immersed in a tank filled with water corresponding to the domain of study. The use of water as working fluid is to neglect the thermal radiation. The numerical use of the tank domain as the external channel environment is to overcome the pressure boundary conditions at the inlet and the outlet of the channel. Indeed, the

<footnote><sup>a,b</sup> Laboratoire de génie climatique, Département de génie climatique, Faculté des Sciences de la Technologie. Université Frères Mentouri – Constantine 1. Campus Ahmed Hamani, 25000 Constantine, Algeria</footnote>

<footnote><sup>b</sup> GRESPI/Thermomécanique, UFR sciences Exactes et Naturelle, Campus du Moulin de la Housse-BP 1039, 51687 Reims, France</footnote>

<footnote><sup>c</sup>Laboratoire d'Energétique Appliquée et de Pollution, Université Constantine 1, Constantine 25000, Algeria</footnote>

<footnote><sup>d</sup>Faculté des Hydrocarbures et Energies renouvelables et Sciences de la terre et de l'univers, Université Kasdi Merbah, 30000 Ouargla, Algérie</footnote>

<footnote><sup>e</sup>Laboratoire de Mécanique, Physique et Modélisation Mathématique, Université de Médéa, Médéa, Algeria</footnote>

<footnote>* Corresponding Author. Email: asais83@yahoo.fr</footnote>
study can be made complicated in the case of natural convection, if the channel is only the domain of the study.

## 2. CHANNEL DESCRIPTION

The channel consists of two parallel planar vertical walls of 376 mm high separated by a gap b (b = 36 mm). One wall is composed of a heated central part (height A = 188 mm) surrounding by two unheated parts (length A/2) while the opposite wall remains entirely unheated (Fig. 1). In addition, for a better control of the flow conditions at the entrance, a quarter circle (R = 36mm) is added at the bottom of each wall. The channel is immersed in a vertical tank (500 x 500 x 1000 mm³) made of 20 mm thick Plexiglas® plates and filled with water.

The second wall is an adiabatic wall with an opening of size equal respectively to b/3, 2b/3, b and 4b/3 and located at the bottom, middle and top of the heated area. The channel heating system consists of a heating fabric (20 x 30 cm²) and covered by an aluminum plate of 3 mm thick, which provides a uniform and controllable heat flux density. The rear face of the heated wall is made adiabatic by adding a layer of polystyrene of 140 mm thick, to minimize the heat loss and to prevent the occurrence of convection in the external environment.

<figure><img src="image_5.png" /><figcaption>Fig. 1 (a) Channel geometry, (b) Boundary conditions</figcaption></figure>

## 3. PHYSICAL MODEL

Natural convection is governed by the continuity equation (1), the momentum equations (2, 3) and the conservation equation of energy (4).

$$ \frac{\partial \rho u}{\partial x} + \frac{\partial \rho v}{\partial y} = 0, \quad \left( \frac{d\rho}{dt} \right) = 0 \qquad (1) $$

$$ \left( \frac{\partial u}{\partial t} + \rho u \frac{\partial u}{\partial x} + \rho v \frac{\partial u}{\partial y} \right) = - \frac{\partial P^*}{\partial x} + \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{\partial u}{\partial x} \frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial x} \frac{\partial \mu}{\partial y} \qquad (2) $$

$$ \left( \frac{\partial v}{\partial t} + \rho u \frac{\partial v}{\partial x} + \rho v \frac{\partial v}{\partial y} \right) = - \frac{\partial P^*}{\partial y} + \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{\partial u}{\partial y} \frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial y} \frac{\partial \mu}{\partial y} - \rho g \qquad (2) $$

$$ \rho C_p \left( \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} \right) = \frac{\partial}{\partial x} \left( k \frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial y} \left( k \frac{\partial T}{\partial y} \right) \qquad (4) $$

To simplify the model that describes the natural convection in the vertical channel, some simplifying assumptions have been used. In this study, it is assumed that the flow is laminar, two-dimensional and steady. The fluid is also assumed Newtonian, incompressible and the physical properties of the fluid depend only on the temperature as described in Abidi-Saad et al. (2017) and Hemmer et al. (2016). The simulations were conducted for a modified Rayleigh number (Ra*) equal to 4.5×10⁶. The modified Rayleigh number is based on the flux density (q") and the width of the channel (b) and takes into account the aspect ratio of the heated portion (A/b) through Equation (5).

$$ Ra^* = \frac{g\beta q''b^4}{kv^2} \frac{b}{A} Pr \qquad (3) $$

## 4. NUMERICAL METHOD

Transport equations of mass, momentum and energy (1) – (4) are solved numerically using the finite volume method (Patankar, 1980). This method is based on the spatial integration of transport equations on control volumes. In the solver configuration the pressure-based solver was used. This solver can solve sequentially transport equations. The coupling between velocity and pressure is achieved with the Coupled algorithm. It solves the equations of continuity and momentum simultaneously and gives an advantage to treat the flows with a strong interdependence between dynamic and thermal fields. Central-Differenced Schemes and Second-Order Upwind Schemes are used for the spatial discretization for diffusive terms and convective terms respectively. The time advance uses the Second-Order Implicit formulation. The convergence criteria were based on the absolute residuals resulting from the integration of the conservation equations over finite control volumes. For all simulations performed in this study, converged solutions were achieved after a residuals decrease larger than 10⁻⁴ for all governing equations. 2D numerical simulations are performed with ANSYS Fluent® CFD commercial software.

The mesh is uniform in the channel then it is non-uniform in the rest of the computational domain as depicted in Fig. 2.

We performed several grid refinement experiments inside the channel (33000, 25000, 27000, 28000 and 30000) relative to the mass flow rate to determine the suitable grid size required for the numerical simulation, for the case of Ra* = 4.5×10⁶ and a size of the opening equal to 2b/3. At the end of the preliminary calculations, a grid size 27000 was judged sufficient to get results with a satisfactory precision, which corresponds to 42072 cells in the entire calculation domain. This grid was used to carry out the final computations. Because there is no much improvement in precision, successive mesh refinement leads to a change in values occurring in the fourth digit after the decimal point, with a relative error in the mass flow rate less than 2%.

<figure><img src="image_6.png" /><figcaption>Fig. 2 Mesh distribution in the fluid domain</figcaption></figure>

## 5. RESULTS

In the following section, the results of the influence of the position and the size of an opening in the adiabatic wall of an asymmetrically heated channel on the flow dynamics and heat transfer are described in details. The opening represents a secondary ventilation window. The results were obtained for a modified Rayleigh number Ra* = 4.5×10⁶, an aspect ratio A/b = 5.2, a heat flux density of 510 W.m⁻² for several opening size (b/3, 2b/3, b and 4b/3) and for three locations; i.e. at the bottom, middle and top of the heated area. The numerical results were stopped after 30 minutes of heating because it has been experimentally shown that the flow structure becomes steady from about 30 min of heating (Ospir et al., 2012). We used the configuration of the channel without opening studied by Popa et al. (2012) as the basic configuration.
<figure><img src="image_7.png" /><figcaption>**Fig. 3** Streamline patterns in the channel for different opening sizes located at the (a) bottom, (b) middle and (c) top of the channel</figcaption></figure>

An established regime is found and an upward flow of boundary layer type is developed along the heated wall fed by fluid coming from the

entrance of the channel with a V-shaped recirculation zone nearby the unheated wall.

For the aim of validating our simulation algorithm, the present numerical streamlines were compared with those found experimentally. This comparison is well described in Hemmer et al. (2016) and it is not necessary to present it here.

## 5.1 The influence of the opening on the flow structure and the temperature field

The results allow to observe the effect of the opening position of various sizes (b/3, 2b/3, b and 4b/3) on the flow structure and the temperature field within the channel. Indeed, it can be seen in Figures (3 a, b and c) that the position of the opening plays an important role on the depth of the recirculation zone. If the opening is located in the lower part of the channel, the recirculation zone is almost similar to that observed in the unopened configuration. Whereas, if the opening is placed in the middle of the channel, the boundary layer is also fed by the cold fluid coming from the opening, leading to the decrease of the recirculation zone size.

However, the presence of the opening at the top of the channel leads to the feeding of only the recirculation zone by cold fluid and not directly the boundary layer flow.

Figures (4 a, b and c) show the temperature field within the channel. It can be seen that the presence and position of the opening on the unheated wall has no effect on the thermal boundary layer which is the motor of the flow. The presence of the opening in the bottom of the channel increases the heat transfer more than the other positions do.

## 5.2 The influence of the opening on the axial velocity in inlet and outlet of the channel

The influence of the opening size on the axial velocity component in the entrance and the exit of the channel is plotted in Figs. (5 a, b and c) and Figs. (6 a, b and c) respectively, for the three different positions of the opening along the channel height i.e. bottom, middle and top of the channel.

From Figs. (5 a, b and c) it can be seen that, at the entrance of the heated zone, the velocity profiles are approaching the laminar flow profiles for all studied cases, with a slight asymmetric velocity peak highlighted by the acceleration of the fluid in the vicinity of the heated wall. This is due to the presence of the opening on the unheated wall. For the case of the bottom position (Fig. 5 a), the maximum velocity magnitude is found for the small opening size (b/3). Because the reversed flow zone is very small only a small quantity of cold fluid enters from the opening, which lets more space for the fluid coming from the channel inlet. For the largest opening size (4b/3), the velocity profile showed a minimal velocity magnitude with negative velocity near the adiabatic wall, due to the blockage caused by the significant quantity of fluid entering through the opening and the presence of small vortices in the channel inlet (Fig. 3 a).

Figs. 6 (a) – (c) respectively show the vertical velocity profiles at the exit of the channel for the three opening positions (bottom, middle and top of the channel) and for various opening sizes. The velocity profiles are of boundary layer type with a strong acceleration close to the heated wall for all studied cases. Negative velocities are highlighted close to the adiabatic wall for all cases indicating the reversed flow presence. For the velocity profiles at the channel inlet and outlet, the maximal velocity magnitude varied with the variation of the opening size and position due to the variation of the reversed flow size and the variation of the quantity of fluid that enters from the opening.

## 5.3 The influence of the opening on the mass flow rate

Tables 1, 2 and 3 respectively depict the effect of the opening sizes on the mass flow rate for the three considered positions i.e. bottom, middle and top. The mass flow rates are accounted at the channel inlet, the opening and the channel outlet.
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Isotherm patterns in the channel for different opening sizes located at the (a) bottom of the channel, (b) middle of the channel and (c) top of the channel</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5** vertical velocity profile at the channel inlet for an opening of different sizes located at the (a) bottom of the channel, (b) middle of the channel and (c) top of the channel</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 6** vertical velocity profile at the channel outlet for the different opening sizes located at the (a) bottom of the channel, (b) middle of the channel and (c) top of the channel</figcaption></figure>

One can see from the three tables that the mass flow rate is perfectly conserved for each case. Indeed, if the opening is placed at the bottom of the channel (Table 1), it is noticed that increasing the size of the opening causes a reduction of the mass flow rate at the entrance of the channel. Therefore, a large opening placed in the bottom of the channel causes a sharp decrease in the flow rate at the inlet of the channel. Note also that an opening of size equal to the channel width b has almost no effect on the flow rate. As it can be seen from streamlines patterns in Fig.

3 that a significant quantity of fluid coming from the channel inlet blocks the opening.

For the middle position (Table 2), the flow rate through the channel entrance takes a “sinusoidal pattern” of amplitude equal to 0.006 kg/s, then, a maximum flow rate value equal to 0.024 kg/s for opening sizes 2b/3 and 4b/3. Moreover, a minimal value equal to 0.018 kg/s for a size equal to b/3 and b (the width of the channel). In addition, it can be seen that the mass flow rate at the exit of the channel increases with the increase in opening size. The maximal amount of mass flow rate at the channel exit is found for the largest opening size (4b/3) and for the middle position which is equal to 0.053 (kg/s) correspond to flow rates equal to 0.024(kg/s), 0.029(kg/s) at the channel inlet and opening, respectively.

The presence of the opening in the top of the channel (Table 3) causes a reduction of the mass flow rate when the size of the opening is greater than 2b/3. Therefore, increasing the size of the secondary ventilation placed in the upper part of the channel more than 2b/3 leads to the decrease of the flow rate values at the inlet of the channel. Due to the direct feeding of the reversed flow.

<table><caption>**Table 1:** mass flow rate variation as a function of the opening size – opening at the bottom of the channel</caption><thead><tr><th>Mass flow rate (kg/s)</th><th>b/3</th><th>2b/3</th><th>b</th><th>4b/3</th></tr></thead><tbody><tr><th>Channel inlet</th><td>0.037</td><td>0.014</td><td>0.0302</td><td>0.006</td></tr><tr><th>Opening</th><td>0.011</td><td>0.029</td><td>0.0007</td><td>0.0427</td></tr><tr><th>Channel outlet</th><td>0.048</td><td>0.043</td><td>0.0308</td><td>0.049</td></tr></tbody></table>

<table><caption>**Table 2:** mass flow rate variation as a function of the opening size – opening at the middle of the channel</caption><thead><tr><th>Mass flow rate (kg/s)</th><th>b/3</th><th>2b/3</th><th>b</th><th>4b/3</th></tr></thead><tbody><tr><th>Channel inlet</th><td>0.018</td><td>0.024</td><td>0.018</td><td>0.024</td></tr><tr><th>Opening</th><td>0.012</td><td>0.010</td><td>0.019</td><td>0.029</td></tr><tr><th>Channel outlet</th><td>0.031</td><td>0.034</td><td>0.038</td><td>0.053</td></tr></tbody></table>

<table><caption>**Table 3:** mass flow rate variation as a function of the opening size – opening at the top of the channel</caption><thead><tr><th>Mass flow rate (kg/s)</th><th>b/3</th><th>2b/3</th><th>b</th><th>4b/3</th></tr></thead><tbody><tr><th>Channel inlet</th><td>0.029</td><td>0.032</td><td>0.029</td><td>0.019</td></tr><tr><th>Opening</th><td>0.016</td><td>0.018</td><td>0.005</td><td>0.018</td></tr><tr><th>Channel outlet</th><td>0.045</td><td>0.050</td><td>0.034</td><td>0.038</td></tr></tbody></table>

## 6. CONCLUSIONS

The effect of an opening (secondary ventilation) on the heat transfer and fluid flow in an asymmetrically heated channel has been numerically analyzed for a modified Rayleigh number equal to $4.5×10^6$, channel aspect ratio equal to 5.2 and for four sizes of the opening and three different positions along the channel height. Overall, the thermal fields and flow characteristics are strongly affected by the opening sizes and positions. The recirculation zone was always observed at the channel outlet whatever the opening size and position.

* The recirculation length gradually decreases with the increase of the opening size for the middle position and the trend is reversed for the top position.

* The flow coming from the opening feeds directly the boundary layer flow for the bottom and middle positions. Whereas it feeds only the recirculation zone for the top position.

* When the opening is in the middle of the channel, it allows to cool more effectively the channel. In addition, the larger the size of the opening is, the more remarkable the cooling becomes.
* For both top and bottom opening positions the mass flow rate at the channel inlet decreases with the increase of the opening size. Whereas, for the middle position the mass flow rate has a sinusoidal form.

* The maximal mass flow rate at the channel outlet is found for the middle position at the largest opening size 4b/3.

## NOMENCLATURE

* b channel wall spacing (m)

* C<sub>p</sub> specific heat capacity (J/kg.K)

* g acceleration of gravity (m/s<sup>2</sup>)

* k thermal conductivity (W/m.K)

* L heated zone length (m)

* P* driving pressure (Pa)

* Pr Prandtl number

* q'' heat flux density (W/m<sup>2</sup>)

* Ra Rayleigh number

* Ra* modified Rayleigh number

* t time (s)

* T temperature (K)

* u,v velocity (m/s)

* x,y coordinate (m)

### Greek Symbols

* β volume expansion coefficient (1/K)

* φ heat flux density (W/m<sup>2</sup>)

* ρ density (kg/m<sup>3</sup>)

* ν kinematic viscosity (m<sup>2</sup>/s)

## REFERENCES

* Elenbaas W., 1942, "Heat Dissipation of Parallel Plates by Free Convection," Physica 9(9), 1-23.
http://dx.doi.org/10.1016/S0031-8914(42)90053-3

* Aung W., Fletcher L.S., Sernas V., 1972, "Developing Laminar Free Convection Between Vertical Flat Plates with Asymmetric Heating," International Journal Of Heat and Mass Transfer, 15, 2293-2308.
http://dx.doi.org/10.1016/0017-9310(72)90048-8

* Azevedo, L.F., and Sparrow, E.M., 1986, "Natural Convection in a Vertical Channel Vented to the Ambient Through an Aperture in the

Channel Wall," International Journal of Heat and Mass Transfer, 29, 819-830.
http://dx.doi.org/10.1016/0017-9310(86)90178-X

Weeb, B.W., Hill, D.P., 1989, "High Rayleigh Number Laminar Natural Convection in an Asymmetrical Heated Vertical Channel," ASME J. Heat transfer, 111, 649-656.
http://dx.doi.org/10.1115/1.3250732

Ospir, D., Popa, C., Chereches, C., Polidori, G., and Fohanno, S., 2012 "Flow Visualization of Natural Convection in a Vertical Channel with Asymmetric Heating," International Communications in Heat and Mass Transfer, 39(4), 486-493.
http://dx.doi.org/10.1016/j.icheatmasstransfer.2012.02.005

Popa, C., Ospir, D., Fohanno, S., and Chereches, C., 2012, "Numerical Simulation of Dynamics Aspects of Natural Convection Flow in a Double-skin Façade," Energy and Buildings, 50, 229-233.
http://dx.doi.org/10.1016/j.enbuild.2012.03.042

Li, R., Bousetta, M., Chenier, E., and Lauriat, G., 2013, "Effect of Surface Radiation on Natural Convective Flows and Onset of Flow Reversal in Asymmetrically Heated Channels," International Journal of Thermal Sciences, 65, 9-27.
http://dx.doi.org/10.1016/j.ijthermalsci.2012.10.023

Polidori, G., Fatnassi, S., Ben Maad, R., Fohanno, S., and Beaumont, F., 2015, "Early-Stage Dynamics in The Onset of Free-convective Reversal Flow in an Open-ended Channel Asymmetrically Heated," International Journal of Thermal Sciences, 88, 40-46.
http://dx.doi.org/10.1016/j.ijthermalsci.2014.09.011

Abidi-Saad, A., Kadja, M., Popa, C., and Polidori, G., 2017, "Effect of Adiabatic Square Ribs on Natural Convection in an Asymmetrically Heated Channel," Heat Mass Transfer, 53(2), 743-752.
http://dx.doi.org/10.1007/s00231-016-1853-y

Hemmer, C., Popa, C.V., Sergent, A., Polidori, G., 2016, "Heat and Fluid Flow in an Uneven Heated Chimney," International Journal of Thermal Sciences, 107, 220-229.
http://dx.doi.org/10.1016/j.ijthermalsci.2016.04.015

Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere Publication Co., NY.