0 引言

有限空间内流体流动的非线性特性研究有重要的理论价值，正确认识流动由稳态层流经过一系列分岔过渡到混沌流的物理机制对认识湍流发生机理有所帮助.国内外很多学者对此问题进行了研究.Angirasa^[1]对一侧为恒温壁面, 具有单一送风口和出风口的二维方腔进行了数值研究，分析了不同理查德森数Ri对腔内换热强度的影响，并提出了在较大Ri时无法获得稳态解.Prasad等^[2]采用有限体积法研究了方腔的宽深比、格拉晓夫数Gr对混合对流的影响，发现在宽深比为2.0，Gr=10⁵时，总动能出现了霍夫型分岔.赵明等^[3]采用有限差分法对顶部设置4个送风口的具有对称结构的方腔内混合对流换热问题进行了数值研究, 计算得到了非对称的解.Kang等^[4]对瑞利数Ra=10⁷时，内置发热圆的封闭方腔内自然对流的分岔现象进行了数值模拟研究，结果发现发热圆在方腔水平中线和对角线上变化位置时，流动在稳态和非稳态之间转变.杨茉^[5]等采用有限差分法对圆内开缝圆非定常自然对流换热进行了数值研究，发现当瑞利数Ra较大时流动开始出现振荡.

上述研究为采用传统的数值方法得到的结果.近年来发展的格子Boltzmann方法由于其物理意义清晰、边界处理容易、演化过程清晰、程序易于实现、并行性好等优点已经发展成为流体流动模拟和物理建模的一种独特和可信的数值方法，不少学者采用该方法在非线性理论研究方面取得了一些成果.Anupindi等^[6]用格子Boltzmann方法研究三维顶盖驱动的腔体内振荡的不稳定性，得出第一霍夫分岔雷诺数与方腔深度负相关.Li等^[7]采用多松弛格子Boltzmann模型模拟不同宽深比下顶盖驱动方腔内的振荡情况，发现随着方腔深度的增加，第一霍夫分岔雷诺数减小.卞恩杰等^[8]用格子Boltzmann方法对二维Rayleigh-Benard流进行模拟与非线性分析，在Ra=10⁶时得到了振荡解，预测出了动态分岔.

由上述可知，非线性问题的研究主要集中于自然对流、顶盖驱动流和单介质的混合对流，而对于双扩散混合对流振荡特性的研究很少.本文研究的内置发热圆的方腔内双扩散混合对流问题，是典型的污染物处理、热舒适性研究以及烟气控制等理论模型，本文采用格子Boltzmann方法探讨有限空间内双扩散混合对流的非线性问题，研究理查德森数Ri和出口位置对流动稳定性的影响.

1 计算模型

本文的双扩散混合对流的物理模型如图 1所示，A为监测点，方腔内置一发热圆，发热圆位于方腔中心.方腔的长宽均为L，发热圆的直径为d(d=0.4 L).方腔设有出入口，出入口的大小均为0.1L.流体从左侧底部流入，可分别从顶部右侧位置a、中间位置b和左侧位置c流出，流体入口速度为u_i，温度为T_i，浓度为C_i，发热圆表面的温度为T_h，浓度为C_h.方腔内其余壁面绝热且浓度梯度为零，重力加速度为g，其中T_h>T_i，C_h>C_i.方腔内流体的普朗特数Pr为0.7.

假设流体不可压且不考虑粘性耗散，采用Boussinesq近似^[9]，即假设除与体积力有关的项外，流体的密度ρ作为常数，并且认为体积力项中的密度

$ \rho {\rm{ = }}{\rho _0}\left[ {1 - {\beta _{\rm{T}}}\left( {T - {T_0}} \right) - {\beta _{\rm{C}}}\left( {C - {C_0}} \right)} \right], $

(1)

其中ρ₀、T₀和C₀分别是参考密度、温度和浓度，β_T为温度引起的体积膨胀系数，K^-1: ${\beta _{\rm{T}}} = - \frac{1}{{{\rho _0}}}\left( {\frac{{\partial \rho }}{{\partial T}}} \right)$, β_C为浓度引起的体积膨胀系数，m³·kg^-1：${\beta _{\rm{C}}} = \frac{1}{{{\rho _0}}}\left( {\frac{{\partial \rho }}{{\partial C}}} \right)$.

当温度和浓度变化不大时，可以假设流体的物性(粘性系数υ、热扩散系数α和质量扩散系数D)为常数.在这些假设条件下，双扩散混合对流的控制方程为

