不可压热流体的离散统一气体动理学简化算法

1、不可压热流体的离散统一气体动理学简化算法摘要：离散统一气体动理学算法是一种适用于连续流和稀薄流的统一数值算法，具有良好的数值 稳定性和较高的数值精度，但算法结构较为复杂，计算效率有待提升&针对不可压缩热流体流动, 在原始的离散统一气体动理学算法的基础上，使用梯形法则代替中点法则求解对流项，提出一个 简化的离散气体动理学算法，用二维自然对流进行数值模拟&理论分析和计算结果表明，对于二 维热对流问题的求解，和原始算法相比，改进算法在保障稳定性和数值精度下，计算效率提升了 30%左右&关键词：离散气体动理学算法;介观流体力学方法热流动数值模拟。引言流动与热交换广泛存在于自然界及工程领域，表现形式多种

2、多样&从自然界中的大气循环、水循 环、风霜雨雪的形成，到工业中微电子的散热冷却、航空航天飞行器的航行活动，都与流动和传热过程密 切相关&传统的热流体力学研究从宏观角度出发，通过求解Navier-Stokes方程来研究热流动问题，逐 渐发展成一套完善的解决热流过程问题的可靠方法1-2。但是，宏观方法难以克服连续性假设的制约， 在稀薄流及微流动领域的局限性日益突出&以格子Boltzmann方法(Lattice Boltzmann Method，LBM) 为代表的介观流体力学方法从分子底层分析流动与传热机理，为热流体的数值研究提供了新的思 路传统的LBM在可压缩及稀薄流领域的研究也遇到难以突破的瓶颈

3、，拓展一种广泛适用于连续 流和稀薄流的统一数值方法仍是一个尚未彻底解决的科学问题& Guo等)*提出的离散统一气体动理学 方法！ Discrete Unified Gas Kinetic Scheme，DUGKS)D发展跨越微观-宏观尺度的统一计算流体力学 方法提供了新的可行性方案& DUGKS是一种基于有限体积的新方法，结合LBM和统一气体动理学 (Gas Kinetic Scheme，GKS)方法的优点，成功应用于连续流和稀薄流的相关研究，并不断被发展完 善句。Wang等)5*拓展了热流耦合的DUGKS方法，通过动力学边界条件分别独立求解流场内流体的 速度和温度，成功模拟热流动现象& Ya

4、ng等),*使用DUGKS求解器，将相场方法应用于两相流的模拟 中，有效地模拟了界面严重变形的几种情况，并且精准捕获了许多微妙的细节& Wu等逆开发了一个守 恒的DUGKS方法用于解决流场多尺度问题，在非结构化粒子速度空间中，宏观量由介观气体分布函数 通量来更新，大大提高了效率&与LBM相比，DUGKS适用于任意努森数的流动，网格适用灵活，其数 值精度和算法稳定性均有所提升10-11，但在一个时间步长内需要演化多个分布函数，算法结构和计算 效率需要进一步优化和提升&本文针对不可压缩热流体，通过优化演化过程中分布函数微通量的计算 方法，提出一个简化的DUGKS，在不改变数值精度和稳定性的前提下，

