高速可压反应流的二维简化离散玻尔兹曼模型

林传栋

(中山大学 中法核工程与技术学院, 珠海 519082)

0 引 言

化学反应流不仅广泛存在于自然界中，更在人类日常生活、工业生产、航空航天和国防建设等诸多领域发挥着越来越重要的作用[1-3]。在化学反应流中，各类流体动力学行为和化学反应过程相互耦合，具有典型的流体力学非平衡、热力学非平衡和化学非平衡效应，表现出丰富多变的物理化学现象[1-3]。对于此类非定常、非平衡的复杂物理化学系统，不能使用定常流场的解析方法或数值模型进行研究，而需要使用包含多物理场耦合和时空演化过程的物理模型。尤其是高超声速飞行器[2-4]还同时面临稀薄气体效应、化学反应效应和壁面热流效应等外界因素的影响。这些极端复杂条件对如何有效模拟和预测高超声速飞行器的运行过程提出挑战。因此，高速反应流的数值建模与模拟研究在航空航天等领域具有特别重要的实际应用背景，是进一步提升航空航天工程技术的关键问题之一。

相对于高超声速飞行器，航天飞行器再入过程跨越了更加宽泛的时空尺度。具体来讲，航天飞行器以（超）高马赫数依次经历稀薄自由分子流区、过渡流区、滑移流区以及连续流区，这一经历伴随着复杂的变密度、变温度、变速度、电离和化学反应等真实飞行环境。在各个流体区域，气体物性表现出显著差异，其绕流状态具有典型的跨尺度效应。同时，几何尺寸差异明显的外形部件也会使飞行器不同区域表现出多流区共存的多尺度混合流动现象。在高超声速再入过程中，还具有强黏性、高摩阻、高熵层、强气动加热、热力学与化学非平衡相互耦合等典型特征[5-6]。为此，需要构建准确高效、稳定可靠的跨尺度模型，用以研究分析高超声速飞行器跨流域飞行的气动力热机理，明晰复杂绕流的时空变化特点以及非平衡热化学行为表征。

然而，传统的宏观数值研究大都基于连续性假设，侧重于模拟研究系统宏观演化过程，忽略了实时伴随流体宏观行为的丰富的热力学非平衡效应，导致其应用范围不够广泛、数值仿真结果不够精确。目前，基于Boltzmann 方程的动理学模型具有微介观尺度的优势，相关方法的应用研究已经成为国际国内的热点课题，包括格子气自动机（lattice gas automaton，LGA）、格子玻尔兹曼方法（lattice Boltzmann method，LBM） 、 离 散 玻 尔 兹 曼 方 法（ discrete Boltzmann method，DBM）、离散速度方法、离散统一气体动理论方法、气体动理论统一算法和矩方法等。

作为一种元胞自动机，LGA 不仅是解方程（恢复方程）的一种全新的数值方法，更是研究各类复杂系统的一种基本工具[7]。最初，LGA 由法国的Hardy、Pomeau、De Pazzis 于1973 年首次提出，目标是能够模拟Navier-Stokes（N-S）方程组描述的流体系统，但是由于该模型旋转不变性的缺失，导致它具有高度各向异性[8]。后来，Frisch、Hasslacher、Pomeau 在1986 年优化了LGA 模型，其格子矢量具有四阶张量对称性[9]。时至今日，LGA 已经被成功应用于流体力学、化学、电磁学、热声学等领域。

LGA 的基本思想是：不同的微观行为可以导致相似的宏观表征，虚拟粒子在网格上迁移和碰撞的过程保持物理量守恒，如质量、动量、能量守恒[7]。为了提高计算效率，这种虚拟世界的微观动力学应该尽可能简单[7]。此外，LGA 可看作是LBM 和DBM 的鼻祖。广义上讲，LGA 及其后代动力学模型被认为具有许多相似优点：1）设计方案简单易行，因此便于编写代码和检验程序。2）碰撞规则具有局部性，因此非常适合大规模并行计算。3）可通过选取合适的边界条件来模拟具有复杂几何结构的系统。然而，LGA具有统计噪声等问题，而LBM 和DBM 能够有效避免该问题。

