The onset of double-diffusive convection in a Brinkman porous layer with convective thermal boundary conditions

Since the commencement of the study of Rayleigh-Bénard problem in the field of fluid mechanics, several analyses have been done considering this physical set up and the study has been extended under various other physical situations. Horton and Rogers¹ and Lapwood² led the analysis in the mid of the twentieth century by considering the flow through a porous layer. Further studies in this area considered the effect of the boundaries with different thermal conditions such as insulated, perfectly conducting or imperfectly conducting boundaries. The basic idea behind the consideration of such kinds of boundaries is just to maintain a purely vertical temperature gradient across the two boundaries, which in turn causes the onset of natural convection. Whenever, some other heating source apart from the two boundaries with different temperatures, is used, the instability in the system is caused due to the combined effect of natural convection and forced convection. All these phenomena in the porous media flow are discussed in great details by Nield and Bejan.³ 

Numerous studies can be found which used the combination of a perfectly insulating boundary and a perfectly conducting boundary across the horizontal porous layer. Some of them are Barletta et al.,^4,5 Barletta and Nield,⁶ Barletta and Celli,⁷ who used such boundaries and carried out the analysis with the combination of Neumann and Dirichlet thermal boundaries. While others used the combination of two Dirichlet thermal conditions at the two boundaries maintained at a basic temperature difference. Walker and Homsy,⁸ Rees,⁹ Nield et al.,¹⁰ Shivkumara et al.,¹¹ Narayana et al.,¹² Roy and Murthy¹³ are some of them who considered the perfectly conducting boundaries. Among them, Rees,⁹ Nield et al.,¹⁰ Shivkumara et al.¹¹ and Roy and Murthy¹³ considered the sparsely packed porous layer and used the Brinkman extended Darcy’s model for their investigation.

Over the past few decades, researchers have used the combination of two imperfectly conducting boundaries for investigating the Darcy-Bénard convection. They used the convective thermal boundary conditions which are given by the Robin-type boundary conditions. The onset of mixed convection in a porous layer with viscous dissipation was discussed by Storesletten and Barletta.¹⁴ Their analysis concluded that the transverse rolls are the most unstable convective rolls and the critical values of R(= GePe⁴) are dependent on the coefficient of external heat transfer at the lower boundary. Barletta and Storesletten^15–17 investigated the effect of the external heating at the two horizontal boundaries for a circular porous duct, a rectangular porous channel and a vertical cylindrical porous channel, respectively. Barletta et al.¹⁸ investigated the effect of free surface and convective boundary condition on the instability of the flow. They considered the case of local thermal non-equilibrium between the two phases of a porous layer. Braga et al.¹⁹ studied the effect of internal heating by considering the external heating at the upper boundary and a constant heat flux condition at the lower boundary. Again, Barletta et al.²⁰ investigated the onset of mixed convection for the vertical porous channel.

The present work is based on the model discussed by Dubey and Murthy²¹ considering the different forms of thermal and solutal boundary conditions and the persistence of the Soret effect, in addition. Most general convective thermal boundary conditions (Alves and Barletta²²) are considered, where the coefficients of external heating at the two boundaries are different. It also includes the effect of the Soret parameter on the onset of double-diffusive convective instability. The analysis is carried out for a range of permeability values. The thermal buoyancy is responsible for driving the flow through the medium. The aim of this investigation is to study the effect of viscous dissipation, Soret parameter and externally supplied heat at the two boundaries, on the instability of the flow through a high permeability porous layer having a horizontal throughflow. The analysis also includes the change in the behavioral pattern of the convective rolls, when the different combinations of the perfectly conducting boundary and the boundary with constant heat flux are used. The applications to this work can be found in the field of geophysics and applied mathematics.

