各向异性的双孔隙多孔介质中Oldroyd-B流体双扩散对流不稳定性

The Double-Diffusive Convection Instability in an Anisotropic Bidispersive Layer of an Oldroyd-B Fluid

摘要: 本文研究了饱和Oldroyd-B流体的各向异性双孔隙多孔层内的双扩散对流不稳定性。建立了相关控制方程和边界条件的理论模型。通过小扰动法和常规展开法推导了Orr-Sommerfeld方程。此外，通过理论分析和数值模拟，得到了中性稳定曲线。结果表明，体系稳定性受到Lewis数Le、浮力比N、大孔与微孔渗透率比K_r、大孔隙各向异性参数K_fr、微孔隙各向异性参数K_pr以及弛豫时间λ₁和滞后时间λ₂的综合影响。各参数在临界Darcy-Rayleigh数与中性稳定曲线的变化中均起关键作用。具体而言，增大λ₁提高流体的弹性储能能力，从而削弱扰动增长、延迟失稳发生；增大λ₂使应力响应具有更强的时间滞后性，进一步提高系统稳定阈值。随着浮力比N的增加，溶质浮力削弱热浮力的驱动作用，导致体系更难发生对流。大孔与微孔隙各向异性参数K_fr、K_pr的增大均会使临界Darcy-Rayleigh数上升，表明孔隙各向异性对系统具有稳定化效应。总体而言，Oldroyd-B流体的粘弹性效应与孔隙结构各向异性共同主导了双扩散对流的不稳定性行为，增强任一因素均能有效推迟失稳的发生，为非牛顿流体在复杂多孔介质中的传热传质调控提供了理论依据。

Abstract: This paper investigates the double-diffusive convection instability in an anisotropic bidispersive layer of an Oldroyd-B fluid. A theoretical model with the relevant control equations and boundary conditions was established. The Orr-Sommerfeld equation was derived using the method of small perturbations and regular expansion. Further, through theoretical analysis and numerical simulation, the neutral stability curve was obtained. The results show that the system stability is jointly influenced by the Lewis number Le, buoyancy ratio N, permeability ratio between macropores and micropores K_r, anisotropy parameters of the macroporous and microporous structures K_fr and K_pr, as well as the relaxation time λ₁ and retardation time λ₂. All these parameters play crucial roles in determining the critical Darcy-Rayleigh number and the shape of the neutral stability curves. Specifically, increasing λ₁ enhances the elastic energy storage capacity of the fluid, suppressing disturbance growth and delaying the onset of convection. A larger λ₂ introduces stronger temporal lag in the stress response, further raising the stability threshold. With the increase of N, the solutal buoyancy counteracts the thermal buoyancy, making the system more resistant to convection. In addition, higher anisotropy parameters K_fr and K_pr of the macro- and micropores lead to an increase in the critical Darcy-Rayleigh number, indicating that pore anisotropy exerts a stabilizing influence on the system. Overall, the viscoelastic effects of the Oldroyd-B fluid and the anisotropic pore structure jointly govern the double-diffusive convection instability, and the enhancement of either factor effectively delays the onset of instability, providing theoretical insights for the control of heat and mass transfer of non-Newtonian fluids in complex porous media.

1. 引言

双扩散对流是热量与溶质(或质量组分)在流体中以不同扩散率共同作用而引起的一类重要流动现象[1]，广泛存在于CO₂地质封存[2]、含盐地热田[3]、污染物治理[4]等领域。当热扩散率与质量扩散率显著不同步时，系统内部会产生复杂的对流模式和不稳定性结构，从而影响物质与能量的输运效率。近年来，随着多孔介质理论的发展，双孔隙多孔介质在化学工程和物理学领域得到广泛研究。大孔隙与微孔隙之间的耦合传递效应对流体流动与传热过程具有显著影响，这使得双扩散对流问题在多孔体系中具有更为复杂的动力学特征。

在许多工程和自然系统中，流体并非理想的牛顿型流体，而呈现出粘弹性或非线性流变行为。Oldroyd-B模型作为典型的粘弹性流体本构关系，能够同时描述流体的粘性与弹性响应，这种效应主要在化学和石油工业[5] [6]以及生物流体力学[7]等领域被观察到。基于适用于粘弹性流体的达西定律的修正形式，Tan和Masuoka [8]将Stokes第一问题扩展到了Oldroyd-B流体在多孔半空间中的情况，从而为后续关于粘弹性流体在多孔结构中的流动研究奠定了基础。在他们的开创性贡献之后，Prema等人[9]对饱和有Oldroyd-B流体的多孔介质中的对流换热进行了全面综述，总结了主要的理论进展，并指出了该领域的开放性挑战。最近，Ren等人[10]研究了垂直多孔层中Oldroyd-B流体的双扩散对流的不稳定性，从而扩展了对粘弹性多孔系统中耦合溶质和热效应的理解。

相较于牛顿流体，Oldroyd-B流体在受热和浓度梯度驱动下的对流不稳定性更为复杂，其应力弛豫和弹性效应会显著改变临界Rayleigh数及流动模式。研究各向异性双孔隙介质中Oldroyd-B流体的双扩散对流特性，对于理解非牛顿流体在复杂介质中的输运规律具有重要的理论与应用价值。

自然与工程多孔介质通常具有明显的各向异性特征，例如渗透率、热导率及扩散系数在不同方向上存在差异。各向异性不仅改变了热量和溶质的传输路径，也对流体不稳定性的发展产生显著影响。Straughan [11]研究了具有水平各向异性的多孔介质中双扩散对流不稳定性。最近，Jia和Jian [12]将研究范围扩展到了Oldroyd-B流体，揭示了粘弹性如何改变具有大宏观孔隙的双分散多孔介质中的热对流不稳定性。

