Heat and nanofluid transfer in baffled channels of different outlet models

Page:

21-28

DOI:

https://doi.org/10.18280/mmep.060103

OPEN ACCESS

Abstract:

The paper is concerned with the effects of baffled obstacles on steady turbulent Al2O3-H2O nanofluid flow and heat transfer characteristics through channels in different outlet models. The first channel has an outlet as its entrance (case A). The second (case B), third (case C), and fourth (case D) channels have narrow, upper, lower, and central exits, with 45 per cent of their entrance, respectively. These effects are investigated with the help of CFD in a 2D model. The numerical data show improvements in the heat transfer rate of about 45.071, 58.404, 82.413, and 92.433 per cent for cases A, B, C, and D compared to the smooth channel using the same solid volume fraction of Al2O3 nanoparticle, respectively. Among the most effective channels on heat transfer is case D, about 37.658, 21.356, and 9.348 per cent compared to cases A, B, and C, respectively for the maximum value of Reynolds number.

Keywords:

*nanofluid, forced convection, turbulent flow, fluid mechanics, baffle, channel*

1. Introduction

The field of nanofluid flows can be employed to enhance thermal transfer in many engineering applications such as heat exchangers, solar energy collectors, etc. Recently, several studies on the modeling of convective heat transfer in nanofluids are published. Izadi et al. [1] presented the mixed convection heat transfer and entropy generation of a nanofluid containing carbon nanotubes, flowing in a 3D rectangular channel, subjected to opposed buoyant forces. Chamkha et al. [2] investigated the entropy generation and natural convection inside a C-shaped cavity filled with CuO-water nanofluid and subjected to a uniform magnetic field. Chamkha et al. [3] performed an analysis to study non-Darcy free convection boundary-layer flow over a permeable vertical cone embedded in a porous medium saturated with a nanofluid in the presence of uniform wall transpiration. Chamkha et al. [4] considered the problem of mixed convection MHD flow of a Nanofluid past a stretching permeable surface. The effects of heat generation or absorption, buoyancy, thermophoresis, Brownian motion and wall suction or injection were investigated. Sadripour and Chamkha [5] investigated the effects of different morphology of supported nanoparticles, including copper (Cu), silver (Ag), aluminum dioxide (Al_{2}O_{3}), boehmite alumina (γ-AlOOH), molybdenum disulfide (MoS_{2}) and silicon dioxide (SiO_{2}), on heat transfer and entropy generation in compression with each other in case of a water-based heat-sink solar collector located in Isfahan city, Iran. Chamkha et al. [6] presented a boundary layer analysis for the mixed convection past a vertical wedge in a porous medium saturated with a non-Newtonian nano fluid. Chamkha and Selimefendigil [7] numerically examined the mixed convection of CuO-water nanofluid due to a rotating inner hot circular cylinder in a 3D cubic enclosure with phase change material (PCM) attached to its vertical surface. Chamkha et al. [8] investigated the effect of uniform lateral mass flux on non-Darcy natural convection of non-Newtonian fluid along a vertical cone embedded in a porous medium filled with a nanofluid. Chamkha and Selimefendigil [9] numerically investigated the efficiency of a PV/T module with SiO_{2}-water nanofluid for nanoparticle properties (shape, volume fraction) and for different operating conditions. Chamkha et al. [10] numerically analyzed the natural convection of a magnetohydrodynamic nanofluid in an enclosure under the effects of thermal radiation and the shape factor of nanoparticles using the control-volume-based finite element method (CVFEM). Chamkha et al. [11] presented a boundary layer analysis for the natural convection past a sphere embedded in a porous medium saturated with a nanofluid. Chamkha et al. [12] studied the problem of steady, laminar, mixed convection boundary-layer flow over an isothermal vertical wedge embedded in a porous medium saturated with a nanofluid, in the presence of thermal radiation. Other similar works can be found in the literature as Minea and Moldoveanu [13], Akhatov et al. [14], Dehaj and Mohiabadi [15], Kaya et al. [16], Mirzaei Kaya et al. [17].

Due to the importance of this new kind of fluid in various industrial applications, i.e., separation processes in chemical industries, petroleum industries, building thermal insulation, geothermal reservoirs, storage of nuclear waste materials, transpiration cooling, and solar heating systems, due to its high thermal conductivity, and due to its distinct physical properties, in this study we highlight the importance of these fluids on fluid flow and heat transfer characteristics inside rectangular channels with different exit configurations and upper and lower wall-attached baffle plates, using the finite volume method by means of the Commercial CFD software FLUENT in a 2D model.