$ \nabla \cdot \mathit{\boldsymbol{u}} = 0, $

(2)

$ \frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \nabla \cdot \left( {\mathit{\boldsymbol{uu}}} \right) = - \frac{1}{{{\rho _0}}}\nabla p + \upsilon {\nabla ^2}\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{F}}, $

(3)

$ \frac{{\partial T}}{{\partial t}} + \nabla \cdot \left( {\mathit{\boldsymbol{u}}T} \right) = \alpha {\nabla ^2}T, $

(4)

$ \frac{{\partial C}}{{\partial t}} + \nabla \cdot \left( {\mathit{\boldsymbol{u}}C} \right) = D{\nabla ^2}C, $

(5)

其中F=gβ_T(T-T₀)j+gβ_C(C-C₀)j为浮力项，g为重力加速度常数, j为y方向的单位矢量，u=(u, v)^T为速度矢量，t、P、T和C分别为时间、压力、温度和浓度.双扩散混合对流流动的特征可以用5个无量纲数描述：普朗克数Pr=υ/α，它反映流体的动量传递与热量传递的相对大小；雷诺数Re=u_iL/υ，它表征由入口速度引起的强制对流的强度；理查德森数Ri=gβ_T(T_h-T_i)L/u_i²，它表示由浮升力引起的自然对流相对于由入口速度引起的强制对流的重要程度；路易斯数Le=α/D，它反映热量传递与质量传递的相对大小；浮升力比Br=β_c(C_h-C_i)/β_T(T_h-T_i)，它表示由浓度引起的浮升力与由温度引起的浮升力的相对大小.

利用Guo等^[10]提出的耦合格子BGK模型，并在此基础上进一步耦合了浓度分布函数的模型.该模型与Guo提出的模型^[11]有所不同，虽然也是用三个独立的LBGK方程分别模拟速度场、温度场和浓度场，并采用Boussinesq假设将这三个方程耦合，但是本文的温度和浓度分布函数采用的是9个离散速度，较已有模型^[11]中采用的4个离散速度更能保证复杂边界条件下数值模拟的精度.

采用三个独立的分布函数模拟速度场、温度场和浓度场，其演化方程分别为

$ {f_i}\left( {\mathit{\boldsymbol{r}} + {\mathit{\boldsymbol{c}}_i}\Delta t,t + \Delta t} \right) - {f_i}\left( {\mathit{\boldsymbol{r}},t} \right) = - \frac{1}{{{\tau _f}}}\left[ {{f_i}\left( {\mathit{\boldsymbol{r}},t} \right) - f_i^{{\rm{eq}}}\left( {\mathit{\boldsymbol{r}},t} \right)} \right] + \Delta t{F_i}, $

(6)

$ {T_i}\left( {\mathit{\boldsymbol{r}} + {\mathit{\boldsymbol{c}}_i}\Delta t,t + \Delta t} \right) - {T_i}\left( {\mathit{\boldsymbol{r}},t} \right) = - \frac{1}{{{\tau _{\rm{T}}}}}\left[ {{T_i}\left( {\mathit{\boldsymbol{r}},t} \right) - T_i^{{\rm{eq}}}\left( {\mathit{\boldsymbol{r}},t} \right)} \right], $

(7)

$ {C_i}\left( {\mathit{\boldsymbol{r}} + {\mathit{\boldsymbol{c}}_i}\Delta t,t + \Delta t} \right) - {C_i}\left( {\mathit{\boldsymbol{r}},t} \right) = - \frac{1}{{{\tau _{\rm{C}}}}}\left[ {{C_i}\left( {\mathit{\boldsymbol{r}},t} \right) - C_i^{{\rm{eq}}}\left( {\mathit{\boldsymbol{r}},t} \right)} \right], $

(8)

