封闭方腔自然对流的格子-Boltzmann方法动态模拟

4.504 4.519 4.510 4.510 0.199%
8.767 8.800 8.806 8.805 0.056%
从表 1 中可以发现，采用本文所介绍的不可压缩双分布函数 TLBM 模型进行数值计算，得 到了比较精确的结果。相对误差
5. 方腔内自然对流的动态模拟
封闭方腔自然对流是热流耦合的经典问题，通过对其进行数值模拟而获得不同 Ra 情况
2. 物理模型
本文所计算的封闭方腔自然对流的物理模型如图 1 所示。封闭方腔高为 H ，上、下壁
1
本课题得到国家杰出青年科学基金资助项目（50425620）及高等学校博士学科点专项科研基金资助项目 -1（20050698036）资助。
Th + Tc ⎞ 面绝热，腔内充满 ρ = 3 ， Pr = 0.71 ，温度 T = ⎛ ⎜ ⎟ 的均质 ⎝ 2 ⎠
p i x + ei dt , t + dt − p x, t = −
(
) ( )
dtτ p Fi dt p i − pieq + τ p + 0.5dt τ p + 0.5dt
(
)
（6）
g i x + ei dt , t + dt − g x, t = −
(
) ( )
p dt dt g i − gieq − Z i 2i τ g + 0.5dt τ g + 0.5dt cs
(
)
（7）
图 2. D2Q9 模型
。 其中 τ p ，τ g 分别为运动和热方程的松弛时间； cs 为声速（ cs = 1/ 3 ） 流体的宏观参量（包括压力，速度，温度及热流等）可按下列各式计算：
p x, t = ∑ pi
i
( )
p0 u = ∑ p i eix +
i
ρ Gcs2δ t
2
⎞ τg ⎟ ⎟ τ + 0.5dt ⎠ g
（a） t =400
（b） t = 2.5 × 103
（c） t = 104
（d） t = 3 × 105
（e） 等压线分布
（f） 等温线分布
图 3. Ra = 103 时方腔内自然对流随时间的动态演化过程（a,b,c,d）及稳定状态下的压强、温度分布
对流的形成使流体在左上方和右下方积蓄而形成高压中心， 与此相反， 在方腔的左下方和右 上方由于流体不断流出而成为低压中心。 于是在水平方向产生了由压差导致的强制流动。 同 时，左上方和左下方之间也存在压差，加之粒子与壁面碰撞产生的涡量，使得在竖壁自然对 流层的内侧由于流体回流而形成旋涡；出于相同原因，图 3b 示出了腔内右侧所形成的一个 相同方向的旋涡。我们把这两个涡称为原始涡。由图 3c 可见，随着时间的步进，压差逐渐 增大，同时旋涡强度亦逐渐增强，某时刻，左、右涡相遇，由于各自涡量的方向相同，两涡 合并，形成一个合成涡。合并初期由于两原始涡的涡核相距较远，故而合成涡呈椭圆形。随 着流场的均匀化，原始涡核消失，合成涡逐渐变成近似圆形。图 3d 所示为充分发展以后的 流场情况。 从以上分析中可以看出， 封闭方腔内自然对流是温差、 压差和方腔有限空间限制三项因 素综合作用的结果。在 Ra 较小时，两原始涡经过一定时间的演化后，最终合成为一个近似 圆形的合成涡。 当 Ra = 104 时流场的演化情况与 Ra = 103 时基本相似。如图 4a 所示，充分发展以后流线 分布呈椭圆形，并且椭圆长轴向高压中心倾斜。这与 Ra 的增大而引起竖壁自然对流强度的 增加有关。 Ra 的增加使浮升力增强，流体速度增加，在左上方和右下方由于粒子与壁面碰 撞而导致的涡的卷起更为强烈， 因而在方腔中产生的合成涡的长轴偏离方腔中线， 而偏向两 个高压中心，使旋涡椭圆化程度加大。从而，如图 4c 所示，方腔内热交换得到增强，等温 线发生较大程度的扭曲。从图中可以看出，方腔内两个低压中心的等温线分布最密，换热最
( )
( )
2
2 ⎤ u ⎥ − 1.5 2 c ⎥ ⎥ ⎦
（5）
-2-

g
eq 5,6,7,8
⎡ ei iu ei iu ρε ⎢ 3 + 6 2 + 4.5 = 4 36 ⎢ ⎢ ⎣ c c
( )
( )
2
− 1.5
2 ⎤ u ⎥ c2 ⎥ ⎥ ⎦
其中， ei 为离散速度，如图 2 所示； u 为宏观速度； ρ ， ε （二维情况下 ε = RT ）分别为当 地密度和内能。 TLBM 的粒子演化方程为：

封闭方腔自然对流的格子-Boltzmann 方法动态模拟
童长青 何雅玲* 王勇 刘迎文
（西安交通大学能源与动力工程学院 动力工程多相流国家重点实验室，710049，西安）
Email: yalinghe@ 摘要：本文构建的以压力分布函数和内能密度分布函数为基本变量的不可压缩双分布函数热
均温气体。初始时刻方腔左、右 ， Tc （ = 1 ） ，并保 壁面分别从 T0 （ = 1 ）被加温至 Th （ = 20 ） 持恒定。 Ra 是衡量自然对流强弱的重要无量纲数，其定义为：
β g ΔTH 3 （1） υa 随着 Ra 的增大，自然对流作用增强， 并由层流转化为湍流。
Ra =
图 1. 封闭方腔自然对流物理模型
（8）
ρ e = ∑ gi −
i
pZ dt ∑ i i 2 i cs2
⎛ ei pi Z i dt q = ⎜ ∑ ei gi − ρ eu − ∑ 2 ⎜ i 2 cS i ⎝
（9）
粘性系数与热扩散系数可表示为松弛因子 τ p ， τ g 的关系式：
υ = τ p RT 其中 R 为气体常数。
a = 2τ g RT
-3-