同时，双孔隙结构提供了两种尺度的流动通道：大孔隙主导整体渗流，而微孔隙则通过相互交换影响局部流场。这种多尺度耦合使得双扩散对流问题更加复杂，必须在理论上综合考虑流体粘弹性、双孔隙交换机制以及各向异性效应的共同作用。已有研究主要集中在牛顿流体在单孔隙或各向同性多孔介质中的双扩散对流不稳定性，对非牛顿流体尤其是Oldroyd-B型流体的系统研究相对较少。最近，Jia和Jian [13]研究了各向异性的双孔隙多孔介质内Oldroyd-B流体对热对流不稳定性的影响，但他们忽略了溶质扩散的影响，从而难以全面揭示复杂多孔系统中热质耦合传输的内在规律。

基于上述背景，本文旨在建立一个描述各向异性双孔隙多孔介质中Oldroyd-B流体双扩散对流不稳定性的数学模型。建立相关控制方程和边界条件，通过线性稳定性分析、数值模拟，讨论流体弹性参数、孔隙耦合系数、各向异性参数以及热、溶质相关参数对系统临界不稳定性的影响。研究结果有助于加深对粘弹性流体在复杂多孔介质中热–质输运机理的理解，为非牛顿流体的多场耦合调控与工程应用提供理论依据。

2. 数学模型

我们考虑一个厚度为d的无限延伸水平多孔层，该多孔层由大孔隙和微孔隙组成，并充满不可压缩的Oldroyd-B流体，如图1所示。建立笛卡尔坐标系(x^*, y^*, z^*)，其原点位于多孔层的下边界处，其中x^* – y^*平面为水平面，z^*轴垂直向上。上下边界z^* = 0和z^* = d假设为等温边界，分别维持在温度$T_1^{\ast }$ 和$T_2^{\ast }$ ，以及溶质浓度$C_1^{\ast }$ 和$C_2^{\ast }$ 。重力加速度g的方向与z^*轴相反。依据已有研究文献[8] [14]中的模型建模方法，采用适用于双分散多孔介质中Oldroyd-B流体的控制方程。

Figure 1. Model diagram of the horizontal bidispersive porous layer

图1. 水平双向分散多孔层的模型图

${\nabla }^{\ast }\cdot u_f^{\ast }=0$ , (1)

${\nabla }^{\ast }\cdot u_p^{\ast }=0$ , (2)

${\nabla }^{\ast }{p}_f^{\ast }+\zeta \left(u_f^{\ast }-u_p^{\ast }\right)+{\rho }^{\ast }g=r_f$ , (3)

${\nabla }^{\ast }{p}_p^{\ast }+\zeta \left(u_p^{\ast }-u_f^{\ast }\right)+{\rho }^{\ast }g=r_p$ . (4)

在这里，速度矢量$u_f^{\ast }=\left(u_f^{\ast },v_f^{\ast },w_f^{\ast }\right)$ 和$u_p^{\ast }=\left(u_p^{\ast },v_p^{\ast },w_p^{\ast }\right)$ 分别表示大孔隙和微孔隙中的渗流速度。ζ表示考虑到大孔隙和微孔隙之间动量交换的相互作用系数。大孔隙和微孔隙内的压力分别用$p_f^{\ast }$ 和$p_p^{\ast }$ 表示。为了简化流动模型，采用Oberbeck-Boussinesq近似[15]。在该近似下，假设流体密度${\rho }^{\ast }$ 随温度和溶质浓度线性变化。

${\rho }^{\ast }={\rho }_0\left[1-{\alpha }_T\left(T^{\ast }-T_2^{\ast }\right)+{\alpha }_C\left(C^{\ast }-C_2^{\ast }\right)\right]$ . (5)

其中，参考温度、溶质浓度和密度分别用$T_2^{\ast }$ 、$C_2^{\ast }$ 和ρ₀表示。变量T^*和C^*分别表示局部温度和溶质浓度。系数α_T和α_C分别表示与温度和溶质浓度相关的体积膨胀系数。固体基质的流动阻力通过r_f和r_p来表征，如文献[8] [14]中的方程(1)和(2)所定义。

$\left(1+{\lambda }_1^{\ast }\frac{\partial }{\partial t^{\ast }}\right)r_f=-K_f^{-1}\mu \left(1+{\lambda }_2^{\ast }\frac{\partial }{\partial t^{\ast }}\right)u_f^{\ast }$ , (6)

$\left(1+{\lambda }_1^{\ast }\frac{\partial }{\partial t^{\ast }}\right)r_p=-{{\rm K}}_p^{-1}\mu \left(1+{\lambda }_2^{\ast }\frac{\partial }{\partial t^{\ast }}\right)u_p^{\ast }$ . (7)

在该模型中，${{\rm K}}_f^1={{\rm K}}_{fx}^{-1}i+{{\rm K}}_{fy}^{-1}j+{{\rm K}}_{fz}^{-1}k$ 和${{\rm K}}_p^1={{\rm K}}_{px}^{-1}i+{{\rm K}}_{py}^{-1}j+{{\rm K}}_{pz}^{-1}k$ 分别表示与大孔隙和微孔隙的渗透率张量。我们假设水平方向的渗透率是相同的，则有${{\rm K}}_{fx}={{\rm K}}_{fy}$ 和${{\rm K}}_{px}={{\rm K}}_{py}$ 。μ表示流体粘性。参数λ₁和λ₂分别代表Oldroyd-B流体的松弛时间和滞后时间。能量平衡控制方程(即温度场方程)基于Nield和Kuznetsov [16]的理论推导，形式如下：

${\left(\rho c\right)}_m\frac{\partial T^{\ast }}{\partial t^{\ast }}+{\left(\rho c\right)}_f\left(u_f^{\ast }+u_p^{\ast }\right)\cdot {\nabla }^{\ast }{T}^{\ast }={\kappa }_m{\nabla }^{\ast 2}{T}^{\ast }$ . (8)