A horizontal porous layer saturated with a newtonian fluid is considered. It is bounded by two impermeable, rigid and imperfectly conducting boundaries situated at $z¯=0$ and $z¯=L$ as shown in Fig. 1. A basic temperature difference, $ΔT¯$⁠, and a basic solute concentration difference, $ΔC¯$⁠, are maintained across these boundaries. The external heat is supplied at both the boundaries, where the coefficients of external heat transfer at the two boundaries are different. A horizontal throughflow is considered inside the medium. The viscous dissipation and the Soret effect are also assumed to persist inside the medium. The fluid is considered to obey the linear Oberbeck-Boussinesq approximation. There is a local thermal equilibrium between the two phases of the medium. The Brinkman extension to the Darcy’s law is adopted for the momentum balance and the energy generated due to viscous dissipation inside the medium. Using all these assumptions and approximations, the set of equations governing the flow inside the medium is given by

subject to the boundary conditions

where $ΔT¯=q0Lk$⁠, q₀ is the constant heat flux at the two boundaries, k is the conductivity of the medium, $u¯=(ū,v¯,w¯)$ is the velocity along $x¯=(x¯,ȳ,z¯)$ direction, $T0¯$ and $C0¯$ are the reference temperature and reference concentration, respectively, σ is the heat capacity ratio, α is the thermal diffusivity, ρ is the reference density at temperature T₀, c is the specific heat per unit mass of the fluid, β_T and β_C are the thermal and solutal expansion coefficients, respectively, h_l and h_u are the coefficients of external heat transfer at the lower and upper boundaries, respectively, μ is the dynamic viscosity, μ′ is the effective dynamic viscosity, D_m is the solutal diffusivity, D_CT is the Soret coefficient, and K is the permeability of the porous medium. Here, overbar represents the dimensional quantity.

The scaling variables used to parameterize the equations governing flow in the porous medium are as follows:

The dimensionless governing equations are given as

subject to the dimensionless boundary conditions

where Ge is the Gebhart number, Ra and Sa are the thermal and solutal Rayleigh numbers, respectively, Sr is the Soret parameter, Le is the Lewis number, B₀ and B₁ are the Biot numbers for the lower and upper boundaries, respectively, and Da is the effective Darcy number. They are given by

It makes some sense to take the limiting cases of the thermal boundary conditions at the two boundaries. Hence, by taking the different combinations of the limiting values of B₀ and B₁, it is possible to get the four submodels out of a single model. The thermal boundaries of the flour submodels are as follows

• Model 1: B₀ → 0 and B₁ → 0 (boundaries with constant heat flux)z = 0, 1: $∂T∂z=−1$⁠.

• Model 2: B₀ → 0 and B₁ → ∞ (lower boundary with constant heat flux and isothermal upper boundary)z = 0: $∂T∂z=−1$⁠,z = 1: T = 0.

• Model 3: B₀ → ∞ and B₁ → 0 (isothermal lower boundary and upper boundary with constant heat flux)z = 0: T = 1,z = 1: $∂T∂z=−1$⁠.

• Model 4: B₀ → ∞ and B₁ → ∞ (isothermal boundaries)

z = 0: T = 1,z = 1: T = 0.

It is assumed that there is an average mass flow in the horizontal xy- plane along the direction of a unit vector, s, inclined at an angle, γ to the x-axis. The magnitude of the average mass flow is given by

where Pe is the Péclet number, u_B is the dimensionless basic velocity such that

where u_B is the magnitude of the basic velocity and

In order to maintain the equilibrium inside the system, the basic temperature and concentration are assumed to be of the form, such that their gradients are purely vertical. Hence, the governing equations for u_B, T_B and C_B can be reduced from Eqs. (7)–(10) and are given as

subject to the following boundary conditions

On solving the Eqs. (15)–(18) and using Eqn. (13), the basic velocity, temperature and concentration profiles are obtained as

and

