粘弹性流体在开口多孔介质中的热对流启动

舒晓晖，穆朗枫，牛 骏，丁 健

(1.中国石油化工股份有限公司，北京100083;2.密苏里科技大学，美国罗拉密苏里州65401;3.中国石油大学，北京102249;4.中国海油 伊拉克有限公司，北京100010)

粘性和非线性的存在使得多孔介质中流体热对流研究极为复杂，但是相关研究成果在石油工程、生物工程等方面均有着广泛的应用前景［1－4］。因此，目前国内外有大量学者正在该领域开展研究工作。现有成果对于多孔介质内牛顿流体的流动与传热已经进行了广泛深入分析［5－9］，例如 Kubitschek［7］等通过引入毕奥数对3类底部加热边界条件下多孔介质中牛顿流体的热对流启动进行了深入研究，并总结了各参数对于热对流启动点的影响。但是对于非牛顿流体在多孔介质内的流动与传热研究还相对匮乏。

Kim等［10］对上、下两边均处于封闭且等温加热边界条件下的多孔介质内粘弹性流体热对流进行了线性和非线性稳定性分析，结果热对流启动时总是振荡对流，振荡对流存在超临界和次临界分叉现象。Jun等［11］对封闭边界多孔介质内黏弹性流体在不同底部加热边界条件下的热对流启动进行了线性稳定性分析，结果等热流加热边界条件下热对流启动最早。但目前极少有针对非封闭多孔介质内黏弹性流体热对流启动特征的研究。

文中通过引入修正后的Darcy定律，对上端开口、底部等温加热条件下多孔介质内粘弹性流体热对流启动进行了线性稳定新分析，得出了热对流启动时的临界条件及黏弹性特征对于启动点的影响。

1 数学模型

考虑一个高度为d，长度为a的二维矩形多孔介质。垂直两侧面绝热，并且流体不能透过。多孔介质上下边界维持恒温T1，T2，且T2－T1=ΔT。上边界开口，并与外界连通，压强为p0。假设多孔介质的渗透率为K，并且其中充满了不可压粘弹性流体，该流体动力学粘度为μ，热膨胀系数为β，密度为ρ，整个系统的热扩散系数为κ。通过运用修正的Darcy定律和Boussinesq近似，可以得到粘弹性流体在多孔介质中运动的控制方程组:

该模型的边界条件为:

其中边条件(5b)为等压边界条件，用于描述多孔介质上端开口边界。将其带入(2)得

为了方便计算，对上述方程组中的变量进行了无量纲化。这里选择的无量纲化参量为:长度d，时间d2/κ，速度 κ/d，压强 κμ/K，温度 ΔT。

2 线性稳定性分析

将控制方程组和边界条件无量纲化后，引入温度扰动量θ和垂直速度扰动量w得:

其中，λ，ε为无量纲化后的松弛时间和延迟时间，Ra=Kρ1βgdΔT/μκ 为瑞利数。边界条件为:

根据流动模式分析，温度、速度扰动量在水平方向上均存在周期性:

其中，α=mπ/a，m为正整数，Θ(z)为关于z的函数。

将方程(16)代入(12)得垂直速度扰动量

式中:D=∂/∂z为微分算子。

将方程(16)，(17)代入方程(12)—方程(14)得关于z的常微分方程:

分析该微分方程形式易知其通解为:

这里，Aj(j=1，2，3，4)是待定常数。将(19)代回原方程(18)得kj满足方程:

由边界条件易得A2=－A1，A4=－A3。

由于方程(20)中仅含有k的偶数次幂，因此k的4个解满足k2=－k1，k4=－k3。进一步求解(20)可得:

由边界条件式(13)和式(14)知 A1，A3和 k1，k3满足下列方程组:

当多孔介质内开始发生对流时，内部流体的速度分布不能处处为零。所以线性方程组(17)和(18)必须有非零解，即方程组系数矩阵的行列式必须为零:

上式中当σ=0时代表的是稳态对流。而对于粘弹性流体，热对流发生时为振荡对流，此时σ应该是一个纯虚数iω，这里ω定义为振荡频率。在此情况下，kj(j=1，2，3，4)与 σ 耦合在一起，需要同时进行求解，因此无法得出解析解形式。我们使用了数值方法对振荡对流情况进行编程求解。

