基于LTNE模型的双层Oldroyd⁃B流体双扩散对流不稳定性分析

Instability Analysis of Double Diffusive Convection in a Two⁃Layer Oldroyd⁃B Fluid Based on the LTNE Model

摘要

基于局部热不平衡模型（LTNE）， 本文研究了流体-多孔介质双层系统中Oldroyd-B流体双扩散对流的不稳定性问题， 利用线性稳定性分析和Chebyshev配点法， 得到了振荡对流下的临界瑞利数和临界波数以及中性曲线。文中主要分析了厚度比、 溶质Rayleigh数、 流固热传导系数比率、 流固界面换热系数和粘弹性参数对系统稳定性的影响。经研究发现： 较大的厚度比会破坏双层系统的稳定性， 且会导致系统双峰特性消失； 增加流固热传导系数比率时， 系统会变得更不稳定， 而增加流固界面换热系数则相反， 且局部热不平衡对系统稳定性的影响发生在小波数区域； 溶质Rayleigh数增强了系统的稳定性， 同时促进了振荡对流的发生。

Abstract

Based on the local thermal non-equilibrium （LTNE） model， this paper investigated the instability of double-diffusive convection for an Oldroyd-B fluid in a fluid-porous system. The linear stability analysis and the Chebyshev collocation method were used， and the critical Rayleigh number， the critical wavenumber， and the neutral curves under oscillatory convection were obtained. The effects of the thickness ratio， the solute Rayleigh number， the fluid-solid thermal conductivity ratio， the fluid-solid interfacial heat transfer coefficient， and the viscoelastic parameters on the system’s stability were analyzed. The results show that a greater thickness ratio destabilizes the two-layer system and results in the disappearance of the system’s bimodal characteristics. Increasing the fluid-solid thermal conductivity ratio makes the system more unstable， while increasing the fluid-solid interfacial heat transfer coefficient has the opposite effect. Local thermal non-equilibrium （LTNE） mainly influences the system’s instability in the lower wavenumber region. The solute Rayleigh number enhances the system’s stability and promotes the occurrence of oscillatory convection.

Graphical abstract

关键词

局部热不平衡 / 双扩散对流 / Oldroyd-B流体 / 热不稳定性 / Chebyshev配点法

Key words

local thermal non-equilibrium / double diffusive convection / Oldroyd-B fluid / thermal instability / Chebyshev collocation method

0 引　言

流体层和多孔层组成的多层系统中的对流不稳定性问题在很多工业领域有着广泛应用， 例如薄膜蒸发、 填充床太阳能储存、 地下水流动和事故后反应堆冷却等， 因此被许多研究人员注意并研究。Nield^［1］最早研究了可变形上表面的流体-多孔双层系统的线性稳定性问题； Chen等^［2］研究得出Rayleigh数与波数具有双谷性， 即有两个局部最小值； Chang^［3］对辐射选择性吸收引起的两层系统的热对流进行了研究； Avramenko等^［4］在垂直通流的影响下对双层系统中旋流微生物悬浮液中的生物对流进行了线性分析； Yin等^［5］研究了双层系统中非牛顿Oldroyd-B流体的热不稳定性； Yin等^［6］还研究了具有可变形表面的双层系统的Rayleigh-Marangoni-Bénard对流稳定性。在很多地质系统中， 流体在多孔介质中流动时， 一部分溶质因温度、 压强等发生溶解或者沉淀进而影响流体的传热， 这种情况下， 双扩散对流就会发生。Nield^［7］最早给出了多孔介质中双扩散系统的线性稳定性分析， 得到了稳态对流与振荡对流的临界Rayleigh数； Rudraiah等^［8］最早开始用非线性微扰理论研究了水平多孔层中的双扩散对流； Poulikakos^［9］利用Darcy-Brinkman模型对热溶质对流进行了线性稳定性分析， Rudraiah等^［10］还分析了多孔介质中存在交叉扩散效应时的双扩散对流， 另外， 流体饱和多孔层中的双扩散对流问题后来被许多学者^{［11‐15］}所研究。为方便起见， 前人在研究时多考虑固相和液相局部热平衡， 但在很多现实领域中， 这两者却是局部热不平衡的。关于这一主题的研究非常广泛， Chen等^［15］研究了基于单一多孔层热不平衡的双扩散对流， 用解析的方法分析了Lewis数等无量纲数对稳定性的影响； Malashetty等^［16］分别用线性稳定性分析和非线性稳定性分析研究了局部热不平衡对流的稳定性； Rees等^{［17‐19］}对处于热不平衡状态多孔介质中的自然对流进行了一系列研究； Straughan等^［20］用Darcy-Brinkman模型研究了多孔层中具有Soret效应的双扩散对流， 并给出了先验解； Gangadharaiah^［21］用正则模态方法研究了在施加自下而上通流的双层热不平衡系统的稳定性， 但没有考虑时间项； Sumithra等^［22］ 在考虑热不平衡和热梯度的情况下， 研究了在磁场作用下双层系统的热溶质对流； Shivakumara等^［23］利用热不平衡模型研究了密度和黏度变化对Darcy-Bénard对流的影响。Hema等^［24］使用热不平衡模型研究了热流对双扩散对流的影响。上述研究中， 对于热不平衡双扩散对流一般都是或者只考虑单一多孔层， 或者没有考虑时间项的影响。

本文在上述研究的基础上， 考虑了处于局部热不平衡状态的双扩散双层Oldroyd-B流体的对流稳定性， 通过线性化分析并使用Chebyshev 配点法， 着重考虑了在振荡对流下， 溶质Rayleigh数、 流固界面换热系数、 流固热传导系数比率以及粘弹性参数对Rayleigh-Bénard对流稳定性的影响。

