蒋励剑，赵文文，陈伟芳，尧少波

浙江大学 航空航天学院，杭州 310027

钱学森[1]根据来流气体稀薄程度将气体流域分为连续流域、滑移流域、过渡流域和自由分子流域，除连续流域外的其他3大流域又被合称为稀薄非平衡流域。在连续流域，气体满足连续介质假设，因而在宏观连续介质流体力学的基础下，形成了以Navier-Stokes(N-S)方程为控制方程的连续流数值模拟计算方法，结合物面滑移边界条件，可将N-S方程的适用领域扩展到滑移流域。随着气体稀薄程度的增加，连续介质假设在过渡流域和自由分子流域失效。当特征尺度逼近分子平均自由程时，在微观气体动理论的基础上形成了直接模拟蒙特卡洛(Direct Simulation of Monto Carlo, DSMC)方法[2]，以微观粒子的运动和碰撞直接模拟宏观流动。DSMC能够在稀薄流域给出与真实飞行试验较为吻合的流场与气动数据，然而其统计特点导致计算需要消耗大量的内存和机时，限制了该方法在工程领域近连续和多尺度流动中的应用。以Boltzmann方程为控制方程的统一气体动理论格式(Unified Gas-Kinetic Scheme, UGKS)方法是一套介于宏观和微观之间的介观方法[3]。然而UGKS在求解过程中依赖速度分布函数在速度空间的离散，同样给存储和计算负荷带来了很大的压力。针对稀薄非平衡流域计算效率与计算精度的冲突，构建同时具备高效率与高精度的数值方法在工程领域显得尤为重要。

近年来，随着大数据科学的兴起、大规模并行计算及GPU设备的普及，机器学习方法在计算机视觉、自然语言处理等各个领域迅速得到应用与发展。人工智能也被认为是21世纪3大尖端技术(基因工程、纳米科学、人工智能)之一[4]。在流体力学领域，机器学习也表现出巨大的潜力。在流动控制领域，诞生了面向流动控制基于机器学习的系统辨识与降阶模型、基于遗传规划的主动流动控制、基于人工神经网络与深度强化学习的主动流动控制等方法[5]。在改善RANS模型方面，机器学习方法的应用思路可以分为两类[6]：一是采用增加源项的方式改变模型的控制方程形式，如Tracey等[7]利用机器学习重建RANS模型方程的源项，提出了开发湍流模型函数形式的新方法；二是构造偏差函数，再与原RANS结果相加作为雷诺应力预示结果，如Wang等[8]用随机森林方法学习RANS模型偏差量，提高了RANS模型的准确性；Xiao等[9]针对RANS和DNS结果的差异，利用机器学习方法并引进了物理约束，对RANS结果进行了修正。另一方面，也有学者选择不再修正现有模型，而是直接寻找流场物理量与湍流之间的关系，如Zhu等[10]构建黑箱模型，采用径向基函数神经网络，建立了基于机器学习的高维数据驱动网络模型；Ling等[11]利用张量不变基将伽利略不变性嵌入到预测的各向异性张量中，证明满足伽利略不变性的神经网络在部分条件下具有更好的预测精度。而在稀薄流域，李廷伟等[12]采用极端随机树模型对非线性本构关系的复杂高维非线性回归问题进行了建模，对基于数据驱动非线性本构方程数值计算方法(DNCR)进行了初步验证。Zhang等[13]基于DSMC数据利用稀疏回归方法推导宏观控制方程，证明数据驱动与分子模拟结合具有广阔的发展前景。由于研究刚刚起步，这些适用于稀薄非平衡流动问题的计算方法目前仅应用于简单的流动问题，尚未扩展到较全面和复杂的多尺度稀薄非平衡流动。