2. Model under Analysis

The purpose of this study is to carry out a numerical study on the dynamic and thermal behavior of a nanofluid with a constant property and flowing turbulently through two-dimensional horizontal rectangular channels in different outlet models. The first channel has an outlet as its entrance (case A). The second (case B), third (case C), and fourth (case D) channels have narrow, upper, lower, and central exits, with 45 per cent of their entrance, respectively. The upper surfaces were put in a constant temperature condition, while the lower surfaces were thermally insulated. Two transverse, solid, flat baffled obstacles were inserted into the channels and fixed to the top and bottom walls, in a periodically staggered manner to force vortices to improve the mixing, and consequently the heat transfer.

**Figure 1.** Baffled channels with various outlet models: (a) case A, (b) case B, (c) case C, and (d) case D. Al_{2}O_{3}/H_{2}O Nanofluid flow is from left to right

Figure 1 depicts a schematic representation of the physical model. Demartini et al. [18] conducted an experimental analysis which served as the basis for the detailed structural parameters used. Al_{2}O_{3}-H_{2}O is the working fluid used, and the Reynolds numbers considered range from 12,000 to 32,000.

The following assumptions were made to develop the numerical model for the fluid flow and heat transfer in the computational domain under consideration: the flow is steady-state, turbulent, and two-dimensional; nanofluid in single phase, Newtonian, and incompressible; the physical properties of the fluid and solid are kept constant; and the body forces, viscous dissipation and radiation heat transfer are not considered. The thermal-physical properties of the Al_{2}O_{3}-H_{2}O nanofluid with a volume fraction 2% at 300 K are given in Ref [19]. The Reynolds-averaged Navier-Stokes equations, with the standard k-epsilon turbulence model and the energy equation are the governing equations to be considered in this study as reported by Ajeel et al. [20].

A uniform one-dimensional velocity profile (U_{in}) was used as the aeraulic boundary condition at the intake of the computational domain. The temperature (T_{in}) of the Al_{2}O_{3}-H_{2}O working fluid was set equal to 300 K at the inlet of the channel. Turbulence intensity (TI_{in}) equal to 2 per cent was selected for the intake height. The thermal boundary condition consisted of the constant temperature (T_{w}) of 375 K which was applied to the upper wall of the computational domain. The bottom surface of the computational domain was taken as adiabatic. Moreover, it was decided to impose the no-slip and impermeability boundary conditions at all the solid walls. Note that the atmospheric pressure (P_{atm}) is prescribed at the channel outlet.

The detail on governing parameters i.e., hydraulic diameter, Reynolds number, friction factor, and average Nusselt number can be found in similar work by Menni et al. [21]. The nanofluid density, thermal conductivity, viscosity and its heat capacity can be calculated as reported by Chamkha et al. [22].

3. Numerical Analysis

The fluid is considered Newtonian, incompressible with constant properties. The Reynolds averaged Navier-Stokes equations, along with the standard k-epsilon turbulence model [23] and the energy equation, are used to control the channel flow model. The finite volume method [24] is used to integrate all the equations in two-dimensions; the commercial CFD software FLUENT along with the SIMPLE-algorithm [24] is used for pressure-velocity coupling. Various values of the Reynolds number and different channel outlet models were selected to perform the numerical runs, using Al_{2}O_{3}-H_{2}O as the working medium. Then, the solutions are said to be converging when the normalized residuals are smaller than

10^{-12} and 10^{-9} for the energy equation and the other variables, respectively. Figure 2 compares the numerical and experimental data of Demartini et al. [18] with those of the present CFD simulation on the dynamic pressure coefficient in the case A with the air conventional fluid, for the value Re = 8.73×10^{4}. The numerical and experimental profiles of the dynamic pressure coefficient were calculated along the height of the channel at positions x = 0.223 m. This comparison indicates that the two results related to the dynamic coefficient of pressure show qualitative agreement and are in very good concordance.

**Figure 2. **Numerical validation with reported data [18]

4. Results Analysis

**4.1 Dynamic pressure **