近30 年来，LBM 被成功应用到诸多复杂流体系统[10-12]，并且在早期就被拓展到化学反应流领域。早在1997 年，欧洲院士Succi 等[13]首次提出燃烧系统的LBM，其限制条件包括：化学反应极快、流场不可压、反应热对流场没有影响。2015 年，日本的Yamamoto 等[14]使用LBM 模拟了三维空隙结构系统中的燃烧现象。2018 年，德国Hosseini 等的课题组[15]构建了多组分系统的LBM，并解决了扩散方程误差引起的质量不守恒问题。同时，国内学者也在相关领域做出了突出贡献。华中科技大学的陈胜等[16]把LBM 应用到二维和三维低马赫数燃烧系统，将流体密度和温度场直接耦合。西安交通大学的陶文铨院士、何雅玲院士等的课题组[17-18]使用LBM 研究包含多相、多组分、化学反应的能源转换系统。另外，国内外还有许多其他优秀科研团队将LBM 应用到包含溶解、沉淀、流固耦合等复杂因素的反应流系统[19-21]。

尽管传统LBM 的应用研究已经取得了巨大进步，但基本都是还原传统宏观控制方程的解法器，基本没有聚焦于详细的热力学非平衡信息，也没有实现真正意义上的化学反应和可压缩流体的耦合。作为一种新型的微介观动理学方法，近几年发展起来的DBM 能够解决物理精度与计算效率之间的矛盾。DBM 是玻尔兹曼方程在速度空间的特殊离散化形式，可以看作传统LBM 的进一步发展或变异[22-36]。值得说明的是，DBM 继承了玻尔兹曼方程描述非平衡物理系统的主要功能。可以证明，一方面DBM 在连续性极限条件下能够恢复宏观流体控制方程组；另一方面DBM 还包含宏观流动方程组之外的热力学非平衡效应[22-36]。因此，DBM 可用于处理宏观模型失效、非平衡效应较为显著的微尺度物理系统、稀薄气体系统和滑移层流体区域等相关问题[31]。

目前，DBM 在化学反应流方面已经得到初步应用。2013 年，DBM 已被初步发展到爆轰领域，并被用于模拟分析爆轰波附近的非平衡效应[22]。2014 年，林传栋等[23]构建并使用极坐标系DBM 模拟研究了内爆和外爆过程的非平衡强度，并从动理学角度分析了粒子速度分布函数的主要特征。2015 年，许爱国等[24]构建并使用多松弛DBM 模拟研究了黏性和热传导在爆轰系统中对局域和全域非平衡效应的影响。2016 年，张玉东等[26]使用DBM 模拟研究了负温度系数对爆轰系统的影响机理。之后，林传栋等构建了多组分DBM 用于模拟预混、非预混和部分预混反应流系统[27]，又使用高效、准确的多松弛DBM研究了不稳定反应热流演化过程的非平衡行为[30]。

为了在连续性极限条件下恢复N-S 方程组或者更高阶的流体力学方程组，并且为了包含更准确和详实的热力学非平衡信息，要求DBM 包含足够多的动理学矩关系。现有的DBM 模型至少能够恢复N-S方程组，与之对应，至少包含了7 个张量形式的动理学矩关系。例如，对于二维DBM，当包含额外自由度时，需要构建16 速度的离散速度模型[30]，与之对应，需要包含16 组形式统一的DBM 演化方程。然而，在实际的物理系统中，有时（或有些区域）物理场蕴含的非平衡效应可能较弱，不需要特精确和太费时的模拟工具，而需要使用较为简洁和高效的模型。

为此，本文首次提出并构建了简化的二维DBM。该模型只包含9 个独立的动理学矩关系，并且在连续性极限条件下能够恢复带有化学反应的Euler 方程组。该DBM 包含9 个离散速度，基于9 组形式统一的DBM 演化方程。相对于之前构建的16 速度模型[30]或者24 速度模型[24]，该DBM 具有更高的运算效率。本项目的开展有助于推动动理学方法在化学反应流的发展应用，并进一步促进其在复杂多尺度物理系统的模拟研究。

