# Analytical and Numerical Study of Soret and Dufour Effects on Thermosolutal Convection in a Horizontal Brinkman Porous Layer with a Stress-Free Upper Boundary

Authors: Ismail Filahi; Mohamed Bourich; Mohammed Hasnaoui; Abdelkhalek Amahmid
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 Figure1, is a two-dimensional horizontal porous cavity saturated with a binary mixture. The height of the cavity is H′, its length is L′ (A=L′/H′≫1), 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, q′, and mass, j′. 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 T′ and the concentration S′ as follows:(1)ρ=ρ01−βTT′−T0′−βSS′−S0′,where ρ0 is the binary fluid density at temperature and concentration T0′,S0′ and βT and βS are the thermal and solutal expansion coefficients, respectively.Figure 1
Schematic diagram of the studied configuration.The dimensional equations governing this problem are written as(2)∇⋅V′⟶=0,(3)V′⟶=−Kμ∇⟶P′+∇2KμeffμV′⟶+ρg⟶,(4)ρCP∂T′∂t+ρCfV′⟶∇T′=∇⋅λ∇T′+ρCfKTS∇S′,(5)ε′∂S′∂t+V′⟶∇S′=∇⋅D∇S′+KST∇T′.The reference scales for time, length, velocity, and pressure areσH′2/α, H′, α/H′, and αμ/K, respectively. The parameters σ and α are defined by σ=ρCP/ρCf and α=λ/ρCf.The dimensionless temperature and mass fraction are given byT=T′−T0′/ΔT′ and S=S′−S0′/ΔS′ where ΔT′=q′H′/λ and ΔS′=j′H′/D.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:(6)η∂ξ∂t+ξ=Da∇2ξ+RT∂T∂x+N∂S∂x,(7)∂T∂t+u∂T∂x+v∂T∂y=∇2T+Du∇2S,(8)ε∂S∂t+u∂S∂x+v∂S∂y=1Le∇2S+Sr∇2T,(9)∇2ψ=−ξ,(10)u=∂ψ∂y,v=−∂ψ∂x.The boundary conditions associated to the present problem are(11a)ψ=∂ψ∂x=∂T∂x=∂S∂x=0,forx=±A/2,(11b)ψ=∂2ψ∂y2=0,∂T∂y=Du−11−Du Sr,∂S∂y=Sr−11−Du Sr,fory=12,(11c)ψ=∂ψ∂y=0,∂T∂y=Du−11−Du Sr,∂S∂y=Sr−11−Du Sr,fory=−12.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, Da, the Dufour number, Du, the Lewis number, Le, the thermal Rayleigh number, RT, and the Soret number, Sr. These parameters are defined as follows:(12)A=L′H′,Da=μeffKμH′2,Du=KTSΔS′αΔT′,Le=αD,RT=gβTKΔT′H′αv,Sr=KSTΔT′DΔS′.The Nusselt and Sherwood numbers are, respectively, given by the following expressions:(13)Nu=1ΔT+DuΔS,Sh=1ΔS+SrΔT,where ΔT=T0,1/2−T0,−1/2 and ΔS=S0,1/2−S0,−1/2 are the temperature and concentration differences, evaluated at x=0.


## 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 thresholdA=12. Thereby, the numerical results reported here were obtained with A=12 and a grid of 201×81. 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 201×81 differ by less than 0.2% from those corresponding to the finest grid 261×121. In addition, the grid 201×81 reproduces the analytical results with relative errors lower than 0.3%.Table 1
Effect of the grid size forRT=100, N=1, Sr=0.5, Du=0.5, and Le=1.1.
Mesh
ψ
0
Nu
Sh
161
×
61
−
4.273
2.710
2.755
201
×
81
−
4.269
2.703
2.748
261
×
121
−
4.269
2.701
2.745
Analytical results
−
4.268
2.697
2.741


## 4. Results and Discussion

### 4.1. Parallel Flow Analysis