下流动和换热情况的文献已屡见不鲜，然而对其形成和演化过程的动态模拟很少见文献报 道。本文采用所建立的不可压双分布函数 TLBM 模型对其进行了动态模拟，尝试从涡运动 和演化[12]的角度对其运动及演化进行分析。 图 3 为 Ra = 103 时流场的演化情况， 图中 t 表示粒子演化步数， 相当于非稳态问题的时间。 从图 3a 可见，流场演化初期，由于左壁壁温 Th 高于壁面附近流体温度 T ，因而壁面附近的 流体在浮升力的作用下向上运动。 同时由于右壁壁温低于壁面附近流体温度而在右壁附近产 生向下流动。这种壁面附近的流动表现为典型的竖壁自然对流[13]。从图 3e 可见，竖壁自然
（10）
4. 模型考核
下面用封闭方腔内稳态自然对流问题来考核上述所提计算模型的可靠性。 文献[5-10]分别使用不同的数值方法模拟了封闭方腔稳态自然对流问题，其中，文献[8] 所得的数值解是被引用较多的数值解。为考查模型的可靠性，分别计算了 Ra 为 103 ， 104 ，
105 ， 106 情况下壁面的平均努塞尔数（ Nu =
1. 引 言
格子 Boltzmann 方法（LBM）[1-2,6]是近十几年来得以迅速发展和逐步健全的一种新型数 值方法。它起源于格子气自动机（LGA） ，继承了 LGA 实时跟踪并记录粒子时间演化的特 点。 模拟与热相关问题的 LBM 方法称为热格子 Boltzmann 方法 （TLBM） 。 He Xiaoyi[3]于 1998 提出的双分布函数模型是稳定性较好、 适用范围较广的一种。 该模型在密度分布函数的基础 上引入独立的内能密度分布函数，从而使流动与换热分离计算，即：通过密度分布函数的演 化获得流场情况， 通过内能密度分布函数的时间演化得到温度分布情况。 这样处理不但提高 [3] 了 TLBM 方法的稳定性，而且从根本上解决了粘性耗散问题 。然而，He 在其论文中指出， 该模型同样存在可压缩效应， 为了实现更为精确地模拟不可压缩问题， 他在论文附录中提出 了一种改善平衡态分布函数的方法计算，即用压力和平均密度来计算密度平衡态分布函数， 并以此密度分布函数进行演化。本文借鉴这个思想，并参考文献[4]的不可压缩 LBM 模型，重 新构建了以压力分布函数和内能密度分布函数为基本变量的双分布函数模型。 该模型具有精 度高、概念清晰、易于实现等特点，可以应用于模拟各类不可压缩热流耦合问题。封闭方腔 自然对流是一个由于流体内部密度差而导致流动的可压缩热流耦合问题， 但是对此经典问题 所进行的理论或数值模拟[5-10]都是基于流体不可压的 Boussinesq 近似，例如，De Val Davis[8] 使用流函数-涡量差分方法得到了一个比较精确的数值解。文献[9]在此基础上，采用有限容 积法模拟了更宽 Ra 范围内的自然对流情况，并采用非均匀网格以提高计算的精确度。 文献[5-10]都发现随着 Ra 的增加，方腔内自然对流作用增强，流场和温场变得更为复 杂。文献[10]从边界层理论和量纲分析的角度，分析了一充满等温流体的封闭方腔，一侧被 突然升温至某一恒定高温， 另一侧被降温至某一恒定低温的非稳态自然对流的演变情况， 发 现封闭方腔内自然对流是热边界层和运动边界层共同作用的结果。本文从 LBM 方法可以实 时地跟踪并记录粒子运动的特点出发，利用 LBM 方法动态模拟了这一非稳态封闭方腔自然 对流问题， 并从涡运动和演化的角度对此进行分析， 发现这一非稳态问题中蕴含着丰富的涡 运动学和涡动力学内容。
3. 不可压双分布函数 TLBM 方法
对自然对流的理论和数值分析都是基于流体不可压的 Boussinesq 近似，而文献[3]的 TLBM 模型具有低马赫数下的弱可压缩效应，其可压缩效应源于所采用的运动模型[1]。为 克服可压缩效应，文献[3]附录中介绍了用压力和平均密度来计算密度平衡态分布函数，然 后以改善后的密度分布函数进行演化，从而克服模型的可压缩效应。另外，文献[4]提出的 以压力分布函数为基本演化变量的 LBM 模型，可以完全克服这一可压缩效应，但并没有讨 论过与热模型的耦合问题。为改进 TLBM 方法模拟不可压缩问题的适用性，本文在文献[3] 和文献[4]的启发下，构造了以压力分布函数和内能密度分布函数为基本演化变量的双分布 函数 TLBM 模型。该模型适用于低马赫数流动。 与文献[3]双分布函数 TLBM 相似，为了保持模型的粘性系数和 N-S 方程的粘性系数的 一致性，并保持方法显性，首先有必要引入了新的分布函数： 0.5dt dt （2） p i = pi + ( pi − pieq ) − Fi τf 2