Abstract
In this paper, thermo-diffusion (Soret effect) and diffusion-thermo (Dufour effect) effects on double-diffusive natural convection induced in a horizontal Brinkman porous layer with a stress-free upper boundary are investigated. The cavity is filled with a binary fluid and subjected to uniform fluxes of heat and mass on its long sides. An analytical solution based on the parallel flow approximation is developed for the problem considered in order to allow prompt determination of the thresholds of stationary and finite amplitude solutions and also heat and mass transfer characteristics. The analytical solution is validated numerically by using a finite difference method. The combined effects of the Soret and Dufour parameters, the thermal Rayleigh number, the buoyancy ratio, and the Darcy number on the flow intensity and heat and mass transfer are illustrated graphically, and some particular behaviors observed are discussed. The analytical solution proves the existence of different regions in the buoyancy ratio-Dufour parameter plane, corresponding to different parallel flow behaviors. The number, the location, and the extent of these regions, which are impossible to predict numerically, depend strongly on Soret and Dufour parameters. The effect of thermo-diffusion and diffusion-thermo on flow intensity and heat and mass transfer is found to be important.
1. Introduction
Natural convection involving binary mixtures in the porous and fluid media is still a relevant field of investigation since it occurs in wide ranges of applications in many engineering problems and natural fields such as geophysics, oil reservoirs, storage of nuclear wastes, operation of solar ponds, chemical reactors, migration of mixtures in fibrous insulation, and metal manufacturing processes. This phenomenon occurs due to temperature differences, concentration differences, or by combination of these two differences [1]. Convection in which the buoyancy forces are due to both temperature and chemical concentration gradients is referred to as thermosolutal or double-diffusive convection.
The fluid flow behavior driven by temperature and concentration gradients has attracted the interest of researchers worldwide through the decades [2]. Different aspects of the problems involving double-diffusive natural convection have been addressed previously [3–10] in the absence of Soret and Dufour effects and neglected also in many other studies since they are of smaller order of magnitude compared to the effects described by Fourier’s and Fick’s laws. Nevertheless, the literature review shows the existence of experiments mainly focused on the determination of the Soret coefficients for mixtures [11–13] or on theoretical aspects to determine the optimal combinations of the governing parameters leading to the separation of species under the effect of thermo-diffusion [14–16]. Other aspects of double-diffusive convection problems in the presence of the Soret effect have been conducted in shallow horizontal cavities by Bourich et al. [17–19]. The results reported in these investigations are concerned with convection thresholds [17], reversal of the horizontal concentration gradient [18], and a comparison study between the cases where the enclosure is filled with a clear binary mixture or a saturated porous medium [19]. More recently, Hasnaoui et al. [20] studied numerically double-diffusive natural convection in an inclined enclosure with heat generation and Soret effect using a hybrid lattice Boltzmann-finite difference method. They show that the negative Soret parameter combined with high internal heat generation and a relatively high inclination is important when the objective is to maintain the fluid at a high concentration of species.
Despite the fact that Dufour and Soret effects are of a smaller order of magnitude, they however may play a significant role in double-diffusive flow processes. The list of published works taking into account both these two effects is limited. Ojjela et al. [21] studied numerically the influence of an external magnetic field and radiation in free convective Jeffrey fluid flow between parallel porous plates in the presence of Soret and Dufour effects. The governing equations were solved numerically using the shooting method with the 4th Runge–Kutta scheme. The numerical study conducted by Zaho et al. [22] focused on heat and mass transfer generated by natural convection in a porous medium saturated by a viscoelastic fluid and subjected to a magnetic field. Kefayati [23] studied the impact of Soret and Dufour effects on entropy generation due to double diffusion in a square cavity filled with non-Newtonian power-law fluid by using a finite-difference lattice Boltzmann method. He found that the Dufour parameter results in the enhancement of the total irreversibility. A similar study was performed by the same author [24, 25] in the case of inclined cavities. Wang et al. [26] investigated numerically the behavior change of solutions accompanying the increase of the buoyancy ratio. Ren and Chan [27] studied transient double-diffusive convection under the effect of numerous parameters of control including thermo-diffusion and diffusion-thermo. Balla and Naikoti [28] analyzed numerically double-diffusive free convection heat and solute transfer in an inclined square porous cavity saturated with a fluid in the presence of Soret and Dufour effects. Both these effects were considered by Al-Mudhaf et al. [29] who studied numerically unsteady double-diffusive natural convection inside trapezoidal inclined enclosures filled with an isotropic porous medium and submitted from one side to nonconstant distributions of temperature and concentration. Double-diffusive natural convection of non-Newtonian power-law fluids in an open cavity subjected to a horizontal magnetic field in the presence of Soret and Dufour effects was analyzed by Kefayati [30, 31]. Among others, he showed that the rise of Soret and Dufour parameters enhances the entropy generations due to heat transfer and fluid friction.
The present study is devoted to analytical and numerical investigation of Soret and Dufour effects within a horizontal Brinkman layer with a stress-free upper boundary. The main objective is to bring out the combined effects of Soret and Dufour parameters on thresholds of convection, fluid flow, and heat and mass transfer characteristics. For this shallow enclosure, the parallel flow approximation allows to predict the critical Rayleigh numbers for the onset of supercritical and subcritical convection. In addition, the simplifications allowed by this approximation leads to identify different regions with various parallel flow behaviors. The adopted Brinkman model permits covering the limits of Darcy and the pure fluid media. The study focuses on the discussion of combined effects of Rayleigh number, Soret parameter, Dufour parameter, buoyancy ratio, and Darcy number.
2. Mathematical Formulation of the Problem
The domain under study, sketched in Figure 1, is a two-dimensional horizontal porous cavity saturated with a binary mixture. The height of the cavity is , its length is (), and its upper horizontal surface is nondeformable and stress-free, while the remaining boundaries are assumed rigid. The short boundaries (vertical walls) of the cavity are maintained adiabatic and impermeable to mass transfer, while its long horizontal walls are subjected to uniform fluxes of heat, , and mass, . The porous matrix is assumed isotropic and homogeneous, and the Brinkman–Hazen–Darcy model is adopted. The binary fluid that saturates the porous matrix is assumed to be Newtonian, and it is modeled as a Boussinesq incompressible fluid whose density varies linearly with the temperature and the concentration as follows:where is the binary fluid density at temperature and concentration and and are the thermal and solutal expansion coefficients, respectively.
The dimensional equations governing this problem are written as
The reference scales for time, length, velocity, and pressure are , , , and , respectively. The parameters σ and α are defined by and .
The dimensionless temperature and mass fraction are given by and where and .
The governing equations that are solved numerically are based on the vorticity-stream function formulation. In their dimensionless forms [26], they are obtained as follows:
The boundary conditions associated to the present problem are
The examination of the dimensionless governing equations, equations (6)–(10), and their associated boundary conditions, equations 11a–11c, shows that the control parameters are the aspect ratio of the enclosure, A, Darcy number, , the Dufour number, , the Lewis number, , the thermal Rayleigh number, , and the Soret number, . These parameters are defined as follows:
The Nusselt and Sherwood numbers are, respectively, given by the following expressions:where and are the temperature and concentration differences, evaluated at .
3. Numerical Solution
The numerical solution of the governing equations together with the associated boundary conditions, equations (6) to (11a)–(11c), is obtained using a finite-difference method. The vorticity equation, equation (6), the energy equation, equation (7), and the equation of conservation of species, equation (8), are written in their transient forms since the iterative procedure is performed using the alternate direction implicit (ADI) method for these equations. The stream function field is derived from the resolution of equation (9) using the point successive over-relaxation method (PSOR). The numerical results reported in this paper were obtained with a globally nonuniform grid; the latter is uniform and tight near the confining rigid walls and the free surface to capture the flow details near these boundaries and uniform but coarser elsewhere. The vorticity values on the rigid boundaries are calculated using Wood’s relation [32], while zero is attributed to the vorticity at the free surface. Further details on the numerical method and its validation were given in a previous study by Amahmid et al. [5].
Numerous preliminary numerical tests have been performed to determine the minimum aspect ratio from which the assumption of the parallel flow is corroborated. Within the ranges of variations of the parameters considered in this investigation, it was found that the numerical results are independent of the aspect ratio from the threshold . Thereby, the numerical results reported here were obtained with and a grid of . The effect of the grid size on the results is illustrated in Table 1 (analytical results are also included in this table). It can be seen from this table that the results obtained with the grid differ by less than 0.2% from those corresponding to the finest grid . In addition, the grid reproduces the analytical results with relative errors lower than 0.3%.
4. Results and Discussion
4.1. Parallel Flow Analysis
Based on the remarks (verified numerically) allowing simplifications in the case of shallow cavities (), the flow is parallel to the long walls of the cavity and the temperature and concentration fields are characterized by a linear and horizontal stratification if we exclude the edges’ effects. Thus, the stream function, the temperature, and the concentration can be approximated as follows:the parameters and are the unknown constant temperature and concentration gradients in the horizontal direction, respectively.
The investigation by Cormack et al. [33] counts among the early works developing these approximations based on the theory of asymptotic developments for . The steady form of the governing equations (6)–(8) while using these approximations can be reduced to the following set of ordinary differential equations:where .
The corresponding boundary conditions on the horizontal walls become
The analytical integration of equations (15)–(17), together with the associated boundary conditions (18) and (19), leads to the following solutions, respectively, in terms of stream function, horizontal velocity, temperature, and concentration:
The parameter is the stream function at the center of the cavity, obtained analytically aswith
The parameters in the previous equation are obtained as
The parameter is the normalized Rayleigh number, and , , and are the functions, respectively, given by the following analytical expressions:
The expressions of are given by
The parameter and the functions and are dependent on the Darcy number as follows:
According to equation (21), the velocity cannot be zero in the range , which means that only unicellular flows are possible. In addition, the constants and are obtained by considering mass and thermal balances across any transversal section of the porous layer [34].
Substituting equations (22) and (23) into equation (13), the respective analytical expressions of Nusselt and Sherwood numbers are obtained as
By substituting the solutions given by equations (21)–(23) into equations (31) and (32), the final analytical expressions of and are obtained as follows:with and .
The analysis of equation (24) indicates that up to five different steady-state solutions are possible including the diffusive regime . The four other convective solutions depend on the signs ± inside and outside the brackets. In fact, the sign ± outside the brackets refers to counter-clockwise/clockwise circulation, while it refers to convective stable/unstable solutions within the brackets. It is to specify that the “unstable” solutions refer to the analytical solutions that could not be obtained numerically. In addition, from a mathematical point of view, equation (24) shows the existence of supercritical and subcritical bifurcations [19]. The supercritical bifurcation is characterized by the transition from the quiescent state to the convective regime through zero flow amplitude , while the subcritical bifurcation is characterized by the onset of motion through finite amplitude .
These parallel flow solutions exist only when the two following conditions are satisfied: and . The satisfaction of these conditions allows the subdivision of the plane into different regions with specific behaviors, depending on the sign of the parameter P defined by equation (26).
4.1.1. Analytical Delineation of the Different Regions
On the basis of the analytical solution, up to six regions with specific behaviors are identified in the plane. These regions are illustrated in Figures 2–4 for , , and , respectively.
The region (1) is characterized by the rest state, which means that the parallel flow solution is not possible in this region whatever the value of . Different conditions must be verified to delineate this region in the plane, depending on the sign of .(i)For (ii)For (iii)For
For the region (2), only the supercritical flow is possible. The corresponding supercritical Rayleigh number, characterizing the transition from the rest state to the convective regime through the zero flow amplitude, is given by
The conditions that should be verified to delineate this region in the plane are dependent on the sign of .(i)For with(ii)For (iii)For
For the region (3), the characteristics are similar to those in region (2). The stable parallel flow solution exists from the supercritical Rayleigh number given by equation (38). Even if the regions (2) and (3) correspond to the same type of convection (stationary convection), they are however characterized by a difference in terms of asymptotic evolutions at large values of the Rayleigh number.
The conditions that should be verified for this region are dependent on the sign of as follows:(i)For with and (ii)For (iii)For
For the region (4), both stable and unstable solutions are existing from the subcritical Rayleigh number, but the unstable branch disappears at some threshold of . In this region, the resulting solutions correspond to a subcritical bifurcation for which the onset of motion occurs through a finite amplitude.
The critical Rayleigh number, , above which the parallel flow exists, is obtained from the condition and its expression is given bywhere .
The conditions that should be verified for the region (4) in the plane are dependent on the sign of as follows:(i)For (ii)For (iii)For
For the region (5), the behaviors are similar to those of region (4) but the existence range of the unstable branches of the former is extended towards infinite .(i)For (ii)For (iii)For
Finally, the region (6) is characterized by the same behavior as that of region (4); the only difference is observed in terms of their asymptotic behaviors. The Soret number appears explicitly in the conditions delineating this region as follows:
The expression of was already defined for region (2).
4.2. Effect of
The effect of on the flow intensity, , the Nusselt number, , and the Sherwood number, , is illustrated in Figures 5(a)–5(c) for , , , and different combinations of corresponding to different regions in this plane. The choice of allows to cover the six identified regions. Moreover, it can be observed from these figures that the stable analytical solution (solid lines) is in a very good agreement with the numerical results depicted by full circles and obtained by solving the full governing equations.
(a)
(b)
(c)
For the combination corresponding to the region (2), only the stable branch exists. The supercritical convection starts from the rest state (characterized by , , and ) at . The evolutions of , , and are characterized by asymptotic behaviors, and their respective asymptotic limits are 2.31, 3.45, and 1.64, reached for sufficiently large values of . Note that for this combination of , the asymptotic value of indicates that diffusion plays an important role in the mass transfer. At large values of , the analytical expressions of the asymptotic expressions of , , and corresponding to the stable solution are obtained by the following expressions:with and .
For the combination , illustrating the region (3), the behavior is also supercritical and the convection starts from the rest state at . The evolution of vs. is characterized by a monotonous increase with different rates that depend on the combination . For , the analytical expression of is reduced to
Equation (58) shows that varies as , while the evolutions of and are characterized by asymptotic behaviors towards the same limit, which is 4.49 (see Figures 5(b) and 5(c)). This value is deduced from the asymptotic expression given by equation (59) that is common for both and at large :
The region (4) is illustrated by the combination in Figures 5(a)–5(c). Unlike the previous combinations of , for which only the stable branches exist, both stable and unstable solutions exist for this combination in the range . This means that the parallel flow solution starts at , and the nascent flow corresponds to a state clearly different from the purely diffusive regime characterized by and . For the unstable branch, /( and ) decreases notably/(weakly) by increasing and disappears for .
For the stable branch, the evolution of versus is very similar to that described for the combination but characterized by a delay in terms of . At large values of the last parameter, the expression of is the same as that established for the region (3) and given by equation (58). For the region (4) also, the evolutions of and are characterized by asymptotic behaviors at large , leading to the common asymptotic value 4.49, obtained with equation (59) for both and .
The results corresponding to the region (5) are illustrated in Figures 5(a)–5(c) by the combination . For this region, both stable and unstable solutions exist for for the considered combination of . In addition, the evolution of the stable branch is similar to that described for region (4). By increasing , the unstable branches are characterized by slow decreases towards asymptotic limits for , , and .
The region (6) is illustrated by the combination and shows a behavior different from those already described. More precisely, both stable and unstable solutions for this region exist in the very short range . In addition, at large values of , the quantities , , and tend, respectively, to the asymptotic limits 9.70, 5.827, and 4.155. It should be outlined that for the region (6) the behaviors are similar to those already described for the regions (4) and (2), respectively, at low and large values of . Consequently, the analytical expressions of the asymptotic values of , , and are given by equations (55)–(57), respectively. Besides, for the region (6), the Nusselt number becomes negative in the range . This change in the sign is exclusively attributed to the negative value of the Dufour number and the small temperature differences induced in some ranges of .
4.3. Combined Effects of Dufour and Soret Numbers
This section is devoted to assess the simultaneous effects of thermal-diffusion (Soret effect) and diffusion-thermo (Dufour effect) on fluid flow and heat and mass transfer characteristics. Soret and Dufour effects are defined as the mass flux caused by a temperature difference and the energy flux engendered by concentration difference, respectively [35]. The Soret and Dufour parameters can be widely varied by changing the mean solute concentration or mean temperature of the system.
The effect of the Dufour number on the flow intensity, the Nusselt number, and the Sherwood number is illustrated in Figures 6(a)–6(c) for = 500, , , and various values of . These figures show that, regardless of the value of , there exists a critical value of above which convection disappears. This critical value increases with ; it rises from 0.124 to 0.897 when goes from −0.9 to 0.9.
(a)
(b)
(c)
The effect of on , , and may be drastically affected by the change of . More specifically, we can observe two different behaviors depending on whether the value of is above or below critical ranges, which are so narrow that they appear as nodes in Figures 6(a)–6(c). In these critical ranges of , the effect of thermo-diffusion becomes insignificant. Another mystery of these critical ranges of lies in the fact that they are not the same for , , and (compare [0.6, 0.732] for , [0.016, 0.042] for , and [−2.29, −2.25] for ). Due to their very narrowness, the critical ranges will be termed as critical nodes in the following discussion. The critical node corresponding to the curves of is seen to be largely shifted towards negative values of . The interaction between convection, themodiffusion, and diffusion-thermo is so complex that, without analytical calculations, it was not possible to predict such specific behaviors attributed to the combined effects of the governing parameters. It is worth mentioning that these behaviors were obtained analytically and validated numerically. By focusing on the stable solution (validated numerically), it can be seen that below the critical nodes, the quantities , , and are increasing functions of . Above the critical node corresponding to , the latter decreases by increasing . However, immediately near the critical nodes, and exhibit a behavior similar to that of and their evolutions are characterized by complex variations far above these critical nodes. Another important difference observed by comparing heat and mass transfer from Figures 6(b) and 6(c) is the fact that decreases by increasing , while increases/decreases by increasing for . Furthermore, the most important variations observed for when is varied are located below/above the critical nodes. Quantitatively, at for instance, by increasing from to 0.9, is multiplied by a factor of 4.1 while is increased by only 11.44%. For , by increasing from to 0.9, is divided by 3.2 while undergoes a relative decrease of about 14.82%.
The combined effects of and on the temperature and concentration fields are illustrated in Figure 7 for RT = 500, N = −0.2, Le = 1.1, , , and . For the combination , it can be seen from Figure 7 that the isotherms are very similar to the isoconcentrations. Even the quantitative comparisons (results not presented) showed that the temperature field is not very different from the concentration one. The temperature and concentration profiles obtained at midwidth of the cavity and plotted in Figure 8 show that the temperature/concentration gradients are important in the vicinity of the horizontal wall, compared to the those (i.e., gradients) around the horizontal median. By changing the combination from to , we can observe an accentuation/attenuation of the isoconcentration/isotherm distortions. The isotherms are nearly straight lines tilted with respect to the horizontal direction, which means that the temperature field in the cavity (if we except the vicinity of the short walls where the end effect is important) is quasi-linear. The temperature profile at midwidth of the cavity (Figure 8) confirms this behavior. It should be noted that the quasi-linearity exhibited by the temperature field for cannot be explained by the dominance of the diffusive regime. In fact, the diffusive regime leads to linear profiles with the isotherms being parallel to the horizontal walls. However, the isotherms obtained for are strongly tilted with respect to the horizontal direction. This means that this behavior is due to the interaction between convection, conduction, and diffusion-thermo (i.e., Dufour effect). The changes observed in the temperature/concentration field when passes from to are similar to those observed in the concentration/temperature field when passes from to . Hence, we can observe a quasi-linear behavior for the concentration field at . The combination leads to the most important variations for the temperature and the concentration among the four combinations considered in this comparison. Quantitatively, by changing from to , the temperature difference between the horizontal wall at a fixed cross section passes from 0.1769 to . Thus, the temperature difference obtained with the combination is at least 8.6 times higher than that obtained for the other combinations. A result with the same order is also registered for the concentration difference.
(a)
(b)
4.4. Effect of Darcy Number
The Darcy number, , measures the relative importance of the permeability of the porous medium as it is proportional to the latter. It is very small for the well-packed porous media and relatively large for sparsely packed ones.
The effect of on , , and is illustrated in Figures 9(a)–9(c) for , , and different combinations corresponding to different regions. In these figures, it can be seen that the numerical results are in excellent agreement with the analytical ones, corresponding to the stable branches. In addition, regardless of the region considered, the parallel flow solution vanishes when exceeds a critical value, , that depends on the combination . For the stable branches, the quantities , , and decrease by increasing the permeability of the porous medium. Similar behaviors were reported in previous studies [5, 36]. Note that the variations versus of the flow intensity corresponding to the unstable branches are remarkably strong in specific ranges of the latter parameter. However, the effect of on and corresponding to the unstable branches is clearly much less important. For the combinations (0.6, −0.5, −2) and (0.6, 0.7, 2) corresponding to the regions 2 and 3 for which only the stable solution exists, the corresponding values of are 0.441 and 0.428, respectively. For these regions (2 and 3) can be computed analytically by solving the following equation:
(a)
(b)
(c)
For the combinations , (0.9, 0.5, 2), and , corresponding, respectively, to regions 4, 5, and 6, both stable and unstable branches exist for specific ranges of and the corresponding critical Darcy numbers, , marking the transition from the parallel flow towards the rest state are 0.4072, 0.1673, and 0.2938, respectively. Note that the bifurcation around is supercritical for regions 2 and 3 and subcritical for regions 4, 5, and 6 (results not reported here). For regions 4, 5, and 6, can be computed analytically by solving the following equation:with
For the combination , a singularity is registered in the curve of at the critical Darcy number , as can be seen in Figure 9(b). For the latter combination, the term that represents the inverse of the Nusselt number is zero for , leading to an infinite value for . This means that, for , the Dufour heat flux induced by the concentration gradient compensates the ordinary heat flux induced by the temperature gradient. Such a compensation is not observed for the mass transfer as and . For the particular combination , induced by the unstable branch (in its existence range) is higher than that induced by the stable branch (see behavior in the vicinity of ). Note that and induced by the stable branches are larger than those corresponding to unstable branches for all the combinations considered in Figures 9(a)–9(b).
5. Conclusion
Combined effects of Soret and Dufour parameters on thermosolutal convection in shallow horizontal Brinkman porous enclosures with stress-free upper surfaces is evaluated analytically and numerically. The analytical results obtained using the parallel flow assumption are validated numerically by solving the full governing equations. The thresholds marking the onset of stationary and finite amplitude convection are determined analytically according to the control parameters. The analytical study shows that the plane can be divided into up to six regions corresponding to different parallel flow regimes. The number, the location, and the extent of these regions are strongly controlled by Dufour and Soret numbers. Depending on the value of the Soret number, the Dufour effect may affect considerably the fluid flow and heat and mass transfer characteristics. By analyzing the effect of on , , and for various , we observe the existence of critical nodes (small ranges of ) for which these quantities are quasi-independent of . These nodes did not correspond to the same range of for , , and .
Nomenclature
| A: | Aspect ratio of the enclosure, |
| D: | Diffusion coefficient |
| Da: | Darcy number, |
| Du: | Dufour number, |
| : | Gravitational acceleration |
| : | Height of the enclosure |
| : | Constant mass flux per unit area |
| K: | Permeability of the porous medium |
| : | Soret coefficient |
| : | Dufour coefficient |
| : | Length of the enclosure |
| Le: | Lewis number, |
| N: | Buoyancy ratio, |
| Nu: | Thermal Nusselt number |
| : | Constant heat flux per unit area |
| : | Thermal Darcy–Rayleigh number, |
| S: | Dimensionless solute concentration, |
| : | Concentration |
| : | Sherwood number |
| : | Soret number, |
| t: | Dimensionless time, |
| T: | Dimensionless temperature, |
| : | Temperature |
| u: | Dimensionless horizontal velocity, |
| : | Dimensionless vertical velocity, |
| x: | Dimensionless distance along the x axis, |
| y: | Dimensionless distance along the y axis, |
| : | Concentration solute difference, |
| : | Temperature difference, . |
| : | Solute expansion coefficient |
| : | Thermal expansion coefficient () |
| α: | Thermal diffusivity |
| : | Kinematic viscosity of the fluid |
| ψ: | Dimensionless stream function, |
| ε: | Normalized porosity of the porous medium, |
| : | Porosity of the porous medium |
| σ: | Heat capacity ratio, . |
| c: | Critical value |
| 0: | Refers to the value taken at the center of the cavity. |
| ′: | Dimensional variables |
| Sub: | Subcritical |
| Sup: | Supercritical. |
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.