1 二维简化DBM

1.1 离散玻尔兹曼方程

图1 离散速度示意图Fig. 1 Sketches for the discrete velocity model

1.2 动理学矩关系

1.3 化学反应进程

为了简单起见，本文忽略电离、热辐射、活性自由基和可逆反应等效应[25]；同时，考虑较为一般的情况，即化学反应的时间尺度大于热力学弛豫时间，且小于流体系统的特征尺度[25]。那么，由此可以进一步假设：化学反应引起的粒子速度分布函数的改变率近似等于平衡态速度分布函数的改变率，且在化学反应的时间尺度内，局部流体的密度和速度尚未来得及改变[25]。

1.4 流体力学模型和非平衡物理量

可以看出，在物理描述功能方面，简化二维DBM除了具有Euler 模型描述流体系统演化的功能之外，还具有描述一定热力学非平衡行为的功能。因此，该DBM 动理学模型具有超越Euler 模型物理描述范围的优势。也就是说，在该简化二维DBM 中，DBM等价于带有化学反应的Euler 方程组和包含少量非平衡效应的粗粒化模型。

需要说明的是，对于二维DBM，当对应包含额外自由度的N-S 方程组时，需要构建16 速度[30]（或24 速度[24]）的离散速度模型，与之对应地需要包含16 组[30]（或24 组[24]）形式统一的DBM 演化方程；当对应Burnett 方程组时，至少需要26 个离散速度和26 组演化方程[36]，计算量相应增大。然而，在实际的物理系统中，有时（或有些区域）物理场蕴含的非平衡效应可能较弱，不需要特精确和太费时的模拟手段，反而使用较为简洁和高效的模型更为可行。此时，仅有9 个离散速度和9 组演化方程的简化DBM 更加可取。总之，简化DBM 的优势为：离散速度数量较少、计算效率高。

另外，在已有的二维DBM 中，平衡态分布函数和反应项都至少各自满足16 个[30]（或24 个[24]）独立的矩关系，DBM 等价于带有化学反应的N-S 方程组和包含较多非平衡效应的粗粒化模型。尤其对于文献[36]，平衡态分布函数满足26 个独立的矩关系，该DBM 等价于Burnett 方程组和包含更多非平衡效应的粗粒化模型。在本文构建的简化模型中，平衡态分布函数和反应项都仅满足9 个独立的矩关系。相比于已有的DBM，该DBM 的局限性为：满足的动理学矩关系较少、物理精度较低，仅适用于热力学非平衡程度较弱的物理系统。

1.5 数值离散格式

DBM 使用离散玻尔兹曼方程描述化学反应流的演化过程，该方程形式简单、易于编程。同时，该方程具有时空局域性特点，易于实现并行化，方便在高性能计算集群运行。另外，对于离散玻尔兹曼方程式（1）的时间和空间偏微分，可以借用已有的数值格式处理。本文使用的是二阶精度的Runge-Kutta 格式和二阶精度的无波动、无自由参数的耗散有限差分格式（nonoscillatory and nonfree-parameter dissipation finite difference scheme，NND）[38]分别用于处理方程式（1）的时间和空间偏微分。

1.5.1 Runge-Kutta 格式

本文采用Runge-Kutta 格式处理时间离散化问题，具体计算流程见图2。

图2 二阶Runge-Kutta 格式流程图Fig. 2 Flow chart of the second-order Runge-Kutta scheme

1）程序的开头是初始化，根据物理场的初始条件，输入参数；

2）根据给定的时间步数，即程序循环周期，判断是否终止运行；

3）根据给定的时间间隔，判断是否输出数值结果；

需要说明的是，上述NND 格式具有二阶空间精度，是二阶中心格式和二阶迎风格式的结合，是具有负系数的四阶耗散项。因此，该NND 格式可以抑制奇偶失联震荡，能够有效捕捉冲击波和爆轰波[38]。