1 数学模型

考虑一个水平的双层系统， 由厚度为d_f的流体层覆盖在厚度为d_p的多孔层上， 多孔层的下面保持恒定温度T_l， 而流体层的上层保持在温度T_u， 流体层与多孔层之间的界面位于z=\mathrm{0}， 交界面处温度为T_{\mathrm{0}}， 如图 1 所示， 流体层和多孔介质层都是饱和Oldroyd-B流体， 其本构方程为^{［25‐26］}

\mathit{\tau }+\overline{\lambda }\left[\frac{\partial \mathit{\tau }}{\partial t}+\left(\mathit{V}\cdot \nabla \right)\mathit{\tau }-{\left(\nabla \mathit{V}\right)}^{\mathrm{T}}\cdot \mathit{\tau }-\mathit{\tau }\cdot \left(\nabla \mathit{V}\right)\right]=\mu \left\{\mathit{A}+\overline{\epsilon }\left[\frac{\partial \mathit{A}}{\partial t}+\left(\mathit{V}\cdot \nabla \right)\mathit{A}-{\left(\nabla \mathit{V}\right)}^{\mathrm{T}}\cdot \mathit{A}-\mathit{A}\cdot \left(\nabla \mathit{V}\right)\right]\right\}，（1）

式中： \mathit{\tau }为应力张量；\mathit{V}为速度； \mathit{A}为应变率张量； \mu为动力黏性系数； \overline{\lambda }为应力松弛时间； \overline{\epsilon }为应变弛豫时间\left(\overline{\lambda }\ge \overline{\epsilon }\right)。同时， 在本问题中， 使用了Boussinesq近似的一种特殊情况， 微观热传递的两相模型， 即固相和液相是热不平衡的， 可以更好地描述多孔层中固相和液相的传热差别， 表示为

\rho ={\rho }_{\mathrm{0}}\left[\mathrm{1}-{\alpha }_T\left(T-T_{\mathrm{0}}\right)+{\alpha }_c\left(C-C_{\mathrm{0}}\right)\right]，（2）

式中： \rho和{\rho }_{\mathrm{0}}分别为温度为T和T_{\mathrm{0}}时的密度； {\alpha }_T为热膨胀系数； {\alpha }_c为由浓度引起的膨胀系数。

1.1 控制方程

流体层在{\mathbf{R}}^{\mathrm{2}}\times \left(\mathrm{0},d_f\right)区域的控制方程为

\nabla \cdot {\mathit{V}}_f=\mathrm{0}，（3）

\begin{array}{c}{\rho }_{\mathrm{0}}\left(\mathrm{1}+\overline{\lambda }\frac{\partial }{\partial t}\right)\left[\frac{\partial {\mathit{V}}_f}{\partial t}+({\mathit{V}}_f\cdot \nabla ){\mathit{V}}_f\right]=\\ -\left(\mathrm{1}+\overline{\lambda }\frac{\partial }{\partial t}\right)\nabla p_f-\left(\mathrm{1}+\overline{\lambda }\frac{\partial }{\partial t}\right)\rho g\mathit{z}+\\ \mu \left(\mathrm{1}+\overline{\epsilon }\frac{\partial }{\partial t}\right){\nabla }^{\mathrm{2}}{\mathit{V}}_f,\end{array}（4）

\frac{\partial T_f}{\partial t}+(\mathit{V}\cdot \nabla )T_f={\kappa }_f{\nabla }^{\mathrm{2}}{T}_f，（5）

\frac{\partial C_f}{\partial t}+({\mathit{V}}_f\cdot \nabla )C_f={\kappa }_c{\nabla }^{\mathrm{2}}{C}_f，（6）

式中： p_f为压力； \mathit{z}为z轴正方向的单位向量； g为重力加速度； T_f为温度； C_f为溶质浓度； {\kappa }_f和{\kappa }_c分别为流体层的热扩散率和溶质热扩散率。

多孔层在{\mathbf{R}}^{\mathrm{2}}\times \left(-d_p,\mathrm{0}\right)上的控制方程为

\nabla \cdot {\mathit{V}}_p=\mathrm{0}，（7）

\begin{array}{c}\frac{{\rho }_{\mathrm{0}}}{\varphi }\left(\mathrm{1}+\overline{\lambda }\frac{\partial }{\partial t}\right)\frac{\partial {\mathit{V}}_p}{\partial t}=-\left(\mathrm{1}+\overline{\lambda }\frac{\partial }{\partial t}\right)\nabla p_p-\\ \left(\mathrm{1}+\overline{\lambda }\frac{\partial }{\partial t}\right)\rho g\mathit{z}-\frac{\mu }{K}\left(\mathrm{1}+\overline{\epsilon }\frac{\partial }{\partial t}\right){\mathit{V}}_P\end{array}，（8）

\begin{array}{c}\varphi {\left({\rho }_{\mathrm{0}}{c}_p\right)}_f\frac{\partial T_{fp}}{\partial t}+{\left({\rho }_{\mathrm{0}}{c}_p\right)}_f({\mathit{V}}_p\cdot \nabla )T_{fp}=\\ \varphi k_f{\nabla }^{\mathrm{2}}{T}_{fp}+h\left(T_{sp}-T_{fp}\right)\end{array}，（9）

\begin{array}{c}\left(\mathrm{1}-\varphi \right){\left({\rho }_{\mathrm{0}}{c}_p\right)}_s\frac{\partial T_{sp}}{\partial t}=\\ \left(\mathrm{1}-\varphi \right)k_s{\nabla }^{\mathrm{2}}{T}_{sp}-h\left(T_{sp}-T_{fp}\right)\end{array}，（10）