其中，$u_f^{\ast }=\phi V_f$ 和$u_p^{\ast }=\gamma \left(1-\phi \right)V_p$ 分别表示大孔隙和微孔隙内的平均渗流速度，并且

${\left(\rho c\right)}_m=\left(1-\gamma \right)\left(1-\phi \right){\left(\rho c\right)}_s+\phi {\left(\rho c\right)}_f+\gamma \left(1-\phi \right){\left(\rho c\right)}_p$ , (9)

${\kappa }_m=\left(1-\gamma \right)\left(1-\phi \right){\kappa }_s+\phi {\kappa }_f+\gamma \left(1-\phi \right){\kappa }_p$ . (10)

其中，下标s、f和p分别表示固体基质、大孔隙和微孔隙。符号κ_m表示热扩散率。V_f和V_p分别表示大孔隙区和微孔隙区中的流体速度。ρc表示材料密度与定压比热的乘积。φ表示大孔隙所占总体积的比例，γ表示微孔隙在多孔介质中所占总体积的比例。由于大孔隙和微孔隙均被相同的流体饱和，因此假设 (ρc)_f 与(ρc)_p相同[17]。令$M={\left(\rho c\right)}_f/{\left(\rho c\right)}_m$ 和${\kappa }_T={\kappa }_m/{\left(\rho c\right)}_f$ ，则方程(8)可改写为：

$\frac{1}{M}\frac{\partial T^{\ast }}{\partial t^{\ast }}+\left(u_f^{\ast }+u_p^{\ast }\right)\cdot {\nabla }^{\ast }{T}^{\ast }={\kappa }_T{\nabla }^{\ast 2}{T}^{\ast }$ . (11)

此外，溶质浓度场的控制方程由以下给出

$\epsilon \frac{\partial C^{\ast }}{\partial t^{\ast }}+\left(u_f^{\ast }+u_p^{\ast }\right)\cdot {\nabla }^{\ast }{C}^{\ast }={\kappa }_C{\nabla }^{\ast 2}{C}^{\ast }$ . (12)

其中

$\epsilon =\phi +\gamma \left(1-\phi \right)$ , (13)

${\kappa }_c=\phi {\kappa }_{cf}+\gamma \left(1-\phi \right){\kappa }_{cp}$ . (14)

其中，κ_cf和κ_cp分别表示大孔隙区和微孔隙区流体的溶质扩散率。在大孔隙区和微孔隙区边界处施加无渗透边界条件，同时保持温度和溶质浓度恒定：

$T^{\ast }=T_1^{\ast }$ , $C^{\ast }=C_1^{\ast }$ , $w_f^{\ast }=0$ , $w_p^{\ast }=0$ , 在$z^{\ast }=0$ 处 (15)

并且

$T^{\ast }=T_2^{\ast }$ , $C^{\ast }=C_2^{\ast }$ , $w_f^{\ast }=0$ , $w_p^{\ast }=0$ , 在$z^{\ast }=d$ 处 (16)

以下这些特征尺度被引入：

$\left(x,y,z\right)=\frac{\left(x^{*},y^{*},z^{*}\right)}{d}$ ，$\left(u_f,u_p\right)=\frac{\left(u_f^{*},u_p^{*}\right)}{{\kappa }_T/d}$ ，$t=\frac{t^{\ast }}{d^2/{\kappa }_T}$ ，$\left({\lambda }_1,{\lambda }_2\right)=\frac{\left({\lambda }_1^{*},{\lambda }_2^{*}\right)}{{\kappa }_T/d^2}$ ，$T=\frac{T^{*}-T_2^{*}}{T_1^{*}-T_2^{*}}$ ，

$C=\frac{C^{*}-C_2^{*}}{C_1^{*}-C_2^{*}}$ ，$Le=\frac{{\kappa }_T}{{\kappa }_C}$ ，${\sigma }_f=\frac{\varsigma K_{fx}}{\mu }$ ，$N=\frac{{\alpha }_C\left(C_1^{*}-C_2^{*}\right)}{{\alpha }_T\left(T_1^{*}-T_2^{*}\right)}$ ，$Da=\frac{K_{fx}}{d^2}$ ，$\left(p_f,p_p\right)=\frac{\left(p_f^{*},p_p^{*}\right)}{\mu {\kappa }_T/K_{fx}}$ ，

$Ra=\frac{g{\alpha }_T\left(T_1^{*}-T_2^{*}\right)d^3}{\nu {\kappa }_T}$ ，$R{a}_D=RaDa$ ，$K_r=\frac{K_f}{K_p}$ ，$K_{pr}=\frac{K_{px}}{K_{pz}}$ ，$K_{fr}=\frac{K_{fx}}{K_{fz}}$ . (17)

在此，Ra_D、N、Le、K_r、K_fr、K_pr分别代表Darcy-Rayleigh数、浮力比以及Lewis数、大孔与微孔渗透率之比、大孔隙的各向异性参数、微孔隙的各向异性参数。

对方程和边界条件进行无量纲化处理，得到无量纲控制方程：

$\nabla \cdot u_f=0$ , (18)

$\nabla \cdot u_p=0$ , (19)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial p_f}{\partial x}+{\sigma }_f\left(u_f-u_p\right)\right]+\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)u_f=0$ , (20)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial p_f}{\partial y}+{\sigma }_f\left(v_f-v_p\right)\right]+\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)v_f=0$ , (21)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial p_f}{\partial z}+{\sigma }_f\left(w_f-w_p\right)-R{a}_{D}T+R{a}_{D}NC\right]+K_{fr}\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)w_f=0$ , (22)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial p_p}{\partial x}+{\sigma }_f\left(u_p-u_f\right)\right]+K_r\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)u_p=0$ , (23)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial p_p}{\partial x}+{\sigma }_f\left(v_p-v_f\right)\right]+K_r\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)v_p=0$ , (24)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial p_p}{\partial z}+{\sigma }_f\left(w_f-w_p\right)-R{a}_{D}T+R{a}_{D}NC\right]+K_r{K}_{pr}\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)w_p=0$ , (25)

