Soret Effect on the Instability of Double-Diffusive Convection in a Saturated Vertical Brinkman Porous Layer of Oldroyd-B Fluid

1. Introduction

A solid with micropores are called porous media, and these micropores are often interconnected and usually filled with fluid. Porous media encompasses natural substances, such as biological tissues, wood, zeolites, soils, and rocks. Due to its numerous fundamental and industrial applications, double diffusive convection (DDC) in porous media is a phenomenon that occurs in a variety of systems and has attracted a lot of attention in recent decades. These applications include producing high-quality crystals, storing liquid gases, migrating moisture in fiber insulation, transporting contaminants in saturated soils, solidifying molten alloys, and heating lakes and magmas through geothermal processes. Ingham and Pop [1], Nield and Bejan [2], and Vafai [3] have provided comprehensive reviews on this subject.

The stability or instability of fluid flow is a crucial aspect of fluid mechanics in porous media. Research on how buoyancy and shearing forces affect fluid flow stability in a vertical layer of porous media has received significant attention. Gill [4] performed an initial analysis of the stability of natural convection in a vertical porous layer whose boundary is impermeable at different temperatures, using Darcy's law as a framework. He concluded that the system always maintains linear stability. Subsequently, Rees [5] added a temporal derivative of the velocity to the momentum equation and determined that the flow remains linearly stable, which is the same as Gill's conclusion. For porous media with high porosity, Lundgren [6] experimentally confirmed that the extension of Darcy's law by Brinkman is more valid. Shankar et al. [7] investigated the effect of inertial terms on the stability of natural convection in a vertical layer of a porous medium at different temperatures using the Brinkman model. Using the classical linear stability theory, they discovered that instability is caused by the inertia effect.

Viscoelastic fluids naturally exhibit convection when flowing through porous media, which has important implications for reservoir engineering, bioengineering and geophysics. Gözüm and Arpaci [8] initially investigated the stability of natural convection in a vertical layer at various temperatures for a viscoelastic Maxwell fluid. Takashima [9] investigated the same problem with Oldroyd-B fluid. Khuzhayorov et al. [10] proposed a revised Darcy's law for analyzing the viscoelastic properties of saturated porous materials. Using the homogenization method, they discovered a general filtration law that describes the flow of a linear viscoelastic fluid in a porous material. Kim et al. [11] used a modified Darcy model to analyze thermal instability in a horizontal porous layer with saturated Oldroyd-B fluid. The conventional linear stability theory provided the critical conditions for the initiation of convective motion. Using a modified Darcy-Brinkman-Oldroyd model, Zhang et al. [12] investigated the convection of saturated Oldroyd-B fluid in a horizontal porous layer heated from below. They calculated the critical number of stationary and oscillating convection. Sun et al. [13] studied the lower horizontal plate that was heated in a time-periodic pattern. The system changed from a steady convection to a periodic or chaotic state. Barletta and Alves [14] investigated the Gill stability problem for power-law fluids and discovered that Gill's findings are also applicable to the power-law fluid. The Gill problem for Oldroyd-B fluid was addressed by Shankar and Shivakumara [15]. In contrast to the results for Newtonian and power-law fluids, the flow of Oldroyd-B fluid is unstable.

Besides, they [16] considered the situation of local thermal non-equilibrium (LTNE) and instability that existed as well. Then, Shanker and Shivakumara [17] considered the case of a porous layer containing an internal heat source. He established that internal heating and relaxation parameters contribute to system instability. Newtonian fluids flow steady, while Oldroyd-B fluids flow unstable. For double-diffusive convection instead of thermal convection, Wang and Tan [18] used a modified Darcy model to study the stability of Maxwell fluids in porous media with two parallel planes.

Malashetty and Biradar [19] investigated the cross-diffusion effect on DDC. They conducted linear and weakly nonlinear stability analyses using a modified Darcy model with a time derivative term as the momentum equation and found the critical Rayleigh number. Malashetty et al. [20,21] then studied the stability of the Oldroyd-B fluid and the formation of DDC saturated isotropic and anisotropic horizontal porous layers, respectively. As a result, the thermal anisotropy parameter has a dual role in the stability of the flow. Kumar and Bhadauria [22] explored the situation of LTNE, adding a time derivative term to the momentum equation. They discovered that when the interphase heat transfer coefficient was big or small, the system behaved similarly to the local thermal equilibrium (LTE) model. Subsequently, Malashetty et al. [23] investigated the situation of an anisotropic rotating porous layer with LTE and found that the thermal anisotropy parameter has the opposite effect on the onset of convection compared with the no-rotation condition. The onset of DDC in a horizontal porous layer was recently studied by Swamy et al. [24] based on the Darcy-Brinkman-Oldroyd model.