$ f_i^{{\rm{eq}}} = \left\{ {\begin{array}{*{20}{l}} {{\rho _0} - \frac{5}{3}\frac{\rho }{{{c^2}}}{\rho _0}{s_i}\left( u \right),}&{i = 0;}\\ {\frac{1}{3}\frac{\rho }{{{c^2}}}{\rho _0}{s_i}\left( u \right),}&{i = 1,2,3,4;}\\ {\frac{1}{{12}}\frac{\rho }{{{c^2}}}{\rho _0}{s_i}\left( u \right),}&{i = 5,6,7,8;} \end{array}} \right. $

(9)

$ {s_i}\left( u \right) = {\omega _i}\left[ {3\frac{{{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}}}{{{c^2}}} + 4.5\frac{{{{\left( {{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}} \right)}^2}}}{{{c^4}}} - 1.5\frac{{{u^2}}}{{{c^2}}}} \right], $

(10)

其中权重系数ω₀=4/9, ω_i=1/9(i=1~4)和ω_i=1/36(i=5~8)，(6)式中的外力F_i为

$ {F_i} = \left( {1 - \frac{1}{{2{\tau _{\rm{f}}}}}} \right){\omega _i}\left[ {3\frac{{{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}}}{{{c^2}}} + 9\frac{{{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}}}{{{c^4}}}{\mathit{\boldsymbol{c}}_i}} \right] \cdot \mathit{\boldsymbol{F}} $

(11)

$ {\mathit{\boldsymbol{c}}_i}\left\{ {\begin{array}{*{20}{l}} {\left( {0,0} \right).}&{i = 0;}\\ {\left( {\cos \left[ {\left( {i - 1} \right)\frac{\pi }{2}} \right],\sin \left[ {\left( {i - 1} \right)\frac{\pi }{2}} \right]} \right)c,}&{i = 1,2,3,4;}\\ {\sqrt 2 \left( {\cos \left[ {\left( {2i - 1} \right)\frac{\pi }{4}} \right],\sin \left[ {\left( {2i - 1} \right)\frac{\pi }{4}} \right]} \right)c,}&{i = 5,6,7,8;} \end{array}} \right. $

(12)

$ T_i^{{\rm{eq}}} = {\omega _i}T\left[ {1 + 3\frac{{{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}}}{{{c^2}}} + 4.5\frac{{{{\left( {{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}} \right)}^2}}}{{{c^4}}} - 1.5\frac{{{u^2}}}{{{c^2}}}} \right], $

(13)

$ C_i^{{\rm{eq}}} = {\omega _i}C\left[ {1 + 3\frac{{{\mathit{\boldsymbol{c}}_i} \cdot \mathit{\boldsymbol{u}}}}{{{c^2}}}} \right]. $

(14)

流体的宏观速度、压力、温度和浓度定义为

$ u = \sum\nolimits_{i = 1}^8 {{c_i}{f_i}} ,\;\;\;\;p = \frac{5}{3}{c^2}\left[ {\sum\nolimits_{i = 1}^8 {{f_i} + {s_0}\left( u \right)} } \right],\;\;\;\;T = \sum\nolimits_{i = 1}^8 {{T_i}} ,\;\;\;\;C = \sum\nolimits_{i = 1}^8 {{C_i}} . $

(15)

流体的粘性系数、热扩散系数和质扩散系数分别为

$ \nu = \frac{1}{3}\left( {{\tau _{\rm{f}}} - \frac{1}{2}} \right){c^2}\Delta t,\;\;\;\;\alpha = \frac{1}{3}\left( {{\tau _{\rm{T}}} - \frac{1}{2}} \right){c^2}\Delta t\;\;\;\;D = \frac{1}{3}\left( {{\tau _{\rm{C}}} - \frac{1}{2}} \right){c^2}\Delta t. $

(16)

这里的i=0, 1, 2, …, 8.其中f_i(r, t)表示密度分布函数，T_i(r, t)表示温度分布函数，C_i(r, t)表示浓度分布函数，f_i^eq(r, t), T_i^eq(r, t)和C_i^eq(r, t)分别为速度、温度和浓度的平衡态分布函数，τ_f, τ_T和τ_c分别为其相应的松弛时间.

发热圆表面的平均努赛尔数Nu_av和平均舍伍德数Sh_av由下式确定

$ \begin{array}{l} N{u_{{\rm{av}}}} = - \frac{1}{{2\pi }}{\int_0^{2\pi } {\left( {\frac{{\partial T}}{{\partial \mathit{\boldsymbol{n}}}}} \right)} _{圆表面}}{\rm{d}}\varphi ,\\ S{h_{{\rm{av}}}} = - \frac{1}{{2\pi }}{\int_0^{2\pi } {\left( {\frac{{\partial C}}{{\partial \mathit{\boldsymbol{n}}}}} \right)} _{圆表面}}{\rm{d}}\varphi , \end{array} $

(17)

其中n为圆表面的法线方向，φ为圆角弧度值, Nu_av表征发热圆表面与方腔内流体之间的对流换热强度，Sh_av表征发热圆表面与方腔内流体之间的质扩散能力大小.

发热圆曲面边界处理格式采用Bouzidi等^[12]提出的插值格式，基本思想是在弹回方向上用空间插值来获得由边界进入流体的分布函数.其余边界的处理格式采用Guo等^[13]提出的具有二阶数值精度的非平衡态外推格式，其基本思想是将边界节点上的分布函数分解为平衡态和非平衡态两部分，平衡态部分由边界条件获得，而非平衡态部分则用非平衡外推确定.

模型中网格数采用200×200, 速度和温度及浓度最大绝对残差均小于10^-6.

2 程序验证

为了验证本文所采用的格子Boltzmann方法的准确性，本文与传统的有限容积法的数值结果^[14]进行对比，等温度线图、流线图、等浓度线图对比结果如图 2所示.

数值方法得出的平均努赛尔数Nu_av和平均舍伍德数Sh_av与文献[14]的结果对比见表 1, 由表可知，误差均小于1%.同时，为了验证本文采用的数值方法对于研究非线性问题的准确性，本文与文献[15]采用有限容积法的模拟结果进行对比，结果如图 3所示.从对比结果可知，本文研究非线性问题的方法是可信的.

3 结果分析

研究流体振荡现象的方法有时间历程分析法、相空间轨迹法和功率谱法^[16].本文采用时间历程分析法和相空间轨迹法，分析不同Ri和出口位置下的监测点的速度随时间的变化以及监测点速度相图的变化规律.本文数值模拟中Le=10.0，Br=1.0.

图 4、5、6给出了当Ri分别等于0.001, 0.005和1.0时，不同出口位置下监测点的速度随时间的变化.从图可知，在Ri=0.005和1.0时，三种出口位置情况下的流体流动随时间变化很快趋于稳定，速度值瞬间由零增大，之后快速收敛，达到稳定状态，保持不变.然而在Ri=0.001, 出口位置为a时(见图 4(a)), 速度达到一定值后，出现振荡情况.这是因为随着Ri的减小，强制对流将逐渐占据主导地位，流动会变得更加剧烈，从而导致流动状态从稳态转变为非稳态.并由图 4(b)Ri=0.001时的速度随时间变化曲线局部放大图可知，该振荡呈现周期性变化规律.而出口位置为b和c的方腔内监测点在Ri=0.001时，流动出现振荡现象，但从图中已分辨不出其周期性，呈现无规律振荡, 这是因为出口位置设置在方腔顶部中间位置b和顶部左侧c时，与出口位置为a的情况相比，更有利于空间内部流体的流动，使得在相同的Ri下方腔内流体流动更加剧烈，流动呈非稳态非周期性振荡.

为了更好地确定方腔内双扩散混合对流的振荡形态，图 7、8、9分别给出了Ri等于0.001, 0.005和1.0时，不同出口位置的方腔内监测点速度相空间轨迹图.当出口位置为方腔顶部右侧a，Ri为0.001时，图 7(a)显示了非稳态周期性振荡解在相空间中的轨迹，其形态为一个封闭环，其吸引子为极限环吸引子，它将周围的轨道吸引到这个周期性的循环之中，反映了流动随时间作周期性运动.而在Ri=0.005和1.0时，图 7(b)、7(c)给出了稳定流场的解在相空间中的轨迹，其吸引子为定常吸引子，其形态最终表现为一个定点.当出口位置为方腔顶部中间b和顶部左侧c时, 在Ri=0.001时, 图 8(a)和9(a)给出了非周期振荡解在相空间中的轨迹，形态表现没有规律，流动呈现混沌流特征.而在Ri=0.005和1.0时，定常解在相空间中的轨迹稳定后最终为一个点，它将周围的轨道吸引过来.

4 结论

通过格子Boltzmann方法对方腔内双扩散混合对流的非线性特性进行研究，随着理查德森数Ri的变化，在不同的出口位置下，方程的解呈现出多样的、不稳定的非线性特征，归纳为三种结果，即定常解、周期性振荡解和非周期性振荡解.在Ri=0.005和Ri=1.0时三种出口位置下均为定常解，表现为流体流动趋于稳定，相轨迹最终为一个点.在Ri=0.001, 出口位于方腔顶部右侧时，为周期性振荡解，表现为流动参数随时间呈周期性变化，速度相轨迹归为一个封闭环.在Ri=0.001，出口位于方腔顶部中间和左侧时，为非周期性振荡解，表现为流动参数随时间做无规律的变化，相轨迹毫无规则，流动进入混沌区.