Based on the remarks (verified numerically) allowing simplifications in the case of shallow cavities (A≫1), 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:(14)ψx,y≈ψy,Tx,y≈CTx+θTy,Sx,y≈CSx+θSy,the parameters CT and CS 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 A≫1. The steady form of the governing equations (6)–(8) while using these approximations can be reduced to the following set of ordinary differential equations:(15)Dad4ψdy4−d2ψdy2+R=0,(16)d2θTdy2=dψdyCT−Le DuCS1−Du Sr,(17)d2θSdy2=dψdyLeCS−SrCT1−Du Sr,where R=−RTCT+NCS.The corresponding boundary conditions on the horizontal walls become(18)ψ=∂2ψ∂y2=0,∂θT∂y=Du−11−Du Sr,∂θS∂y=Sr−11−Du Sr,fory=12,(19)ψ=∂ψ∂y=0,∂θT∂y=Du−11−Du Sr,∂θS∂y=Sr−11−Du Sr,fory=−12.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:(20)ψy=ψ0ZyZ0,(21)uy=ψ0HyZ0,(22)Tx,y=CTx+Du−11−Du Sry−CT−Le DuCS1−Du Srψ0SyZ0,(23)Sx,y=CSx+Sr−11−Du Sry−LeCS−SrCT1−Du Srψ0SyZ0.The parameterψ0 is the stream function at the center of the cavity, obtained analytically as(24)ψ0=±Z02Le2ΓDa−b±b2−4Le2c1/2,with(25)b=P−RT0F1,c=Q−RT0F2.The parameters in the previous equation are obtained as(26)P=1+Le2+2Sr Du Le,Q=1−Du Sr2,(27)F1=Le2+LeN,F2=NDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.The parameterRT0=RT/Rsup=RT/ADa is the normalized Rayleigh number, and Zy, Hy, and Sy are the functions, respectively, given by the following analytical expressions:(28)Zy=−χ1ζ2coshζy−χ2ζ2sinhζy+y2−χ3y−χ4,Hy=−χ1ζsinhζy−χ2ζcoshζy+2y−χ3,Sy=χ1ζ3sinhζy+χ2ζ3coshζy−y33+χ32y2+χ4y.The expressions ofχii=1,2,3,4 are given by(29)χ1=ζ2sinhζ/2+2ζcoshζ/2−4sinhζ/22ζcosh2ζ/2−ζ−2sinhζ/2coshζ/2,χ2=ζ2sinhζ/2+2ζcoshζ/2−4sinhζ/22ζcosh2ζ/2−ζ−2sinhζ/2coshζ/2,χ3=2Ωsinhζ/2−ζ2coshζ/22Ωcosh2ζ/2−ζ−2sinhζ/2coshζ/2,χ4=−χ32+ζ2−84ζ2.The parameterζ=Da−1/2 and the functions ADa and ΓDa are dependent on the Darcy number as follows:(30)ADa=χ1sinhζ/2ζ3−coshζ/22ζ2+112,ΓDa=χ12−3coshζ/2sinhζ/24ζ5+cosh2ζ/24ζ4+18ζ4+χ1−2sinhζ/2ζ5+coshζ/2ζ4−coshζ/212ζ4+χ222sinh2ζ/2ζ6+3coshζ/2sinhζ/24ζ4+sinh2ζ/212ζ4−18ζ4+1120.According to equation (21), the velocity cannot be zero in the range 0<y<1/2, which means that only unicellular flows are possible. In addition, the constants CT and CS are obtained by considering mass and thermal balances across any transversal section of the porous layer [34].(31)∫−1/21/2uT−∂T∂x−Du∂S∂xdy=0,(32)∫−1/21/2uS−1Le∂S∂x+Sr∂T∂xdy=0.Substituting equations (22) and (23) into equation (13), the respective analytical expressions of Nusselt and Sherwood numbers are obtained as(33)Nu=11+ADaRCT,Sh=11+ADaRLeCS.By substituting the solutions given by equations (21)–(23) into equations (31) and (32), the final analytical expressions of CT and CS are obtained as follows:(34)CT=−ADaRLe1−SrLe DuR2ΓDa−MDu−ADaR1−DuM+Le2R2ΓDaM+R2ΓDaM+Le2R2ΓDa−Le DuR2ΓDa−MDuLe SrR2ΓDa−SrM,CS=−ADaRLe1−SrM+R2ΓDa−ADaR1−DuLe SrR2ΓDa−SrMM+R2ΓDaM+Le2R2ΓDa−Le DuR2ΓDa−MDuLe SrR2ΓDa−SrM,with M=1−Du Sr and R=2ψ0/Z0.The analysis of equation (24) indicates that up to five different steady-state solutions are possible including the diffusive regime ψ0=0. 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 ψ0=0, while the subcritical bifurcation is characterized by the onset of motion through finite amplitude ψ0=±Z0−b/2Le2ΓDa.These parallel flow solutions exist only when the two following conditions are satisfied:−b±b2−4Le2c>0 and b2−4Le2c>0. The satisfaction of these conditions allows the subdivision of the N,Du 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 theN,Du plane. These regions are illustrated in Figures 2–4 for Sr>0, Sr<0, and Sr=0, respectively.Figure 2
Delimitation of the regions corresponding to different behaviors forSr=0.6.Figure 3
Delimitation of the regions corresponding to different behaviors forSr=−0.6.Figure 4
Delimitation of the regions corresponding to different behaviors forSr=0.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 ofRT. Different conditions must be verified to delineate this region in the N,Du plane, depending on the sign of Sr.(i)
ForSr>0(35)N<−Le,for Du≥1+LeSr−1,N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,forSr1+Le−LeSr≤Du≤1+LeSr−1.(ii)
ForSr<0(36)N<−Le,for Du≥1+LeSr−1,N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≤1+LeSr−1.(iii)
ForSr=0(37)N<−Le,for Du≥1−Le,N<Du1+LeLe−1Le,for Du≤1−Le.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(38)RTCsup=Rsup1−Du Sr2NDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.The conditions that should be verified to delineate this region in theN,Du plane are dependent on the sign of Sr.(i)
ForSr>0(39)N<fDu,for Du<−1−Le22Le Sr,N<−Le,for−1−Le22Le Sr≤Du≤Sr1+Le−LeSr,Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤−Le,forSr1+Le−LeSr≤Du≤1+LeSr−1,with(40)fDu=Du22Sr Le1+Le−Sr2Le2+Du1+Le+Le2+Le3−3Sr Le−Le3Sr−1LeSr2Du2+DuSr+3SrLe2−2Sr2Le2−2Sr2Le2+2Le Sr+Le3−SrLe3−Sr−Sr Le−SrLe2.(ii)
ForSr<0(41)Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤−Le,for Du≤1+LeSr−1.(iii)
ForSr=0(42)N<−Le,for Du≥1−Le,N>Du1+LeLe−1Le,for Du≤1−Le.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 ofSr as follows:(i)
ForSr>0(43)−Le≤N<fDu,for−1−Le22Le Sr<Du<Du1,N>−Le,forDu1<Du<1+LeSr−1,N>fDu,for Du>1+LeSr−1,with(44)Du1=g+1+10Le2−16SrLe3+4Sr3Le4−8SrLe2+5Le4−8SrLe4+12Le3+4Sr2Le2+8Sr2Le3+4Le2Sr Le,Du2=g−1+10Le2−16SrLe3+4Sr3Le4−8SrLe2+5Le4−8SrLe4+12Le3+4Sr2Le2+8Sr2Le3+4Le2Sr Le,andg=1+3Le2+2Le−2SrLe2−2Sr Le(ii)
ForSr<0(45)N>−Le,for Du≤1+LeSr−1,N>fDu,for1+LeSr−1≤Du<Du2.(iii)
ForSr=0(46)N>−Le,for Du<1−Le,N>maxDu1+Le+Le2+Le3Le3−1Le3,Du1+LeLe−1Le,for Du>1−Le.For the region (4), both stable and unstable solutions are existing from the subcritical Rayleigh number, but the unstable branch disappears at some threshold ofRT. 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,RTCsub, above which the parallel flow exists, is obtained from the condition b2−4Le2c=0 and its expression is given by(47)RTCsub=G+2Le2Le2F22+1−Du Sr2Le2+LeN−F2Le2+LeN1+Le2+2Sr Du LeLe2+LeN2Rsup,where G=Le4−Le2+NLe−Le3+2NSr DuLe3−Le2+2NSr+2DuLe2+Le3.The conditions that should be verified for the region (4) in theN,Du plane are dependent on the sign of Sr as follows:(i)
ForSr>0(48)−Le<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du<−1−Le22Le Sr,fDu<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for−1−Le22Le Sr≤Du≤Sr1+Le−LeSr,N>fDu,forSr1+Le−LeSr<Du<Du1,Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤fDu,for Du>1+LeSr−1.(ii)
ForSr<0(49)Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le<N<fDu,for1+LeSr−1≤Du<Du2,N>Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,forDu2≤Du≤Sr1+Le−LeSr.(iii)
ForSr=0(50)N>−Le,for Du≤1−Le,Du1+LeLe−1Le<N<Du1+Le+Le2+Le3Le3−1Le3,for Du>1−Le.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 infiniteRT.(i)
ForSr>0(51)N>Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≤Sr1+Le−LeSr,−Le<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≥1+LeSr−1.(ii)
ForSr<0(52)−Le≤N≤Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≥1+LeSr−1.(iii)
ForSr=0(53)−Le≤N≤Du1+LeLe−1Le,for Du≥1−Le.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:(54)fDu≤N≤−Le,for Du≤Sr1+Le−LeSr.The expression offDu was already defined for region (2).