$\frac{1}{M}\frac{\partial T}{\partial t}+\left(u_f+u_p\right)\cdot \nabla T=\Delta T$ , (26)

$\frac{1}{M}\frac{\partial C}{\partial t}+\left(u_f+u_p\right)\cdot \nabla C=\Delta C$ . (27)

无量纲边界条件：

$T=1$ , $C=1$ , $w_f=0$ , $w_p=0$ , 在$z=0$ 处 (28)

并且

$T=0$ , $C=0$ , $w_f=0$ , $w_p=0$ , 在$z=1$ 处. (29)

本研究在建立各向异性的双孔隙多孔介质中Oldroyd-B流体的双扩散对流模型时，基于若干常见的物理近似与合理假设。首先，采用了Oberbeck-Boussinesq近似，即假设流体密度的变化仅在浮力项中保留，而在其余方程中视为常数。这种处理使方程体系显著简化，适用于温度与浓度差异较小、对流强度较弱的情形。其次，流体的渗流行为遵循针对粘弹性流体的修正Darcy定律，该假设能够较好地反映Oldroyd-B流体在多孔介质中的弹性响应，但当孔隙结构复杂或流速较高时，仍需进一步考虑Brinkman或Forchheimer修正以提高模型的适用性。此外，模型中假定流体的粘度、扩散率及热导率等物性参数不随温度或浓度变化，从而突出了粘弹性特征与孔隙结构对系统稳定性的主导作用。然而，这一常物性假设也在一定程度上忽略了热质耦合引起的物性变化对不稳定性的次级影响。总体而言，上述假设确保了模型具有良好的理论一致性和解析可行性，其结果可为研究更复杂情况下的非牛顿流体多场耦合问题提供基础参考。

假设在基本流动(未受干扰)状态下，速度、温度和浓度仅取决于垂直坐标z，即：$u_{fb}=0$ , $u_{pb}=0$ , $T_b=T_b\left(z\right)$ 和$C_b=C_b\left(z\right)$ 。我们通过求解得到：

$T_b\left(z\right)=1-z$ , $C_b\left(z\right)=1-z$ . (30)

基本流动受到微小扰动，其形式可表示为：

$\left(u_f,v_f,w_f,u_p,v_p,w_p\right)=\delta \left(u_f^{*},v_f^{*},w_f^{*},u_p^{*},v_p^{*},w_p^{*}\right)$ ，$T=T_b\left(z\right)+\delta T^{\prime }$ ，$C=C_b\left(z\right)+\delta C^{\prime }$ ，

$p_f=p_{fb}\left(z\right)+\delta {p^{\prime }}_f$ ，$p_p=p_{pb}\left(z\right)+\delta {p^{\prime }}_p$ . (31)

这里，δ是一个很小的参数。通过将方程(31)代入方程(18)至(29)，并进行线性化处理，得到了以下线性化的控制方程和边界条件：

$\nabla \cdot {u^{\prime }}_f=0$ , (32)

$\nabla \cdot {u^{\prime }}_p=0$ , (33)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial {p^{\prime }}_f}{\partial x}+{\sigma }_f\left({u^{\prime }}_f-{u^{\prime }}_p\right)\right]+\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){u^{\prime }}_f=0$ , (34)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial {p^{\prime }}_f}{\partial y}+{\sigma }_f\left({v^{\prime }}_f-{v^{\prime }}_p\right)\right]+\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){v^{\prime }}_f=0$ , (35)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial {p^{\prime }}_f}{\partial z}+{\sigma }_f\left({w^{\prime }}_f-{w^{\prime }}_p\right)-R{a}_{D}T+R{a}_{D}NC\right]+K_{fr}\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){w^{\prime }}_f=0$ , (36)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial {p^{\prime }}_p}{\partial x}+{\sigma }_f\left({u^{\prime }}_p-{u^{\prime }}_f\right)\right]+K_r\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){u^{\prime }}_p=0$ , (37)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial {p^{\prime }}_p}{\partial x}+{\sigma }_f\left({v^{\prime }}_p-{v^{\prime }}_f\right)\right]+K_r\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){v^{\prime }}_p=0$ , (38)

$\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[\frac{\partial {p^{\prime }}_p}{\partial z}+{\sigma }_f\left({w^{\prime }}_f-{w^{\prime }}_p\right)-R{a}_D{T}^{\prime }+R{a}_{D}N{C}^{\prime }\right]+K_r{K}_{pr}\left(1+{\lambda }_2\frac{\partial }{\partial t}\right){w^{\prime }}_p=0$ , (39)

$\frac{1}{M}\frac{\partial T^{\prime }}{\partial t}+\left({w^{\prime }}_f+{w^{\prime }}_p\right)\frac{\partial T_b}{\partial z}=\Delta T^{\prime }$ , (40)

$\frac{1}{M}\frac{\partial C^{\prime }}{\partial t}+\left({w^{\prime }}_f+{w^{\prime }}_p\right)\frac{\partial C_b}{\partial z}=\Delta C^{\prime }$ . (41)

$T^{\prime }=0$ , $C^{\prime }=0$ , ${w^{\prime }}_f=0$ , ${w^{\prime }}_p=0$ , 在$z=0,1$ 处 (42)

根据Squire理论[18] [19]，将三维模型简化为二维问题，然后引入流函数：

${u^{\prime }}_f\left(x,z,t\right)=\frac{\partial {{\psi }^{\prime }}_f\left(x,z,t\right)}{\partial z}$ ，${w^{\prime }}_f\left(x,z,t\right)=-\frac{\partial {{\psi }^{\prime }}_f\left(x,z,t\right)}{\partial x}$ ，