where $ξ=12Da$ and is associated to the permeability of the medium. It can be referred as a switching parameter or the Brinkman parameter. Here, ξ → ∞ refers to the Darcy flow regime and ξ → 0 refers to the clear fluid regime. However, the case of ξ → 0 is avoided in the present investigation, for the viscous dissipation term (in Eqs. (8) and (9)) corresponding to ξ → 0 is not exactly the same as the viscous dissipation term in the Navier-Stokes equation for the clear fluid flow. This fact is supported by the results produced in the analysis by Barletta et al.⁵ 

The basic velocity and concentration fields remain unaltered in the limiting submodels. However, the basic temperature field changes in the four limiting submodels such that the temperature gradients in the four submodels are given as follows

• Model 1: B₀ → 0 and B₁ → 0

• Model 2: B₀ → 0 and B₁ → ∞

• Model 3: B₀ → ∞ and B₁ → 0

• Model 4: B₀ → ∞ and B₁ → ∞

The basic temperature gradients given by Eqs. (22)–(25) are obtained by taking the limits of the gradient of the basic temperature field given by Eqn. (20) at the corresponding limiting values of B₀ and B₁ mentioned in Models 1–4, respectively. If the set of Eqs. (15)–(17) subjected to the four set of boundary conditions for the four submodels mentioned in Sec. II A are solved, the corresponding temperature gradients obtained for the Models 2–4 are same as the respective gradients given by Eqs. (23)–(25). However, the temperature gradient given by Eqn. (22) is a mere approximation obtained to discuss the limiting case with the constant heat flux at both the boundaries. The actual solution of Eqs. (15)–(17) subjected to the boundary condition for Model 1. gives us no solution, since for some branches of the general solution, the given boundary conditions lead to an empty solution.

Small disturbances are superimposed on the basic flow in order to investigate the onset of convective instability. Hence, the velocity, temperature and concentration fields are given as

where U = (U, V, W), θ and χ are the disturbance functions for velocity, temperature and concentration, respectively and ε is a very small positive quantity.

Considering only the O(ε) terms, the linearized governing equations are given as

subject to the boundary conditions

The disturbances in the flow are assumed in the form of two dimensional oblique structures inclined to the base flow at an angle, $γ(0≤γ≤π2)$⁠. The values of γ at the two extreme inclinations correspond to the transverse rolls (γ = 0) and the longitudinal rolls $(γ=π2)$⁠. Hence, the disturbance functions can be considered as the functions of x, z and t and are given as

Introducing the stream function formulation, the velocity components can be written as

Thus, the new set of equations for the oblique structures is given as

Wave like solutions of ψ, θ, and χ are assumed for the normal mode analysis and are expressed as

where λ = λ_R + iλ_I is the growth rate, and λ_R and λ_I are the real and imaginary parts of λ, respectively and a is the wave number. Here, f(z), g(z) and h(z) are complex valued functions.

The imaginary part of the growth rate, λ_I is taken as zero, following which the oblique structures neither grow nor decay. The real part of the growth rate, λ_R is considered as non-zero, due to which the parallel structures oscillate with time. In that case, the final set of governing equations for the general oblique structures is given as

subject to the boundary conditions

Here, ω = λ_R, which in general, can be non-zero when λ_I vanishes. For both the cases, when the base flow gets destabilized for ω ≠ 0, and for ω = 0, the corresponding structures are described by Eqn. (37). The present investigation is mainly focussed on the non-decaying disturbances that are non-oscillatory in time. So, the principle of exchange of stabilities is considered to hold, for which the problem must be governed by the self-adjoint system of differential equations. Thus, for the principle of exchange of stabilities to hold, ω should be equal to the term au_Bcos(γ) for the general oblique rolls. Following we have, ω = 0 for the longitudinal rolls (⁠$γ=π2$⁠) and ω = au_B for the transverse rolls (γ = 0).

The set of Eqs. (38)–(40) along with the boundary conditions given by Eqn. (41) admit a non-zero solution only for certain values of ω and Ra. Thus, it could be considered as an eigenvalue problem for the general oblique rolls, where (Ra, ω) is the set of eigenvalues. The eigenvalue Ra is a function of the wave number, a. As usual, the eigenfunctions are defined except for a normalization condition. For the numerical computation, a normalization conditions is defined which is given by