Due to numerous applications in geothermal systems, energy storage devices, thermal insulation, drying technologies, catalytic reactors, and nuclear waste repositories, theoretical and experimental research has focused on heat and mass transmission in porous systems. Density gradients in fluid-saturated media can cause heat and mass transmission to be coupled because of heat and mass inhomogeneities when examining heat and mass transfer processes in porous media channels. The cross-diffusion effect is simultaneously brought on by mass and heat fluxes. The Soret effect transfers mass via a temperature gradient, while the Dufour effect transfers heat via a concentration gradient. The Dufour coefficient has a negligible energy flux in liquids and is an order of magnitude smaller than the Soret coefficient (see Straughan and Hutter [25]). Therefore, while discussing liquid flow, we can ignore the Dufour term. The majority of recent research on DDC in porous layers that takes the Soret effect into account has concentrated on horizontal porous layers. With the Darcy model, Bahloul et al. [26] examined the beginning of natural convection in a horizontal porous layer. They calculated the critical values of finite amplitude, oscillatory, and monotonic convective instability for DDC and Soret convection using linear stability analysis. Subsequently, the case of anisotropic horizontal porous layer based on the previous paper was studied by Gaikwad et al. [27] Using linear and nonlinear stability analysis, they analyzed the critical Rayleigh number, wave number, and oscillation frequency of steady and oscillatory modes. Utilizing linear and weakly nonlinear stability studies, Gaikwad and Dhanraj [28] initially used the Darcy-Brinkman model to study the DDC with the Soret effect in a binary viscoelastic fluid-saturated horizontally porous layer. They discovered that whereas a positive Soret value accelerates the onset of DDC in stationary mode and has the opposite impact in oscillatory and finite amplitude modes, a negative Soret parameter increases system stability. Bouachir et al. [29] investigated DDC flow in a vertical porous chamber filled with a binary mixture exhibiting Soret and Dufour effects. The Darcy-Brinkman model and the Oberbeck-Boussinesq approximation are used to study the impact of the Soret and Dufour effects on convective stability. Overall, the thresholds of oscillatory, overstable, and stationary convection were greatly impacted by the Soret and Dufour effects.

Due to the significant differences in molecular diffusivity, DDC leads to complex flow structures and, consequently, the transport of heat and solute concentrations with different time and length scales. The majority of prior investigations used pure DDC for natural convection in a vertical porous material influenced by horizontal heat and solute concentration differences. Therefore, the Soret effect is introduced in this paper to describe the mass flux induced by the temperature gradient. Until now, there has been no study on the instability of the Soret effect on DDC in a saturated vertical porous layer of Oldroyd-B fluid. Therefore, the present work aims to investigate the linear stability of DDC of Oldroyd-B fluid in a vertical porous layer with the Soret effect. The manuscript is structured as follows. Section 2 presents the mathematical model that outlines the governing equations. The linear stability analysis and numerical procedures are given in Section 3 and Section 4, respectively. In Section 5, we present the results and discussion, while the last section contains the conclusions that we have drawn.

2. Mathematical Model

2.1. Mathematical Description of the Model

Figure 1 depicts an infinitely extended vertical porous layer surrounded by two impermeable plates and saturated with Oldroyd-B fluid in the region −d ≤ x ≤ d, with vertically downward gravity acting on it. The vertical limits of the porous layer are maintained at constant but distinct temperatures T₁ and T₂ (> T₁) and solute concentrations S₁ and S₂ (>S₁). We begin with a Cartesian coordinate frame, in which the z-axis points vertically and the x-axis extends horizontally over the breadth.

Given a temperature T and a solute concentration S, the density ρ of the mixture fluid is expected to change linearly in the following form:

$${\rho }^{*}={\rho }_0\left[1-\alpha \left(T^{*}-T_0\right)+\beta \left(S^{*}-S_0\right)\right],$$

(1)

The volumetric thermal expansion coefficient α and volumetric solutal expansion coefficient β of the fluid are represented by the variables ρ₀ at T₀=(T₁+T₂)/2 and S₀=(S₁+S₂)/2.

It is assumed that the fluid and medium have constant properties and that few chemical reactions occur in the homogeneous, isotropic media. Furthermore, the porous media exhibits a local thermal equilibrium. After adopting the modified Darcy-Brinkman-Oldroyd model [24], which incorporates the Oberbeck-Boussinesq approximation and a time derivative term, the system governing equations are as follows

$${\nabla }^{*}\cdot u^{*}=0,$$

(2)

$$\left(1+{\lambda }_1^{*}{\displaystyle \frac{\partial }{\partial t^{*}}}\right)\left({\displaystyle \frac{{\rho }_0}{\epsilon }}{\displaystyle \frac{\partial u^{*}}{\partial t^{*}}}+{\nabla }^{*}{p}^{*}-{\rho }^{*}g\right)=\left(1+{\lambda }_2^{*}{\displaystyle \frac{\partial }{\partial t^{*}}}\right)\left(-{\displaystyle \frac{\mu }{K}}{u}^{*}+{\mu }_e{\nabla }^{*2}{u}^{*}\right),$$

(3)

$$\gamma {\displaystyle \frac{\partial T^{*}}{\partial t^{*}}}+\left(u^{*}\cdot {\nabla }^{*}\right)T^{*}={\kappa }_T{\nabla }^{*2}{T}^{*},$$

(4)

$$\epsilon {\displaystyle \frac{\partial S^{*}}{\partial t^{*}}}+\left(u^{*}\cdot {\nabla }^{*}\right)S^{*}={\kappa }_S{\nabla }^{*2}{S}^{*}+{\kappa }_{ST}{\nabla }^{*2}{T}^{*},$$

(5)

where t^*, p^*, and T^* stand for time, pressure, and temperature in the porous media, respectively, and u^*= (u^*, v^*, w^*) is the velocity vector. The symbols for dynamic viscosity and effective viscosity are μ and μ_e. The medium's porosity and permeability are represented by K and ε. The strain retardation time is represented by λ₂^*, the stress relaxation time by λ₁^*, and the gravity acceleration vector is given by g= (0, 0, -g). κ_ST is the Soret coefficient expressed as a constant, and κ_T and κ_S stand for thermal and solute diffusivity, respectively. The ratio of specific heats is represented by γ=(ρc)_m /(ρc)_f, where (ρc)_f represents the volumetric heat capacity of the fluid and (ρc)_m =ε(ρc)_f + (1-ε)(ρc)_s represents the saturated medium's overall volumetric heat capacity. The subscripts f and s stand for the fluid and solid matrix of the porous medium.

The boundaries are impermeable, and the applicable boundary conditions are

$$u^{*}=0,T^{*}=T_1,S^{*}=S_1\mathrm{at} x^{*}=-d,$$

(6)