\varphi \frac{\partial C_p}{\partial t}+({\mathit{V}}_p\cdot \nabla )C_p=\varphi {\kappa }_c{\nabla }^{\mathrm{2}}{C}_p，（11）

式中： p_p为多孔层的压力； \varphi为孔隙率； T_{fp}和T_{sp}分别为多孔层中液相和固相的温度； C_p为多孔层溶质浓度； k_f和k_s分别为多孔层中液相和固相的热导率； {\left({\rho }_{\mathrm{0}}{c}_p\right)}_f和{\left({\rho }_{\mathrm{0}}{c}_p\right)}_s分别为多孔层中液相和固相的单位热容，多孔层整体单位热容为{\left({\rho }_{\mathrm{0}}{c}_p\right)}_p=\varphi {\left({\rho }_{\mathrm{0}}{c}_p\right)}_f+\left(\mathrm{1}-\varphi \right){\left({\rho }_{\mathrm{0}}{c}_p\right)}_s， 方程（9）~（10）分别为LTNE模型中固相和液相的方程。

基本状态方程为

\begin{array}{c}{\mathit{V}}_{fb}=\mathrm{0},T_{fb}\left(z_f\right)=T_{\mathrm{0}}-\left(T_{\mathrm{0}}-T_u\right)\frac{z_f}{d_f},\hspace{1em}\\ C_{fb}\left(z_f\right)=C_{\mathrm{0}}-\left(C_{\mathrm{0}}-C_u\right)\frac{z_f}{d_f}\hspace{1em}\hspace{1em}\left(\mathrm{0}\le z_f\le d_f\right),\\ {\mathit{V}}_{pb}=\mathrm{0},T_{pb}\left(z_p\right)=T_{\mathrm{0}}-\left(T_l-T_{\mathrm{0}}\right)\frac{z_p}{d_p},\hspace{1em}\\ C_{pb}\left(z_p\right)=C_{\mathrm{0}}-\left(C_l-C_{\mathrm{0}}\right)\frac{z_p}{d_p}\hspace{1em}\hspace{1em}\left(-d_p\le z_p\le \mathrm{0}\right)\mathrm{。}\end{array}（12）

由交界面处通量的连续性， 得到

T_{\mathrm{0}}=\frac{k_p{d}_f{T}_l+k_f{d}_p{T}_u}{k_p{d}_f+k_f{d}_p},C_{\mathrm{0}}=\frac{\varphi {\kappa }_c{d}_f{C}_l+{\kappa }_c{d}_p{C}_u}{\varphi {\kappa }_c{d}_f+{\kappa }_c{d}_p}\mathrm{。}（13）

为了分析其流动稳定性， 在基础态上加入小扰动

M=M_{fb}+{M_f}^{\mathrm{\text{'}}},M_p=M_{pb}+{M_p}^{\mathrm{\text{'}}},（14）

式中： M=V,p,T,C； M_p=V_p,p_p,T_{fp},C_p,T_{sp}。将方程\left(\mathrm{14}\right)代入方程（3）~（11）， 得到线性化控制方程， 对动量方程进行两次旋度变换来消除压力项\nabla p_f和\nabla p_p， 为了方便， 省略扰动的上标。引入下面的尺度， 包括长度、 速度、 时间、 温度、 溶质浓度， 对方程进行无量纲化， 流体层和多孔层的尺度分别为d_f,\frac{{\kappa }_f}{d_f},\frac{d_f^{\mathrm{2}}}{{\kappa }_f},T_{\mathrm{0}}-T_u,C_{\mathrm{0}}-C_u和d_p,\frac{\varphi {\kappa }_p}{d_p},\frac{d_p^{\mathrm{2}}}{{\kappa }_p},T_l-T_{\mathrm{0}},C_l-C_{\mathrm{0}}， 经过上述线性化和无量纲化处理过后的方程为

\left[\frac{\mathrm{1}}{P{r}_f}\left(\mathrm{1}+{\lambda }_f\frac{\partial }{\partial t}\right)\frac{\partial }{\partial t}-\left(\mathrm{1}+{\epsilon }_f\frac{\partial }{\partial t}\right){\nabla }_{}^{\mathrm{2}}\right]{\nabla }_{}^{\mathrm{2}}{w}_f=R{a}_f\left(\mathrm{1}+{\lambda }_f\frac{\partial }{\partial t}\right){\nabla }_{\mathrm{1}}^{\mathrm{2}}{T}_f-R{a}_c\left(\mathrm{1}+{\lambda }_f\frac{\partial }{\partial t}\right){\nabla }_{\mathrm{1}}^{\mathrm{2}}{C}_f,（15）

\frac{\partial T_f}{\partial t}-w_f\mathit{z}={\nabla }_{}^{\mathrm{2}}{T}_f,（16）

\frac{\partial C_f}{\partial t}-w_f=\frac{\mathrm{1}}{L{e}_f}{\nabla }_{}^{\mathrm{2}}{C}_f,（17）

\frac{Da}{\varphi P{r}_p}\left(\mathrm{1}+{\lambda }_p\frac{\partial }{\partial t}\right)\frac{\partial }{\partial t}{\nabla }_{}^{\mathrm{2}}{w}_p+\left(\mathrm{1}+{\epsilon }_p\frac{\partial }{\partial t}\right){\nabla }_{}^{\mathrm{2}}{w}_p=R{a}_p\left(\mathrm{1}+{\lambda }_p\frac{\partial }{\partial t}\right){\nabla }_{\mathrm{1}}^{\mathrm{2}}{T}_{fp}-R{a}_{cp}\left(\mathrm{1}+{\lambda }_p\frac{\partial }{\partial t}\right){\nabla }_{\mathrm{1}}^{\mathrm{2}}{C}_p,（18）

