二维 Rayleigh-Bénard 对流的插值格子-Boltzmann 方.doc

精品论文推荐二维 rayleigh-bnard 对流的插值格子-boltzmann 方法模拟研究1王勇，何雅玲，童长青，刘迎文 西安交通大学动力工程多相流国家重点实验室，西安 (710049) e-mail: 摘要：本文采用插值格子-boltzmann 方法对较大 ra 数范围下的二维 rayleigh-bnard 对流进行了模拟研究。模拟中针对不同的 ra 数，采用不同的插值比。模拟获得了系统的最大垂 直方向速度分量随时间的变化规律、系统的流线以及等温线分布、nu 数与 ra 数的变化规律以及壁面附近水平截面平均温度的分布。模拟结果与相关文献数据做了对比，相互吻合良好。 关键词：格子-boltzmann 方法； 插值； rayleigh-bnard 对流中图分类号：tk124文献标识码:a1. 引言格子-boltzmann 方法（lbm）是由格子气自动机（lga）发展而来的一种新的流体数 值模拟方法，被认为是一种介观的动理论格式。与传统计算流体力学（cfd）及数值传热学（nht）方法相比，lbm 具有许多独特的优点，并因此吸引了众多领域学者的研究兴趣1-3。 从理论上讲，lbm可以模拟从介观到宏观的一系列流动与传热问题。但对于宏观大尺度下（如大re数、大ra数）的流动与换热问题的模拟，因需要大量的计算网格，由此，对 计算资源提出了很高的要求。不少学者致力于这方面的研究，例如清华大学杨帆等在lbm 中引入smagorinsky亚格子概念4，并通过大re数下二维顶盖驱动流的模拟，证实了方法的 有效性。本文采用一种基于差值的 lbm 模型5，以二维 rayleigh-bnard 对流为例，证实其对大ra 数问题模拟的可靠性。在较大的 ra 数范围内，本文模拟结果与相关文献结果吻合很好。2. 插值热 lbm 模拟方法的实施对于热 lbm 模型，本文采用 he 等提出的双分布函数（ddf）模型6。该模型是在密- 5 -度分布函数 fi的基础上引入了内能密度分布函数 gi ，并用 fi演化得到速度场，用 gi 演化得到温度场。演化方程、平衡态分布函数及宏观物理量的定义详见文献6。插值 lbm 的实质是在常规的 lbm 基础上引入适体网格的概念。下面以二维物理问题 为例，简述插值 lbm 的实施过程2, 3, 5：（1）在物理平面上选区适当的曲线坐标系 ( ,)，使得 ( ,)上网格为正方形。取坐标变换 x = x( ,)， y = y( ,)得到网格节点的直角坐标 r( ,)。 ( x, y)为计算平面。对于本文所研究的物理问题，则 ( x, y) = ( , ) r ， r 为插值比；（2）给定网格节点上的宏观量，由平衡态分布函数计算初始时刻的 fi (r, 0)、gi (r, 0) ；（ 3 ）按照 lbm 演化 方程 在计 算平 面上 进 行 演 化 ， 得 到fi (r + ci t, t + t) 、gi (r + ci t, t + t ) ；（4）坐标变换 x = x( , )、y = y( , )对应逆变换 x = ( ,)、 y = ( ,)。在 ( ,)1 本课题得到国家杰出青年科学基金（50425620），高等学校博士学科点专项科研基金（20050698036）资 助，教育部重点项目（306014）。i i i平面上寻找 r + c t 相对应的坐标 ( x , y )，以及最接近的网格节点 (x , y ) 。取id = x x ， d = y y ，按如下公式插值：i2 2x yfi (r, t + t ) = am, k bn,l fi (m + k md ,n +l nd , t + t )（1）k =0 l =0其中， md = sign(1, d ) ， nd = sign(1, d ) 。本文采用二次抛物线插值，详见文献2；（5）反复循环第三、四步，直到程序收敛，并输出物理平面上的结果。3. 数值模拟及结果分析3.1 模拟细节rayleigh-bnard 对流的物理模型如图 1 所示，流体位于一 x 方向上无限长的狭窄通道内，3通道上壁面温度t1 ,下壁面温度t0（t0 t1 ），通道高 ly 。定义 ra 数为 ra = (t0 t1 ) gly / a ，其中 为热膨胀系数， g 为重力加速度， 和 a 分别表示动力粘度系数和热扩散系数。线性 稳定性分析表明，在 boussinesq 假设下的流体运动中，当流动 ra 数大于临界值 rac = 1707.76 时，流动将丧失稳定性，换热机制由导热转变为对流换热1, 6。 lx t1yt0g lyx图 1 物理模型模拟中，取 lx ly = 2ly ly =(rn x ) (rn y ) ， n x 、 n y 分别为 x、y 方向上的网格数。x方向上采用周期性边界条件，y 方向上采用非平衡外推方法。取 pr = 0.71 ，g = g (t tm ) j ，tm = (t0 + t1 ) 2 。此外，定义特征速度 uc =gly ，特征时间 tc = lyvc =ly g 。在下面的模拟中，网格数固定为 n x n y = 200 100 ，通过改变插值比 r 的大小来改变模拟区域的大 小，从而达到调节 ra 数的目的。3.2 模拟结果及分析初始时刻，模拟区域内均给定一固定的温度和压强的小扰动6。此外，当 ra 10000 时， r 取为 1，此时插值 lbm 退化为标准的 lbm；当 ra 10000 时，r 随 ra 数的增大而增大。 具体地，当 ra = 10000 、20000 时，r 取为 2；当 ra = 50000 时，r 取为 4；当 ra = 100000 时， r 取为 8。图 2 给出了不同 ra 数下，模拟区域内的最大垂直方向速度分量 uy ,max 随时间的变化规律。由图可知，当 ra 数小于 rac 时（ ra = 1000 ），黏性力和扩散作用抑制了系统内的 小扰动，uy ,max 保持零不变，系统处于静止状态；当 ra 数大于 rac 时（ ra = 2000100000 ）， 初始时刻的扰动得到放大，uy ,max 迅速增大并最终保持为一恒定值。当 uy ,max 保持不变时，系 统达到稳定状态。此外，比较不同 ra 数可发现， ra 数越大，uy ,max 最终达到的恒定值也越 大，即系统内的对流越强烈。图 2 最大垂直方向速度分量随时间的变化图 3 为不同 ra 数下系统稳定后的流线、等温线分布情况。由图 3（a）可见，当 ra 数 超过 rac 时，上下壁面温差引起的浮升力的作用超过了黏性力和扩散作用，系统丧失稳定性， 出现对流流动，y 方向上的温度不再呈线性分布，开始有涡的产生，此时涡的形状近似圆形。 随着 ra 数增加（图 3（b）、（c），壁面附近的温度变化剧烈，梯度逐渐增大；等温线 分布变得复杂，局部区域出现冷热流体混合；涡的形状也逐渐变为扁平。流线等温线(a)ra =2000(b)ra =20000(c)ra =100000图 3 不同 ra 数下的流线、等温线分布当系统达到稳定后，统计其 nu 数，以作为程序正确性的定量考核。 nu 数定义为6nu = 1 + lya(t0 t1 )（2）式中 表示系统平均。图 4 给出了本文获得的 nu 数随 ra 数的变化规律，并与文献 7采用 galerkin 方法以及文献8采用有限差分法所获得的结果做了对比。由图可见，当 ra 数较小（ ra 2000 ）时，本文结果介于文献7与8所获得的数据之间。总的说来，本文模 拟结果与文献吻合良好。lbm : r =1 5 lbm : r =2 lbm : r =4 4 lbm : r =8 nuclever & busse73 from m 8211000 10000100000ra 图 4 nu 数随 ra 数的变化进一步，我们统计了流体的水平截面平均温度ta ( y) = t ( x, y)dx dx 。由图 5 可见，在下壁面附近，流体的温度梯度随 ra 数的增大而增大，也即边界层厚度随着 ra 数的增大 而减小。因此，为了更详尽地分析壁面附近的物理场景，需要在壁面附近划分较细的网格。 此外，当 ra 数继续增大时，流体的运动呈现出三维特性，并逐渐向湍流过渡。这些都是本 文下一步将要开展的工作。0.5(t -t )/(t -t )a m0 ra=3000ra=5000ra=10000ra=20000ra=500000.10.00.0 0.1 0.2 0.3 0.4 0.5y/l4. 结论y图 5 不同 ra 数下水平截面平均温度随空间位置的变化本文采用一种插值格子-boltzmann 方法对较大 ra 数范围下的二维 rayleigh-bnard 对流 进行了模拟研究。模拟中针对不同的 ra 数，采用不同的插值比。模拟获得了系统最大垂直 方向速度分量随时间的变化规律，系统的流线以及等温线分布，nu 数与 ra 数的关系以及壁 面附近水平截面平均温度的分布。模拟结果与相关文献数据做了对比，相互吻合良好。由此 证明了该方法的有效性，从而为下一步将 lbm 应用于诸如热声制冷机复杂系统等宏观问题 的数值研究奠定了基础。参考文献1shan xiaowen. simulation of rayleigh-benard convection using a lattice boltzmann method. phys. rev. e,1997, 35(3): 2780-27882李明秀，陶文铨，何雅玲等. 格子-boltzmann 方法非均匀网格的实施. 工程热物理学报, 2003, 24(1):73-753程永光. 基于插值的 lattice boltzmann 方法非均匀网格算法. 武汉水利电力大学学报, 2000, 33(5):26-314杨帆，刘树红，唐学林，吴玉林. 格子 boltzmann 亚格子模型的研究. 工程热物理学报, 25(增刊): 43-465he xiaoyi, luo lishi and dembo m. some progress in lattice boltzmann method. part i. nonuniform meshgrids. j. comput. phys., 1996, 129(2), 357-3636he xiaoyi, chen shiyi and doolen d g. a novel thermal model for the lattice boltzmann method inincompressible limit. j. comput. phys., 1998, 146(1): 282-3007clever r m and busse f h. transition to time-dependent convection. j. fluid mech., 1974, 65: 625-6458fromm j e. numerical solutions of the nonlinear equations for a heated fluid layer. phys. fluids, 1965, 8(10):1757-1769numerical studied of two-dimensional rayleigh-bnardconvection using lattice-boltzmann methodbased on interpolationwang yong, he yaling, tong changqing, liu yingwenstate key laboratory of multiphase flow in power engineering, xian jiaotong university, xian (710049)abstractnumerical studies are presented for two-dimensional rayleigh-bnard convection in a large scope ofra number using a lattice boltzmann method based on interpolation. in the simulation, various ratios of interpolation are adopted for corresponding ra numbers. numerical results such as the maximal time dependent vertical velocities of the system, steady-state streamlines and isotherms