$$u^{*}=0,T^{*}=T_2,S^{*}=S_2\mathrm{at} x^{*}=d.$$

(7)

Dimensionless quantities are introduced through the scaling

$$(x^{\ast },y^{\ast },z^{\ast })=(x,y,z)d{t}^{*}=(\gamma d^2/{\kappa }_T)t\hspace{0.33em}{u}^{*}=({\kappa }_T/d)u{p}^{*}=\left(\mu {\kappa }_T/K\right)p$$

(8)

$${\lambda }_1^{*}=(\gamma d^2/{\kappa }_T){\lambda }_1{\lambda }_2^{*}=(\gamma d^2/{\kappa }_T){\lambda }_2{T}^{\ast }-T_0=(\Delta T)T{S}^{\ast }-S_0=(\Delta S)S$$

to non-dimensionalize Eqs. (2) -(7) in the form,

$$\nabla \cdot u=0,$$

(9)

$$\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left({\displaystyle \frac{1}{P{r}_D}}{\displaystyle \frac{\partial u}{\partial t}}+\nabla p-R{a}_{T}T{\widehat{e}}_z+R{a}_{S}S{\widehat{e}}_z\right)=\left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\left(-u+Da{\nabla }^{2}u\right),$$

(10)

$${\displaystyle \frac{\partial T}{\partial t}}+\left(u\cdot \nabla \right)T={\nabla }^{2}T$$

(11)

$$\eta {\displaystyle \frac{\partial S}{\partial t}}+\left(u\cdot \nabla \right)S={\displaystyle \frac{1}{Le}}{\nabla }^{2}S+S_r{\nabla }^{2}T$$

(12)

here, Pr_D=εγd²v/κ_TK is the Darcy–Prandtl number, Ra_T=αg∆TdK/vκ_T is the Darcy–Rayleigh number, Ra_S=βg∆SdK/vκ_T is solute Rayleigh number, Le=κ_T/κ_S is the Lewis number and η=ε/γ is normalized porosity. Da=μ_eK/μd² is the modified Darcy number. S_r=κ_ST∆T /κ_T∆S is the Soret parameter. There are two possible values for the Soret parameter S_r: positive and negative. The solute diffuses toward the cooler plate when the Soret parameter is positive; the opposite is true when the parameter is negative.

The boundary conditions are

$$u=0,\mathrm{at}x=\pm 1$$

(13)

$$T=\pm {\displaystyle \frac{1}{2}},S=\pm {\displaystyle \frac{1}{2}},\mathrm{at}x=\pm 1$$

(14)

2.2. Basic Flow

It is assumed that the base flow is stable, unidirectional, and fully developed. In this condition, the physical variables are provided by

$$u_b=(0,0,w_b(x)),T=T_b(x),S=S_b(x),p=p_b(x)$$

(15)

Applying the preceding assumptions to the governing equations (9)-(14) reduces the system of ordinary differential equations and easily obtains the basic flow equations:

$$Da{D}^2{w}_b-w_b+R{a}_T{T}_b-R{a}_S{S}_b={\displaystyle \frac{\partial p_b}{\partial z}},$$

(16)

$$D^2{T}_b=0,$$

(17)

$${\displaystyle \frac{D^2{S}_b}{Le}}+S_r{\displaystyle \frac{{\partial }^2{T}_b}{\partial x^2}}=0,$$

(18)

where D=d/dx and the suffix b serves to denote the basic flow. The associated boundary conditions are

$$w_b=0, \mathrm{at} x=\pm 1,$$

(19)

$$T_b=\pm {\displaystyle \frac{1}{2}},S_b=\pm {\displaystyle \frac{1}{2}}, \mathrm{at} x=\pm 1.$$

(20)

The solution to the basic flow is found to be

$$w_b(x)={\displaystyle \frac{R{a}_S-R{a}_T}{2}}\left[{\displaystyle \frac{\sinh (D{a}^{-\frac{1}{2}}x)}{\sinh D{a}^{-\frac{1}{2}}}}-x\right],T_b={\displaystyle \frac{1}{2}}x,S_b={\displaystyle \frac{1}{2}}x,p_b(z)=\mathrm{const}.$$

(21a,b,c,d)

Ra_T, Ra_S, and Da all have an impact on the fundamental velocity profile. When Ra_S →0, it is possible to quickly verify that the basic solution provided by Eqs. (21a, b) coincides with that of Shankar et al. [31]. The base flow velocity profiles are shown in Figure 2(a)–(c). Figure 2(a) shows scaled basic velocity profiles w_b(x) for various Da values. The no-slip condition on the wall and the buoyant term's acceleration of the fluid near the border cause a point of inflection in the velocity profile, as seen in the image. With decreasing Da, it is found that the point of inflection moves closer to the walls. The appearance of an inflection point indicates a potential basis for flow instability. Ra_T and Ra_S have different effects on the base flow, as seen in Figure 2(b) and (c), respectively. As Ra_T or Ra_S rises, the velocity value rises as well. The velocity profiles are significantly asymmetric when compared to the vertical line at x = 0.

3. Linear Stability Analysis

In order to examine linear stability, we must first develop the perturbation equations for the infinitesimal disturbance. Small-amplitude disturbances of the following form are applied to the basic state:

$$u=u^{\prime },v=v^{\prime },w=w_b+w^{\prime },p=p_b+p^{\prime },T=T_b+T^{\prime },S=S_b+S^{\prime },$$

(22)

where the primed amounts denote the perturbations field across the basic state.

The three-dimensional governing linear stability equations are obtained by applying Eq. (22) in Eqs. (9)-(14) and by following the standard process of linear stability theory [30]:

$${\displaystyle \frac{\partial u^{\prime }}{\partial x}}+{\displaystyle \frac{\partial v^{\prime }}{\partial y}}+{\displaystyle \frac{\partial w^{\prime }}{\partial z}}=0,$$