在各种机器学习模型中，决策树是一种基本的分类与回归方法[14]。决策树算法的概念最初可以追溯到1966年提出的CLS学习系统，而可用于回归预测的CART算法在1984年由Breiman等提出。2001年，Breiman[15]继续在此基础上提出随机森林算法，其本质相当于分类树的组合，通过几个模型降低泛化误差，针对数值属性的特征，将不同特征的样本分配到对应的分支。由于随机森林的计算效率、精度等综合性能较好，在如环境保护[16]、推荐系统[17]、道路交通情况预测[18]、结构裂纹监测[19]、动力装置可靠性分析[20]等许多领域都有应用。对比随机森林，极端随机树[21]则在随机提取样本的基础上随机提取特征，因而随机性与抗干扰能力更强。

综上所述，现有DNCR方法采用了极端随机树方法，直接提取流场物理量与物理量梯度作为特征值与标记值，然而这些物理量大部分天然不具备伽利略不变性，使得所训练模型的预测范围严重依赖于训练集坐标系的固有方向，限制了模型的泛化性能。针对上述问题，本文在原始DNCR方法的基础上提出了运用旋转不变量进行训练与学习的改进DNCR方法。改进DNCR方法仍以N-S方程与UGKS方法作为理论基础与数据来源，采用旋转不变量代替特征值与标记值中包含空间坐标方向信息的流场物理量，使得所得训练模型能够适用于任意旋转后的坐标系框架，增强了模型的泛化能力和迁移学习能力。为了验证改进方法的有效性，本文以不同稀薄来流条件和不同几何模型下的顶盖方腔驱动流和二维圆柱/椭圆柱绕流为例，开展了特征参数选取、模型参数调优、泛化性能评估等工作，并评估改进方法与原有DNCR方法的计算精度提升情况。计算结果表明：采用旋转不变量进行训练和学习后，改进DNCR方法在模型的泛化性能和流场物理量预示精度上明显优于采用极端随机树直接学习流场物理量的原始DNCR方法。

1 DNCR计算方法与旋转不变性修正

数据驱动非线性本构方程数值计算方法的基本原理是将N-S方程和UGKS方法的计算结果作为流场样本训练数据，采用机器学习方法构建热流/应力项与流场特征参数的高维复杂非线性回归关系模型，对N-S方程本构关系进行非线性修正，最后通过求解耦合了数据驱动非线性回归关系模型的宏观守恒方程得到待预测状态的稀薄非平衡流数值解[12]。三维N-S流动控制方程在笛卡尔直角坐标系下的具体表达式为

(1)

式中：Q为守恒量；E、F、G为无黏通量；Ev、Fv、Gv为黏性通量。其黏性应力张量和热流项的表达式为

(2)

随着气体稀薄程度变高，努森数增大，基于连续介质假设的N-S流动控制方程逐渐失效。在稀薄非平衡流域，流动问题主要通过Boltzmann方程来描述，该方程运用概率速度分布函数来描述稀薄气体分子的运动状态，是气体动理学的基础控制方程[22]。UGKS采用由BGK[23]简化的Boltzmann模型方程计算气体分子的迁移和碰撞过程：

e-(t-t′)/τdt′+e-t/τf(x-ut,u,0)

(3)

式中：分布函数f为坐标x、分子速度u和时间t的函数，用来表示t时刻以速度u到达坐标x的分子数密度；g为平衡态分布函数，依赖于宏观物理量；τ为松弛时间，量级和分子平均碰撞时间相同。然而，UGKS将分子运动和碰撞过程耦合在一起的同时，也使用了额外的速度离散空间，使得在面对复杂高速流动问题时收敛较慢，计算效率较低。

以二维问题为例，对热流/应力项进行如下非线性修正：

(4)

将式(4)代入质量、动量与能量守恒方程中进行迭代计算至收敛，最终获得DNCR方法的流场预示结果。由于非线性修正式(4)的约束，最终收敛得到的流场解将趋近于训练模型给出的应力、热流与流场梯度。

在原始DNCR方法中，由于特征值和标记值均直接选取流场的物理量与物理量梯度，不具备天然的伽利略不变性。在面对不同的待预测状态与几何模型时，只有当待预测状态的坐标系方向与训练集坐标系方向相同时，才能得到较为理想的预测结果。而当坐标系发生诸如旋转变化时，预测结果可能无法收敛，导致训练模型适用范围受限，泛化性能亟待提升。