\frac{\partial T_{fp}}{\partial t}-w_p={\nabla }_{}^{\mathrm{2}}{T}_{fp}+H\left(T_{sp}-T_{fp}\right),（19）

{ϵ}_s\frac{\partial T_{sp}}{\partial t}={\nabla }_{}^{\mathrm{2}}{T}_{sp}-\gamma H\left(T_{sp}-T_{fp}\right),（20）

\frac{\partial C_p}{\partial t}-w_p=\frac{\mathrm{1}}{L{e}_p}{\nabla }_p^{\mathrm{2}}{C}_p,（21）

式中： {\kappa }_f=\frac{k_f}{{\left({\rho }_{\mathrm{0}}{c}_p\right)}_f}， {\kappa }_p=\frac{k_p}{{\left({\rho }_{\mathrm{0}}{c}_p\right)}_f}， {\kappa }_s=\frac{k_s}{{\left({\rho }_{\mathrm{0}}{c}_p\right)}_s}；P{r}_f=\frac{\mu }{{\rho }_{\mathrm{0}}{\kappa }_f}，P{r}_p=\frac{\mu }{{\rho }_{\mathrm{0}}{\kappa }_{fp}}为Prandtl数； Da={\delta }^{\mathrm{2}}=\frac{K}{d_p^{\mathrm{2}}}为达西数； H=\frac{h{d}_p^{\mathrm{2}}}{\varphi k_{fp}}为无量纲流固界面换热系数； \gamma =\frac{\varphi k_f}{\left(\mathrm{1}-\varphi \right)k_s}为流固热传导系数比率； {ϵ}_T=\frac{{\kappa }_f}{{\kappa }_p}， {ϵ}_s=\frac{{\kappa }_p}{{\kappa }_s}为热扩散率比； {\chi }_s=\frac{k_s}{k_p}为热导率比； L{e}_f=\frac{{\kappa }_f}{{\kappa }_c},L{e}_p=\frac{{\kappa }_p}{{\kappa }_c}为Lewis数； {\lambda }_f=\frac{\overline{\lambda }{\kappa }_f}{d_f^{\mathrm{2}}},{\lambda }_p=\frac{\overline{\lambda }{\kappa }_p}{d_p^{\mathrm{2}}}为无量纲应力松弛时间； {\epsilon }_f=\frac{\overline{\epsilon }{\kappa }_f}{d_f^{\mathrm{2}}}，{\epsilon }_p=\frac{\overline{\epsilon }{\kappa }_p}{d_p^{\mathrm{2}}}为无量纲应变弛豫时间； R{a}_f=\frac{{\rho }_{\mathrm{0}}g{\alpha }_T{d}_f^{\mathrm{3}}\left(T_{\mathrm{0}}-T_u\right)}{\mu {\kappa }_f}，R{a}_p=\frac{{\rho }_{\mathrm{0}}g{\alpha }_T{d}_{p}K\left(T_l-T_{\mathrm{0}}\right)}{\varphi \mu {\kappa }_p}和R{a}_c=\frac{{\rho }_{\mathrm{0}}g{\alpha }_c{d}_f^{\mathrm{3}}\left(C_{\mathrm{0}}-C_u\right)}{\mu {\kappa }_f}，R{a}_{cp}=\frac{{\rho }_{\mathrm{0}}g{\alpha }_c{d}_{p}K\left(C_l-C_{\mathrm{0}}\right)}{\varphi \mu {\kappa }_p}分别为热Rayleigh数和溶质Rayleigh数； {\nabla }_{\mathrm{1}}^{\mathrm{2}}=\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}为拉普拉斯算子； 下标f表示流体层， 下标p表示多孔层， \nabla为哈密顿算子。

1.2 边界条件

在z=d_f和z=d_p处的边界条件为

z=d_f:w_f=\mathrm{0},T_f=\mathrm{0},\frac{\partial w_f}{\partial z}=\mathrm{0},C_f=\mathrm{0},（22）

z=d_p:w_p=\mathrm{0},T_{fp}=\mathrm{0},T_{sp}=\mathrm{0},C_p=\mathrm{0}\mathrm{。}（23）

在交界面z=\mathrm{0}处的边界条件为

{ϵ}_T{w}_f=\varphi \widehat{d}{w}_p,\widehat{d}{T}_f={ϵ}_T{T}_{fp},\widehat{d}{T}_{fp}={ϵ}_T{T}_{sp},

\frac{\mathrm{d}{T}_f}{\mathrm{d}{z}_f}=\varphi {\mathrm{ϵ}}_T\frac{\mathrm{d}{T}_{fp}}{\mathrm{d}{z}_p}+\left(\mathrm{1}-\varphi \right){\chi }_s\frac{\mathrm{d}{T}_{sp}}{\mathrm{d}{z}_p},

\widehat{d}\varphi C_f=C_p,\varphi \frac{\mathrm{d}{C}_f}{\mathrm{d}{z}_f}=\frac{\mathrm{d}{C}_p}{\mathrm{d}{z}_p},（24）

\frac{{\widehat{d}}^{\mathrm{4}}}{{ϵ}_{T}P{r}_p}\left(\mathrm{1}+{\lambda }_p\frac{\partial }{\partial t}\right)\frac{\partial }{\partial t}\frac{\partial w_p}{\partial z_p}+\frac{\varphi {\widehat{d}}^{\mathrm{4}}}{{ϵ}_T{\delta }^{\mathrm{2}}}\left(\mathrm{1}+{\epsilon }_p\frac{\partial }{\partial t}\right)\frac{\partial w_p}{\partial z_p}=