The contour plots of dynamic pressure fields are represented in Figure 3 in the various studied cases. Four models to exit the fluid from the channel. The first outlet is trapped between the upper and lower surfaces of the channel (see Fig. 3a). The second exit is located at the top of the channel in contact with the hot surface (see Fig. 3b). The third exit is located on the lower side of the channel next to the thermally insulated surface (see Fig. 3c).

**Figure 3.** Dynamic pressure fields over baffle plates in channels with different outlet models, Re = 12,000

The final exit is located in the middle of the canal with an estimated 45% of the flow area (see Fig. 3d). The values of these fields were analyzed for a fixed value of the number of Reynolds, of 12,000.

In the first case, the pressure values rise from the top edge of the front of the fin (first obstacle) to the baffle (second obstacle). These values are very important next to the upper edge of the front of the last obstacle and also in contact with the upper surface of the channel to the exit. This increase in pressure is due to the decrease in flow area due to the presence of obstacles. Pressure values in the front and back areas of each obstacle are reduced by the current split into two main streams from left to right and reverse in the recycling cells. In the rest of the studied cases, the pressure is also raised at the outlet of the channel at the bottom, top or center of the channel due to the increase in flow area and decrease in the exit area. The pressure values rise to 86.888 Pa in the second case. While this value decreases by 1.929, 0.920, and 4.788 % in cases A, C, and D, respectively.

**4.2 Streamlines**

In Figure 4, current lines are shown in all studied channels. The fluid enters regularly and then is disturbed by the fin attached to the hot wall, where small recycling cells are formed in the front of the fin due to the decrease in pressure in this area, while the main current moves downward, passing under the same fin, due to a decrease in flow area, due to the presence of the obstacles. The main current at the top edge of the fin is divided into two currents, a reverse in the recycling rings, behind the same fin, as a result of the lower pressure in this area, and the main towards the second obstacle. Part of this current is in the lower front of the baffle, due to the decrease in pressure, while the bulk is directed towards the top of the channel, above the same baffle, to the exit, due to the rise in pressure, because of the recycling rings behind the fin and due to the presence of the baffle.

**Figure 4.** Streamlines for various baffled channel outlet models, Re = 12,000

In the A and C cases, the current comes out in contact with the hot wall, with a decrease in pressure in the second side next to the bottom channel surface and the right side of the baffle. In the second case, the current is diverted from the top of the channel, behind the baffle, towards the bottom of the channel. Pressure decreases in the area between the hot upper wall, the top of the outlet, and the main stream, creating weak cells for recycling, while forming a strong recycling cell next to the right side of the baffle on the bottom of the channel. The same phenomenon is observed when the channel outlet is at its middle end, with the formation of cells for recycling in front and behind each obstacle. In this figure, the stream function value is estimated at about 19.710 Kg/s in the case A. This value increases by 14.057, 8.447, and 14.541 % in the cases B, C, and D, respectively.

**4.3 Contours of velocity magnitude**

Figure 5 shows the magnitude of the velocity field. As shown by the figure, the velocity of the fluid is constant at the channel entrance and is equal in all cases. This speed changes in the rest of the channel areas, where it decreases next to the obstacle, especially in the rear, due to the presence of reverse currents due to the decrease in pressure, while the velocity values are high near the upper side of each obstacle especially near the hot surface due to the augmentation in pressure due to the reduction in the flow area by obstacles, in all studied cases. In addition, the speed is very high at the outlet of the channel, especially in the last three cases, due to the high pressure, due to lack of output.

**Figure 5.** Contour plots of mean velocity fields for various investigated cases, Re = 12,000

Its maximum value is detected in the second case, of 0.406 m/s, while in the other cases; it decreases by 0.969 % in the first case, 0.461 % in the third case, and 2.449 % in the fourth case.

**4.4 Contours of axial velocity **

As expected, the velocity values are very low next to the obstacles, Figure 6. The negative values indicate the presence of recycling cells in the back areas of each obstacle. While the speed is important in the areas between the tip of the obstacles and the inner walls of the channel. Also, the speed is very important from the upper left side of the baffle to the channel exit next to the hot wall. The speed values vary according to the position of the channel outlet, which is the greatest in the third case, estimated at 0.404 m/s, while this value decreases by 0.508, 4.717, and 2.480 % in cases A, B, and D, respectively. In term of the intensity of recirculation, it is very important in the latter case that carries an outlet from the center of the channel end.