为了使训练得到的机器学习模型能够在任意旋转坐标系下适用，提高模型的泛化性能，拓展模型的适用范围。本文对原方法中的特征值和标记值进行了优化与旋转不变性处理，其中DNCR方法的计算流程如图1所示。

图1 DNCR方法计算流程图

(5)

由于应力Δτij为对称张量，可以将其表示为Δτij=VΛVT形式，其中Λ为特征值λ1、λ2、λ3构成的对角阵，代表应力Δτij的主应力方向，不随坐标变换而变化，由此得到3个标记值λ1、λ2、λ3；V为特征值对应的特征向量v1、v2、v3构成的矩阵，特征值对应的特征向量两两正交，并且在此问题下表现为旋转矩阵，即|V|=-1，V记录了当前坐标系的方位信息。由于矩阵无法直接作为标记值参与训练和预测，将旋转矩阵变换为四元数(具体变换过程见附录A)，因而得到4个标记值q0、q1、q2、q3，代表应力Δτij的坐标方向信息。速度梯度的合成量采取与应力同样的处理方式。

在预测过程中，对N-S计算收敛得到的KnP、Knρ、udp/dx等23个旋转不变量作为机器学习模型输入；将UGKS与N-S的应力、热流、速度梯度、温度梯度的张量特征值差量λi和由N-S特征向量变换为UGKS特征向量的旋转矩阵Vq(Vq由四元数qj表示，具体变换过程详见附录A)作为模型输出。以应力张量为例，对DNCR续算所需的Δτij，有

(6)

2 机器学习模型选择与参数调优

改进DNCR方法使用极端随机树作为学习模型，极端随机树本质上是多个决策树的集合，决策树是一种基本的分类与回归方法。决策树的生成就是递归地构建二叉决策树的过程，回归树采用平方误差最小化准则。每一棵回归树对应着输入空间(即特征空间)的一个划分以及在划分的单元上的输出值。假设已将输入空间划分为M个单元，并且在每个单元 上有一个固定的输出值，于是回归树模型可表示为[14]

(7)

在输入空间的划分确定后，一般采用平方误差来表示回归树对训练集的预测误差，用使平方误差最小的准则来求解每个单元的最优输出。通常使用的最小二乘回归树算法，在训练数据集的输入空间中，递归的将每个区域分为两个子区域并决定其输出值，来构建二叉决策树：

1) 选择最优切分变量j与切分点s，求解

(8)

遍历变量，对固定的切分变量 扫描切分点，选择使式(8)达到最小值的对。

2) 用选定的划分区域并决定相应的输出值：

(9)

(10)

3) 继续对两个字区域调用步骤1)和步骤2)，直至满足停止条件。

4) 将输入空间划分为M个区域R1,R2,…,RM，生成决策树：

(11)

极端随机树作为二叉决策树的一种集成算法，对所有训练样本直接随机选择分叉特征得到结果，实现简单，不容易产生过拟合，易找到对结果影响较大的特征，并且可并行化、速度快。

本文对基学习器数目以及最大特征数进行调优，通过可决系数的大小对两个超参数进行寻优，可决系数公式为

(12)

图2 基学习器数目调优

图3 最大特征数调优

3 算例与分析

UGKS方法的数据作为训练集中最重要数据来源之一，需要确保能够对稀薄非平衡流场进行较为准确的描述，因此首先与文献中DSMC[25]数据进行对比，验证UGKS方法的准确性。计算采用单原子氩气圆柱绕流，来流条件Ma=5，Kn=0.1，来流分子数密度为n=1.294 4×1021/m3，黏性系数μ∞=2.117×10-5Pa·s，壁面温度为Tw=273 K，壁面与DSMC保持一致均设置为漫反射边界条件。UGKS方法物理空间网格剖分为50×64，速度空间离散网格剖分为93×93。