\frac{\mathrm{1}}{P{r}_f}\left(\mathrm{1}+{\lambda }_f\frac{\partial }{\partial t}\right)\frac{\partial }{\partial t}\frac{\partial w_f}{\partial z_f}-

\left(\mathrm{1}+{\epsilon }_f\frac{\partial }{\partial t}\right)\left({\nabla }^{\mathrm{2}}\frac{\partial w_f}{\partial z_f}+\mathrm{2}{\nabla }_{\mathrm{1}}^{\mathrm{2}}\frac{\partial w_f}{\partial z_f}\right),（25）

\left(\mathrm{1}+{\epsilon }_f\frac{\partial }{\partial t}\right)\frac{{\partial }^{\mathrm{2}}{w}_f}{\partial {z_f}^{\mathrm{2}}}=

\frac{\alpha \widehat{d}}{\delta }\left(\mathrm{1}+{\lambda }_f\frac{\partial }{\partial t}\right)\left(\frac{\partial w_f}{\partial z_f}-\frac{\varphi {\widehat{d}}^{\mathrm{2}}}{{ϵ}_T}\frac{\partial w_p}{\partial z_p}\right),（26）

式中： \widehat{d}为流体层与多孔层的厚度比。

无量纲参数关系为

P{r}_f=\frac{P{r}_p}{{\mathrm{ϵ}}_T},R{a}_f=\frac{\varphi {\widehat{d}}^{\mathrm{4}}}{{ϵ}_T^{\mathrm{2}}{\delta }^{\mathrm{2}}}R{a}_p,

R{a}_c=\frac{{\varphi }^{\mathrm{2}}{\widehat{d}}^{\mathrm{4}}}{{ϵ}_T{\delta }^{\mathrm{2}}}R{a}_{cp},{\lambda }_f=\frac{{ϵ}_T}{{\widehat{d}}^{\mathrm{2}}}{\lambda }_p,{\epsilon }_f=\frac{{ϵ}_T}{{\widehat{d}}^{\mathrm{2}}}{\epsilon }_p,

L{e}_f={ϵ}_{T}L{e}_p\mathrm{。}（27）

2 简正模态

利用简正模态的方法， 设置扰动量为

\left(w_f,T_f,C_f\right)=\left(W_f\left(z_f\right),{\Theta }_f\left(z_f\right),C_f\left(z_f\right)\right)f\left(x,y\right)e^{{\sigma }_{f}t},（28）

\left(w_p,T_{fp},T_{sp},C_p\right)=\left(W_p\left(z_p\right),{\Theta }_p\left(z_p\right),{\Phi }_s\left(z_p\right),C_p\left(z_p\right)\right)f\left(x,y\right)e^{{\sigma }_{p}t},（29）

式中： {\mathrm{\sigma }}_f={\mathrm{\sigma }}_f^r+\mathrm{i}{\mathrm{\sigma }}_f^i,\hspace{1em}{\mathrm{\sigma }}_p={\mathrm{\sigma }}_p^r+\mathrm{i}{\mathrm{\sigma }}_p^i， 两个实部{\sigma }_f^r,{\sigma }_p^r代表扰动随时间的增长率， 而虚部{\sigma }_f^i,{\sigma }_p^i是时间频率， 显然， 实部{\sigma }^r=\mathrm{0}是稳定的临界状态， 此时若满足虚部{\sigma }^i=\mathrm{0}则为稳态对流， {\sigma }^i\ne \mathrm{0}则振荡对流， {\nabla }_f^{\mathrm{2}}{f}_f+a_f^{\mathrm{2}}{f}_f=\mathrm{0},\hspace{1em}{\nabla }_p^{\mathrm{2}}{f}_p+a_p^{\mathrm{2}}{f}_p=\mathrm{0}， a_f,a_p是无量纲水平波数。由于无量纲尺度的关系， 上述无量纲物理量需满足

a_f=\widehat{d}{a}_p{\sigma }_f=\frac{{\widehat{d}}^{\mathrm{2}}}{{ϵ}_T}{\sigma }_p\mathrm{。}（30）

将方程（28）~（29）代入方程（15）~（21）得

\left[\left(\mathrm{1}+{\mathrm{\epsilon }}_f{\sigma }_f\right)\left(D_f^{\mathrm{2}}-a_f^{\mathrm{2}}\right)-\frac{{\mathrm{\sigma }}_f}{P{r}_f}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right)\right]\left(D_f^{\mathrm{2}}-a_f^{\mathrm{2}}\right)W_f=R{a}_f{a}_f^{\mathrm{2}}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right){\Theta }_f-R{a}_c{a}_f^{\mathrm{2}}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right)C_f,（31）

\left(D_f^{\mathrm{2}}-a_f^{\mathrm{2}}-{\sigma }_f\right){\Theta }_f=-W_f,（32）

\frac{\mathrm{1}}{L{e}_f}\left(D_f^{\mathrm{2}}-a_f^{\mathrm{2}}\right)C_f=-W_f+{\mathrm{\sigma }}_f{C}_f,（33）

\left(\mathrm{1}+{\epsilon }_p{\sigma }_p\right)\left(D_p^{\mathrm{2}}-a_p^{\mathrm{2}}\right)W_p+{\sigma }_p\frac{Da}{\varphi P{r}_p}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right)\left(D_p^{\mathrm{2}}-a_p^{\mathrm{2}}\right)W_p=-R{a}_p{a}_p^{\mathrm{2}}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right){\Theta }_p+R{a}_{cp}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right)a_p^{\mathrm{2}}{C}_p,（34）

\left(D_p^{\mathrm{2}}-a_p^{\mathrm{2}}\right){\Theta }_p+H\left({\Phi }_s-{\Theta }_p\right)+W_p={\sigma }_p{\Theta }_p,（35）