3 结果与讨论

为了研究λ，ε对多孔介质中的粘弹性流体对流产生和流动模式的影响，选取下面四组粘弹性参数值作为算例:(1)λ=0.003，ε=0.001;(2)λ =0.3，ε=0.2;(3)λ =0.2，ε =0.1;(4)λ =0.3，ε =0.1。

3.1 粘弹性对热对流启动的影响

算例(1)中的松弛时间和延迟时间值均很小，其热对流启动时的临界瑞利数Rac随多孔介质边长变化如图1中所示。随着边长的增加，Rac值振荡衰减。其全局最小瑞利数值为29.298，该值与前人得出的开口多孔介质中牛顿流体热对流启动时的临界瑞利数29.3非常接近。由于牛顿流体热对流启动时总是稳态对流，因此可以推断当粘弹性参数都很小时，微弱的粘弹性效应无法引起振荡对流，首先出现的是稳态对流。

图1 算例(1)中临界瑞利数随多孔介质边长变化Fig.1 Variation of the critical Rayleigh number as a function of the medium geometry for case(1)

图2 算例(2)，(3)和(4)中临界瑞利数随多孔介质边长变化Fig.2 Variation of the critical Rayleigh number as a function of the medium geometry for cases(2)，(3)and(4)respectively

图2a—c显示算例(2)，(3)和(4)对应的临界瑞利数随多孔介质边长的变化，图中也标出了对应不同阶段最佳模态m的取值。3个算例的全局最小临界瑞利数为28.953，27.567和19.112，均小于稳态对流启动时的临界瑞利数29.298，计算结果也表明对流启动时均为振荡对流。由此得出，当流体的粘弹性参数足够大时，振荡对流优先启动，并且其启动点早于稳态对流。此外，从图2中还可以发现，随着多孔介质长度的增加，速度场和温度分布中的最佳模态在不断增大。

算例(2)和(4)具有相同的松弛时间0.3s，延迟时间不同。对比两例的计算结果，图2a和图2c显示，对每一个固定边长a，临界瑞利数Rac随延迟时间ε值的增大而增大，说明较大的延迟时间能够推迟振荡对流的发生。两幅子图中Rac的两个相邻极小值之间的距离相同，这个距离随着ε值的减小而缩短。

算例(3)和(4)中的延迟时间相同，均为0.1，松弛时间不同，对比图2b和图2c，对于任意一个固定的边长a，Rac值随着λ值的增大而减小。因此，较大的松弛时间值能够促进粘弹性流体在多孔介质内振荡对流的发生。上述两个算例对应的临界瑞利数值不尽相同，但是两条曲线的变化趋势非常相似，且对应的两幅图中的两个相邻极值间的距离一致。由此可得出结论，此时松弛时间λ的取值对于最优模态的变化没有影响。

3.2 粘弹性对温度分布的影响

取边长为3.2，时间 t=2nπ/ω，图 3 a，b 和 c分别为算例(1)，(2)和(4)中热对流启动时的温度扰动量等值线图，其流动模式分别为m=2，3，4。其中亮斑表示温度较高区域，暗斑表示温度较低区域。可以看出，亮斑与暗斑交替出现，亮、暗斑总数与流动模式有关。在上下两边界处的温度等值线均与边界平行，这与上下两边的等温边界条件一致。算例(2)，(4)对应的振荡对流温度等值线图明显比算例(1)对应的稳态对流温度等值线图复杂。由于松弛时间λ的取值对于最优模态的变化没有影响，因此算例(3)对应的热对流启动时的温度等值线图形式应与算例(4)一致。

3.3 与闭口多孔介质内计算结果的对比

图4中显示的是在开口和闭口边条件下对应算例(2)的临界瑞利数随多孔介质边长变化曲线，其中闭口情况对于曲线一直位于开口情况对应曲线上方，二者的全局最小临界瑞利数分别为28.953和32.643。余下3个算例(1)，(3)和(4)在开口和闭口条件下的临界瑞利数计算结果也呈现出相同的规律。上述结果说明开口边界条件对于多孔介质中的流动起到扰动作用，有利于热对流的启动。