### 4.2. Effect ofRT

The effect ofRT on the flow intensity, ψ0=ψ0,0, the Nusselt number, Nu, and the Sherwood number, Sh, is illustrated in Figures 5(a)–5(c) for Sr=0.6, Da=0.01, Le=1.1, and different combinations of N,Du corresponding to different regions in this plane. The choice of Sr=0.6 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.Figure 5
(a) Variations ofψ0 vs. RT for Sr=0.6, Le=1.1, and different regions delineated in the N,Du plane. (b) Variations of Nu vs. RT for Sr=0.6, Le=1.1, and different regions delineated in the N,Du plane. (c) Variations of Sh vs. RT for Sr=0.6, Le=1.1 and different regions delineated in the N,Du plane.
(a)
(b)
(c)For the combinationN,Du=−2,−0.5 corresponding to the region (2), only the stable branch exists. The supercritical convection starts from the rest state (characterized by ψ0=ψ0,0=0, Nu=1, and Sh=1) at RTCSup=11.12. The evolutions of ψ0, Nu, and Sh are characterized by asymptotic behaviors, and their respective asymptotic limits are 2.31, 3.45, and 1.64, reached for sufficiently large values of RT. Note that for this combination of N,Du, the asymptotic value of Sh indicates that diffusion plays an important role in the mass transfer. At large values of RT, the analytical expressions of the asymptotic expressions of ψ0, Nu, and Sh corresponding to the stable solution are obtained by the following expressions:(55)ψ0asy=±Z02NDu Sr−Sr1+Le+Le−Du1+Le−Sr Le+1ΓDaLeLe+N,(56)Nuasy=hh−ADa2χ2Le1−SrLe Duχ2ΓDa−MDu−ADa2χ21−DuM+Le2χ2ΓDa,(57)Shasy=hh−χ2ADa2Le21−SrM+χ2ΓDa−Le1−DuLe Srχ2ΓDa−SrM,with h=M+χ2ΓDaM+Le2χ2ΓDa−Le Duχ2ΓDa−MDuLe Srχ2ΓDa−SrM and χ=NDu Sr−Sr1+Le+Le−Du1+Le−Sr Le+1/ΓDaLeLe+N.For the combinationN,Du=2,0.7, illustrating the region (3), the behavior is also supercritical and the convection starts from the rest state at RTCSup=11.42. The evolution of ψ0 vs. RT is characterized by a monotonous increase with different rates that depend on the combination N,Du. For RT⟶+∞, the analytical expression of ψ0 is reduced to(58)ψ0=±Z02N+LeADaΓDaLeRT1/2.Equation (58) shows that ψ0 varies as RT1/2, while the evolutions of Nu and Sh 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 Nu and Sh at large RT:(59)Nu=Sh=ΓDaΓDa−ADa2.The region (4) is illustrated by the combinationN,Du=0.3,0.7 in Figures 5(a)–5(c). Unlike the previous combinations of N,Du, for which only the stable branches exist, both stable and unstable solutions exist for this combination in the range RTCSub=45.12≤RT<RTCSup=83.50. This means that the parallel flow solution starts at RTCSub=45.12, and the nascent flow corresponds to a state clearly different from the purely diffusive regime characterized by ψ0=0 and Nu=Sh=1. For the unstable branch, ψ0/(Nu and Sh) decreases notably/(weakly) by increasing RT and disappears for RT≥RTCSup=83.50.For the stable branch, the evolution ofψ0 versus RT is very similar to that described for the combination N,Du=2,0.7 but characterized by a delay in terms of RT. At large values of the last parameter, the expression of ψ0 is the same as that established for the region (3) and given by equation (58). For the region (4) also, the evolutions of Nu and Sh are characterized by asymptotic behaviors at large RT, leading to the common asymptotic value 4.49, obtained with equation (59) for both Nu and Sh.The results corresponding to the region (5) are illustrated in Figures5(a)–5(c) by the combination N,Du=0.2,0.9. For this region, both stable and unstable solutions exist for RT≥RTCSub=69.16 for the considered combination of N,Du. In addition, the evolution of the stable branch is similar to that described for region (4). By increasing RT, the unstable branches are characterized by slow decreases towards asymptotic limits for ψ0, Nu, and Sh.The region (6) is illustrated by the combinationN,Du=−1.2,−1.8 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 RTCSub=14.793≤RT<RTCSup=14.797. In addition, at large values of RT, the quantities ψ0, Nu, and Sh 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 RT. Consequently, the analytical expressions of the asymptotic values of ψ0, Nu, and Sh are given by equations (55)–(57), respectively. Besides, for the region (6), the Nusselt number becomes negative in the range 27.34≤RT≤453.99. 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 RT.