\left(D_p^{\mathrm{2}}-a_p^{\mathrm{2}}\right){\Phi }_s-\gamma H\left({\Phi }_s-{\Theta }_p\right)={ϵ}_s{\mathrm{\sigma }}_p{\Phi }_s,（36）

\frac{\mathrm{1}}{L{e}_p}\left(D_p^{\mathrm{2}}-a_p^{\mathrm{2}}\right)C_p=-W_p+{\sigma }_p{C}_p\mathrm{。}（37）

边界条件经过该处理后为

z=\mathrm{1},W_f=\mathrm{0},{\Theta }_f=\mathrm{0},C_f=\mathrm{0},D_f{W}_f=\mathrm{0},（38）

z=-\mathrm{1},W_p=\mathrm{0},{\Theta }_p=\mathrm{0},{\Phi }_s=\mathrm{0},C_p=\mathrm{0}\mathrm{。}（39）

z=\mathrm{0}的边界条件为

\begin{array}{c}{ϵ}_T{W}_f=\varphi \widehat{d}{W}_p,\widehat{d}{\Theta }_f={ϵ}_T{\Theta }_p,\widehat{d}{\mathrm{\Theta }}_f={ϵ}_T{\Phi }_s,\\ D_f{\Theta }_f=\varphi {ϵ}_T{D}_p{\Theta }_p+\left(\mathrm{1}-\varphi \right){\chi }_p{D}_p{\Phi }_s,\\ \widehat{d}\varphi C_f=C_p,\varphi D_f{C}_f=D_p{C}_p，\end{array}（40）

\frac{{\widehat{d}}^{\mathrm{4}}{\sigma }_p}{{ϵ}_{T}P{r}_p}\left(\mathrm{1}+{\lambda }_p{\sigma }_p\right)D_p{W}_p+\frac{\varphi {\widehat{d}}^{\mathrm{2}}}{{ϵ}_T{\delta }^{\mathrm{2}}}\left(\mathrm{1}+{\epsilon }_p{\sigma }_p\right)D_p{W}_p=\frac{\mathrm{1}}{P{r}_f}\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right)D_f{W}_f-\left(\mathrm{1}+{\epsilon }_f{\sigma }_f\right)\left(D_f^{\mathrm{3}}{W}_f-\mathrm{3}{a_f}^{\mathrm{2}}{D}_f{W}_f\right)，（41）

\begin{array}{c}\left(\mathrm{1}+{\epsilon }_f{\sigma }_f\right)D_f^{\mathrm{2}}{W}_f=\\ \frac{\alpha \widehat{d}}{\delta }\left(\mathrm{1}+{\lambda }_f{\sigma }_f\right)\left(D_f{W}_f-\frac{\varphi {\widehat{d}}^{\mathrm{2}}}{{ϵ}_T}{D}_p{W}_p\right)，\end{array}（42）

式中： D_f=\frac{\mathrm{d}}{\mathrm{d}{z}_f},D_p=\frac{\mathrm{d}}{\mathrm{d}{z}_p}为流体层和多孔层的微分算子。

3 数值方法

上述过程通过线性化和无量纲化将偏微分方程转化为了一个十六阶常微分方程组， 然后使用简正模态理论， 将小扰动分解为各种模态， 每一模态都满足线性系统， 对于一个模态如果实部大于0， 则不稳定， 等于0是中性稳定， 小于0则是稳定。本文要求解出每个波数下， 保持所有扰动模态都稳定的临界瑞利数， 并通过每次调整一个物理参数， 得到对应瑞利数与波数的关系曲线来分析该参数对稳定性的影响。数值方法采用Chebyshev配点法^［27］， 其基本思路如下：

变量W_f展开的Chebyshev级数为

W_f={\displaystyle \sum _{n=\mathrm{0}}^N}{W}_n{T}_n\left({\widehat{z}}_f\right)\mathrm{。}（43）

其他的变量如{\Theta }_f,C_f,W_p,{\Theta }_p,{\Phi }_s,C_p也类似， 由此， 问题转化为了广义特征值问题\mathit{A}\mathit{X}=R{a}_p\mathit{B}\mathit{X}， 其中\mathit{A}和\mathit{B}均为\mathrm{8}\times (N+\mathrm{1})阶方阵， \mathit{X}=（W_f，{\Theta }_f，Q，C_f，W_p，{\Theta }_p，{\Phi }_s，C_p{）}^{\mathrm{T}}， Q=\left(D_f^{\mathrm{2}}-a_f^{\mathrm{2}}\right)W_f。对比不同N值下的稳态对流和振荡对流的临界波数a_{mc}、 临界Rayleigh数R{a}_{mc}， 分别取基函数个数为N=\mathrm{5,10,15,20}， 通过对比发现N=\mathrm{15}之后的数值结果已经稳定， 为了同时兼顾数值精度和计算速度， 本文选择取N=\mathrm{15}进行数值计算。

4 结果与分析

本文研究的是非牛顿流体的Rayleigh-Bénard对流， 其本质是从下方加热流体会产生热不稳定性， 这种不稳定性达到一定程度就会产生热对流， 通过上述过程， 方程转化为了一个广义特征值问题， 特征值数R{a}_p是临界瑞利数， 代表不稳定化的浮力和稳定化的扩散力之比， 当大于这个值时， 就会出现不稳定的对流， 小于这个值则代表在当前波数下所有扰动模态都稳定。通过数值算法得到关于临界瑞利数和波数之间的中性曲线， 中性曲线描述的是临界瑞利数与波数的关系， 也是稳定区域与不稳定区域的分界线， 曲线下方为稳定区域， 上方为不稳定区域， 利用中性曲线可以分析出非牛顿流体的对流稳定性。