2 数值验证

本部分通过3 个典型算例对DBM 进行数值验证。首先，通过均匀化学反应系统验证本模型能够将化学反应与物理场自然耦合，模型既适用于吸热反应，也适用于放热反应，并且反应流系统的比热比可调；其次，通过Sod 激波管验证本模型适用于可压缩流体系统，并且能够同时捕捉稀疏波、物质界面和冲击波的时空演化；最后，通过爆轰波算例开展了网格无关性验证，并验证本模型适用于超声速可压缩化学反应流。

2.1 均匀化学反应系统

为了验证本模型能够将化学反应与物理场自然耦合，本节首先考虑均匀化学反应系统：在初始时刻，化学反应物均匀、静止分布在计算区域，初始温度为T0=2。经过一段足够长的时间后，化学反应结束，反应物完全转变成产物，化学能完全释放并转变成系统内能。由于该系统为均匀系统，每处的物理量分布都相同，所以为了提高计算速度，可以只用一个网格点，即NX×NY=1×1。 松 弛 时 间 τ=4×10−6，空间步长Δx=Δy=1×10−4，时间步长 Δt=2×10−6，离散参数va=−3.7、vb=−2、vc=−1.5 ， ηa=4、 ηb=0、 ηc=0。另外，对计算区域四周都采用周期边界条件。

图3 展示了化学反应后不同化学能对应的系统温度，其关系见式（37），此处设置比热比 γ=1.4。符号代表DBM 模拟结果，实线代表理论解：

可以看出，DBM 的模拟结果与理论解完全吻合。而且，图3 中化学能的数值涵盖负值到正值的范围。其中负值区域代表化学反应吸热，正值区域代表化学反应放热。也就是说，该DBM 既适用于化学吸热系统，也适用于化学放热系统。

图3 化学反应后不同化学能对应的系统温度Fig. 3 Temperature after chemical reaction under various chemical energies

下面考虑固定化学能不变Q=1，调节比热比数值的情况。图4 展示了化学反应后不同比热比对应的系统温度。符号代表DBM 模拟结果，实线代表理论解（式（37））。同样，DBM 给出的模拟结果与理论解完全吻合。因此，数值模拟证明该模型能够适用于不同比热比的化学反应流系统。

图4 化学反应后不同比热比对应的系统温度Fig. 4 Temperature after chemical reaction under various specific heat ratios

2.2 激波管

图5 展示了在t= 0.02 时刻激波管中的物理量分布，图5（a～d）依次是密度、温度、压强和水平速度。符号代表DBM 模拟结果，实线代表理论解。可以看出：在最左侧的是稀疏波（向左传播），在中间的是物质分界面，在最右侧的是冲击波（向右传播）。其中，稀疏波在水平方向的跨度最大，各物理量都有清晰的空间梯度；在物质界面两侧，密度和温度具有明显空间梯度，而压强和速度均匀分布；在冲击波波阵面，各物理量梯度最大，具有剧烈的空间变化。在以上演化过程中，DBM 与理论解具有很高的匹配度。这充分说明，DBM 能够很好地捕捉稀疏波、物质界面和冲击波的空间分布。

图5 物理量在激波管中的分布Fig. 5 Profiles of physical quantities in the shock tube

2.3 爆轰波

爆轰是一种特殊的化学反应流现象，在其化学反应区前沿是以超声速传播的冲击波[1]。在爆轰演化过程中，密度、温度和压强急剧上升，流场产生剧烈变化。物理量时空变化如此剧烈，在数值模拟中容易产生数值振荡，所以对爆轰波的数值模拟方法要具有较好的鲁棒性。本节将验证我们的DBM 模型能够有效模拟爆轰波的传播过程。

下面考虑沿水平方向传播的爆轰波。在一个水平放置的长度为l0=0.06的直通道内，初始物理场设置如下：