The set of higher order ordinary differential equations is further converted into a set of first order ordinary differential equations. On separating the real and imaginary parts of the equations along with the boundary conditions, an eigenvalue problem is obtained, consisting of sixteen first order ordinary differential equations with eighteen variables subject to eighteen boundary conditions. Finally, the eigenvalue problem is solved numerically using bvp4c in Matlab R2016a for each set of assigned values of a, γ, Pe, Ge, Sa, Le, Sr, B₀, B₁ and ξ. The relative tolerance has been taken as RelTol = 10⁻⁶, and the absolute tolerance as AbsTol = 10⁻⁹ to obtain all the data values displayed throughout this article. The results are obtained which consist of the eigenfunctions f(z), g(z) and h(z) along with the parameter Ra(a). The critical value of the eigenvalue is obtained by taking the minimum of the function Ra(a) at some critical wave number, a_C.

The validity of the numerical technique used is established by comparing the present results with the results obtained by other authors using some well established numerical techniques. The limiting case of the Darcy flow regime (ξ → ∞) with both the boundaries as isothermal (B₀, B₁ → ∞), no viscous dissipation (Ge = 0) and no throughflow (Pe = 0) is considered for this comparison. The presence of solute is also ignored and hence the Soret effect is also considered to be absent. In this case, Ra_c = 39.4784 is obtained at a_c = 3.14159 which matches with the corresponding results of Barletta and Nield²³ as well as Kaloni and Qiao²⁴ for the case of longitudinal rolls with no horizontal thermal gradient (Ra_H = 0). Barletta and Nield²³ used the Shooting technique with the fourth order explicit Runge-Kutta method with adaptive step size, whereas, Kaloni and Qiao²⁴ used the compound matrix method to get the results. It can even be compared with the results by Guo and Kaloni²⁵ (by using compound matrix method), for different values of Sa, as shown in Table I. It is observed that the present results are in excellent agreement with the existing ones.

Tables II and III illustrate the effect of viscous dissipation on the longitudinal rolls in different flow regimes. ξ = 1.58114, 0.5 and 0.158114 (Da = 10⁻¹, 1 and 10, respectively) are taken to represent different flow regimes under the effect of external heating at the boundaries. The other physical parameters are fixed as Pe = 10, Sa = 1 and Le = 1 with varying B₀ and B₁. The Soret parameter is suppressed for this analysis, in order to study the sole effect of viscous dissipation on the instability of the flow. The values of B₀ and B₁ are taken as 0.01, 0.1, 1, 10 and 100 (Alves and Barletta²²). It is observed that the viscous dissipation has a stabilizing effect on the flow as long as the amount of external heat supplied at the bottom boundary is higher than is supplied at the upper boundary, in any flow regime. The increase in the coefficient of external heating at the upper boundary over that at the lower boundary opposes the direction of thermal buoyancy, which is responsible for the onset of convective instability. Viscous dissipation has a destabilizing effect (very small though) on the flow, when the amount of heat supplied at the two boundaries are equal. There is no qualitative change in this behavior with the changing flow regime.

It is found important to discuss about the inclination angle at which the two dimensional convective rolls are most sound to instability. Different combinations of the values of B₀ and B₁ are considered in order to see the behavior of the rolls with the change in the coefficient of external heating at the two boundaries. The values of the other physical parameters are taken as Ge = 0.1, Sa = 1, Le = 1, Sr = 0.1 and Pe = 10, while considering in the different flow regimes (ξ = 0.05, 0.5, 5 and ∞). Fig. 2 clearly indicates that the longitudinal rolls at the inclination, $γ=π2$⁠, are the preferred mode of instability in all the flow regimes represented by the different values of the Brinkman parameter, ξ. However, in the Darcy flow regime (ξ → ∞), all the oblique rolls at different inclinations to the base flow are almost equally unstable. It can also be observed that the oblique rolls at different inclinations have qualitatively similar response towards the instability, with the increase in the coefficient of external heating at the two boundaries. However, the quantitative response varies a bit in the Brinkman flow regime, depending upon the combination of the values of B₀ and B₁. It can analyzed from the curves in Figs. 2(a) and 2(b) that all the oblique rolls at any arbitrary inclination, in the Brinkman flow regime or transition flow regime, tend to become stable with the increasing values of the function, F(B₀, B₁) = B₀ + B₁ + B₀B₁.