${u^{\prime }}_p\left(x,z,t\right)=\frac{\partial {{\psi }^{\prime }}_p\left(x,z,t\right)}{\partial z}$ ，${w^{\prime }}_p\left(x,z,t\right)=-\frac{\partial {{\psi }^{\prime }}_p\left(x,z,t\right)}{\partial x}$ . (43)

并消除压力项，得到方程：

$\begin{array}{l}\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)\left(\frac{{\partial }^2{{\psi }^{\prime }}_f}{\partial z^2}+K_{fr}\frac{{\partial }^2{{\psi }^{\prime }}_f}{\partial x^2}\right)\\ +\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[{\sigma }_f\left(\frac{{\partial }^2{{\psi }^{\prime }}_f}{\partial z^2}-\frac{{\partial }^2{{\psi }^{\prime }}_p}{\partial z^2}-\frac{{\partial }^2{{\psi }^{\prime }}_p}{\partial x^2}+\frac{{\partial }^2{{\psi }^{\prime }}_f}{\partial x^2}\right)-R{a}_D\frac{\partial T^{\prime }}{\partial x}+R{a}_{D}N\frac{\partial C^{\prime }}{\partial x}\right]=0,\end{array}$ (44)

$\begin{array}{l}{K}_r\left(1+{\lambda }_2\frac{\partial }{\partial t}\right)\left(\frac{{\partial }^2{{\psi }^{\prime }}_p}{\partial z^2}+K_{pr}\frac{{\partial }^2{{\psi }^{\prime }}_p}{\partial x^2}\right)\\ +\left(1+{\lambda }_1\frac{\partial }{\partial t}\right)\left[{\sigma }_f\left(\frac{{\partial }^2{{\psi }^{\prime }}_p}{\partial z^2}-\frac{{\partial }^2{{\psi }^{\prime }}_f}{\partial z^2}-\frac{{\partial }^2{{\psi }^{\prime }}_f}{\partial x^2}+\frac{{\partial }^2{{\psi }^{\prime }}_p}{\partial x^2}\right)-R{a}_D\frac{\partial T^{\prime }}{\partial x}+R{a}_{D}N\frac{\partial C^{\prime }}{\partial x}\right]=0,\end{array}$ (45)

$\frac{1}{M}\frac{\partial T^{\prime }}{\partial t}+\left(-\frac{\partial {{\psi }^{\prime }}_f}{\partial x}-\frac{\partial {{\psi }^{\prime }}_p}{\partial x}\right)\frac{\partial T_b}{\partial z}=\frac{{\partial }^2{T}^{\prime }}{\partial x^2}+\frac{{\partial }^2{T}^{\prime }}{\partial z^2}$ , (46)

$\frac{1}{M}\frac{\partial C^{\prime }}{\partial t}+\left(-\frac{\partial {{\psi }^{\prime }}_f}{\partial x}-\frac{\partial {{\psi }^{\prime }}_p}{\partial x}\right)\frac{\partial C_b}{\partial z}=\frac{{\partial }^2{C}^{\prime }}{\partial x^2}+\frac{{\partial }^2{C}^{\prime }}{\partial z^2}$ . (47)

边界条件：

$T^{\prime }=0$ , $C^{\prime }=0$ , $\frac{\partial {{\psi }^{\prime }}_f}{\partial x}=0$ , $\frac{\partial {{\psi }^{\prime }}_p}{\partial x}=0$ , 在$z=0,1$ 处. (48)

引入以下正则模型来代替扰动流函数、温度和浓度：

$\left({{\psi }^{\prime }}_f,{{\psi }^{\prime }}_p,T^{\prime },C^{\prime }\right)=\left({\Psi }_f\left(z\right),{\Psi }_p\left(z\right),\Theta \left(z\right),\Omega \left(z\right)\right){\text{e}}^{iax+wt}$ , (49)

将其代入到方程(44)~(48)，得到常微分特征值问题：

$\left(1+{\lambda }_2\omega \right)\left(D^2-K_{fr}{a}^2\right){\Psi }_f+\left(1+{\lambda }_1\omega \right)\left[{\sigma }_f\left(D^2-a^2\right)\left({\Psi }_f-{\Psi }_p\right)-R{a}_{D}ia\Theta +R{a}_{D}Nia\Omega \right]=0$ , (50)

$K_r\left(1+{\lambda }_2\omega \right)\left(D^2-K_{pr}{a}^2\right){\Psi }_p+\left(1+{\lambda }_1\omega \right)\left[{\sigma }_f\left(D^2-a^2\right)\left({\Psi }_p-{\Psi }_f\right)-R{a}_{D}ia\Theta +R{a}_{D}Nia\Omega \right]=0$ , (51)

$\frac{\omega }{M}\Theta -ia\left({\Psi }_f+{\Psi }_p\right)-\left(D^2-a^2\right)\Theta =0$ , (52)

$\epsilon \omega \Omega -ia\left({\Psi }_f+{\Psi }_p\right)-\frac{1}{Le}\left(D^2-a^2\right)\Omega =0$ , (53)

$\frac{\partial {\Psi }_f}{\partial x}=0$ , $\frac{\partial {\Psi }_p}{\partial x}=0$ , $\Theta =0$ , $\Omega =0$ , at $z=0,1$ , (54)

其中D = d/dz。

3. 数值模拟

通过切比雪夫节点法求解了线性化扰动方程。第n个切比雪夫多项式定义为

${\xi }_n\left(z\right)=\cos \left(n{\cos }^{-1}z\right)$ , $z_j=\cos \left(\frac{\pi j}{H}\right)$ , $j=0,\cdots ,H$ . (55)

在此，z_j表示切比雪夫节点，而H为任意正整数。每个依赖变量均按照切比雪夫多项式的形式展开如下：

