化学反应流不仅广泛存在于自然界中,更在人类日常生活、工业生产、航空航天和国防建设等诸多领域发挥着越来越重要的作用[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 离散玻尔兹曼方程化学反应流的离散玻尔兹曼方程如下:
$ \frac{{\partial {f_i}}}{{\partial t}} + {{\boldsymbol{v}}_i} \cdot \nabla {f_i} = - \frac{1}{\tau }\left( {{f_i} - f_i^{{\rm{eq}}}} \right) + {R_i} $ | (1) |
式中:
本文构建了二维九速模型,其离散速度
$ {{\boldsymbol{v}}_i} = \left\{ \begin{gathered} {v_a}\left( {\cos \frac{{2{\text{π}} i - 2{\text{π}} }}{3},\sin \frac{{2{\text{π}} i - 2{\text{π}} }}{3}} \right),\;\;\;\;{\text{1}} \leqslant i \leqslant 3 \hfill \\ {v_b}\left( {\cos \frac{{2{\text{π}} i - 2{\text{π}} }}{3},\sin \frac{{2{\text{π}} i - 2{\text{π}} }}{3}} \right),\;\;\;\;{\text{4}} \leqslant i \leqslant 6 \hfill \\ {v_c}\left( {\cos \frac{{2{\text{π}} i - 2{\text{π}} }}{3},\sin \frac{{2{\text{π}} i - 2{\text{π}} }}{3}} \right),\;\;\;\;{\text{7}} \leqslant i \leqslant 9 \hfill \\ \end{gathered} \right. $ | (2) |
式中:
图1给出了两种情况下的离散速度示意图,其中图1(a)对应
同时,为了描述分子转动和振动对应的内能,引入额外自由度
$ {\eta _i} = \left\{ \begin{gathered} {\eta _a},\;\;\;\; 1 \leqslant i \leqslant 3 \hfill \\ {\eta _b},\;\;\;\; 4 \leqslant i \leqslant 6 \hfill \\ {\eta _c},\;\;\;\; 7 \leqslant i \leqslant 9 \hfill \\ \end{gathered} \right. $ | (3) |
与式(2)中的离散速度类似,控制参数
$ \gamma = \frac{{D + I + 2}}{{D + I}} $ | (4) |
其中
流体的密度
$ \rho = \sum\nolimits {{f_i}} $ | (5) |
$ \rho {\boldsymbol{u}} = \sum\nolimits {{f_i}{{\boldsymbol{v}}_i}} $ | (6) |
$ E = \frac{1}{2}\sum\nolimits {{f_i}\left( {{{\left| {{{\boldsymbol{v}}_i}} \right|}^2} + \eta _i^2} \right)} $ | (7) |
通过式(5)~式(7),由离散分布函数可以计算得到系统的密度、速度和能量。根据温度与能量的关系:
$ E = \frac{1}{2}\rho {\left| {\boldsymbol{u}} \right|^2} + \frac{1}{{\gamma - 1}}\rho T $ | (8) |
可得温度:
$ T = \left( {\gamma - 1} \right)\left( {\frac{1}{\rho }E - \frac{1}{2}{{\left| {\boldsymbol{u}} \right|}^2}} \right) $ | (9) |
另外,平衡态离散速度分布函数
$ \rho = \sum\nolimits {f_i^{{\rm{eq}}}} $ | (10) |
$ \rho {u_{\alpha} } = \sum\nolimits {f_i^{{\rm{eq}}}{v_{i\alpha }}} $ | (11) |
$ \rho \left[ {{u^2} + \left( {D + I} \right)T} \right] = \sum\nolimits {f_i^{{\rm{eq}}}\left( {v_i^2 + \eta _i^2} \right)} $ | (12) |
$ \rho \left( {{\delta _{\alpha \beta }}T + {u_{\alpha} }{u_{\beta} }} \right) = \sum\nolimits {f_i^{{\rm{eq}}}{v_{i\alpha }}{v_{i\beta }}} $ | (13) |
$ \rho {u_{\alpha} }\left[ {\left( {D + I + 2} \right)T + {u^2}} \right] = \sum\nolimits {f_i^{{\rm{eq}}}\left( {v_i^2 + \eta _i^2} \right){v_{i\alpha }}} $ | (14) |
其中
$ {{\boldsymbol{\hat f}}^{{{\rm{eq}}}}} = {\boldsymbol{M}}{{\boldsymbol{f}}^{{{\rm{eq}}}}} $ | (15) |
上式中
$ {{\boldsymbol{f}}^{{\rm{eq}}}} = {{\boldsymbol{M}}^{ - 1}}{{\boldsymbol{\hat f}}^{{\rm{eq}}}} $ | (16) |
其中
对于化学反应进程的描述,可以使用(半)经验模型、简化或者详细的化学反应机理。本文采用两步化学反应模型[37],即使用变量
$ \xi ' = H{k_I}\exp \left[ {{E_I}\left( {T_s^{ - 1} - {T^{ - 1}}} \right)} \right] $ | (17) |
$ \lambda ' = \left( {1 - H} \right){k_R}\left( {1 - \lambda } \right)\exp \left( { - {E_R}{T^{ - 1}}} \right) $ | (18) |
其中
为了简单起见,本文忽略电离、热辐射、活性自由基和可逆反应等效应[25];同时,考虑较为一般的情况,即化学反应的时间尺度大于热力学弛豫时间,且小于流体系统的特征尺度[25]。那么,由此可以进一步假设:化学反应引起的粒子速度分布函数的改变率近似等于平衡态速度分布函数的改变率,且在化学反应的时间尺度内,局部流体的密度和速度尚未来得及改变[25]。
下面推导化学反应项的数学表达式。首先引入以下物理量:单位质量的反应物释放的化学能为
$ R = \frac{{ - \left( {1 + D} \right)IT + I{{\left| {{\boldsymbol{v}} - {\boldsymbol{u}}} \right|}^2} + {\eta ^2}}}{{2I{T^2}}}{f^{{\rm{eq}}}}T' $ | (19) |
离散形式的化学反应项
$ \iint {R\varPsi {{\rm{d}}} {\boldsymbol{v}}{{\rm{d}}} \eta } = \sum\nolimits {{R_i}} {\varPsi _i} $ | (20) |
其中,
$ {\boldsymbol{\hat R}} = {\boldsymbol{MR}} $ | (21) |
上式中
$ {\boldsymbol{R}} = {{\boldsymbol{M}}^{ - 1}}{\boldsymbol{\hat R}} $ | (22) |
该式即为化学反应项的计算方式。
1.4 流体力学模型和非平衡物理量理想气体状态方程
通过Chapman-Enskog多尺度分析可以证明,在连续性极限条件下,该DBM模型可以恢复包含化学反应的Euler方程组,即:
$ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho {u_{\alpha} }}}{{\partial {r_{\alpha} }}} = 0 $ | (23) |
$ \frac{{\partial \rho {u_{\alpha} }}}{{\partial t}} + \frac{\partial }{{\partial {r_{\alpha} }}}\left( {{\delta _{\alpha \beta }}p + \rho {u_{\alpha} }{u_{\beta} }} \right) = 0 $ | (24) |
$ \frac{{\partial E}}{{\partial t}} + \frac{\partial }{{\partial {r_{\alpha} }}}\left( {p{u_{\alpha} } + E{u_{\alpha} }} \right) = \rho \lambda 'Q $ | (25) |
式(23)~式(25)分别描述了质量守恒、动量守恒和能量守恒定律。特别是对于能量方程式(25),在化学反应过程中,化学能转化为系统其他形式的能量(转化速率为
需要说明的是,动理学关系式(10)~式(14)是从DBM获得Euler方程组的前提条件。对于前三个子式(10)~式(12),将方程中的
$ \Delta _{\alpha \beta }^* = \sum\nolimits {\left( {{f_i} - f_i^{{\rm{eq}}}} \right)v_{i\alpha }^*v_{i\beta }^*} $ | (26) |
$ \Delta _{\alpha} ^* = \sum\nolimits {\left( {{f_i} - f_i^{{\rm{eq}}}} \right)\left( {{v_i^*}^2 + \eta _i^2} \right)v_{i\alpha }^*} $ | (27) |
其中
可以看出,在物理描述功能方面,简化二维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。
1)程序的开头是初始化,根据物理场的初始条件,输入参数;
2)根据给定的时间步数,即程序循环周期,判断是否终止运行;
3)根据给定的时间间隔,判断是否输出数值结果;
4)根据离散速度分布函数
5)根据宏观物理量(
6)根据宏观物理量(
7)考虑对流项、碰撞项和反应项对分布函数的改变,得到离散速度分布函数
8)根据离散速度分布函数
9)根据宏观物理量(
10)根据宏观物理量(
11)考虑对流项、碰撞项和反应项对分布函数的改变,得到离散速度分布函数
需要补充说明的是第4~7步依次与第8~11步相似。
1.5.2 NND格式NND有限差分格式用于处理空间偏微分,具体如下:
$ {v_{ir}}\frac{{\partial {f_i}}}{{\partial r}} = \frac{1}{{\Delta r}}\Bigg[H\Bigg(ir + \frac{1}{2}\Bigg) - H\Bigg(ir - \frac{1}{2}\Bigg)\Bigg] $ | (28) |
$ H\Bigg(ir + \frac{1}{2}\Bigg) = {H_L}\Bigg(ir + \frac{1}{2}\Bigg) + {H_R}\Bigg(ir + \frac{1}{2}\Bigg) $ | (29) |
$\begin{split} & {H_L}\Bigg(ir + \frac{1}{2}\Bigg) = f_i^ + (ir) + \\ & \qquad\frac{1}{2}\min \,\bmod \,\Bigg[\Delta f_i^ + \Bigg(ir + \frac{1}{2}\Bigg),\Delta f_i^ + \Bigg(ir - \frac{1}{2}\Bigg)\Bigg] \end{split} $ | (30) |
$ \begin{split} &{H_R}\Bigg(ir + \frac{1}{2}\Bigg) = f_i^ - (ir + 1)- \\ & \qquad \frac{1}{2}\min \,\bmod \,\Bigg[\Delta f_i^ - \Bigg(ir + \frac{1}{2}\Bigg),\Delta f_i^ - \Bigg(ir + \frac{3}{2}\Bigg)\Bigg] \end{split} $ | (31) |
$ \Delta f_i^ + \Bigg(ir + \frac{1}{2}\Bigg) = f_i^ + (ir + 1) - f_i^ + (ir) $ | (32) |
$ \Delta f_i^ - \Bigg(ir + \frac{1}{2}\Bigg) = f_i^ - (ir + 1) - f_i^ - (ir) $ | (33) |
$ f_i^ + (ir) = \frac{1}{2}({v_{ir}} + \left| {{v_{ir}}} \right|){f_i} $ | (34) |
$ f_i^ - (ir) = \frac{1}{2}({v_{ir}} - \left| {{v_{ir}}} \right|){f_i} $ | (35) |
在式(28)~式(35)中符号
$ \begin{split}\min\mathrm{mod}[X,Y] = \left\{\begin{array}{l}X,Y中绝对值小的那个值\text{,}若XY\geqslant 0\\ 0\text{,}若XY < 0\end{array}\right. \\[-9pt]\end{split}$ | (36) |
需要说明的是,上述NND格式具有二阶空间精度,是二阶中心格式和二阶迎风格式的结合,是具有负系数的四阶耗散项。因此,该NND格式可以抑制奇偶失联震荡,能够有效捕捉冲击波和爆轰波[38]。
2 数值验证本部分通过3个典型算例对DBM进行数值验证。首先,通过均匀化学反应系统验证本模型能够将化学反应与物理场自然耦合,模型既适用于吸热反应,也适用于放热反应,并且反应流系统的比热比可调;其次,通过Sod激波管验证本模型适用于可压缩流体系统,并且能够同时捕捉稀疏波、物质界面和冲击波的时空演化;最后,通过爆轰波算例开展了网格无关性验证,并验证本模型适用于超声速可压缩化学反应流。
2.1 均匀化学反应系统为了验证本模型能够将化学反应与物理场自然耦合,本节首先考虑均匀化学反应系统:在初始时刻,化学反应物均匀、静止分布在计算区域,初始温度为
图3展示了化学反应后不同化学能对应的系统温度,其关系见式(37),此处设置比热比
$ T = {T_0} + \left( {\gamma - 1} \right)Q $ | (37) |
可以看出,DBM的模拟结果与理论解完全吻合。而且,图3中化学能的数值涵盖负值到正值的范围。其中负值区域代表化学反应吸热,正值区域代表化学反应放热。也就是说,该DBM既适用于化学吸热系统,也适用于化学放热系统。
下面考虑固定化学能不变
为了验证该DBM能够适用于可压缩流体系统,本节考虑Sod激波管的模拟。初始物理场设置如下:
$ \left\{ {\begin{array}{*{20}{l}} {{{\left( {\rho ,T,{u_x}} \right)}_L} = \left( {1,1,0} \right)} \\ {{{\left( {\rho ,T,{u_x}} \right)}_R} = \left( {{\text{0}}{\text{.125}},{\text{0}}{\text{.8}},0} \right)} \end{array}} \right. $ |
下标L表示左侧区域
图5展示了在t = 0.02时刻激波管中的物理量分布,图5(a~d)依次是密度、温度、压强和水平速度。符号代表DBM模拟结果,实线代表理论解。可以看出:在最左侧的是稀疏波(向左传播),在中间的是物质分界面,在最右侧的是冲击波(向右传播)。其中,稀疏波在水平方向的跨度最大,各物理量都有清晰的空间梯度;在物质界面两侧,密度和温度具有明显空间梯度,而压强和速度均匀分布;在冲击波波阵面,各物理量梯度最大,具有剧烈的空间变化。在以上演化过程中,DBM与理论解具有很高的匹配度。这充分说明,DBM能够很好地捕捉稀疏波、物质界面和冲击波的空间分布。
爆轰是一种特殊的化学反应流现象,在其化学反应区前沿是以超声速传播的冲击波[1]。在爆轰演化过程中,密度、温度和压强急剧上升,流场产生剧烈变化。物理量时空变化如此剧烈,在数值模拟中容易产生数值振荡,所以对爆轰波的数值模拟方法要具有较好的鲁棒性。本节将验证我们的DBM模型能够有效模拟爆轰波的传播过程。
下面考虑沿水平方向传播的爆轰波。在一个水平放置的长度为
$ \left\{ {\begin{array}{*{20}{l}} {{{\left( {\rho ,p,{u_x}} \right)}_L} = \left( {1.48043,3.05433, - 1.69953} \right)} \\ {{{\left( {\rho ,p,{u_x}} \right)}_R} = \left( {1,1, - 2.51603} \right)} \end{array}} \right. $ |
下标L表示左侧区域
首先进行网格无关性验证。图6展示了不同空间步长的密度分布图,即
图7展示了在空间步长
基于动理学理论,本文提出了适用于超声速可压缩化学反应流的二维简化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的优势为离散速度数量较少、计算效率较高;局限性为物理精度较低,适用于热力学非平衡程度较弱的物理系统。
[1] |
LAW C K. Combustion physics[M]. Cambridge: Cambridge University Press, 2010.
|
[2] |
陈坚强, 张益荣, 郭勇, 等. 高超声速流动数值模拟方法及应用[M]. 北京: 科学出版社, 2019.
|
[3] |
艾邦成. 临近空间高超声速飞行器计算空气动力学[M]. 北京: 科学出版社, 2020.
|
[4] |
毛枚良, 陈亮中, 万钊, 等. 高超声速流动多物理效应对美国航天飞机气动力影响研究的回顾[J]. 空气动力学学报, 2017, 35(1): 1-12. MAO M L, CHEN L Z, WAN Z, et al. Review on influence of the multi-physical effects in hypersonic flow on aerodynamic characteristics of the United States space shuttle[J]. Acta Aerodynamica Sinica, 2017, 35(1): 1-12. DOI:10.7638/kqdlxxb-2015.0159 (in Chinese) |
[5] |
沈青. 稀薄气体动力学[M]. 北京: 国防工业出版社, 2003.
|
[6] |
樊菁. 稀薄气体动力学: 进展与应用[J]. 力学进展, 2013, 43(2): 185-201. FAN J. Rarefied gas dynamics: Advances and applications[J]. Advances in Mechanics, 2013, 43(2): 185-201. DOI:10.6052/1000-0992-13-018 (in Chinese) |
[7] |
WOLF-GLADROW D A. Lattice-gas cellular automata and lattice Boltzmann models: An introduction[M]. Springer, 2004. https://epic.awi.de/id/eprint/3739/1/Wol2000c.pdf
|
[8] |
HARDY J, POMEAU Y, DE PAZZIS O. Time evolution of a two-dimensional model system. I. Invariant states and time correlation functions[J]. Journal of Mathematical Physics, 1973, 14(12): 1746-1759. DOI:10.1063/1.1666248 |
[9] |
FRISCH U, HASSLACHER B, POMEAU Y. Lattice-gas automata for the Navier-Stokes equation[J]. Physical Review Letters, 1986, 56(14): 1505-1508. DOI:10.1103/physrevlett.56.1505 |
[10] |
SUCCI S. The lattice Boltzmann equation for fluid dynamics and beyond[M]. Oxford : Clarendon Press ; New York : Oxford University Press, 2001.
|
[11] |
何雅玲, 王勇, 李庆. 格子Boltzmann方法的理论及应用[M]. 北京: 科学出版社, 2009.
|
[12] |
郭照立, 郑楚光. 格子Boltzmann方法的原理及应用[M]. 北京: 科学出版社, 2009.
|
[13] |
SUCCI S, BELLA G, PAPETTI F. Lattice kinetic theory for numerical combustion[J]. Journal of Scientific Computing, 1997, 12(4): 395-408. DOI:10.1023/A:1025676913034 |
[14] |
YAMAMOTO K, TAKADA N, MISAWA M. Combustion simulation with lattice Boltzmann method in a three-dimensional porous structure[J]. Proceedings of the Combustion Institute, 2005, 30(1): 1509-1515. DOI:10.1016/j.proci.2004.08.030 |
[15] |
HOSSEINI S A, DARABIHA N, THÉVENIN D. Mass-conserving advection-diffusion lattice Boltzmann model for multi-species reacting flows[J]. Physica A:Statistical Mechanics and Its Applications, 2018, 499: 40-57. DOI:10.1016/j.physa.2018.01.034 |
[16] |
CHEN S, LIU Z H, ZHANG C, et al. A novel coupled lattice Boltzmann model for low Mach number combustion simulation[J]. Applied Mathematics and Computation, 2007, 193(1): 266-284. DOI:10.1016/j.amc.2007.03.087 |
[17] |
CHEN L, HE Y L, TAO W Q, et al. Pore-scale study of multiphase reactive transport in fibrous electrodes of vanadium redox flow batteries[J]. Electrochimica Acta, 2017, 248: 425-439. DOI:10.1016/j.electacta.2017.07.086 |
[18] |
CHEN L, WANG M Y, KANG Q J, et al. Pore scale study of multiphase multicomponent reactive transport during CO2 dissolution trapping
[J]. Advances in Water Resources, 2018, 116: 208-218. DOI:10.1016/j.advwatres.2018.02.018 |
[19] |
ZHANG L, WANG M R. Modeling of electrokinetic reactive transport in micropore using a coupled lattice Boltzmann method[J]. Journal of Geophysical Research:Solid Earth, 2015, 120(5): 2877-2890. DOI:10.1002/2014JB011812 |
[20] |
YAN Z F, YANG X F, LI S L, et al. Two-relaxation-time lattice Boltzmann method and its application to advective-diffusive-reactive transport[J]. Advances in Water Resources, 2017, 109: 333-342. DOI:10.1016/j.advwatres.2017.09.003 |
[21] |
ASADOLLAHI A, RASHIDI S, ESFAHANI J A. Simulation of liquid reaction and droplet formation on a moving micro-object by lattice Boltzmann method[J]. Meccanica, 2018, 53(4-5): 803-815. DOI:10.1007/s11012-017-0780-4 |
[22] |
YAN B, XU A G, ZHANG G C, et al. Lattice Boltzmann model for combustion and detonation[J]. Frontiers of Physics, 2013, 8(1): 94-110. DOI:10.1007/s11467-013-0286-z |
[23] |
LIN C D, XU A G, ZHANG G C, et al. Polar coordinate lattice Boltzmann kinetic modeling of detonation phenomena[J]. Communications in Theoretical Physics, 2014, 62(5): 737-748. DOI:10.1088/0253-6102/62/5/18 |
[24] |
XU A G, LIN C D, ZHANG G C, et al. Multiple-relaxation-time lattice Boltzmann kinetic model for combustion[J]. Physical Review E, 2015, 91(4): 043306. DOI:10.1103/physreve.91.043306 |
[25] |
LIN C D, XU A G, ZHANG G C, et al. Double-distribution-function discrete Boltzmann model for combustion[J]. Combustion and Flame, 2016, 164: 137-151. DOI:10.1016/j.combustflame.2015.11.010 |
[26] |
ZHANG Y D, XU A G, ZHANG G C, et al. Kinetic modeling of detonation and effects of negative temperature coefficient[J]. Combustion and Flame, 2016, 173: 483-492. DOI:10.1016/j.combustflame.2016.04.003 |
[27] |
LIN C D, LUO K H, FEI L L, et al. A multi-component discrete Boltzmann model for nonequilibrium reactive flows[J]. Scientific Reports, 2017, 7: 14580. DOI:10.1038/s41598-017-14824-9 |
[28] |
LIN C D, LUO K H. Mesoscopic simulation of nonequilibrium detonation with discrete Boltzmann method[J]. Combustion and Flame, 2018, 198: 356-362. DOI:10.1016/j.combustflame.2018.09.027 |
[29] |
LIN C D, LUO K H. MRT discrete Boltzmann method for compressible exothermic reactive flows[J]. Computers & Fluids, 2018, 166: 176-183. DOI:10.1016/j.compfluid.2018.02.012 |
[30] |
LIN C D, LUO K H. Discrete Boltzmann modeling of unsteady reactive flows with nonequilibrium effects[J]. Physical Review E, 2019, 99: 012142. DOI:10.1103/physreve.99.012142 |
[31] |
ZHANG Y D, XU A G, ZHANG G C, et al. Discrete Boltzmann method with Maxwell-type boundary condition for slip flow[J]. Communications in Theoretical Physics, 2018, 69(1): 77. DOI:10.1088/0253-6102/69/1/77 |
[32] |
ZHANG Y D, XU A G, ZHANG G C, et al. Entropy production in thermal phase separation: a kinetic-theory approach[J]. Soft Matter, 2019, 15(10): 2245-2259. DOI:10.1039/c8sm02637h |
[33] |
GAN Y B, XU A G, ZHANG G C, et al. Lattice Boltzmann study on Kelvin-Helmholtz instability: Roles of velocity and density gradients[J]. Physical Review E, 2011, 83(5): 056704. DOI:10.1103/physreve.83.056704 |
[34] |
YE H Y, LAI H L, LI D M, et al. Knudsen number effects on two-dimensional Rayleigh–Taylor instability in compressible fluid: based on a discrete Boltzmann method[J]. Entropy, 2020, 22(5): 500. DOI:10.3390/e22050500 |
[35] |
LIN C, LUO K H, XU A, et al. Multiple-relaxation-time discrete Boltzmann modeling of multicomponent mixture with nonequilibrium effects[J]. Physical Review E, 2021, 103(1): 013305. DOI:10.1103/physreve.103.013305 |
[36] |
GAN Y, XU A, ZHANG G, et al. Discrete Boltzmann trans-scale modeling of high-speed compressible flows[J]. Physical Review E, 2018, 97(5): 053312. DOI:10.1103/physreve.97.053312 |
[37] |
NG H D, RADULESCU M I, HIGGINS A J, et al. Numerical investigation of the instability for one-dimensional Chapman-Jouguet detonations with chain-branching kinetics[J]. Combustion Theory and Modelling, 2005, 9(3): 385-401. DOI:10.1080/13647830500307758 |
[38] |
ZHANG H X, ZHUANG F G. NND schemes and their applications to numerical simulation of two- and three-dimensional flows[J]. Advances in Applied Mechanics, 1991, 29: 193-256. DOI:10.1016/S0065-2156(08)70165-0 |