(23)

$$\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left({\displaystyle \frac{1}{P{r}_D}}{\displaystyle \frac{\partial u^{\prime }}{\partial t}}+{\displaystyle \frac{\partial p^{\prime }}{\partial x}}\right)=\left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\left[-u^{\prime }+Da\left({\displaystyle \frac{{\partial }^2{u}^{\prime }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{u}^{\prime }}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{u}^{\prime }}{\partial z^2}}\right)\right],$$

(24)

$$\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left({\displaystyle \frac{1}{P{r}_D}}{\displaystyle \frac{\partial v^{\prime }}{\partial t}}+{\displaystyle \frac{\partial p^{\prime }}{\partial y}}\right)=\left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\left[-v^{\prime }+Da\left({\displaystyle \frac{{\partial }^2{v}^{\prime }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{v}^{\prime }}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{v}^{\prime }}{\partial z^2}}\right)\right],$$

(25)

$$\begin{array}{c}\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left({\displaystyle \frac{1}{P{r}_D}}{\displaystyle \frac{\partial w^{\prime }}{\partial t}}+{\displaystyle \frac{\partial p^{\prime }}{\partial z}}-R{a}_T{T}^{\prime }+R{a}_S{S}^{\prime }\right)\\ =\left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\left[-w^{\prime }+D\mathrm{a}\left({\displaystyle \frac{{\partial }^2{w}^{\prime }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{w}^{\prime }}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{w}^{\prime }}{\partial z^2}}\right)\right],\end{array}$$

(26)

$${\displaystyle \frac{\partial T^{\prime }}{\partial t}}+u^{\prime }D{T}_b+w_b{\displaystyle \frac{\partial T^{\prime }}{\partial z}}={\displaystyle \frac{{\partial }^2{T}^{\prime }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{T}^{\prime }}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{T}^{\prime }}{\partial z^2}}$$

(27)

$$\eta {\displaystyle \frac{\partial S^{\prime }}{\partial t}}+u^{\prime }D{S}_b+w_b{\displaystyle \frac{\partial S^{\prime }}{\partial z}}={\displaystyle \frac{1}{Le}}\left({\displaystyle \frac{{\partial }^2{S}^{\prime }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{S}^{\prime }}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{S}^{\prime }}{\partial z^2}}\right)+S_r\left({\displaystyle \frac{{\partial }^2{T}^{\prime }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{T}^{\prime }}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{T}^{\prime }}{\partial z^2}}\right)$$

(28)

The boundary conditions are

$$u^{\prime }=v^{\prime }=w^{\prime }=T^{\prime }=S^{\prime }=0x=\pm 1$$

(29)

To construct linear stability for arbitrary infinitesimal disturbances, superpose normal modes in the following manner:

$$\left\{u^{\prime },v^{\prime },w^{\prime },p^{\prime },T^{\prime },S^{\prime }\right\}=\left\{\stackrel{⌢}{u},\stackrel{⌢}{v},\stackrel{⌢}{w},\stackrel{⌢}{p},\stackrel{⌢}{T},\stackrel{⌢}{S}\right\}\left(x\right)e^{i(az+by-act)}$$

(30)

where a and b represent the wave numbers in the z and y directions, respectively, and c=c_r+ic_i is the complex parameter whose real part expresses the wave speed and whose imaginary part yields the disturbance growth rate. The basic solution threshold for stability and instability is as follows: if c_i < 0, the system is stable; if c_i > 0, the system is unstable; and if c_i = 0, the system is neutrally stable. When we replace Eq. (30) with Eqs. (23)-(29), we get

$$D\stackrel{⌢}{u}+i(b\stackrel{⌢}{v}+a\stackrel{⌢}{w})=0,$$

(31)

$$(1-aci{\lambda }_1)\left(-{\displaystyle \frac{aci\stackrel{⌢}{u}}{P{r}_D}}+D\stackrel{⌢}{p}\right)=(1-aci{\lambda }_2)\left[Da\left(D^2-a^2-b^2\right)-1\right]\stackrel{⌢}{u},$$

(32)

$$(1-aci{\lambda }_1)\left(-{\displaystyle \frac{aci\stackrel{⌢}{v}}{P{r}_D}}\right)+b\left(i+ac{\lambda }_1\right)\stackrel{⌢}{p}=(1-aci{\lambda }_2)\left[Da\left(D^2-a^2-b^2\right)-1\right]\stackrel{⌢}{v},$$

(33)

$$+a\left(i+ac{\lambda }_1\right)\stackrel{⌢}{p}=(1-aci{\lambda }_2)\left[Da\left(D^2-a^2-b^2\right)-1\right]\stackrel{⌢}{w},$$

(34)

$$(1-aci{\lambda }_1)\left(-{\displaystyle \frac{aci\stackrel{⌢}{w}}{P{r}_D}}-R{a}_T\stackrel{⌢}{T}+R{a}_S\stackrel{⌢}{S}\right)$$

$$ia\left(w_b-c\right)\stackrel{⌢}{T}+\stackrel{⌢}{u}D{T}_b=\left(D^2-a^2-b^2\right)\stackrel{⌢}{T},$$

(35)

$$ia\left(w_b-c\eta \right)\stackrel{⌢}{S}+\stackrel{⌢}{u}D{S}_b={\displaystyle \frac{1}{Le}}\left(D^2-a^2-b^2\right)\stackrel{⌢}{S}+S_r\left(D^2-a^2-b^2\right)\stackrel{⌢}{T},$$

(36)

$$\stackrel{⌢}{u}=\stackrel{⌢}{v}=\stackrel{⌢}{w}=\stackrel{⌢}{T}=\stackrel{⌢}{S}=0, \mathrm{at} x=\pm 1.$$