**4.5 Contours of transverse velocity **

As shown in Figure 7, there is a decrease in fluid velocity near the upper left side of the fin, and this is in all studied cases, especially in the fourth case, while there is a significant rise in fluid velocity in contact with the upper left edge of the baffle, especially in the second case.

**4.6 Contours of kinetic energy of turbulence**

The kinetic energy values are very low next to the fin in all cases, and at the bottom part of the first channel outlet model, Figure 8. The values of kinetic energy of turbulence rise in a wide space above the baffle to the exit, especially in the second case of the lower outlet. The energy in this case is 0.017 m^{2}/s^{2}, where it drops by 51.938, 43.563, and 41.393 % in the cases A, C, and D, respectively.

**Figure 6.** Fields of axial velocity for various channel outlet models, Re = 12,000

**4.7 Contours of intensity of turbulence**

The contour plots of intensity of turbulence are shown in the following figure. The disturbance decreases at the inlet, while the flow is disturbed as it approaches the fin. This disturbance increases from the upper front of the same fin to the outlet, where it rises in the upper area adjacent to the two surfaces, the top of the baffle, and the bottom of the hot wall. The intensity of disturbance also augments behind the baffle to the end of the channel in all cases except for the second lower part of the first channel exit model. The turbulence intensity value is 7.527, 10.858, 9.018, and 8.246 % in the treated cases A, B, C, and D, respectively. The second outlet model has high disturbance intensity while the first case shows a decrease in intensity, Figure 9.

**Figure 7.** Fields of transverse velocity for different channel outlet models under study, Re = 12,000

**4.8 Skin friction coefficients**

The profiles of normalized skin friction coefficient along the hot surface for various models of the channel exit at Re = 12,000 are shown in Figure 10. The friction value tends to drop considerably to almost zero when it reaches and passes the fin, due to the change in the flow path from the upper side of the channel to the bottom, due to the presence of this same fin. The friction coefficient rises slightly near the right corner of the fin due to the presence of recycling cells in this region. This value decreases and then is missing at the point of separation between the reverse current behind the fin and the direct field above the baffle. The coefficient of friction gradually increases next to the second obstacle until it reaches its maximum value near the right side of the same obstacle. This augmentation is due to the extreme deviation of the flow field due to the presence of the baffle as well as the pressure from the recycling cells located behind the fin on the direct current near the left and upper sides of the baffle.

Finally, the friction values are reduced near the channel exit due to direct field collision with the hot surface where the current is going down, especially if the outlet is in the middle or lower part of the channel. For more analysis, the range of Reynolds was expanded from 12,000 to 30,000. In all investigated cases, the normalized friction factor decreases as the flow rate augments, Figure 11. For all Reynolds value, the friction factor increases in the case of the outlet next to the hot wall of the channel in the first and third cases, while this factor decreases if the flow field is directed towards the lower side, especially near the bottom surface of the channel. The friction factor reaches its maximum value in the third case, so that the outlet is situated at the top region of the channel while the friction decreases in the second case so that the exit is located at the bottom area of the channel.

**Figure 8.** Fields of kinetic energy of turbulence for various channel outlet cases, Re = 12,000

**Figure 9.** Fields of intensity of turbulence for various channel outlet cases, Re = 12,000. Turbulence intensity values in (× 100 %)

**Figure 10.** Profiles of normalized skin friction coefficient along the hot surface for various channel outlet models, Re = 12,000

**Figure 11.** Variation of normalized friction factor with Reynolds number for various channel exit models

**4.9 Nusselt numbers**

The normalized local Nusselt number values reduce next to the fin due to the gradual decrease in contact between the heat exchange surface of the channel and the working fluid due to the downward trend of the flow field while its value rises next to the baffle due to the rapid current in this area, the increase in the thermal gradient and thus the good heat exchange. The heat transfer rate also improves where recycling cells are present in the upper part in the back of the fin. The heat exchange improved in the last three cases compared to the first case with a wide bottom-up outlet, which is an indication of the effect of the exit size and position in the channel. The heat transfer reaches its high value in the fourth case where the outlet area is located in the center of the channel end while it is reduced in the first case. In the figure, the thermal exchange are also larger the unit, which means that there is an improvement in heat transfer compared to the present smooth channel that does not contain fin and baffle. This is a reflection of the role of recycling cells in improving the heat transfer within the channel, Figure 12.