图3 算例(1)，(2)和(4)中，取边长为3.2、时间为2nπ/ω时热对流启动时的温度扰动量等值线Fig.3 Contours of the disturbed temperature at the onset point for cases(1)，(2)and(4)respectively，with the length of the medium being 3.2 and the time being 2nπ/ω

图4 算例(2)对应开口和闭口边界条件时的临界瑞利数随多孔介质边长变化Fig.4 Variation of the critical Rayleigh number with the length of porous media corresponding the open-top boundary conditions and closed-top boundary conditions for case(2)

4 结论

1)对于上端开口底部受等温加热多孔介质内的粘弹性流体，当粘弹性参数都很小时微弱的粘弹性效应无法引起振荡对流，首先出现的是稳态对流。

2)当流体的黏弹性参数足够大时，振荡对流优先启动，且振荡对流启动时的临界瑞利数随多孔介质边长增大而呈振荡衰减。

3)松弛时间的增大及延迟时间的减小能够使得粘弹性流体在多孔介质内的振荡对流提前启动。但是松弛时间的取值对于最优模态的变化没有影响。

4)相对于闭口边界条件，开口边界条件对多孔介质内流动具有更强的扰动作用，能够提前热对流启动。

符号注释

d:多孔介质高度，m;a:多孔介质长度，m;T1:多孔介质上边界温度，K;T2:多孔介质下边界温度，K;p0:多孔介质上边界压力，Pa;K:多孔介质渗透率，m2;μ:粘弹性流体粘度，Pa·s;β:粘弹性流体热膨胀系数，1/K;ρ:粘弹性流体密度，kg/m3;κ:多孔介质内粘弹性流体热扩散系数多孔介质内的速度矢量，m/s;u*:水平速度分量，m/s;v*:垂直速度分量，m/s;p*:粘弹性流体压力，Pa;g:重力加速度，m/s2:延迟时间，s:松弛时间，s;:垂直方向单位矢量;ρ1:温度在T1时的流体密度，kg/m3;θ:无量纲温度扰动量;w:无量纲垂直速度扰动量;ε:无量纲延迟时间;λ:无量纲松弛时间;Ra:瑞利数;m:任意整数;Θ(z):关于z的函数;D:微分算子;ω:振荡频率。

［1］王补宣.多孔介质传热研究的意义和现状［J］.中国科学基金，1991，1:32 －34.Wang Buxuan.Meaning and of research on heat transfer in porous media［J］.China Science Foundation，1991，1:32 －34.

［2］刘伟.多孔介质传热传质理论与应用［M］.科学出版社，2006.Liu Wei.Theory and application and heat and mass transfer in porous media［M］.Science Press，2006.

［3］Khaled A R A，Vafai K.The role of porous media in modeling flow and heat transfer in biological tissues［J］.International Journal of Heat and Mass Transfer，2003，46(26):4989 －5003.

［4］Kuznetsov A V，Avramenko A A.Stability analysis of bioconvection of gyrotactic motile microorganisms in a fluid saturated porous media［J］.Transport in Porous Media，2003，53(1):95 －104.

［5］Horton C W，Rogers F T.Convection current in porous medium heated from below［J］.Journal of Applied Physics，1945，6:367 －370.

［6］Wang C Y.Onset of convection in a fluid-saturated rectangular box bottom heated by constant flux［J］.Physics of Fluids，1999，11:1673 －1675.

［7］Kubitschek J P，Weidman P D.Stability of a fluid － saturated porous medium heated from below by forced convection［J］.International Journal of Heat and Mass Transfer，2003，46(9):3697 －3705.

［8］Hayat T，Javed T，Abbas Z.Slip flow and heat transfer of a second grade fluid past a stretching sheet through a porous space［J］.International Journal of Heat and Mass Transfer，2008，51(17):4528－4534.

［9］Fetecau C，Hayat T，Fetecau C.Starting solutions for oscillating motions of Oldroyd-B fluids in cylindrical domains［J］.Journal of Non-Newtonian Fluid Mechanics，2006，153(2):191 －201.

［10］Kim MC，et al.Thermal instability of viscoelastic fluids in porous media［J］.International Journal of Heat and Mass Transfer，2003，46(26):5065－5072.