### 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 Figures6(a)–6(c) for RT = 500, N=−0.2, Le=1.1, Da=0.01 and various values of Sr. These figures show that, regardless of the value of Sr, there exists a critical value of Du above which convection disappears. This critical value increases with Sr; it rises from 0.124 to 0.897 when Sr goes from −0.9 to 0.9.Figure 6
(a) Variations ofψ0 vs. Du for RT=500, N=−0.2, Le=1.1, and various values of Sr. (b) Variations of Nu vs. Du for RT=500, N=−0.2, Le=1.1, and various values of Sr. (c) Variations of Sh vs. Du for RT=500, N=−0.2, Le=1.1, and various values of Sr.
(a)
(b)
(c)The effect ofSr on ψ0, Nu, and Sh may be drastically affected by the change of Du. More specifically, we can observe two different behaviors depending on whether the value of Du 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 Du, the effect of thermo-diffusion becomes insignificant. Another mystery of these critical ranges of Du lies in the fact that they are not the same for ψ0, Nu, and Sh (compare [0.6, 0.732] for ψ0, [0.016, 0.042] for Nu, and [−2.29, −2.25] for Sh). 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 Sh is seen to be largely shifted towards negative values of Du. 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 ψ0, Nu, and Sh are increasing functions of Sr. Above the critical node corresponding to Sh, the latter decreases by increasing Sr. However, immediately near the critical nodes, ψ0 and Nu exhibit a behavior similar to that of Sh 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 Nu decreases by increasing Du, while Sh increases/decreases by increasing Du for Sr<−0.1/Sr≥−0.1. Furthermore, the most important variations observed for Nu/Sh when Sr is varied are located below/above the critical nodes. Quantitatively, at Du=−3 for instance, by increasing Sr from −0.9 to 0.9, Nu is multiplied by a factor of 4.1 while Sh is increased by only 11.44%. For Du=2, by increasing Sr from −0.9 to 0.9, Sh is divided by 3.2 while Nu undergoes a relative decrease of about 14.82%.The combined effects ofDu and Sr on the temperature and concentration fields are illustrated in Figure 7 for RT = 500, N = −0.2, Le = 1.1, Da=0.01, Sr=±0.9, and Du=±0.9. For the combination Du,Sr=0.9,0.9, 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 Du,Sr from 0.9,0.9 to 0.9,−0.9, 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 Du,Sr=0.9,−0.9 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 Du,Sr=0.9,−0.9 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 Du,Sr passes from 0.9,0.9 to −0.9,0.9 are similar to those observed in the concentration/temperature field when Du,Sr passes from 0.9,0.9 to 0.9,−0.9. Hence, we can observe a quasi-linear behavior for the concentration field at Du,Sr=−0.9,0.9. The combination Du,Sr=−0.9,−0.9 leads to the most important variations for the temperature and the concentration among the four combinations considered in this comparison. Quantitatively, by changing Du,Sr from 0.9,0.9 to 0.9,−0.9/−0.9,0.9/−0.9,−0.9, the temperature difference between the horizontal wall at a fixed cross section passes from 0.1769 to 0.0796/0.2416/2.08. Thus, the temperature difference obtained with the combination Du,Sr=−0.9,−0.9 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.Figure 7
Isotherms and isosolutes for Le = 1.1,N=−0.2, and RT=500.Figure 8
Temperature (a) and concentration (b) profiles along y at midwidth of the cavity forLe=1.1, Da=0.01, N=−0.2, and RT=500.
(a)
(b)


### 4.4. Effect of Darcy Number

The Darcy number,Da, 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 ofDa on ψ0, Nu, and Sh is illustrated in Figures 9(a)–9(c) for RT=100, Le=1.1, and different combinations Sr,Du,N 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 Da exceeds a critical value, Dacr, that depends on the combination Sr,Du,N. For the stable branches, the quantities ψ0, Nu, and Sh decrease by increasing the permeability of the porous medium. Similar behaviors were reported in previous studies [5, 36]. Note that the variations versus Da of the flow intensity corresponding to the unstable branches are remarkably strong in specific ranges of the latter parameter. However, the effect of Da on Nu and Sh 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 Dacr are 0.441 and 0.428, respectively. For these regions (2 and 3) Dacr can be computed analytically by solving the following equation:(60)ADacr=1−Du Sr2RTNDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.Figure 9
(a) Variations ofψ0 vs. Da for Le=1.1, RT=100 and different regions. (b) Variations of Nu vs. Da for Le=1.1, RT=100 and different regions. (c) Variations of Sh vs. Da for Le=1.1, RT=100 and different regions.
(a)
(b)
(c)For the combinations0.4,−0.3,2, (0.9, 0.5, 2), and 0.9,−2,−2, corresponding, respectively, to regions 4, 5, and 6, both stable and unstable branches exist for specific ranges of Da and the corresponding critical Darcy numbers, Dacr, 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 Dacr 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, Dacr can be computed analytically by solving the following equation:(61)ADacr=f+2Le21−Du Sr2Le2+LeN−F2Le2+LeN1+Le2+2Sr Du Le+Le2F22RTLe2+LeN2,with(62)f=Le4+Le2+NLe−Le3+2NSr DuLe3−Le2+2NSr+2DuLe2+Le3.For the combination0.9,−2,−2, a singularity is registered in the curve of Nu at the critical Darcy number DaNu=∞=0.145, as can be seen in Figure 9(b). For the latter combination, the term ΔT+DuΔS that represents the inverse of the Nusselt number is zero for Da=DaNu=∞, leading to an infinite value for Nu. This means that, for Da=DaNu=∞, 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 Sr≠Du and Le≠1. For the particular combination 0.9,−2,−2, Sh induced by the unstable branch (in its existence range) is higher than that induced by the stable branch (see Sh behavior in the vicinity of Dacr). Note that ψ0 and Nu induced by the stable branches are larger than those corresponding to unstable branches for all the combinations considered in Figures 9(a)–9(b).


## 4.1. Parallel Flow Analysis