图4和图5给出了UGKS和DSMC两种方法沿圆柱表面从驻点到后缘的物面应力θ、热流q结果对比曲线。其中横坐标为圆柱表面位置与圆心连线与y轴夹角角度，纵坐标分别为表面热流与黏性应力值。从计算结果可以看出，UGKS结果与DSMC结果吻合较好，验证了UGKS计算模块的可靠性。

图4 UGKS与DSMC热流对比

图5 UGKS与DSMC应力对比

本文所发展的改进DNCR方法采用旋转不变量作为特征值和标记值，在任意方向的坐标系下均能够得到相同的流场结果，与实际物理现象吻合。原有DNCR方法直接采用流场物理量作为机器学习模型特征值和标记值，对训练集坐标系发生旋转、平移后的待预测状态无法开展准确预测。为了便于方法之间对比，算例仅展示固定坐标系下的对比结果。需要说明的是，改进DNCR方法针对任意坐标系(旋转、平移)后的待预测状态都能得到相同的流场计算结果。

为了反映新的特征参数构造所带来的模型泛化性能提升，选取二维圆柱和二维方腔两种外形进行验证和对比，并对部分待预测状态进行适当的几何外形变化。二维圆柱采用直径0.304 8 m(12 in)下的数据为训练集，分别预测直径 0.304 8 m下的二维圆柱外形，长轴0.609 6 m、短轴0.304 8 m的二维椭圆柱外形两种几何形状在不同Kn数下的流场结果。二维方腔以外形 1 m×1 m的数据为训练集，分别预测1 m×1 m，2 m×1 m两种几何外形方腔在不同Kn数下的流场结果。最终所得结果的精度提升采用式(14)进行评估：

(13)

3.1 二维圆柱氩气高超声速绕流

首先验证二维圆柱外形下的预测结果，计算采用单原子气体氩气，来流马赫数与来流温度分别为Ma=10,T0=198.439 K，气体的物性参数取μ0=5.071 158×10-5Pa·s，γ=1.667，Pr=0.667，R=208.16 J/(kg·K)，逆幂律幂次ε=0.734，UGKS计算模型采用简单硬球模型。模型的训练数据来自Kn=0.73与Kn=1.09这2组状态，网格大小为33×91(周向×径向)，圆柱直径为0.304 8 m，第1层网格间距为1×-3m。

第1组待预测稀薄非平衡状态为Kn=0.91同尺寸二维圆柱绕流，除来流密努森数变化外，其余气体属性相同，计算网格如图6所示。图7与图8给出了不同计算方法温度云图对比,DNCR-RI(DNCR-Rotarional Invariants)为当前采用旋转不变量的改进DNCR方法。为了更好的表现N-S、UGKS、改进DNCR流场结果的差异，截取图中网格对称轴y=-x截线(45°截线)，对比截线上各项流场物理量的结果，如图9～图12所示。

图6 直径0.304 8 m二维圆柱网格示意图

图7 N-S与UGKS圆柱温度对比

图8 DNCR-RI与UGKS圆柱温度对比

由计算结果可以看出，在来流Kn=0.91条件下，圆柱绕流头部弓形激波附近区域表现出明显稀薄非平衡特征，激波厚度由于平均分子自由程增大明显变大。与UGKS与DNCR方法相比，N-S方程预示得到的激波位置更加靠近物面，激波脱体距离更小，表明N-S方程的连续性假设在这个来流条件下部分失效，无法对稀薄非平衡流场进行准确描述，而UGKS与DNCR流场结构较为一致。对比N-S与UGKS计算结果，采用精度提升式(13)分别考察DNCR计算精度，可以得到压力精度提升49.09%，温度精度提升61.02%，x轴方向速度U精度提升43.85%，y轴方向速度V精度提升48.26%，表明现有DNCR方法在相同外形、不同Kn数下有较好的预测表现。

图9 DNCR与UGKS圆柱截线y=-x上压力对比