${\Psi }_f\left(z\right)={\displaystyle \sum _{n=0}^H{\Psi }_{fn}{\xi }_n\left(z\right)}$ , ${\Psi }_p\left(z\right)={\displaystyle \sum _{n=0}^H{\Psi }_{pn}{\xi }_n\left(z\right)}$ , $\Theta \left(z\right)={\displaystyle \sum _{n=0}^H{\Theta }_n{\xi }_n\left(z\right)}$ , $\Omega \left(z\right)={\displaystyle \sum _{n=0}^H{\Omega }_n{\xi }_n\left(z\right)}$ . (56)

方程(50)至(54)通过切比雪夫多项式进行了离散化处理，从而得到

$\begin{array}{l}\left(1+{\lambda }_2\omega \right)\left({\displaystyle \sum _{k=0}^H{B}_{jk}{\Psi }_{fk}}-K_{fr}{a}^2{\Psi }_{fj}\right)+\left(1+{\lambda }_1\omega \right)[{\sigma }_f({\displaystyle \sum _{k=0}^H{B}_{jk}{\Psi }_{fk}}-a^2{\Psi }_{fj}\\ -{\displaystyle \sum _{k=0}^H{B}_{jk}{\Psi }_{pk}}+a^2{\Psi }_{pj})-R{a}_{D}ia{\Theta }_j+R{a}_{D}Nia{\Omega }_j]=0,\left(j=1,\cdots ,H-1\right),\end{array}$ (57)

$\begin{array}{l}{K}_r\left(1+{\lambda }_2\omega \right)\left({\displaystyle \sum _{k=0}^H{B}_{jk}{\Psi }_{pk}}-K_{pr}{a}^2{\Psi }_{pj}\right)+\left(1+{\lambda }_1\omega \right)K_r[{\sigma }_f({\displaystyle \sum _{k=0}^H{B}_{jk}{\Psi }_{pk}}-a^2{\Psi }_{pj}\\ -{\displaystyle \sum _{k=0}^H{B}_{jk}{\Psi }_{fk}}+a^2{\Psi }_{fj})-R{a}_{D}ia{\Theta }_j+R{a}_{D}Nia{\Omega }_j]=0,\left(j=1,\cdots ,H-1\right),\end{array}$ (58)

$\frac{\omega }{M}{\Theta }_j-ia\left({\Psi }_{fj}+{\Psi }_{pj}\right)-{\displaystyle \sum _{k=0}^H{B}_{jk}{\Theta }_k}+a^2{\Theta }_j=0$ , ($j=1,\cdots ,H-1$ ), (59)

$\epsilon \omega {\Omega }_j-ia\left({\Psi }_{fj}+{\Psi }_{pj}\right)-\frac{1}{Le}{\displaystyle \sum _{k=0}^H{B}_{jk}{\Omega }_k}+a^2{\Omega }_j=0$ , ($j=1,\cdots ,H-1$ ), (60)

${\Psi }_{f0}={\Psi }_{fH}=0$ , ${\Psi }_{p0}={\Psi }_{pH}=0$ , ${\Theta }_0={\Theta }_H=0$ , ${\Omega }_0={\Omega }_H=0$ , ($j=0$ 和$H$ ), (61)

其中