Based on the remarks (verified numerically) allowing simplifications in the case of shallow cavities (A≫1), 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:(14)ψx,y≈ψy,Tx,y≈CTx+θTy,Sx,y≈CSx+θSy,the parameters CT and CS 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 A≫1. The steady form of the governing equations (6)–(8) while using these approximations can be reduced to the following set of ordinary differential equations:(15)Dad4ψdy4−d2ψdy2+R=0,(16)d2θTdy2=dψdyCT−Le DuCS1−Du Sr,(17)d2θSdy2=dψdyLeCS−SrCT1−Du Sr,where R=−RTCT+NCS.The corresponding boundary conditions on the horizontal walls become(18)ψ=∂2ψ∂y2=0,∂θT∂y=Du−11−Du Sr,∂θS∂y=Sr−11−Du Sr,fory=12,(19)ψ=∂ψ∂y=0,∂θT∂y=Du−11−Du Sr,∂θS∂y=Sr−11−Du Sr,fory=−12.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:(20)ψy=ψ0ZyZ0,(21)uy=ψ0HyZ0,(22)Tx,y=CTx+Du−11−Du Sry−CT−Le DuCS1−Du Srψ0SyZ0,(23)Sx,y=CSx+Sr−11−Du Sry−LeCS−SrCT1−Du Srψ0SyZ0.The parameterψ0 is the stream function at the center of the cavity, obtained analytically as(24)ψ0=±Z02Le2ΓDa−b±b2−4Le2c1/2,with(25)b=P−RT0F1,c=Q−RT0F2.The parameters in the previous equation are obtained as(26)P=1+Le2+2Sr Du Le,Q=1−Du Sr2,(27)F1=Le2+LeN,F2=NDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.The parameterRT0=RT/Rsup=RT/ADa is the normalized Rayleigh number, and Zy, Hy, and Sy are the functions, respectively, given by the following analytical expressions:(28)Zy=−χ1ζ2coshζy−χ2ζ2sinhζy+y2−χ3y−χ4,Hy=−χ1ζsinhζy−χ2ζcoshζy+2y−χ3,Sy=χ1ζ3sinhζy+χ2ζ3coshζy−y33+χ32y2+χ4y.The expressions ofχii=1,2,3,4 are given by(29)χ1=ζ2sinhζ/2+2ζcoshζ/2−4sinhζ/22ζcosh2ζ/2−ζ−2sinhζ/2coshζ/2,χ2=ζ2sinhζ/2+2ζcoshζ/2−4sinhζ/22ζcosh2ζ/2−ζ−2sinhζ/2coshζ/2,χ3=2Ωsinhζ/2−ζ2coshζ/22Ωcosh2ζ/2−ζ−2sinhζ/2coshζ/2,χ4=−χ32+ζ2−84ζ2.The parameterζ=Da−1/2 and the functions ADa and ΓDa are dependent on the Darcy number as follows:(30)ADa=χ1sinhζ/2ζ3−coshζ/22ζ2+112,ΓDa=χ12−3coshζ/2sinhζ/24ζ5+cosh2ζ/24ζ4+18ζ4+χ1−2sinhζ/2ζ5+coshζ/2ζ4−coshζ/212ζ4+χ222sinh2ζ/2ζ6+3coshζ/2sinhζ/24ζ4+sinh2ζ/212ζ4−18ζ4+1120.According to equation (21), the velocity cannot be zero in the range 0<y<1/2, which means that only unicellular flows are possible. In addition, the constants CT and CS are obtained by considering mass and thermal balances across any transversal section of the porous layer [34].(31)∫−1/21/2uT−∂T∂x−Du∂S∂xdy=0,(32)∫−1/21/2uS−1Le∂S∂x+Sr∂T∂xdy=0.Substituting equations (22) and (23) into equation (13), the respective analytical expressions of Nusselt and Sherwood numbers are obtained as(33)Nu=11+ADaRCT,Sh=11+ADaRLeCS.By substituting the solutions given by equations (21)–(23) into equations (31) and (32), the final analytical expressions of CT and CS are obtained as follows:(34)CT=−ADaRLe1−SrLe DuR2ΓDa−MDu−ADaR1−DuM+Le2R2ΓDaM+R2ΓDaM+Le2R2ΓDa−Le DuR2ΓDa−MDuLe SrR2ΓDa−SrM,CS=−ADaRLe1−SrM+R2ΓDa−ADaR1−DuLe SrR2ΓDa−SrMM+R2ΓDaM+Le2R2ΓDa−Le DuR2ΓDa−MDuLe SrR2ΓDa−SrM,with M=1−Du Sr and R=2ψ0/Z0.The analysis of equation (24) indicates that up to five different steady-state solutions are possible including the diffusive regime ψ0=0. 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 ψ0=0, while the subcritical bifurcation is characterized by the onset of motion through finite amplitude ψ0=±Z0−b/2Le2ΓDa.These parallel flow solutions exist only when the two following conditions are satisfied:−b±b2−4Le2c>0 and b2−4Le2c>0. The satisfaction of these conditions allows the subdivision of the N,Du 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 theN,Du plane. These regions are illustrated in Figures 2–4 for Sr>0, Sr<0, and Sr=0, respectively.Figure 2
Delimitation of the regions corresponding to different behaviors forSr=0.6.Figure 3
Delimitation of the regions corresponding to different behaviors forSr=−0.6.Figure 4
Delimitation of the regions corresponding to different behaviors forSr=0.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 ofRT. Different conditions must be verified to delineate this region in the N,Du plane, depending on the sign of Sr.(i)
ForSr>0(35)N<−Le,for Du≥1+LeSr−1,N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,forSr1+Le−LeSr≤Du≤1+LeSr−1.(ii)
ForSr<0(36)N<−Le,for Du≥1+LeSr−1,N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≤1+LeSr−1.(iii)
ForSr=0(37)N<−Le,for Du≥1−Le,N<Du1+LeLe−1Le,for Du≤1−Le.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(38)RTCsup=Rsup1−Du Sr2NDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.The conditions that should be verified to delineate this region in theN,Du plane are dependent on the sign of Sr.(i)
ForSr>0(39)N<fDu,for Du<−1−Le22Le Sr,N<−Le,for−1−Le22Le Sr≤Du≤Sr1+Le−LeSr,Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤−Le,forSr1+Le−LeSr≤Du≤1+LeSr−1,with(40)fDu=Du22Sr Le1+Le−Sr2Le2+Du1+Le+Le2+Le3−3Sr Le−Le3Sr−1LeSr2Du2+DuSr+3SrLe2−2Sr2Le2−2Sr2Le2+2Le Sr+Le3−SrLe3−Sr−Sr Le−SrLe2.(ii)
ForSr<0(41)Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤−Le,for Du≤1+LeSr−1.(iii)
ForSr=0(42)N<−Le,for Du≥1−Le,N>Du1+LeLe−1Le,for Du≤1−Le.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 ofSr as follows:(i)
ForSr>0(43)−Le≤N<fDu,for−1−Le22Le Sr<Du<Du1,N>−Le,forDu1<Du<1+LeSr−1,N>fDu,for Du>1+LeSr−1,with(44)Du1=g+1+10Le2−16SrLe3+4Sr3Le4−8SrLe2+5Le4−8SrLe4+12Le3+4Sr2Le2+8Sr2Le3+4Le2Sr Le,Du2=g−1+10Le2−16SrLe3+4Sr3Le4−8SrLe2+5Le4−8SrLe4+12Le3+4Sr2Le2+8Sr2Le3+4Le2Sr Le,andg=1+3Le2+2Le−2SrLe2−2Sr Le(ii)
ForSr<0(45)N>−Le,for Du≤1+LeSr−1,N>fDu,for1+LeSr−1≤Du<Du2.(iii)
ForSr=0(46)N>−Le,for Du<1−Le,N>maxDu1+Le+Le2+Le3Le3−1Le3,Du1+LeLe−1Le,for Du>1−Le.For the region (4), both stable and unstable solutions are existing from the subcritical Rayleigh number, but the unstable branch disappears at some threshold ofRT. 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,RTCsub, above which the parallel flow exists, is obtained from the condition b2−4Le2c=0 and its expression is given by(47)RTCsub=G+2Le2Le2F22+1−Du Sr2Le2+LeN−F2Le2+LeN1+Le2+2Sr Du LeLe2+LeN2Rsup,where G=Le4−Le2+NLe−Le3+2NSr DuLe3−Le2+2NSr+2DuLe2+Le3.The conditions that should be verified for the region (4) in theN,Du plane are dependent on the sign of Sr as follows:(i)
ForSr>0(48)−Le<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du<−1−Le22Le Sr,fDu<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for−1−Le22Le Sr≤Du≤Sr1+Le−LeSr,N>fDu,forSr1+Le−LeSr<Du<Du1,Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤fDu,for Du>1+LeSr−1.(ii)
ForSr<0(49)Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le<N<fDu,for1+LeSr−1≤Du<Du2,N>Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,forDu2≤Du≤Sr1+Le−LeSr.(iii)
ForSr=0(50)N>−Le,for Du≤1−Le,Du1+LeLe−1Le<N<Du1+Le+Le2+Le3Le3−1Le3,for Du>1−Le.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 infiniteRT.(i)
ForSr>0(51)N>Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≤Sr1+Le−LeSr,−Le<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≥1+LeSr−1.(ii)
ForSr<0(52)−Le≤N≤Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≥1+LeSr−1.(iii)
ForSr=0(53)−Le≤N≤Du1+LeLe−1Le,for Du≥1−Le.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:(54)fDu≤N≤−Le,for Du≤Sr1+Le−LeSr.The expression offDu was already defined for region (2).