**Figure 12.** Profiles of normalized local Nusselt number along the hot surface for various channel outlet models, Re = 12,000

The variation of the normalized average Nusselt number with Reynolds number values (12,000-30,000) for various channel exit types is shown in Figure 13. It is very clear that the values of average heat transfer rate are proportional to the rise in Reynolds number values. The rise in Reynolds values by increasing flow velocity causes increased disturbance and recycling cells for all studied channels. The fourth case of the narrow-medium output shows an augmentation in heat exchange coefficient, which confirms the previous figure, with a value of 92.433. This value reduces in the first three cases through all suggested Reynolds values. The fourth case shows an increase of about 37.658, 21.356, and 9.348 per cent compared to the first, second, and third cases, respectively, and for the largest number of Reynolds. While the first simple case has the lowest heat exchange value about 26.148, 45.411, and 60.406 per cent compared to the last three cases, respectively. The figure also shows a quantitative improvement of heat in all channels compared to those that do not contain obstacles. This indicates the effectiveness of the fins and baffles, by creating turbulence and forcing the recycling regions for effective friction with the heat exchange surfaces of the channel and thus qualitative and qualitative heat exchange. Also, the size and position of the channel outlet area have a significant impact on heat transfer phenomenon within the channel.

**Figure 13.** Variation of normalized average Nusselt number with Reynolds number for various channel exit models

5. Concluding Remarks

A numerical study has been carried out on the dynamic and thermal behavior of a nanofluid having a constant property and flowing turbulently through a two-dimensional horizontal channel with a rectangular cross-section. The upper wall of the channel was kept at a constant temperature, while it was made sure to maintain the adiabatic condition of the lower surface of the channel. Two obstacles were inserted into the channel; they were fixed to the top and bottom walls of the channel, in a periodically staggered manner to force vortices to improve the mixing, and consequently the heat transfer. The effects of presence of obstacles and nanofluids on the heat transfer phenomenon were investigated in detail. Three different channel outlet models were considered. These effects are investigated with the help of CFD in a steady 2D model. The heat and nanofluid transfer analysis showed a quantitative improvement of heat in all channels compared to those that do not contain obstacles. This indicates the effectiveness of the fins and baffles, by creating turbulence and forcing the recycling regions for effective friction with the heat exchange surfaces of the channel and thus qualitative and qualitative heat exchange. Also, the size and position of the channel outlet area have a significant impact on heat transfer phenomenon within the channel.

References

[1] Izadi M, Pour SMRH, Yasuri AK, Chamkha AJ. (2018). Mixed convection of a nanofluid in a three-dimensional channel: Effect of opposed buoyancy force on hydrodynamic parameters, thermal parameters and entropy generation. Journal of Thermal Analysis and Calorimetry 1-5. https://doi/10.1007/s10973-018-7889-0

[2] Chamkha A, Ismael M, Kasaeipoor A, Armaghani T. (2016). Entropy generation and natural convection of CuO-water nanofluid in C-shaped cavity under magnetic field. Entropy 18: 50-67. https://doi:10.3390/e18020050

[3] Chamkha AJ, Rashad AM, Aly AM. (2012). Non-Darcy natural convection of a nanofluid about a permeable vertical cone embedded in a porous medium. International Journal of Microscale and Nanoscale Thermal 4(2): 99-114.

[4] Chamkha AJ, Aly AM, Al-Mudhaf HF. (2011). Laminar MHD mixed convection flow of a nanofluid along a stretching permeable surface in the presence of heat generation or absorption effects. International Journal of Microscale and Nanoscale Thermal 2(1): 51-70.

[5] Sadripour S, Chamkha AJ. (2018). The effect of nanoparticle morphology on heat transfer and entropy generation of supported nanofluids in a heat sink solar collector. Thermal Science and Engineering Progress 9: 266-280. https://doi.org/10.1016/j.tsep.2018.12.002

[6] Chamkha AJ, Rashad M, Gorla RSR. (2014). Non-similar solutions for mixed convection along a wedge embedded in a porous medium saturated by a non-Newtonian nanofluid. International Journal of Numerical Methods for Heat & Fluid Flow 24(7): 1471-1486.