文中部分参数取{ϵ}_T=0.7， {ϵ}_s=0.3， \delta=0.001， \alpha=0.1， P{r}_p=15， {\lambda }_p=0.5， {\epsilon }_p=0.4， \varphi=0.3， \gamma=2， H=25， \widehat{d}=0.06， R{a}_{cp}=10作为参考， 相比于牛顿流体， 非牛顿流体存在振荡对流， 即考虑{\sigma }^r=\mathrm{0},{\sigma }^i\ne \mathrm{0}， 一般来说振荡对流先于稳态对流出现， 因而其中性曲线一般位于稳态对流的下方， 图 2 所示曲线即表明了在厚度比为\widehat{d}=\mathrm{0.06}的情况下振荡对流和稳态对流的区别。图中振荡对流曲线位于稳态对流曲线的下方， 代表振荡对流稳定性低于稳态对流， 本文主要讨论非牛顿流体的振荡对流。

4.1 流固热传导系数比率\gamma对热不平衡Rayleigh⁃Bénard对流稳定性的影响

流固热传导系数比率\gamma对Rayleigh数曲线的影响如图 3 所示。

从图 3 中可以看出： 不同于单层的流体或者多孔系统， 双层系统的Rayleigh数曲线大多呈现双模态， 曲线呈现双谷性， 即存在两个极小值； 曲线的最小值即临界瑞利数位于短波区域， 说明系统稳定性由多孔层主导； 在波数较小的时候， 随着\gamma的增加Rayleigh数曲线逐渐向下移动， Rayleigh数减小， 说明流固热传导系数比率\gamma降低了系统的稳定性， 这种作用在波数较大的时候， 又逐渐消失； 由此说明流固热传导系数比率的影响不是全区域的， 其对系统稳定性的影响集中在小波数区域。

4.2 流固界面换热系数H对热不平衡Rayleigh⁃Bénard对流稳定性的影响

不同流固界面换热系数H的Rayleigh-Bénard对流如图 4 所示。从图 4 中可以看出： 随着H的增加， 振荡对流的瑞利数逐渐上升， 因而当H越大时， 对流越稳定； 曲线呈现双谷性， 且临界瑞利数位于短波区域； 当波数较大时， Rayleigh数的曲线趋于重合， 说明H对较大波数的Rayleigh-Bénard对流稳定性的影响较小。

4.3 厚度比对热不平衡Rayleigh⁃Bénard对流稳定性的影响

不同厚度比对Rayleigh数的影响如图 5 所示。

从图 5 中可以看出： 随着厚度比的增加， 曲线逐渐向下移动， 临界Rayleigh数不断降低， 说明系统稳定性下降； 小波数区域的下降幅度较小， 而大波数区域下降更快； 曲线的最小值先是在小波数区域， 后转变为在大波数区域， 说明系统稳定性由多孔层主导转向流体层主导； 随着厚度比的增加， 曲线的双谷性逐渐减弱， 这是因为厚度比对大波数区域稳定性的破坏更大， 导致峰值减弱直到消失。

4.4 溶质Rayleigh数对热不平衡Rayleigh⁃Bénard对流稳定性的影响

探讨溶质Rayleigh数对对流稳定性的影响是本文的一个重点， 在其他参数固定的情况下， 溶质Rayleigh数对稳定性影响的曲线如图 6 所示。从图 6 中可以看出， 随着溶质Rayleigh数的增加， 振荡对流的曲线都呈上升趋势， 即溶质Rayleigh数对系统的稳定性起正向影响， 增强了对流稳定性； 曲线呈现双谷性， 且稳定性由多孔层主导； 溶质Rayleigh数越大， 双谷性越往大波数区域偏移。

4.5 不同应力松弛时间对热不平衡Rayleigh⁃Bénard对流稳定性的影响

不同应力松弛时间下的Rayleigh曲线如图 7 所示。从图 7 中可以看出： 随着{\lambda }_p的增加， Rayleigh数曲线逐渐向下移动， 说明应力松弛时间会加强系统的不稳定性， 而这种作用对于较大波数区域相对更明显； 曲线呈现双模态， 在短波区域和长波区域各有一个最小值， 而短波区域最小值小于长波区域， 系统的不稳定性由多孔层控制。

4.6 不同应变弛豫时间对热不平衡Rayleigh⁃Bénard对流稳定性的影响

不同应变弛豫时间的Rayleigh曲线如图 8 所示。从图 8 中可以看出： 随着{\epsilon }_p的增加， Rayleigh数曲线逐渐向上移动， 说明应变弛豫时间会加强系统的稳定性， 这种作用对于较大波数区域相对更明显； 曲线呈现双模态， 在短波区域和长波区域各有一个最小值， 而长波区域最小值小于短波区域， 系统的不稳定性由流体层控制。

5 结　论

本文研究了双层Oldroyd-B流体在局部热不平衡（LTNE）条件下， 多个物理参数对双层Oldroyd-B流体双扩散对流不稳定性的影响， 并对中性曲线进行了数值模拟。结果表明： 溶质Rayleigh数对系统稳定性影响很大， 它不仅增强了系统稳定性， 还促进了振荡对流的发生； 热不平衡对系统稳定性的影响较小且集中在小波数区域， 其中， 流固界面换热系数增强了系统稳定性， 而流固热传导系数比率则对系统稳定性起到破坏作用； 较大的厚度比破坏了系统稳定性， 且会改变对流稳定性的主导层； 粘弹性中， 粘性越大， 弹性越小， 系统越稳定。

参考文献