(37)

We exclusively examine two-dimensional flows due to the validity of Squire's theorem [31]. Thus, employing Squire's transformation in its extended form

$$\begin{array}{c}\tilde{a}=\sqrt{a^2+b^2},\tilde{w}={\displaystyle \frac{a(a\stackrel{⌢}{w}+b\stackrel{⌢}{v})}{{\tilde{a}}^2}},{\tilde{w}}_b={\displaystyle \frac{a}{\tilde{a}}}{w}_b, \tilde{c}={\displaystyle \frac{a}{\tilde{a}}}c, \tilde{T}={\displaystyle \frac{a}{\tilde{a}}}\stackrel{⌢}{T}, \tilde{S}={\displaystyle \frac{a}{\tilde{a}}}\stackrel{⌢}{S}, \tilde{\eta }=\eta , \tilde{L}e=Le,\\ \tilde{p}={\displaystyle \frac{a}{\tilde{a}}}\stackrel{⌢}{p}, \tilde{D}a=Da, \tilde{P}{r}_D=P{r}_D, {\tilde{\lambda }}_1={\lambda }_1, {\tilde{\lambda }}_2={\lambda }_2, \tilde{R}{a}_T={\displaystyle \frac{a}{\tilde{a}}}R{a}_T, \tilde{R}{a}_S={\displaystyle \frac{a}{\tilde{a}}}R{a}_S.\end{array}$$

(38)

replacing Eq. (38) into Eqs. (31)-(37) yields:

$$D\tilde{u}+i\tilde{a}\tilde{w}=0,$$

(39)

$$(1-i\tilde{a}\tilde{c}{\tilde{\lambda }}_1)\left(-{\displaystyle \frac{i\tilde{a}\tilde{c}\tilde{u}}{\tilde{P}{r}_D}}+D\tilde{p}\right)=(1-i\tilde{a}\tilde{c}{\tilde{\lambda }}_2)\left[\tilde{D}a\left(D^2-{\tilde{a}}^2\right)-1\right]\tilde{u},$$

(40)

$$(1-i\tilde{a}\tilde{c}{\tilde{\lambda }}_1)\left(-{\displaystyle \frac{i\tilde{a}\tilde{c}\tilde{w}}{\tilde{P}{r}_D}}-\tilde{R}{a}_T\tilde{T}+\tilde{R}{a}_S\tilde{S}\right)+\tilde{a}\left(i+\tilde{a}\tilde{c}{\tilde{\lambda }}_1\right)\tilde{p}$$

$$=(1-i\tilde{a}\tilde{c}{\tilde{\lambda }}_2)\left[\tilde{D}a\left(D^2-{\tilde{a}}^2\right)-1\right]\tilde{w},$$

(41)

$$i\tilde{a}\left({\tilde{w}}_b-\tilde{c}\right)\tilde{T}+\tilde{u}D{T}_b=\left(D^2-{\tilde{a}}^2\right)\tilde{T},$$

(42)

$$i\tilde{a}\left({\tilde{w}}_b-\tilde{c}\tilde{\eta }\right)\tilde{S}+\tilde{u}D{S}_b={\displaystyle \frac{1}{\tilde{L}e}}\left(D^2-{\tilde{a}}^2\right)\tilde{S}+S_r\left(D^2-{\tilde{a}}^2\right)\tilde{T},$$

(43)

$$\tilde{u}=\tilde{w}=\tilde{T}=\tilde{S}=0, \mathrm{at} x=\pm 1.$$

(44)

Equations (39)-(44) are two-dimensional equations with the same mathematical structure as Eqs. (31)-(37), when v=b=0 (i.e., two-dimensional equations). When b≠0, $\tilde{R}{a}_T$and$\tilde{R}{a}_S$are smaller than Ra_T and Ra_S, respectively. Therefore, it is sufficient to consider only the two-dimensional equations, as two-dimensional perturbations are less favorable for stabilizing the base flow than three-dimensional perturbations. Then, introduce the stream function $\tilde{\psi }$(x, z, t) such that

$$\tilde{u}={\displaystyle \frac{\partial \tilde{\psi }}{\partial z}},\tilde{w}=-{\displaystyle \frac{\partial \tilde{\psi }}{\partial x}},$$

(45)

the following is how the normal mode analysis approach is used:

$$\left\{\tilde{\psi },\tilde{T},\tilde{S}\right\}=\left\{\widehat{\psi },\widehat{T},\widehat{S}\right\}\left(x\right)e^{ia\left(z-ct\right)},$$

(46)

eliminating the pressure term by cross differentiation and subtracting the resulting equations, we obtain (on dropping the wave symbols)

$$\begin{array}{c}\left[D^2-a^2-Da{(D^2-a^2)}^2\right]\widehat{\psi }+R{a}_{T}D\widehat{T}-R{a}_{S}D\widehat{S}\\ +c\{[-{\displaystyle \frac{ia}{P{r}_D}}(D^2-a^2)+ia{\lambda }_{2}Da{(D^2-a^2)}^2-ia{\lambda }_2(D^2-a^2)]\widehat{\psi }\\ -ia{\lambda }_{1}R{a}_{T}D\widehat{T}+ia{\lambda }_{1}R{a}_{S}D\widehat{S}\}+c^2\left[-{\displaystyle \frac{a^2{\lambda }_1}{P{r}_D}}(D^2-a^2)\widehat{\psi }\right]=0\end{array}$$

(47)

$${\displaystyle \frac{1}{2}}ia\widehat{\psi }+\left(ia{w}_b+a^2-D^2\right)\widehat{T}+c\left(-ia\widehat{T}\right)=0$$

(48)