5、大幅提高了计算效率&1数值计算方法本文采用双分布Boltzmann方程描述热流体的输运过程，碰撞项采用Bhatnagar-Gross-Krook(BGK)碰撞模型:(1)竺1 ：fa-fe(1)式中，fo为流体分布函数，当0=1,2时，:分别表示流场分布函数和温度分布函数(为流场内流体粒 子速度为松弛时间，feq为Maxwellian平衡态分布函数:# p/ _ | ( | 2 1(2)一| (2(RT)D/2 exp(2RT ) (2)T/ _ ( |2=9(2(RT)D(2 exp(2RT ) 0 2式中,R为气体常数,D为空间维度,p为宏观量密度,为流体速度,T为温度)通过分布函数得到宏

6、观 量密度、速度以及温度:(3)p ifl* ,4 = 7：i ,T (3)方程(1)的求解采用类似于Richtmyer格式%12&的两步算法，网格节点分别采用实心与空心点表示, f=,fc分别为网格格点和网格中心点上的流体分布函数，其结构如图1所示。格点 中心& 边界好首先,求解半个时间步长后格点q(.处流体分布函数,碰撞项采用梯形法则，将方程(1)格点 中心& 边界好首先,求解半个时间步长后格点q(.处流体分布函数,碰撞项采用梯形法则，将方程(1)离散为:丁)= + =()=.()n /J a n /? |【a /(4)其中=+2 rn+4 n+2J a ,ijJ a ,ij& #a ,i

7、j(=) = (7 :) / : = (xh)dS式(6)中，：为控制体表面的外法向量，S为控制体面积,微通量F *的计算也采用梯形法则，以通过 的控制体右侧Jn的通量为例，其计算过程如下:M L ,n.=+2 打)M+ =( f a H (X+1 h+1 ) + f a , H (X+1 H ) /_J (:) j a, ij ( j ) os72(,+1 H+1+1 h然后，利用分布函数f=，2在格点处的值，采用梯形法则求出通过也的控制体在n + 1/2时刻的微通量:F：，1 (8)= J ()： (QdSSV*进一步使用如下离散格式求得下一时刻的分布函数：(5)(6)(7)(8)$ n+

8、1/、_ $ + ,= /、厂 n+1 /J a ,j (c)= J a j (8) j V O a, j (8)(9)(10) TOC o 1-5 h z $=+1 _c=C 1 Zc=+1 $ + , = rnCn(10)：a, H J a, H2 Xa，H ! ： a，H J a, H2 H在热流问题的求解中，需要考虑传热与流动的耦合，根据Boussmesq假设,流体密度与温度呈线性关系： pbpo*po3(T*To)(11)式中,4o,To分别为流场中的平均密度和平均温度,3为热扩散系数，记重力加速度为g0。在Boussmesq 假设下，重力与外力加速度表示为：(12)#G=4go =

9、 &40*403(TTo)g。 a = go3(T*To(12)(13)在演化过程中，外力项可以合并到碰撞项中，于是将流场分布函数控制方程中的碰撞项修正为： 。三 *1 J*J1*Ji *1 &/-/eq +；(* 碰撞项的修正不改变演化过程，算法将继续演化流场内分布函数，(13)简化DUGKS算法的计算步骤如下：初始化控制中心处的流场宏观量，计算平衡态Jq,并用其初始化控制中心处的分布函数J” (采用插值函数计算得到格点处的分布函数fa, H，根据分布函数之间对应关系* + 一 2, a 2 $+32 reqJ a，H = 2,a + max +1) ej max + m ax e (jma

10、x + 1)13 max ej max )(14)本文采用梯形法则计算微通量，先后演化控制中心兀与格点如 不再需要计算边界中心处的流场 及温度分布函数。一个时间步长的演化次数约为：(15)NNSDUGKS max X j max + max + 1 ) X (jmax + 1)12 max ej max) 理论上，本文算法的计算量相较于原始算法大约下降了 1/3(15)Ra = 106时，分别采用本文算法DUGKS算法进行二维自然对流模拟，得到方腔竖直中线处水平 速度、方腔水平中线处的垂直速度以及温度分布如图6所示#从图6可以看出，在流场和温度场的模拟 中，本文算法与DUGKS算法取得的结果非

11、常吻合0.030.020.010.00-0.01-0.020.10、-K!Mffi8 6 4 2 0 2 4 .OOOOOOO 0.030.020.010.00-0.01-0.020.10、-K!Mffi8 6 4 2 0 2 4 .OOOOOOO O.O.O.1.0X(c)水平方向上的温度T分布图6 水平中线处速度分布和温度分布比较3结束语f/min图3结束语f/min图7不同算法的演化逮度对比为了比较2种算法的演化速度，设置方腔内水平速度残差E为终止条件,(16)E=2 (16)(兀+1 000) | 2 ； II（C=+1(兀+1 000) | 2 ； II式中，诚=（k，G =）为方腔内的水平速度$在本文计算中将E取值为106 $本文算法%DUGKS算法达到相同终止条件 时的演化速度如图7所示$从图7可以看出，本 文算法的收敛速度明显快于DUGKS算法，2种 算法达到终止条件所需时间分别为26. 80 min和 40. 20 min$由此可见，本文算法在不改变DUGKS 算法稳定性的前提下，提高