The change in the instability condition with the changing flow regime is shown by Fig. 3 for different values of the Soret parameter taken as Sr = −0.1, 0 and 0.1 in Figs. 3(a), 3(b) and 3(c), respectively. The considered values of ξ range from 0 to 15, whereas, the values of other parameters are taken as same as in Fig. 2. The analysis is done for the longitudinal rolls only. The parameter which is chosen to measure the onset condition for instability is log₁₀ Ra_C. This choice is justified, since log₁₀ Ra_C monotonically increases with increasing values of Ra_C. This is done just to minimize the scale of the vertical coordinates of the plots. It is observed that the flow is most stable in the Brinkman regime, which is represented by the values of ξ close to 0. The instability gradually increases as the value of ξ increases. This means that the flow destabilizes as one moves towards the Darcy flow regime from the Brinkman flow regime. This property remains unaffected by the presence of the Soret parameter. The magnified view of the curves in each subfigure confirms that the periodic convective rolls stabilize with the increasing value of the function, F(B₀, B₁) = B₀ + B₁ + B₀B₁ in the Brinkman flow regime or the transition flow regime.

Fig. 4 illustrates the effect of coefficients of external heating on the instability of the longitudinal rolls in different flow regimes represented by ξ = 0.05, 0.5, 5 and ∞ in Figs. 4(a)– 4(d), respectively. The values of the physical parameters are considered as Ge = 0.1, Pe = 10, Le = 1, Sa = 1, Sr = 0.1 and B₁ = 0, 0.1, 1, 5 and 10. Each curve in the figure shows the variation of Ra_C against B₀ for a fixed value of B₁. It can be inferred from the figure that the stability increases with the increasing values of the coefficient of external heating at the lower boundary, while maintaining a fixed amount of externally supplied heat at the upper boundary. This property remains unaltered in all the flow regimes. However, the scenario changes, if the variation of Ra_C against B₁ is observed, for a fixed value of B₀. In that case, the system stabilizes with the increasing value of the coefficient of external heating at the upper boundary in the Brinkman flow regime (represented by Figs. 4(a)–4(b)). On the other hand, in the Darcy flow regime (represented by Figs. 4(c)–4(d)), the instability varies non-monotonically with the increasing values of B₁ up to a certain range of values of B₀, and it decreases further with the increasing values of B₁. This irregularity in the onset condition is observed due to the presence of solute concentration gradient in the upward direction. However, the effect of solute concentration gradient is dominated by the effect of viscous stresses present in the Brinkman flow regime. Hence, this irregularity in the onset condition of instability is not seen in the Brinkman flow regime.

The effect of solute concentration gradient on the instability of the longitudinal rolls is illustrated by Figs. 5(a)–5(c) for Le = 0.1, 1 and 10, respectively. The analysis is done for the Brinkman regime (ξ = 0.5) only, since there is no qualitative change in the behavior of the curves in the Darcy regime. Same set of values of rest of the parameters as earlier are considered. The figures depict that under the effect of the Soret parameter and viscous dissipation, the solute concentration gradient has a linearly destabilizing effect on the flow. This destabilizing nature of Sa remains same for all the values of the diffusivity ratio represented by Le < 1, Le = 1 and Le > 1. The Soret parameter also has a linearly stabilizing or destabilizing effect on the flow, depending upon the sign of solute concentration gradient, as evident in Fig. 6. Figs. 6(a) and Fig. 6(b) also have the same plot legends as in Fig. 5. When the solute concentration at the upper boundary is higher than that at the lower boundary, the Soret parameter has a linearly destabilizing effect on the flow, otherwise, it has a linearly stabilizing effect on the flow.