$${\displaystyle \frac{1}{2}}ia\widehat{\psi }+\left(ia{w}_b+{\displaystyle \frac{a^2}{Le}}-{\displaystyle \frac{D^2}{Le}}\right)\widehat{S}+S_r\left(a^2-D^2\right)\widehat{T}+c\left(-ia\eta \widehat{S}\right)=0$$

(49)

The boundary conditions are

$$\widehat{\psi }=D\widehat{\psi }=\widehat{T}=\widehat{S}=0, \mathrm{at}x=\pm 1.$$

(50)

4. Numerical Procedure

In general, stability analysis of eigenvalue problems uses numerical methods. Equations (47)-(50) are solved via the Chebyshev collocation method. The following gives the Chebyshev polynomial of k-th order:

$${\xi }_k\left(x\right)=\cos k\theta ,\theta ={\cos }^{-1}x.$$

(51)

The following gives the Chebyshev collocation points:

$$x_j=\cos \left({\displaystyle \frac{j\pi }{N}}\right),j=0,1,\cdots N,$$

(52)

where j=0 and j=N represent the left and right wall limits, respectively, and N can be any positive number. The field variables $\widehat{\psi }$, $\widehat{T}$, and$\widehat{S}$ are approximated using Chebyshev polynomials:

$$\widehat{\psi }\left(x\right)={\displaystyle \sum _{j=0}^N{\psi }_j{\xi }_j\left(x\right)},\widehat{T}\left(x\right)={\displaystyle \sum _{j=0}^N{T}_j{\xi }_j\left(x\right)},\widehat{S}\left(x\right)={\displaystyle \sum _{j=0}^N{S}_j{\xi }_j\left(x\right)},$$

(53)

The values of ψ_j, T_j, and S_j are constants. Discretizing equations (47)-(50) in terms of Chebyshev polynomials yields:

$$\begin{array}{c}\left[{\displaystyle \sum _{k=0}^N{B}_{jk}}{\psi }_k-a^2{\psi }_j-Da({\displaystyle \sum _{k=0}^N{E}_{jk}}{\psi }_k-2a^2{\displaystyle \sum _{k=0}^N{B}_{jk}}{\psi }_k+a^4{\psi }_j)\right]+R{a}_T{\displaystyle \sum _{k=0}^N{A}_{jk}}{T}_j\\ -R{a}_S{\displaystyle \sum _{k=0}^N{A}_{jk}}{S}_j+c[-{\displaystyle \frac{ia}{{\mathrm{Pr}}_D}}({\displaystyle \sum _{k=0}^N{B}_{jk}}{\psi }_k-a^2{\psi }_j)+ia{\lambda }_{2}Da({\displaystyle \sum _{k=0}^N{E}_{jk}}{\psi }_k-2a^2{\displaystyle \sum _{k=0}^N{B}_{jk}}{\psi }_k+a^4{\psi }_j)\\ -ia{\lambda }_2({\displaystyle \sum _{k=0}^N{B}_{jk}}{\psi }_k-a^2{\psi }_j)]-c\left(ia{\lambda }_{1}R{a}_T{\displaystyle \sum _{k=0}^N{A}_{jk}}{T}_j+ia{\lambda }_{1}R{a}_S{\displaystyle \sum _{k=0}^N{A}_{jk}}{S}_j\right)\\ +c^2\left[-{\displaystyle \frac{a^2{\lambda }_1}{{\mathrm{Pr}}_D}}({\displaystyle \sum _{k=0}^N{B}_{jk}}{\psi }_k-a^2{\psi }_j)\right]=0,j=1,\cdots ,N-1,\end{array}$$

(54)

$${\displaystyle \frac{1}{2}}ia{\psi }_j+\left((ia{w}_b+a^2)T_j-{\displaystyle \sum _{k=0}^N{B}_{jk}}{T}_k\right)+c\left(-ia{T}_j\right)=0,j=1,\cdots ,N-1,$$

(55)

$$\begin{array}{c}{\displaystyle \frac{1}{2}}ia{\psi }_j+\left((ia{w}_b+{\displaystyle \frac{a^2}{Le}})S_j-{\displaystyle \frac{1}{Le}}{\displaystyle \sum _{k=0}^N{B}_{jk}}{S}_k\right)+S_r\left(a^2{T}_j-{\displaystyle \sum _{k=0}^N{B}_{jk}}{T}_k\right)\\ +c\left(-ia\eta S_j\right)=0,j=1,\cdots ,N-1,\end{array}$$

(56)

$${\psi }_0={\psi }_N=0$$

(57)

$${\displaystyle \sum _{k=0}^N{A}_{jk}}{\psi }_k=0j=0$$

(58)

$$T_0=T_N=0$$

(59)

$$S_0=S_N=0$$

(60)

Where