图10 DNCR与UGKS圆柱截线y=-x上温度对比

图11 DNCR与UGKS圆柱截线y=-x上速度U对比

图12 DNCR与UGKS圆柱截线y=-x上速度V对比

图13和图14反映了N-S、UGKS、DNCR 3种计算方式得到的摩擦阻力系数Cf和表面热流对比，可以看到DNCR方法得到的结果与UGKS十分吻合，表明DNCR方法在气动热和摩阻预测方面也呈现了较好的预测能力。

为了测试改进DNCR方法模型的泛化性能，针对外形改变的情况进行了预测和结果对比。第2组待预测状态采用Kn=0.91，外形为长轴0.609 6 m、短轴0.304 8 m的二维椭圆柱，其余各项设置相同，计算网格如图15所示，重点考察进DNCR方法模型在变外形条件下泛化性能提升情况。图16和图17给出了不同计算方法温度云图对比。

图13 圆柱表面沿x方向摩擦阻力系数

图18～图21中DNCR-RV(DNCR-Rotational Variants)为原有采用旋转可变物理量的DNCR方法，可以明显看出改进DNCR方法得到的结果趋势与UGKS结果更吻合。

若采用精度提升式(13)进行量化分析，改进DNCR方法对压力的精度提升为-10.01%，对温度的精度提升为71.19%，速度U的精度提升为63.80%，速度V的精度提升为53.36%。而原始DNCR方法的压力精度提升为-90.52%，温度精度提升为71.14%，速度U精度提升为57.07%，速度V精度提升为-42.76%。通过精度提升对比可以发现，当待预示几何外形发生变化时，原始DNCR方法以含有空间坐标系方向信息的数据作为训练集，得到的预测结果精度提升不如采用基于旋转不变量改进DNCR方法，尤其是对诸如压力与y方向速度V的预示。计算结果初步表明，改进DNCR方法拥有更强的模型泛化性能，针对变几何构型的流场预示能力明显增强。

图14 圆柱表面沿x方向热流分布

图15 长轴0.609 6 m、短轴0.304 8 m二维椭圆柱网格示意图

图16 N-S与UGKS椭圆柱温度对比示意图

图17 DNCR-RI与UGKS椭圆柱温度对比示意图

图18 椭圆柱流场截线y=-x上压力对比

图19 椭圆柱流场y=-x截线温度对比

图20 椭圆柱流场截线y=-x速度对比

图21 椭圆柱流场y=-x截线速度对比

3.2 顶盖方腔驱动流

为进一步验证改进DNCR方法的计算能力与泛化性能，对二维顶盖方腔驱动流进行了计算与分析。计算采用单原子气体氩气，初始流场温度T0=300 K，气体物性参数取μ0=2.272×10-5Pa·s，γ=1.667，Pr=0.667，R=208.16 J/(kg·K)，逆幂律幂次取ε=0.81，UGKS计算模型采用简单硬球模型。模型的训练数据来自Kn=0.003 2、Re=44.0与Kn=0.002 375、Re=59.3这2组状态，物理计算网格为61×61的均匀网格，速度离散空间网格为28×28均匀分布，计算域尺寸为1 m×1 m。

第1组待预测状态Kn=0.002 85、Re=49.4，几何尺寸保持不变，同为1 m×1 m的方腔，其余各项设置相同，计算网格如图22所示。图23～图26分别展示了x=0.5 m截线处压力、温度、速度U、速度V的对比。

图22 1 m×1 m方腔示意图 Fig.22 Schematic diagram of a 1 m×1 m squarecavity

图23 1 m×1 m方腔截线x=0.5 m上压力

图24 1 m×1 m方腔截线x=0.5 m上温度

图25 1 m×1 m方腔截线x=0.5 m上速度U

图26 1 m×1 m方腔截线x=0.5 m上速度V