[1] NIELD D A. Onset of convection in a fluid layer overlying a layer of a porous medium［J］. Journal of Fluid Mechanics， 1977， 81（3）： 513-522.

[2] CHEN F， CHEN C F. Onset of finger convection in a horizontal porous layer underlying a fluid layer［J］. Journal of Heat Transfer， 1988， 110（2）： 403-409.

[3] CHANG M H. Stability of convection induced by selective absorption of radiation in a fluid overlying a porous layer［J］. Physics of Fluids， 2004， 16（10）： 3690-3698.

[4] AVRAMENKO A A， KUZNETSOV A V. The onset of convection in a suspension of gyrotactic microorganisms in superimposed fluid and porous layers： Effect of vertical throughflow［J］. Transport in Porous Media， 2006， 65（2）： 159-176.

[5] YIN C， FU C， TAN W. Stability of thermal convection in a fluid-porous system saturated with an Oldroyd-B fluid heated from below［J］. Transport in Porous Media， 2013， 99（2）： 327-347.

[6] YIN C， LUAN， Z， WANG S. Rayleigh-marangoni-Bénard instability in an oldroyd‐b fluid layer overlying a highly porous layer with a deformable surface［J］. International Journal of Heat and Mass Transfer， 2023， 209： 124148.

[7] NIELD D A. Onset of Thermohaline convection in a porous medium［J］. Water Resources Research， 1968， 4（3）： 553-560.

[8] RUDRAIAH N， SRIMANI P K， FRIEDRICH R. Finite amplitude convection in a two-component fluid saturated porous layer［J］. International Journal of Heat and Mass Transfer， 1982， 25（5）： 715-722.

[9] POULIKAKOS D. Double diffusive convection in a horizontal sparcely packed porous layer［J］. International Communications in Heat and Mass Transfer， 1986， 13（5）： 587-598.

[10] RUDRAIAH N， MALASHETTY M S. The influence of coupled molecular diffusion on double-diffusive convection in a porous medium［J］. ASME Journal of Heat Transfer， 1986， 108（4）： 872-876.

[11] TAUNTON J W， LIGHTFOOT E N， GREEN T. Thermohaline instability and salt fingers in a porous medium［J］. The Physics of Fluids， 1972， 15（5）： 748-753.

[12] TASLIM M E， NARUSAWA U. Binary fluid convection and double-diffusive convection in a porous medium［J］. Journal of Heat and Mass Transfer， 1986， 108（1）： 221-224.

[13] TREVISAN O V， BEJAN A. Mass and heat transfer by natural convection in a vertical slot filled with porous medium［J］. International Journal of Heat and Mass Transfer， 1986， 29（3）： 403-415.

[14] MURRAY B T， CHEN C F. Double-diffusive convection in a porous medium［J］. Journal of Fluid Mechanics， 1989， 201（4）： 147-166.

[15] CHEN X， WANG S， TAO J， et al. Stability analysis of thermosolutal convection in a horizontal porous layer using a thermal non-equilibrium model［J］. International Journal of Heat and Fluid Flow， 2011， 32（1）： 78-87.

[16] MALASHETTY M S， SHIVAKUMARA I S， KULKARNI S. The onset of convection in an anisotropic porous layer using a thermal non-equilibrium model［J］. Transport in Porous Media， 2005， 60（2）： 199-215.

[17] REES D A S， POP I. Free convective stagnation-point flow in a porous medium using a thermal nonequilibrium model［J］. International Communications in Heat and Mass Transfer， 1999， 26（7）： 945-954.

[18] REES D A S， POP I. Vertical free convective boundary-layer flow in a porous medium using a thermal nonequilibrium model［J］. Journal of Porous Media， 2000， 3（1） ： 31-44.

[19] REES D A S， POP I. Vertical free convective boundary layer flow in a porous medium using a thermal non-equilibrium model： Elliptic effects［J］. Zeitschrift für angewandte Mathematik und Physik ZAMP， 2003， 54（3）： 437-448.

[20] STRAUGHAN B， HUTTER K. A priori bounds and structural stability for double-diffusive convection incorporating the soret effect［J］. Proceedings of the Royal Society of London. Series A： Mathematical， Physical and Engineering Sciences， 1999， 455（1983）： 767-777.

[21] GANGADHARAIAH Y H. LTNE effects on two-layer configuration with throughflow［J］. Heat Transfer， 2024， 53（5）： 2294-2310.

[22] SUMITHRA R， VENKATRAMAN S. Outcomes of uniform as well as non-uniform temperature profiles on the onset of double diffusive magneto-darcy-rayleigh-benard convection in a two layer set up in the presence of local thermal non-equilibrium［J］. Journal of Mines， Metals and Fuels， 2022， 70（7A） ： 38-52.

[23] SHIVAKUMARA I S， MAMATHA A L， RAVISHA M. Effects of variable viscosity and density maximum on the onset of Darcy-Bénard convection using a thermal nonequilibrium model［J］. Journal of Porous Media， 2010， 13（7）： 613-622.

[24] HEMA M， SHIVAKUMARA I S， RAVISHA M. Double diffusive LTNE porous convection with cattaneo effects in the solid［J］. Heat Transfer， 2020， 49（6）： 3613-3629.

[25] NIELD D A， BEJAN A. Convection in porous media［M］. New York： Springer， 2006.

[26] BIRD R B， ARMSTRONG R C， HASSAGER O. Dynamics of polymeric liquids： Vol.1， fluid mechanics［M］. New York： Wiley， 1977.

[27] TREFETHEN L N. Spectral methods in matlab［M］. Philadelphia： Society for Industrial and Applied Mathematics， 2000.

基金资助

江苏省青年基金(BK20170775)