$$A_{jk}=\left\{\begin{array}{l}{\displaystyle \frac{c_j{\left(-1\right)}^{k+j}}{c_k\left(x_j-x_k\right)}}\hspace{0.33em}\hspace{1em}j\ne k,\\ {\displaystyle \frac{x_j}{2\left(1-x_j^2\right)}}\hspace{1em}1\le j=k\le N-1,\\ {\displaystyle \frac{2N^2+1}{6}}\hspace{1em}j=k=0,\\ -{\displaystyle \frac{2N^2+1}{6}}\hspace{1em}j=k=N,\end{array}\right.\hspace{0.33em}$$

(61)

$$B_{jk}=A_{jm}\cdot A_{mk}\hspace{0.33em}\mathrm{and} E_{jk}=B_{jm}\cdot B_{mk},$$

(62)

with

$$c_j=\left\{\begin{array}{l}1\hspace{1em}1\le j\le N-1,\\ 2\hspace{1em}j=0,N.\end{array}\right.$$

(63)

A generalized eigenvalue problem can be obtained from the discretization equation:

$$A_{0}X+c{A}_{1}X+c^2{A}_{2}X=0$$

(1.1)

where A₀, A₁, and A₂ are square matrices, c and X represent the complex eigenvalue and eigenfunction, respectively. The linear stability boundary of the fundamental flow is typically described by the eigenvalue issue. We designate the mode with the biggest imaginary portion c of the eigenvalue as the most growing mode among the eigenvalue and eigenfunction. This mode also happens to be the least decaying. The eigenvalues and eigenfunctions of the generalized eigenvalue problem are found using a QZ method. This algorithm is accessible as a built-in function called polyeig in the MatLab software package.

5. Results and Discussion

With the use of the Chebyshev collocation method, the linear stability of the Soret impact on DDC in a vertical layer of Darcy-Brinkman porous media is numerically explored. The Darcy-Rayleigh number Ra_T, the solute Darcy-Rayleigh number Ra_S, the Lewis number Le, the relaxation parameter λ₁, the retardation parameter λ₂, the Soret parameter S_r, the Darcy-Prandtl number Pr_D, the Darcy number Da, and the normalized porosity of the porous medium η are some of the significant non-dimensional parameters that control the flow.

Firstly, the parameter ranges are discussed based on the works by Shanker et al. [7] and Swamy et al. [24] The value of λ₁ must be greater than that of λ₂ (Bird et al. [32]; Hirata et al. [33]). The range of values for the Soret parameter is −1 < Sr < 1 (Bouachir et al. [29]). In fact, the Darcy-Prandtl number Pr_D can be expressed by Prandtl number, Darcy number, porosity, and specific heat ratio, i.e., Pr_D=γεPr/Da. Pr_D depends on the properties of the fluid and the porous matrix. For a sparse porous media, Da ranges from 0.001 to 1, γ changes from 0 to 1, ε is around 0.5, and the typical Prandtl number for viscoelastic fluids is 10. Therefore, Pr_D varies from 5 to 5000.

Table 1 provides a convergence analysis for numerical calculations by altering the order of the Chebyshev polynomials. Employing 30 to 45 configuration points in Eq. (50), it is found that the value of ac_i is independent of the polynomial order N between 35 and 45 and achieves 4-digit point precision. Therefore, all numerical results given subsequently are obtained by taking N = 40 in the Chebyshev expansion.

Figure 3(a)–3(d) show the growth rate for the Soret parameter as well as the neutral stability curves. In the neutral stability plot, it is stable outside of the tongue-shaped region and unstable within it. The minimum value of the neutral stability curve indicates the critical condition for the flow changing from a stable to an unstable state. In this paper, the corner symbol c represents the critical value, and the value of Ra_Tc is the critical value of the flow between stable and unstable; when Ra_T> Ra_Tc, it means that the flow is unstable, and when Ra_T < Ra_Tc, it means that the flow is stable. Figure 3(a) depicts the effect of different Soret parameters S_r on the neutral stability curve when Le = 1. When S_r increases, the minimum value of Ra_Tc moves towards the large value of wave number a, indicating that S_r reduces the width of the cell. In addition, S_r =0 denotes the neutral stability curve without the Soret effect. For S_r <0, we find that the minimum value of Ra_Tc decreases with the increase of the negative Soret parameter, which indicates that it weakens the stability of the system. For S_r >0, the minimum value of Ra_Tc increases with the increase of the positive Soret parameter, which indicates that it enhances flow instability. It shows that positive S_r results in a more stable flow, and negative S_r has the opposite effect. Figure 3(b) shows the growth rate curve for Ra_T = 1750, i.e., the growth rate case at the dotted line in Figure 3(a). In Figure 3(b), the growth rate is plotted using the same parameters. The larger the positive S_r is, the smaller the growth rate is, indicating that the positive S_r increases the stability of the system; the more significant the negative S_r is, the larger the growth rate is, indicating that the negative S_r weakens the stability of the system, and the same conclusion as Figure 3(a) is obtained. Figure 3(c) shows the neutral stability curve when Le=2.

We find that for S_r<0, the minimum value of Ra_Tc increases with the increase of the negative Soret parameter, which indicates that it improves stability. For S_r >0, the minimum value of Ra_Tc decreases with the rise of the positive Soret parameter, which suggests that it reduces stability. It shows that the positive S_r parameter has an unstable effect, and the negative S_r parameter has a stabilizing effect. It is opposite to the conclusion for Le =1. Figure 3(d) shows the growth rate curve when Ra_T=1750 and Le=2. The larger the positive S_r is, the bigger the growth rate is; the larger the negative S_r is, the smaller the growth rate is, indicating that the negative S_r enhances the stability, and the same conclusion as that of Figure 3(c) is obtained. In summary, we find that the positive and negative values of S_r have different effects on the instability, and the size of Le affects the impact of S_r on the instability of the flow.

Figure 4(a)-(d) illustrate the neutral stability curves and growth rate curves for different Lewis numbers Le. Figure 4(a) depicts that the minimum value of Ra_Tc decreases with increasing Le for Le <0.5, indicating that Le promotes flow instability at this point. Figure 4(b) shows the growth rate curve when Ra_T = 1750 under the parameters of Figure 4(a). We find that the larger Le is, the larger the growth rate is and the more unstable the flow is. It indicates that Le plays a weakening effect on the stability of the flow when Le <0.5. Figure 4(c) suggests that Le promotes the stability of the flow when Le > 1. Under this parameter, Figure 4(d) depicts the growth rate when Ra_T = 1750, i.e., the growth rate corresponding to the neutral stabilization point intersecting the dashed line in Figure 4(c), and the growth rate is found to decrease with increasing Le. In summary, we find that Le has a dual effect on the instability of the flow.

Table 2 shows the values of Le_c1 for different S_r. According to the table, Le_c1 decreases with increasing S_r and remains around 0.7. Therefore, Le enhances the flow instability when Le < 0.7 and inhibits the flow instability when Le > 0.7.

Figure 5 illustrates the neutral stability curves of Ra_Tc as Le varies for different Soret parameters S_r. As Le increases from zero, the value of Ra_Tc first decreases and then increases. This suggests that there is a critical value Le_c1. This indicates that Le has both promoting and inhibiting effects on the instability of the flow. The minimum value of Ra_Tc corresponds to a critical Le value Le_c1, which promotes the instability of the flow when Le < Le_c1 and inhibits the instability of the flow when Le > Le_c1. According to the graph, it is found that there is also a critical value of Le, Le_c2=1.5735, when Le < Le_c2, S_r contributes to the instability of the flow, and when Le > Le_c2, S_r inhibits the instability of the flow. Thus, we find that Le and S_r have a dual role in the stability of the flow and that the magnitude of Le influences the role of S_r.

Figure 6(a)-(d) illustrate the neutral stability curves and growth rate curves for different normalized porosity η. Figure 6(a) shows that η promotes flow instability at this point. Figure 6(b) shows the growth rate profile when Ra_T = 1130 under the parameters of Figure 6(a). We find that the larger η is, the larger the growth rate is and the more unstable the flow is. It indicates that η plays a weakening effect on the stability of the flow when 0.2< η < 0.6. Figure 6(c) depicts that the minimum value of Ra_Tc increases with η when 0.8< η < 0.88, indicating that η promotes the stability of the flow. In Figure 6(d), the growth rate decreases with increasing η. In summary, η acts differently on the instability of the flow at different values.

Figure 7(a) demonstrates the neutral stability curve of Ra_Tc as η varies for different relaxation parameters λ₂. As η increases from zero, the value of Ra_Tc decreases and then increases. This suggests that η has a dual effect on the instability. The minimum value of Ra_Tc corresponds to a critical value η_c, which enhances the instability when η < η_c and inhibits the instability when η > η_c. In addition, the minimum value of Ra_Tc increases with λ₂, indicating that λ₂ promotes flow stability. Figure 7(b) shows the growth rate when Ra_T = 1400 for the parameters of Figure 7(a). It shows that when λ₂ increases, the growth rate is smaller, and the flow is more stable.

According to Table 3, η_c increases with increasing λ₂ and stays around 0.8. Thus, we conclude that η enhances the instability of the flow when η<0.8 and suppresses the instability of the flow when η>0.8.

The neutral stability curves for various parameters λ₁ when Pr_D = 50 are shown in Figure 8(a). When λ₁ is larger, the minimum value of Ra_Tc is smaller, showing that the flow is more stable. In Figure 8(b), for Pr_D =200, it is found that the minimum value of Ra_Tc increases for larger λ₁, indicating a more unstable flow. When Pr_D =600, λ₁ promotes the stability of the flow in Figure 8(c). In order to explore the effect of Pr_D on the role of λ₁, the neutral stability curve of (Pr_D, Ra_Tc) is plotted, i.e., Figure 8(d). As Pr_D increases, the value of Ra_Tc decreases, then increases and decreases continuously. This indicates that Pr_D has a dual effect on the stability of the flow. According to the images, we can find the critical values of PrD, PrDc1 and PrDc2. When PrD < PrDc1 or PrD > PrDc2, PrD has an inhibiting effect on the instability, and when PrDc1 < PrD < PrDc2, PrD has a facilitating effect on the instability. According to Figure 8(d), when λ1 = 0.65, the specific critical values PrDc1 is 50.2 and PrDc2 is 398.8. Moreover, the magnitude of PrD affects the stability of λ1. When PrD < 110.8 or PrD > PrDc2, the larger λ1 is, the smaller RaTc is, indicating that λ1 inhibits the stability. Interestingly, when 110.8 < PrD < PrDc2, λ1 suppresses the flow instability. This suggests that λ1 also has a dual effect on the flow stability and is influenced by the size of PrD.

6. Conclusions

The study examines the instability of the Soret Effect on the DDC of an Oldroyd-B fluid in a vertical porous layer. Based on the Darcy-Brinkman-Oldroyd model, an Orr-Sommerfeld eigenvalue problem is obtained using perturbation theory and the Oberbeck-Boussinesq approximation, which is numerically simulated employing the Chebyshev collocation method. Then, the neutral stability curves, as well as growth rates, are obtained. As a result, Pr_D has a dual effect on instability, and it has two critical values, Pr_Dc1 and Pr_Dc2, which depend on the other parameters. When Pr_Dc1 < Pr_D < Pr_Dc2, Pr_D boosts flow stability; when Pr_D < Pr_Dc1 or Pr_D > Pr_Dc2, Pr_D enhances the instability. The impact of λ₁ on the instability depends on the value of Pr_D. When 110.8< Pr_D < 398.8, λ₁ suppresses the flow instability. However, when taking other values of Pr_D, λ₁ promotes flow instability. In addition, it is obtained that Le has a dual effect on stability, and a critical value exists for Le. Le enhances the flow instability when Le < 0.7 and inhibits the flow instability when Le > 0.7. At the same time, S_r has a dual role in the stability of the flow. According to the positive and negative values of S_r, S_r has different effects on the instability of the flow. Positive S_r have stable effects and negative S_r parameters have unstable effects. The magnitude of Le influences the role of S_r. When Le < Le_c2, S_r boosts the instability, and when Le> Le_c2, S_r plays an inhibiting role. Furthermore, there is a critical value for η. When η < 0.8, η enhances the instability; when η >0.8, it makes the flow more stable. Finally, the relaxation parameter λ₂ promotes flow stability. In conclusion, Pr_D, λ₁, Le, S_r, and η have a dual impact on the instability.