On the validity of a numerical model predicting heat and mass transfer in porous square cavities with a bottom thermal and solute source: case of pollutants spreading and fuel leaks

– The present work refers to the study of natural convection into a conﬁned porous medium, driven by cooperating thermal and solutal buoyancy forces. The side walls are maintained at a uniform temperature and concentration, lower than that of a heat and solute source, which located at the center of the bottom wall, the rest of the horizontal walls are kept insulated. The physical model for the momentum conservation equation makes use of the Brinkman extension of the classical Darcy equation, the set of coupled equations is solved using the ﬁnite volume method and the SIMPLER algorithm. To account for the eﬀects of the main parameters such the buoyancy ratio, the Lewis and porous thermal Rayleigh numbers, as well as the source length, heat and mass transfer characteristics are widely inspected and then, new powerful correlations are proposed, which predict within ± 1% the numerical results. Note that the validity of the used code was ascertained by comparing our results with experimental data and numerical ones already available in the literature.


Introduction
Over the past four decades, double-diffusive natural convection analysis into a porous medium has been the subject of a very intense research activity, due to the importance of related industrial and contemporary technological applications such grain storage installation; geothermal energy resources; petroleum reservoirs; pollutant dispersion in aquifers; fibrous insulating materials, electrochemical processes, and some modes of assisted oil recuperation [1].
With both temperature and concentration gradients present to drive the fluid flow, an increased number of transport configurations was possible, with parallel or perpendicular gradients, and the body forces augmenting or opposing [2].
Several experimental [3][4][5][6] and numerical investigations [7][8][9][10] caused by horizontal thermal and solutal gradients have been done. For certain conditions, a multicellular flow was observed experimentally or predicted numerically. In contrast, a little attention was dedicated to the case of vertical gradients [11][12][13][14]. a Corresponding author: ragui-karim@live.fr To predict the heat and mass rates into such configurations, some powerful correlations were proposed: in 1987, Trevisan and Bejan [15] projected the thermal and solutal transfer characteristics into a porous square as a function of pertinent parameters such: Lewis and porous thermal Rayleigh numbers; the cavity aspect ratio and the buoyancy ratio as well. In 1990, Lin et al. [16] proposed new Nusselt and Sherwood correlations as a function of the thermal Grashof number, and that for small values of the buoyancy ratio (|N | < 5). A few years later, Bennacer [17] suggested a general correlation for the mass transfer, which can be used in a wide range of porous thermal Rayleigh number, buoyancy ratio, and Lewis number as well.
In these works and many others, the problem of double-diffusive convection was commonly inspected using the entire wall length of the enclosure as a heat (and solute) source, what made these correlations only available using this condition. Thus, the main purpose of the current work is to complete these cases by investigating the double-diffusive natural convection within a cold (and less concentric) sides porous square, which including a heat and solute source located at the center of its bottom wall. This especial geometry is an interesting work

Subscripts
h Hot c Cold for many industrial applications, especially in predicting pollutants spreading and fuel leaks in the Refining and Petrochemical industry.

Problem statement and mathematical formulation
The studied configuration, shown in Figure 1, consists of a square porous cavity (of impermeable boundaries) saturated by a binary fluid. The side walls are subjected to uniform temperature and concentration lower than that of a heat and solute source located at the center of the bottom wall, when zero heat and mass fluxes are imposed at the rest parts of the horizontal walls. The binary fluid is assumed to be Newtonian and to satisfy the Boussinesq approximation; the flow is incompressible, laminar, bi-dimensional and in the steady state. On the other hand, the porous medium is supposed to be isotropic, homogeneous and in thermodynamic equilibrium with the fluid. The Soret and Dufour effects are assumed to be negligible. The density variations upon temperature and concentration are described by the state equation: where β T and β C are thermal and solutal expansion coefficients.
The dimensionless conservation equations, describing transport phenomena in the square, can be written as: where Da is the Darcy number, Gr T (=Ra T /P r) is the thermal Grashof number, N is the buoyancy ratio, P r and Sc are the Prandtl and Schmidt numbers. The appropriate boundary conditions used to solve the above system (3 to 7) can be resumed as follows: Noted that average heat and mass fluxes across the cavity' side walls, as well as the bottom source, can be expressed in a dimensionless form by the Nusselt and Sherwood numbers such as:

Numerical procedure and code validation
The governing conservation equations are discretized in space using the finite volume approach, when the convection-diffusion terms are treated with a Power-Law scheme. The resulting algebraic equations, with the associated boundary conditions, are then solved using the line by line method. As the momentum equation is formulated in terms of the primitive variables (U , V and P ), the iterative procedure includes a pressure correction calculation method, namely SIMPLER [18] to solve the pressurevelocity coupling. Compared to other velocity-pressure coupling approaches such SIMPLE and SIMPLEC, the  adopted approach is proved to be faster (about 30 to 50% fewer iterations). The convergence criterion for the temperature, the concentration, the pressure, and the velocity as well, is given as:

Present predictions
where both m and n are the numbers of grid points in Xand Y -directions, respectively; ξ is any of the computed field variables and t is the iteration number. The performance of the using code via the doublediffusive natural convection problem into a porous medium is accomplished by comparing predictions with other numerical results and experimental data, and by verifying the grid independence of the present results.
First, the present results are consistent with previous computations namely those of Hadidi et al. [19] which deal with double-diffusive natural convection into a square enclosure formed by two different porous media. By taking into account the same hypotheses, Figure 2 shows the developed streamlines, isotherms and iso-concentration plots into the square. Tables 1 and 2 illustrate the comparison of mean Nusselt and Sherwood numbers computed with various values of the buoyancy ratio N and the porous thermal Rayleigh number Ra*, respectively.  As we can see, the present results and those of Hadidi et al. [19] are in excellent agreement with a maximum discrepancy of about 2%.
Then, to check the numerical code validity with experimental results, those obtained by Weaver and Viskanta [20] for an Ethanol/Nitrogen and Ethanol/Helium binary fluids have been selected. Figures 3 and 4 display the comparison between the experimental data and both the numerical Weaver and Viskanta predictions and the present ones in term of velocity contours. Once again, the numerical results show a good qualitative concordance with the experimental data and a great agreement with the numerical Weaver and Viskanta predictions.
In order to determine a proper grid for the numerical simulations, a grid independence study is conducted. Several mesh distributions ranging from 91 2 to 151 2 are tested; the mean Nusselt and Sherwood numbers of the porous square, previously shown in Figure 1, are presented in Figure 5 as a function of the above uniform grids. Through the latter, a 121 2 uniform grid is found to be adequate for a grid independent solution. However, a fine structured mesh of 131 2 will be used to avoid roundoff error for all other calculations in our investigation.

Results and discussion
In this paper, the predictions are performed for various values of the Lewis number (10 ≤ Le ≤ 300), the buoyancy ratio (0 ≤ N ≤ 50), and the source length (0.20 ≤ W ≤ 0.80). Both; the Darcy and Prandtl numbers are fixed at 10 −3 and 10, respectively. As the square cavity is formed by a solid matrix, our results are presented using a porous thermal Rayleigh number (100 ≤ Ra * ≤ 2000) rather than the thermal one (Ra T ).

Impact of Lewis number and buoyancy ratio
Starting our investigation by a fixed source length equal to 0.60, Figure 6 displays the streamlines, the isotherms and the iso-concentration plots, according to Lewis number and that, for a buoyancy ratio and a porous thermal Rayleigh number equal 10 and 100, respectively.
For these patterns, two symmetrical counter-rotating cells are formed into the porous square, the fluid flows up in the centre on the enclosure and down all along the side walls at a relatively high speed with a low Lewis number, compared to the fluid speed with a high one.
The corresponding isotherms are found to be nearly the same, and barely changed, demonstrating a dominated-conduction heat transfer regime, unlike the isoconcentration plots where the latter becomes more distorted by increasing the Lewis number. Then, mass boundary layers are created all along the active walls and become thinner with increasing Lewis values. Also a solutal plume appears in the center of the enclosure indicating then, the domination of the diffusion mechanism.
For more clarification, the hydrodynamic, thermal and solutal behaviors are summarized in Figure 7 where the vertical velocity component, the temperature and the concentration profiles, along the horizontal mid-plane of the porous enclosure, are presented. In fact, the examination of the magnitude of V -velocity at various Lewis numbers confirms the results previously obtained from analysing the streamlines. The vertical velocity profiles verify the existence of clockwise and anti-clockwise circulating cells inside the porous square. Besides, the decrease in the velocity magnitude with increasing Lewis number is a clear indication of a weak buoyant flow at high values of the latter. In the same way, the mass transfer mechanism is expected to be more pronounced by increasing Lewis number, as illustrated via the φ profiles, whilst the conduction is responsible for the heat transfer as shown by the θ profiles.
Regarding the buoyancy ratio effect, Figure 8 shows the fluid dynamic, thermal and solutal behaviors, in terms of streamlines isotherms and iso-concentrations, for various values of N (the Lewis number is fixed now at 10). Once again two symmetrical counter-rotating cells are formed into the porous square. Its absolute streamfunction value is an increasing function of the buoyancy ratio as the fluid motion becomes stronger. The corresponding isotherms are more pronounced by increasing values of N , indicating the transition from a pseudoconductive regime (at N = 1) to a primary-convective one (see for instance N = 30). As to the concentration field which is very affected by increasing buoyancy ratio, thin mass boundary layers are created all around the active walls, denoting large concentration gradients along these surfaces. The thickness of the latter (mass bound- ary layer), is a decreasing function of enhanced values of N which may explain the domination of the diffusion mechanism compared to the convection one. Referred to Figure 9, the velocity magnitude and φ profiles depend not only to the Lewis values but to the buoyancy ratio ones as well. As we can see, an important buoyant flow and mass transfer is also produced using a high value of buoyancy ratio compared to a lower one. Note, at the end, that the effect of the buoyancy ratio is still barely observed through the θ profiles, as the convection mechanism still less pronounced. Figure 10 summarizes a series of predictions designed to document the impact of Lewis number and buoyancy ratio on the mass transfer rate. Through the latter, and for each value of Lewis number, the mean Sherwood number is found as an increasing function of the increase buoyancy ratio which can be linked to the flow intensity as reported previously (see Fig. 9). For a given N , the mass transfer still increases with Lewis number. Indeed, the Lewis number increases through the Schmidt number, since the Prandtl number is already fixed at 10, what reduces the thickness of the solutal boundary layer and leads to important mass transfer.
About the heat transfer, Figure 10b, the previous description is shown to hold only for the buoyancy ratio. In other terms, the mean Nusselt number is also an increasing function of the latter (N ) as the V -velocity optimum continuously increases; what decreases the thickness of the corresponding boundary layer and leads to the primary-convective regime (see Fig. 9). With Lewis number all the opposite is observed, the heat transfer rate decreases as the thermal boundary layer thickness is a function of Le 1/2 [21].