The streamlines (solid), isotherms (dashed) and isosolutes (dotted) at $γ=π2$ in the Brinkman flow regime for the four limiting Models 1–4 are shown in Figs. 7–10, respectively. The streamlines, isotherm and isosolute patterns are represented for Sr = −0.1, 0 and 0.1 at the onset condition. The range of the z coordinate is taken from 0 to 1, whereas, the range of y coordinate is taken from 0 to $πaC$⁠. Here, a_C is the critical value of the wave number computed for a given set of data values. This means, the higher the range of y coordinate, smaller is the value of a_C. The streamlines pattern in all the four models are concentrated towards the center of the medium due to the boundary effects in the Brinkman flow regime and are symmetric about $y=π2aC$ and z = 1/2. The isotherm patterns are different in all the four models.

For Model 1.(Fig. 7), representing the medium bounded by the two boundaries with constant heat flux condition, the isotherms are the vertical lines, and are symmetric about the vertical and horizontal mid-axes. For Model 2.(Fig. 8), representing the medium bounded by the lower boundary with constant heat flux and isothermal upper boundary, the periodic isotherm patterns are symmetric about $y=π2aC$ and are formed towards the lower boundary which is having the constant heat flux condition. For Model 3.(Fig. 9), representing the medium bounded by isothermal lower boundary and the upper boundary with constant heat flux, the periodic isotherm patterns are symmetric about $y=π2aC$ and are formed towards the upper boundary which is having the constant heat flux condition. For Model 4.(Fig. 10), representing the medium bounded by the two isothermal boundaries, the periodic isotherm patterns are symmetric about $y=π2aC$ and z = 1/2 both and are concentrated towards the center. The isosolutes for all the models are the periodic patterns symmetric about $y=π2aC$⁠. However, there is no net qualitative change in the isosolute patterns with the changing values of the Soret parameters.

The onset of double-diffusive convective instability in a horizontal throughflow along a porous layer confined between two rigid boundaries with a basic temperature difference and solute concentration difference, is investigated. The effect of viscous dissipation, coefficient of external heating at the two boundaries, and the Soret parameter, on the instability of the base flow caused due to the superimposed disturbances in the form of two dimensional convective rolls, are studied. The analysis is done for different flow regimes and for the oblique rolls at different inclinations. The conclusions drawn from the analysis are as follows:

• The viscous dissipation has a stabilizing effect on the flow, as long as the coefficient of external heating at the bottom boundary is higher than that at the upper boundary. It has a destabilizing effect in case of the equal or higher amount of external heating at the upper boundary than at the bottom boundary.

• The increase in the amount of heat supplied at the lower boundary stabilizes the flow. As long as the external heat supplied at the upper and lower boundaries are same, there is no appreciable effect of external heating on the instability of the medium.

• The flow subjected to the external heating at the two boundaries, is most stable in the Brinkman flow regime and the oblique rolls inclined at an angle $π2$ to the base flow are the most unstable ones.

• The solute concentration gradient has a linearly destabilizing effect on the base flow in all the possible situations. On the other hand, the Soret parameter has a linearly destabilizing effect on the flow, when the solute concentration at the upper boundary is higher than that at the lower boundary, otherwise, it has a linearly stabilizing effect on the flow.

• Under a given condition, when the throughflow is considered in the Brinkman flow regime, while considering the effect of viscous dissipation, Soret parameter, Solute concentration gradient and the external heating at the two boundaries, the stability of the flow increases with the increasing values of the function F(B₀, B₁) = B₀ + B₁ + B₀B₁.