## 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 theN,Du plane. These regions are illustrated in Figures 2–4 for Sr>0, Sr<0, and Sr=0, respectively.Figure 2
Delimitation of the regions corresponding to different behaviors forSr=0.6.Figure 3
Delimitation of the regions corresponding to different behaviors forSr=−0.6.Figure 4
Delimitation of the regions corresponding to different behaviors forSr=0.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 ofRT. Different conditions must be verified to delineate this region in the N,Du plane, depending on the sign of Sr.(i)
ForSr>0(35)N<−Le,for Du≥1+LeSr−1,N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,forSr1+Le−LeSr≤Du≤1+LeSr−1.(ii)
ForSr<0(36)N<−Le,for Du≥1+LeSr−1,N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≤1+LeSr−1.(iii)
ForSr=0(37)N<−Le,for Du≥1−Le,N<Du1+LeLe−1Le,for Du≤1−Le.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(38)RTCsup=Rsup1−Du Sr2NDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.The conditions that should be verified to delineate this region in theN,Du plane are dependent on the sign of Sr.(i)
ForSr>0(39)N<fDu,for Du<−1−Le22Le Sr,N<−Le,for−1−Le22Le Sr≤Du≤Sr1+Le−LeSr,Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤−Le,forSr1+Le−LeSr≤Du≤1+LeSr−1,with(40)fDu=Du22Sr Le1+Le−Sr2Le2+Du1+Le+Le2+Le3−3Sr Le−Le3Sr−1LeSr2Du2+DuSr+3SrLe2−2Sr2Le2−2Sr2Le2+2Le Sr+Le3−SrLe3−Sr−Sr Le−SrLe2.(ii)
ForSr<0(41)Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤−Le,for Du≤1+LeSr−1.(iii)
ForSr=0(42)N<−Le,for Du≥1−Le,N>Du1+LeLe−1Le,for Du≤1−Le.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 ofSr as follows:(i)
ForSr>0(43)−Le≤N<fDu,for−1−Le22Le Sr<Du<Du1,N>−Le,forDu1<Du<1+LeSr−1,N>fDu,for Du>1+LeSr−1,with(44)Du1=g+1+10Le2−16SrLe3+4Sr3Le4−8SrLe2+5Le4−8SrLe4+12Le3+4Sr2Le2+8Sr2Le3+4Le2Sr Le,Du2=g−1+10Le2−16SrLe3+4Sr3Le4−8SrLe2+5Le4−8SrLe4+12Le3+4Sr2Le2+8Sr2Le3+4Le2Sr Le,andg=1+3Le2+2Le−2SrLe2−2Sr Le(ii)
ForSr<0(45)N>−Le,for Du≤1+LeSr−1,N>fDu,for1+LeSr−1≤Du<Du2.(iii)
ForSr=0(46)N>−Le,for Du<1−Le,N>maxDu1+Le+Le2+Le3Le3−1Le3,Du1+LeLe−1Le,for Du>1−Le.For the region (4), both stable and unstable solutions are existing from the subcritical Rayleigh number, but the unstable branch disappears at some threshold ofRT. 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,RTCsub, above which the parallel flow exists, is obtained from the condition b2−4Le2c=0 and its expression is given by(47)RTCsub=G+2Le2Le2F22+1−Du Sr2Le2+LeN−F2Le2+LeN1+Le2+2Sr Du LeLe2+LeN2Rsup,where G=Le4−Le2+NLe−Le3+2NSr DuLe3−Le2+2NSr+2DuLe2+Le3.The conditions that should be verified for the region (4) in theN,Du plane are dependent on the sign of Sr as follows:(i)
ForSr>0(48)−Le<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du<−1−Le22Le Sr,fDu<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for−1−Le22Le Sr≤Du≤Sr1+Le−LeSr,N>fDu,forSr1+Le−LeSr<Du<Du1,Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le≤N≤fDu,for Du>1+LeSr−1.(ii)
ForSr<0(49)Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le<N<fDu,for1+LeSr−1≤Du<Du2,N>Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,forDu2≤Du≤Sr1+Le−LeSr.(iii)
ForSr=0(50)N>−Le,for Du≤1−Le,Du1+LeLe−1Le<N<Du1+Le+Le2+Le3Le3−1Le3,for Du>1−Le.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 infiniteRT.(i)
ForSr>0(51)N>Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≤Sr1+Le−LeSr,−Le<N<Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≥1+LeSr−1.(ii)
ForSr<0(52)−Le≤N≤Du1+Le−Le Sr−1Du Sr−Sr1+Le+Le,for Du≥1+LeSr−1.(iii)
ForSr=0(53)−Le≤N≤Du1+LeLe−1Le,for Du≥1−Le.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:(54)fDu≤N≤−Le,for Du≤Sr1+Le−LeSr.The expression offDu was already defined for region (2).