为量化DNCR计算精度提升情况，借助精度提升计算式(13)，得到x=0.5 m处压力精度提升68.15%，温度精度提升55.89%，速度U精度提升31.74%，速度V精度提升-20.38%，速度V截线处绝对值较小，由于特征值和标记值不涉及方向，所以预示结果较差，而绝对值较大的速度U结果明显较好，当地马赫数精度提升了31.66%，受速度V影响较小。图27～图30给出了方腔流场截线y=0.5 m上流场物理量分布曲线。

图27 1 m×1 m方腔截线y=0.5 m上压力

图28 1 m×1 m方腔截线y=0.5 m上温度

图29 1 m×1 m方腔截线y=0.5 m上速度U

如图27～图30所示，y=0.5 m截线上压力精度提升77.47%，温度精度提升58.66%，速度U精度提升69.43%，速度V精度提升44.49%，整体预示结果较好。

图30 1 m×1 m方腔截线y=0.5 m上速度V

为了考察方法的泛化性能，对变外形情形进行了预测和对比。第2组待预测状态为Kn=0.002 85，Re=49.4，计算几何外形为2 m×1 m的方腔，其余各项设置相同，网格如图31所示。截取相对位置与第1组算例相同，分别给出x=1 m、y=0.5 m处的压力、温度、速度U、速度V结果，其中x=1 m截线计算结果如图32～图35所示。

图31 2 m×1 m方腔示意图

图32 2 m×1 m方腔截线x=1 m上压力

图33 2 m×1 m方腔截线x=1 m上温度

图34 2 m×1 m方腔截线x=1 m上速度U

在x=1 m处压力精度提升80.32%，温度精度提升45.12%，速度U精度提升36.57%，速度V精度提升39.65%。y=0.5 m的截线结果如图36～图39所示。

在y=0.5 m处压力精度提升54.63%，温度精度提升4.49%，速度U精度提升42.49%，速度V精度提升3.72%。

图35 2 m×1 m方腔截线x=1 m上速度V

图36 2 m×1 m方腔截线y=0.5 m上压力

图37 2 m×1 m方腔截线y=0.5 m上温度

图38 2 m×1 m方腔截线y=0.5 m上速度U

图39 2 m×1 m方腔截线y=0.5 m上速度V

通过对固定坐标系不同状态下二维圆柱与方腔算例的验证，从流场关键位置截线对比压力、温度、速度分布可以看出：改进DNCR方法对不同来流状态、训练集相同外形、变几何外形的流场参数均能作出较好的预测，且当外形发生改变时，预测效果优于原有DNCR方法，表明改进模型拥有更好的泛化性能。更为重要的是，若坐标系发生旋转时，改进DNCR方法的流场结果与提升效果保持不变，原始DNCR方法则囿于特征值和标记值不具备旋转不变性而无法进行预测，因此未在文中进行量化比较，这也从另一方面反映了改进DNCR特征量具备旋转不变性的重要意义。

4 结 论

本文在原始DNCR方法的基础上发展了基于旋转不变量的改进DNCR方法，对机器学习模型的特征值和标记值进行了重新筛选与处理，使其满足旋转不变性约束，最终使得训练模型的预示能力与范围得到增强，解决了原始DNCR方法不适用于旋转坐标后的流场精确预示这一难题。改进DNCR方法的计算流程与原始方法一致，即将修正应力与热流预测结果代入守恒方程进行迭代并计算至收敛，在模型训练完成后，DNCR总计算时长等于N-S方程计算时长+模型预测时长+DNCR续算时长，预测时长相较N-S与DNCR计算时长可忽略不计，因此总计算时长与N-S方程处于同一量级，但计算精度较N-S方程得到大幅提升。以固定坐标系下的氩气二维高超声速圆柱绕流和顶盖方腔驱动流为例，针对原有外形和变几何外形的待预示状态，将改进的DNCR方法结果与N-S方程、UGKS方法及原有DNCR方法进行了对比，结果表明使用旋转不变量能够显著提升训练模型对外形变化的适应能力，使其表现出更好的泛化性能，进一步提高了DNCR方法的适用场景与应用范围。