Impact of porous thermal Rayleigh number
To complete the previous results, the evolution of heat and mass transfers with porous thermal Rayleigh number (Ra*) is also discussed (Fig. 11). The Lewis number is taken as 10, which represents the case of hydrocarbon fuel [22], when the buoyancy ratio value is ranged between 0 and 30. As for the latter, the increase in the porous thermal Rayleigh number (through Grashof number) improves significantly the transfer rates, especially with high values of the buoyancy ratio. This behavior confirms the fact that the global buoyancy term in the momentum equation, Gr T (θ + N φ), increases with both; the porous thermal Rayleigh number (Da Gr T P r) and the buoyancy ratio (N ), promotes consequently the flow intensity and so the overall transfer. These primary observations allow us to conclude the proportional relation between the thermosolutal transfer, and the governing parameters such Le, N and Ra*.

Source length effect and proposed models
For fixed values of Lewis number and buoyancy ratio, (Le = N = 10), the impact of the source length on the fluid dynamic, thermal and solutal behaviors into the porous square is displayed in Figure 12. As before, the two counter-rotating cells created by the fluid motion are formed inside the square. By increasing the source length the counter-rotating cells get stronger as the fluid motion is more pronounced. On the other hand, the development of thin thermal and solutal boundary layers over the active walls, as shown in Figures 12b and 12c, is a clear indication of important gradients near these surfaces, especially when the source length is about 0.80.
For the entire range of governing parameters, and by taking into account the expression of the global buoyancy term in the momentum equation Gr T (θ + N φ) our Nusselt and Sherwood predictions obtained for various values of the source length can be represented as a function of dimensionless groups such Ra * N/Le and Ra * Le(N + 1) as shown in Figure 13. Through the latter, the proportional relation between the heat and mass transfers and the source length is well approved, what leads us eventu- ally to the following global correlations (Fig. 14): These powerful correlations, found to predict the numerical results within ±1%, may count as a complement to previous researches done in the case of square porous enclosures saturated by a Newtonian fluid, and still available for a Lewis number ranging from 10 to 300, a buoyancy ratio taking between 0 and 50, a porous thermal Rayleigh number between 100 and 2000, and a source length between 0.10 and 0.90.

Conclusion
The analysis of double-diffusive natural convection phenomenon inside a square Darcy-Brinkman porous enclosure having a bottom heat and solute source was realized numerically through this paper. By taking into account the effects of pertinent parameters such the buoyancy ratio, the source length, the Lewis and the porous thermal Rayleigh numbers, new powerful Sherwood and Nusselt correlations which display the heat and mass transfer rates in such geometry were proposed with an accuracy of about 1%. These unique correlations may count as a complement to previous researches done into square porous enclosures and might prove particularly useful in verifying any instability analysis which might be put forth in the future.