下标L表示左侧区域 0 ≤x＜0.99l0，下标R表示右侧区域 0.99l0≤x≤l0，左右两侧的物理量满足Hugoniot关系。设置化学反应热Q=2， 马赫数Ma=2.12643，比热比 γ=1.4 ， 松弛时间 τ=5×10−6， 离散参数va=−3.7、vb=−2、vc=−1.5， ηa=4、 ηb=0、 ηc=0。另外，上下采用周期边界条件，左右采用出口边界条件。

首先进行网格无关性验证。图6 展示了不同空间步长的密度分布图，即 Δx1=5×10−5、 Δx2=1×10−4、Δx3=2×10−4、Δx4=4×10−4；与之对应的水平网格点数分别是NX=1 200、 600、 300、 150。可以看出，随着空间步长的减小（网格数的增加），模拟结果之间的差异越来越小。尤其是在 Δx1和 Δx2下的模拟结果几乎完全重合。因此，考虑到计算效率，可以使用空间步长 Δx2进行数值模拟。同样，可以对不同的时间步长进行数值对比。结果表明，时间步长 Δt=2×10−6和4×10−6之间的模拟结果差异可以忽略。

图6 在不同空间步长下爆轰波附近的密度分布图Fig. 6 Profiles of density around the detonation wave under various spatial steps

图7 展示了在空间步长Δx=1×10−4、 时间步长Δt=2×10−6下爆轰波附近的模拟结果。图7（a～d）依次对应密度、温度、压强和水平速度的空间分布，其中符号代表DBM 模拟结果，实线代表Zeldovich-Neumann-Doering（ZND）理论解[1]。可以看出，DBM模拟结果与ZND 理论解具有非常好的一致性。尤其是当化学反应结束后，在爆轰波后侧形成稳定的“平台区域”，此处的DBM 模拟结果为 ρ=1.48065、T=2.06216、p=3.05334、ux=−1.69974。对比ZND 理论解，相对误差依次为： 0 .015%、 0 .047%、0 .032%、0 .012%，误差非常小。因此该DBM 能够有效模拟爆轰现象。

图7 爆轰波附近的物理量分布Fig. 7 Profiles of physical quantities around the detonation wave

3 结 论

基于动理学理论，本文提出了适用于超声速可压缩化学反应流的二维简化DBM。该模型具有下列特点：

1）构建并使用了仅有9 个离散速度的模型，即离散速度分为3 组，每组含有3 个速度，各组速度的方向均匀分布、大小独立可调。由于离散速度较少，模型具有较高的运算效率。

2）为了描述分子转动和振动对应的额外自由度，引入了三组独立可调的参数用于描述额外自由度部分的内能。由此，该DBM 具备了模拟比热比可调的反应流系统的功能。

3）平衡态离散速度分布函数与化学反应项各自都满足9 个独立的动理学矩关系，并可以通过矩阵求逆的方式获得其数学表达式。该方法具有物理精确、运算高效的特点。

4）通过Chapman-Enskog 多尺度分析可以证明，该DBM 除了在连续性极限条件下可以恢复包含化学反应的Euler 方程组之外，还具有描述一定热力学非平衡行为的功能。

5）使用形式统一的离散玻尔兹曼方程，方法简单、易于编程。同时，离散玻尔兹曼演化方程具有时空局域性特点，易于实现并行化，方便在高性能计算集群运行。

本文开展了3 个数值算例，验证了该DBM 的可靠性。首先，通过均匀反应系统验证模型能够将化学反应与物理场自然耦合，适用于吸热或放热反应，并且反应流系统的比热比可调；其次，通过激波管验证模型适用于可压缩流体系统，并且能够同时捕捉稀疏波、物质界面和冲击波的时空演化；最后，通过爆轰波算例开展了网格无关性验证，并验证模型适用于超声速可压缩化学反应流。

已有的二维DBM 至少包含16 个离散速度、动理学矩关系和演化方程，等价于N-S（或者Burnett）方程组和包含更多非平衡效应的粗粒化模型[24,30,36]。与之相比，简化DBM 的优势为离散速度数量较少、计算效率较高；局限性为物理精度较低，适用于热力学非平衡程度较弱的物理系统。