## 4.2. Effect ofRT

The effect ofRT on the flow intensity, ψ0=ψ0,0, the Nusselt number, Nu, and the Sherwood number, Sh, is illustrated in Figures 5(a)–5(c) for Sr=0.6, Da=0.01, Le=1.1, and different combinations of N,Du corresponding to different regions in this plane. The choice of Sr=0.6 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.Figure 5
(a) Variations ofψ0 vs. RT for Sr=0.6, Le=1.1, and different regions delineated in the N,Du plane. (b) Variations of Nu vs. RT for Sr=0.6, Le=1.1, and different regions delineated in the N,Du plane. (c) Variations of Sh vs. RT for Sr=0.6, Le=1.1 and different regions delineated in the N,Du plane.
(a)
(b)
(c)For the combinationN,Du=−2,−0.5 corresponding to the region (2), only the stable branch exists. The supercritical convection starts from the rest state (characterized by ψ0=ψ0,0=0, Nu=1, and Sh=1) at RTCSup=11.12. The evolutions of ψ0, Nu, and Sh are characterized by asymptotic behaviors, and their respective asymptotic limits are 2.31, 3.45, and 1.64, reached for sufficiently large values of RT. Note that for this combination of N,Du, the asymptotic value of Sh indicates that diffusion plays an important role in the mass transfer. At large values of RT, the analytical expressions of the asymptotic expressions of ψ0, Nu, and Sh corresponding to the stable solution are obtained by the following expressions:(55)ψ0asy=±Z02NDu Sr−Sr1+Le+Le−Du1+Le−Sr Le+1ΓDaLeLe+N,(56)Nuasy=hh−ADa2χ2Le1−SrLe Duχ2ΓDa−MDu−ADa2χ21−DuM+Le2χ2ΓDa,(57)Shasy=hh−χ2ADa2Le21−SrM+χ2ΓDa−Le1−DuLe Srχ2ΓDa−SrM,with h=M+χ2ΓDaM+Le2χ2ΓDa−Le Duχ2ΓDa−MDuLe Srχ2ΓDa−SrM and χ=NDu Sr−Sr1+Le+Le−Du1+Le−Sr Le+1/ΓDaLeLe+N.For the combinationN,Du=2,0.7, illustrating the region (3), the behavior is also supercritical and the convection starts from the rest state at RTCSup=11.42. The evolution of ψ0 vs. RT is characterized by a monotonous increase with different rates that depend on the combination N,Du. For RT⟶+∞, the analytical expression of ψ0 is reduced to(58)ψ0=±Z02N+LeADaΓDaLeRT1/2.Equation (58) shows that ψ0 varies as RT1/2, while the evolutions of Nu and Sh 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 Nu and Sh at large RT:(59)Nu=Sh=ΓDaΓDa−ADa2.The region (4) is illustrated by the combinationN,Du=0.3,0.7 in Figures 5(a)–5(c). Unlike the previous combinations of N,Du, for which only the stable branches exist, both stable and unstable solutions exist for this combination in the range RTCSub=45.12≤RT<RTCSup=83.50. This means that the parallel flow solution starts at RTCSub=45.12, and the nascent flow corresponds to a state clearly different from the purely diffusive regime characterized by ψ0=0 and Nu=Sh=1. For the unstable branch, ψ0/(Nu and Sh) decreases notably/(weakly) by increasing RT and disappears for RT≥RTCSup=83.50.For the stable branch, the evolution ofψ0 versus RT is very similar to that described for the combination N,Du=2,0.7 but characterized by a delay in terms of RT. At large values of the last parameter, the expression of ψ0 is the same as that established for the region (3) and given by equation (58). For the region (4) also, the evolutions of Nu and Sh are characterized by asymptotic behaviors at large RT, leading to the common asymptotic value 4.49, obtained with equation (59) for both Nu and Sh.The results corresponding to the region (5) are illustrated in Figures5(a)–5(c) by the combination N,Du=0.2,0.9. For this region, both stable and unstable solutions exist for RT≥RTCSub=69.16 for the considered combination of N,Du. In addition, the evolution of the stable branch is similar to that described for region (4). By increasing RT, the unstable branches are characterized by slow decreases towards asymptotic limits for ψ0, Nu, and Sh.The region (6) is illustrated by the combinationN,Du=−1.2,−1.8 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 RTCSub=14.793≤RT<RTCSup=14.797. In addition, at large values of RT, the quantities ψ0, Nu, and Sh 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 RT. Consequently, the analytical expressions of the asymptotic values of ψ0, Nu, and Sh are given by equations (55)–(57), respectively. Besides, for the region (6), the Nusselt number becomes negative in the range 27.34≤RT≤453.99. 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 RT.