[7] Chamkha AJ, Selimefendigil F. (2018). MHD mixed convection of nanofluid due to an inner rotating cylinder in a 3D enclosure with a phase change material. International Journal of Numerical Methods for Heat & Fluid Flow. https://doi.org/10.1108/HFF-07-2018-0364

[8] Chamkha A, Abbasbandy S, Rashad AM. (2015). Non-Darcy natural convection flow for non-Newtonian nanofluid over cone saturated in porous medium with uniform heat and volume fraction fluxes. International Journal of Numerical Methods for Heat & Fluid Flow 25(2): 422-437.

[9] Chamkha AJ, Selimefendigil F. (2018). Numerical analysis for thermal performance of a photovoltaic thermal solar collector with SiO2-water nanofluid. Appl. Sci. 8: 2223. https://doi:10.3390/app8112223

[10] Chamkha AJ, Dogonchi AS, Ganji DD. (2018). Magnetohydrodynamic nanofluid natural convection in a cavity under thermal radiation and shape factor of nanoparticles impacts: a numerical study using CVFEM. Appl. Sci. 8: 2396. https://doi.org/10.3390/app8122396

[11] Chamkha A, Gorla RSR, Ghodeswar K. (2011). Non-similar solution for natural convective boundary layer flow over a sphere embedded in a porous medium saturated with a nanofluid, Transp Porous Med 86: 13-22. https://doi.org/10.1007/s11242-010-9601-0

[12] Chamkha AJ, Abbasbandy S, Rashad AM, Vajravelu K. (2012). Radiation effects on mixed convection over a wedge embedded in a porous medium filled with a nanofluid, Transp Porous Med 91: 261-279. https://doi.org/10.1007/s11242-011-9843-5

[13] Minea AA, Moldoveanu MG. (2017). Studies on Al2O3, CuO, and TiO2 water-based nanofluids: a comparative approach in laminar and turbulent flow. Journal of Engineering Thermophysics 26(2): 291-301.

[14] Akhatov JS, Mirzaev SZ, Halimov AS, Telyaev SS, Juraev ET. (2017). Study of the possibilities of thermal performance enhancement of flat plate solar water collectors by using of nanofluids as heat transfer fluid. Applied Solar Energy 53(3): 250-257.

[15] Dehaj MS, Mohiabadi MZ. (2019). Experimental investigation of heat pipe solar collector using MgO nanofluids, Solar Energy Materials and Solar Cells 191: 91-99.

[16] Kaya H, Arslan K, Eltugral N. (2018). Experimental investigation of thermal performance of an evacuated U-Tube solar collector with ZnO/Etylene glycol-pure water nanofluids. Renewable Energy 122: 329-338.

[17] Mirzaei M, Hosseini SMS, Kashkooli AMM. (2018). Assessment of Al2O3 nanoparticles for the optimal operation of the flat plate solar collector. Applied Thermal Engineering 134: 68-77.

[18] Demartini LC, Vielmo HA, Möller SV. (2004). Numeric and experimental analysis of the turbulent flow through a channel with baffle plates. Journal of the Brazilian Society of Mechanical Sciences and Engineering 26(2): 153-159.

[19] Belhadj A, Bouchenafa R, Saim R. (2018). Numerical investigation of forced convection of nanofluid in microchannels heat sinks. Journal of Thermal Engineering, Research Article 4(5): 2263-2273.

[20] Ajeel RK, Salim WS-IW, Hasnan K. (2019). Design characteristics of symmetrical semicircle-corrugated channel on heat transfer enhancement with nanofluid, International Journal of Mechanical Sciences 151: 236-250.

[21] Menni Y, Azzi A, Chamkha AJ, Harmand S. (2018). Effect of wall-mounted V-baffle position in a turbulent flow through a channel: Analysis of best configuration for optimal heat transfer. International Journal of Numerical Methods for Heat & Fluid Flow. https://doi.org/10.1108/HFF-06-2018-0270

[22] Chamkha AJ, Ismael MA. (2014). Natural convection in differentially heated partially porous layered cavities filled with a nanofluid. Numerical Heat Transfer, Part A 65(11): 1089-1113. https://doi.org/10.1080/10407782.2013.851560

[23] Launder BE, Spalding DB. (1974). The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 3: 269-289.

[24] Patankar SV. (1980). Numerical heat transfer and fluid flow. McGraw-Hill, New York, NY.