$A_{jk}=\{\begin{array}{l}\frac{d_j{\left(-1\right)}^{k+j}}{d_k\left(x_j-x_k\right)},j\ne k\hfill \\ \frac{x_j}{2\left(1-x_j^2\right)},1\le j=k\le H-1\hfill \\ \frac{2N^2+1}{6},j=k=0\hfill \\ -\frac{2N^2+1}{6},j=k=H\hfill \end{array}$ , $j=0,1,\cdots ,H$ , (62)

$B_{jk}=A_{jm}\cdot A_{mk}$ , (63)

$d_j=\{\begin{array}{l}2,j=0,H\hfill \\ 1,1\le j\le H-1\hfill \end{array}$ . (64)

方程(57)~(61)可以整理为

$AX=\omega BX$ . (65)

其中，ω和X分别代表与复特征值相关的特征向量，而A和B则是秩为4(H + 1)的四阶复矩阵。

4. 结果和讨论

根据参考文献[12] [16] [17] [20]-[25]，与本研究相符的参数范围选择如下：ε = 0.2，M = 1.5，σ_f = 0.1，λ₁ ~ 0~0.5，λ₂ ~ 0~0.5，Le ~ 2~50，K_r ~ 1~500，K_fr ~0.1~10，K_pr ~ 0.1~10，N ~ 1/30~5/8。

(a) (b)

Figure 2. The numerical neutral stability curves in Ra_D-Le plane when a = π/2, N = 1, K_r = 1, K_pr = 1, K_fr = 1 with different (a) λ₁, (b) λ₂

图2. 当a = π/2，N = 1，K_r = 1，K_pr = 1，K_fr = 1时，Ra_D-Le平面中的数值中性稳定性曲线，分别对应不同的(a) λ₁，(b) λ₂

图2给出了在a = π/2，N = 1，K_r = 1，K_pr = 1，K_fr = 1条件下的数值中性稳定曲线在Ra_D-Le平面内随弛豫时间λ₁ (图2(a))与滞后时间 λ₂(图2(b))的分布特征。两幅子图均呈现明显的峰形结构：随Le从小到大，Ra_D先上升达到极大值(对应的Le记为Le^c)，随后Ra_D随Le继续增大而下降并趋于平缓。由于图中下方为稳定区、上方为不稳定区，峰值处的Ra_D表明在该Le^c条件下系统最不易发生失稳——此处需施加最大的Darcy-Rayleigh数才能越过中性边界进入不稳定区。由图可见，随λ₁或λ₂的增大，整条中性曲线整体上移且峰值位置Le^c向更高的Lewis数方向发生可观测的偏移(各自子图内曲线按参数增大保持有序排列且峰值高度明显增加)，即在各自给定的参数集合内，增大任一粘弹性时间尺度均会提高在所有Le下触发不稳定所需的Ra_D，并使“最不易失稳”的Le^c向更大的Le移动。

图3展示了在Le = 8，N = 1，K_r = 1，K_pr = 1，K_fr = 1条件下，Oldroyd-B流体在双孔隙多孔介质中的数值中性稳定曲线在Ra_D-a平面内的分布规律。图3(a)给出了弛豫时间λ₁的影响，图3(b)则对应滞后时间λ₂。横坐标a为扰动波数，纵坐标Ra_D为Darcy-Rayleigh数；曲线下方为稳定区，上方为不稳定区。曲线的最低点对应临界Darcy-Rayleigh数$R{a}_D^c$ 及其相应波数a^c，表征系统最易发生失稳的模式。

从图3(a)可以看出，随着弛豫时间λ₁的增大，整个中性稳定曲线明显上移，临界$R{a}_D^c$ 由46.72增至68.94，系统稳定性随之增强。同时，曲线两侧逐渐变平，说明波数变化对稳定性的影响趋于减弱。由于λ₁描述了流体应力对形变速率的响应滞后，其增大意味着流体的弹性效应增强，使得扰动所需的临界热驱动力提高，从而推迟对流失稳的发生。

(a) (b)

Figure 3. The numerical neutral stability curves in Ra_D-a plane when Le = 8, N = 1, K_r = 1, K_pr = 1, K_fr = 1 with different (a) λ₁, (b) λ₂

图3. 当 Le = 8，N = 1，K_r = 1，K_pr = 1，K_fr = 1时，不同(a) λ₁，(b) λ₂下Ra_D-a平面内的数值中性稳定性曲线

图3(b)显示，滞后时间λ₂的增大同样使中性稳定曲线上移，临界$R{a}_D^c$ 由44.40升高至98.71，表现出相同的稳定化趋势。这说明无论是应力响应的弛豫效应还是应变速率反馈的滞后效应λ₂，其增强均能抑制双扩散对流的不稳定性。

综合图2与图3的结果可知，图2揭示了粘弹性参数对系统在Ra_D-Le平面内的整体稳定性影响，而图3进一步说明了这种效应在扰动波数空间中的体现形式。两者相互印证，表明增大Oldroyd-B流体的粘弹性时间参数λ₁、λ₂均会延迟对流失稳的发生，并改变最易失稳模态的分布特征。

从力学与能量学角度来看，Oldroyd-B流体的粘弹性对体系稳定性的增强主要来源于弹性应力的储能作用。流体的总应力张量可以表示为τ = 2μD + τ_e，其中D为形变速率张量，τ_e为满足Oldroyd-B型本构关系的弹性应力分量。当弛豫时间λ₁增大时，弹性应力对局部应变速率的响应滞后，使部分扰动动能以弹性能的形式被暂时储存，从而减缓扰动能量的积累速率；而滞后时间λ₂的增加则进一步强化了这种时间延迟效应，使体系对扰动的响应更具粘性阻尼特征，临界Darcy-Rayleigh数因此显著提高。另一方面，Ra_D-Le曲线中出现的峰值可归因于热扩散与溶质扩散时间尺度的竞争。当Lewis数较小时，溶质扩散较快，浮力主要由溶质浓度梯度驱动；随着Lewis数增大，热扩散逐渐滞后，温度梯度增强，导致热浮力在驱动中占据主导，从而使稳定阈值升高；但当Lewis数进一步增大时，浓度场响应减弱，浮力耦合削弱，系统稳定性重新增强，从而形成具有峰值特征的中性曲线。这表明，Oldroyd-B流体的弹性储能效应与热、质扩散之间的竞争共同决定了双扩散对流不稳定性的动力学特征。

图4给出了在不同弛豫时间λ₁与滞后时间λ₂条件下，Darcy-Rayleigh数Ra_D随浮力比N的变化。图中可以清楚地看到，随着N的增加，Ra_D明显上升，表明在所选参数范围内增大浮力比会提高系统的稳定阈值——即需要更大的驱动力(更高的Ra_D)才能触发对流不稳定。由此可推断，在本模型参数设置下，溶质浮力与热浮力的耦合对总体驱动起到削弱作用，从而表现为随N增大系统更难发生失稳。

此外，两类粘弹参数对稳定性具有不同的放大效应。图4(a)显示，增大弛豫时间λ₁会普遍抬高Ra_D，反映出流体的弹性应力能阻碍扰动增长，从而起到稳定作用；而图4(b)中滞后时间λ₂的影响更为显著：随着λ₂增大，不仅整体Ra_D升高，而且在中高N区间出现更陡的上升趋势，表明滞后效应与浮力比耦合后会引发更强的稳定化响应。综上，浮力比N与粘弹参数λ₁，λ₂在双扩散对流的失稳阈值上存在竞争耦合：增大N、λ₁或λ₂均倾向于增强稳定性。

(a) (b)

Figure 4. The numerical neutral stability curves in Ra_D-N plane when a = π/2, Le = 8, K_r = 1, K_pr = 1, K_fr = 1 with different (a) λ₁, (b) λ₂

图4. 当 a = π/2，Le = 8，K_r = 1，K_pr = 1，K_fr = 1时，不同(a) λ₁，(b) λ₂下Ra_D-N平面中的数值中性稳定曲线

图5展示了数值中性稳定曲线在Ra_D-K_r平面随弛豫时间和滞后时间变化的数值结果。图5(a)对应 λ₁，图5(b)对应λ₂。两幅子图均表现出相似的形态学特征：在K_r的低值区间(约K_r ≲ 10) Ra_D随K_r显著增加，曲线斜率较大；随着K_r继续增大，曲线变得趋于平缓并在中高K_r区间呈出平台化趋势，表明K_r对Ra_D的影响存在明显的非线性并趋于饱和。图5(a)中随λ₁增大所得曲线在整个绘制区间内保持有序分布——即较大λ₁对应的曲线位于较小λ₁曲线的上方，曲线间在整个K_r区间均保持可观测的间隔；图5(b)中随λ₂增大的曲线亦呈类似的有序上移，并在各K_r值处保持可辨别的分离。

(a) (b)

Figure 5. The numerical neutral stability curves in Ra_D-K_r plane when a = π/2, Le = 8, N = 1, K_pr = 1, K_fr = 1 with different (a) λ₁, (b) λ₂

图5. 当a = π/2，Le = 8，N = 1，K_pr = 1，K_fr = 1且具有不同的(a) λ₁，(b) λ₂时，在Ra_D-K_r平面内的数值中性稳定性曲线

图6给出了在不同弛豫时间λ₁与滞后时间λ₂条件下，Darcy-Rayleigh数Ra_D随微孔隙各向异性参数K_pr的变化规律。在两幅图中，曲线均呈单调上升趋势，说明随着K_pr的增大，Ra_D持续升高，即系统对微孔隙各向异性增强表现出更强的稳定性。换言之，微孔结构的各向异性越显著，流体在孔隙中的传输阻滞效应越明显，从而需要更高的驱动力(更大的Ra_D)才能触发对流不稳定。

(a) (b)

Figure 6. The numerical neutral stability curves in Ra_D-K_pr plane when a = π/2, Le = 8, N = 1, K_r = 1, K_fr = 1, with different (a) λ₁, (b) λ₂

图6. 当a = π/2，Le = 8，N = 1，K_r = 1，K_fr = 1且具有不同的(a) λ₁，(b) λ₂时，在Ra_D-K_pr平面内的数值中性稳定性曲线

进一步观察图6(a)可见，随弛豫时间λ₁增大，整组中性稳定曲线整体上移，且在各K_pr区间内保持近似平行关系，说明λ₁的增加对系统稳定阈值的提升作用较为均匀。弹性应力的积累使得扰动增长受到抑制，表现为整体稳定性的增强。图6(b)中滞后时间λ₂的作用趋势与λ₁相似：较大的λ₂对应更高的Ra_D，且曲线之间的间隔在整个K_pr区域内保持明显分离，表明滞后效应同样对系统稳定化起积极作用。

综合两幅图可知，微孔隙各向异性与粘弹性时间尺度对系统稳定性的影响具有一致性：增大K_pr、λ₁或λ₂均会提高Ra_D，使双扩散对流更难被激发。可以认为，微孔隙结构各向异性对流体运动的约束，与Oldroyd-B流体中的粘弹性效应相互协同，共同抑制不稳定模态的发展，从而在整体上强化了系统的稳定性。

图7展示了Darcy-Rayleigh数Ra_D随大孔隙各向异性参数K_fr的变化规律，分别考察了弛豫时间λ₁与滞后时间λ₂对系统稳定性的影响。如图7(a)所示，随着K_fr的增加，Ra_D单调上升，表明体系的稳定性增强。物理上，K_fr的增大意味着大孔隙中水平方向的渗透率大于垂直方向的渗透率，从而导致流体在水平方向上的流动占主导地位。这种各向异性效应削弱了流体在垂直方向上的对流能力，从而抑制了双扩散对流的不稳定性。另一方面，在任意固定的K_fr下，增大弛豫时间λ₁会使Ra_D进一步升高，说明Oldroyd-B流体的弹性效应对系统具有稳定作用。较大的λ₁会增强流体的弹性应力，使其对变形的响应更具滞后性，从而抑制扰动能量的增长并延缓对流失稳的发生。

图7(b)进一步显示了滞后时间λ₂的影响趋势。可以看到，Ra_D随K_fr增大而升高的总体趋势与图7(a)一致，表明大孔隙的水平各向异性仍起到稳定作用。而在不同λ₂下，Ra_D均随λ₂的增加显著上升，说明滞后时间同样对系统稳定性具有强化效应。较大的λ₂使流体的黏滞响应具有更强的时间滞后性，从而有效抑制扰动的增长速率，延缓不稳定性的发展。总体而言，弛豫时间和滞后时间的增大均能增强体系的稳定性，而大孔隙各向异性的增强进一步放大了这一稳定效应。

(a) (b)

Figure 7. The numerical neutral stability curves in Ra_D-K_fr plane when a = π/2, Le = 8, N = 1, K_r = 1, K_pr = 1 with different (a) λ₁, (b) λ₂

图7. 当a = π/2，Le = 8，N = 1，K_r = 1，K_pr = 1且具有不同的(a) λ₁ (b) λ₂时，在Ra_D-K_fr平面内的数值中性稳定性曲线

5. 总结

本研究系统探讨了各向异性的双孔隙多孔介质中Oldroyd-B流体的双扩散对流不稳定性特征。结果表明，体系稳定性受Lewis数Le、浮力比N、大孔与微孔渗透率之比K_r、大孔隙的各向异性参数K_fr、微孔隙的各向异性参数K_pr、弛豫时间λ₁与滞后时间λ₂的共同作用显著影响，各参数在对流临界条件及中性稳定曲线形态中均起关键作用。

• 增大λ₁增强流体的弹性储能能力，削弱扰动增长，延迟对流失稳的发生。

• 增大λ₂使应变速率对应力的反馈更具时间滞后性，降低扰动能量积累速率，进一步提高稳定阈值。

• N的增大会提高临界Darcy-Rayleigh数，说明溶质浮力削弱了热浮力的驱动，使体系更难失稳。

• 大孔隙K_fr与微孔隙K_pr的增大均导致Darcy-Rayleigh数升高，表明孔隙各向异性对系统具有稳定作用。

总体而言，Oldroyd-B流体的粘弹性效应与孔隙结构各向异性共同主导了双扩散对流的不稳定特性，增强任一因素均能有效推迟失稳的发生，为非牛顿流体在多孔介质中的传热传质控制提供了理论依据。

基金项目

本研究由国家自然科学基金地区项目(12262026)、内蒙古自治区研究生科研创新项目(KC2024008B)资助。

NOTES

^*通讯作者。