## 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 Figures6(a)–6(c) for RT = 500, N=−0.2, Le=1.1, Da=0.01 and various values of Sr. These figures show that, regardless of the value of Sr, there exists a critical value of Du above which convection disappears. This critical value increases with Sr; it rises from 0.124 to 0.897 when Sr goes from −0.9 to 0.9.Figure 6
(a) Variations ofψ0 vs. Du for RT=500, N=−0.2, Le=1.1, and various values of Sr. (b) Variations of Nu vs. Du for RT=500, N=−0.2, Le=1.1, and various values of Sr. (c) Variations of Sh vs. Du for RT=500, N=−0.2, Le=1.1, and various values of Sr.
(a)
(b)
(c)The effect ofSr on ψ0, Nu, and Sh may be drastically affected by the change of Du. More specifically, we can observe two different behaviors depending on whether the value of Du 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 Du, the effect of thermo-diffusion becomes insignificant. Another mystery of these critical ranges of Du lies in the fact that they are not the same for ψ0, Nu, and Sh (compare [0.6, 0.732] for ψ0, [0.016, 0.042] for Nu, and [−2.29, −2.25] for Sh). 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 Sh is seen to be largely shifted towards negative values of Du. 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 ψ0, Nu, and Sh are increasing functions of Sr. Above the critical node corresponding to Sh, the latter decreases by increasing Sr. However, immediately near the critical nodes, ψ0 and Nu exhibit a behavior similar to that of Sh 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 Nu decreases by increasing Du, while Sh increases/decreases by increasing Du for Sr<−0.1/Sr≥−0.1. Furthermore, the most important variations observed for Nu/Sh when Sr is varied are located below/above the critical nodes. Quantitatively, at Du=−3 for instance, by increasing Sr from −0.9 to 0.9, Nu is multiplied by a factor of 4.1 while Sh is increased by only 11.44%. For Du=2, by increasing Sr from −0.9 to 0.9, Sh is divided by 3.2 while Nu undergoes a relative decrease of about 14.82%.The combined effects ofDu and Sr on the temperature and concentration fields are illustrated in Figure 7 for RT = 500, N = −0.2, Le = 1.1, Da=0.01, Sr=±0.9, and Du=±0.9. For the combination Du,Sr=0.9,0.9, 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 Du,Sr from 0.9,0.9 to 0.9,−0.9, 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 Du,Sr=0.9,−0.9 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 Du,Sr=0.9,−0.9 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 Du,Sr passes from 0.9,0.9 to −0.9,0.9 are similar to those observed in the concentration/temperature field when Du,Sr passes from 0.9,0.9 to 0.9,−0.9. Hence, we can observe a quasi-linear behavior for the concentration field at Du,Sr=−0.9,0.9. The combination Du,Sr=−0.9,−0.9 leads to the most important variations for the temperature and the concentration among the four combinations considered in this comparison. Quantitatively, by changing Du,Sr from 0.9,0.9 to 0.9,−0.9/−0.9,0.9/−0.9,−0.9, the temperature difference between the horizontal wall at a fixed cross section passes from 0.1769 to 0.0796/0.2416/2.08. Thus, the temperature difference obtained with the combination Du,Sr=−0.9,−0.9 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.Figure 7
Isotherms and isosolutes for Le = 1.1,N=−0.2, and RT=500.Figure 8
Temperature (a) and concentration (b) profiles along y at midwidth of the cavity forLe=1.1, Da=0.01, N=−0.2, and RT=500.
(a)
(b)


## 4.4. Effect of Darcy Number

The Darcy number,Da, 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 ofDa on ψ0, Nu, and Sh is illustrated in Figures 9(a)–9(c) for RT=100, Le=1.1, and different combinations Sr,Du,N 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 Da exceeds a critical value, Dacr, that depends on the combination Sr,Du,N. For the stable branches, the quantities ψ0, Nu, and Sh decrease by increasing the permeability of the porous medium. Similar behaviors were reported in previous studies [5, 36]. Note that the variations versus Da of the flow intensity corresponding to the unstable branches are remarkably strong in specific ranges of the latter parameter. However, the effect of Da on Nu and Sh 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 Dacr are 0.441 and 0.428, respectively. For these regions (2 and 3) Dacr can be computed analytically by solving the following equation:(60)ADacr=1−Du Sr2RTNDu Sr−Sr1+Le+Le−Du1+Le−Le Sr+1.Figure 9
(a) Variations ofψ0 vs. Da for Le=1.1, RT=100 and different regions. (b) Variations of Nu vs. Da for Le=1.1, RT=100 and different regions. (c) Variations of Sh vs. Da for Le=1.1, RT=100 and different regions.
(a)
(b)
(c)For the combinations0.4,−0.3,2, (0.9, 0.5, 2), and 0.9,−2,−2, corresponding, respectively, to regions 4, 5, and 6, both stable and unstable branches exist for specific ranges of Da and the corresponding critical Darcy numbers, Dacr, 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 Dacr 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, Dacr can be computed analytically by solving the following equation:(61)ADacr=f+2Le21−Du Sr2Le2+LeN−F2Le2+LeN1+Le2+2Sr Du Le+Le2F22RTLe2+LeN2,with(62)f=Le4+Le2+NLe−Le3+2NSr DuLe3−Le2+2NSr+2DuLe2+Le3.For the combination0.9,−2,−2, a singularity is registered in the curve of Nu at the critical Darcy number DaNu=∞=0.145, as can be seen in Figure 9(b). For the latter combination, the term ΔT+DuΔS that represents the inverse of the Nusselt number is zero for Da=DaNu=∞, leading to an infinite value for Nu. This means that, for Da=DaNu=∞, 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 Sr≠Du and Le≠1. For the particular combination 0.9,−2,−2, Sh induced by the unstable branch (in its existence range) is higher than that induced by the stable branch (see Sh behavior in the vicinity of Dacr). Note that ψ0 and Nu 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 theN,Du 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 Du on ψ0, Nu, and Sh for various Sr, we observe the existence of critical nodes (small ranges of Du) for which these quantities are quasi-independent of Sr. These nodes did not correspond to the same range of Du for ψ0, Nu, and Sh.



---

*Source: 4046570-2020-03-31.xml*