基于深度神经网络的高阶非线性激波判别式

2023-09-02陈树生高正红向星皓

空气动力学学报 2023年7期

杨华，陈树生，2，*，李猛，庞超，高正红，向星皓

(1.西北工业大学航空学院，西安 710072；2.中航工业第一飞机设计研究院，西安 710089；3.空气动力学国家重点实验室，绵阳 621000)

0 引言

激波是广泛存在于可压缩流动中的重要现象之一，对飞行器的气动特性影响巨大。计算流体力学（Computational Fluid Dynamics，CFD）作为可压缩流动问题研究的主要手段之一[1]，近几十年来，其数值方法主要是围绕如何准确捕捉激波以及如何能够在光滑区域满足高精度的要求中不断发展的[2-3]。

对存在激波、接触间断等现象的流场进行数值计算时，为了保证数值格式的非线性、稳定性，同时保持较高分辨率，可以采用一种高分辨率混合格式[4]，该混合方法将非线性激波捕捉格式与针对光滑流动区域的高分辨率线性格式结合，能准确且高效地识别激波位置，在包含间断的流场中清晰地捕捉激波，最终得到稳定的计算结果。为了提高混合格式的计算效率和求解精度，有必要发展更加精准的激波识别方法。传统的激波识别方法均为人工构造的判别式且多含有经验参数[5]，因此在实际应用中存在一定局限性。近年来，机器学习在偏微分方程（partial differential equation,PDE）、反问题等领域的应用日渐增多。借助机器学习方法构建新型间断识别器可以有效避免对人工经验参数的依赖，具有巨大潜力。

Ray 等[6]构建了基于多层感知机（multilayer perceptron,MLP）的间断识别器并将其作为黑箱，在采用Runge-Kutta 间断 Galerkin（Runge-Kutta discontinuous Galerkin,RKDG）格式框架下，将其与minmod 限制器和基于minmod 的TVB 限制器进行比较，数值结果表明MLP 识别器具有更好的鲁棒性，并且不依赖于问题相关的参数。随后Ray 等[7]在此基础上，将该方法扩展到二维非结构网格中检测间断位置，并缩减了人工神经网络的结构，使其具有更高的计算效率。Abgrall 等[8]利用RKDG 方法和高阶限制器的模拟数据训练神经网络，并将其作为激波识别器；数值结果表明其在笛卡尔网格以及非结构网格上中均有较高的精度。Sun 等[9]构建了基于卷积神经网络的通用间断检测器，在间断处采用高阶WENO 格式，光滑处采用高阶中心差分格式，求解一维和二维问题的数值结果表明该检测器能够更加清晰地识别出间断位置。

在前人的基础上，本文基于深度前馈神经网络模型构建了高阶非线性激波判别式，与Huynh[10]提出的二阶无振荡格式（centered non-oscillatory scheme of second order,Cn2）中的激波判别方法进行了对比，并将其应用于二维超声速无黏圆柱绕流算例中。

1 训练及测试数据集的生成

本文构建的基于神经网络的高阶非线性激波判别式流程如图1 所示，分别是生成训练集和验证集、训练模型以及验证。

图1 基于神经网络的激波判别式建模流程Fig.1 Process of shock discriminant modeling based on Neural Network

生成训练数据集与测试数据集的步骤如下：

Step1：选择一个函数u(x)，在区间 [a,b]内包含其规律性的完整信息，如选择正弦函数、阶跃函数等。

Step2：选取任意点xi∈[a,b]及网格大小h，确保a≤xi-5/2·h＜xi＜xi+5/2·h≤b。以xi为中心，在区间 [a,b]内创建5 个网格单元，如图2 所示。

图2 在区间[a,b]内随机抽取网格点示意图Fig.2 Random extraction of grid points within an interval [a,b]

Step4：对应于输入向量X的实际输出取决于选取的函数性质及规律。理想情况下，如果网格单元Ii处的解不连续，则该网格单元会被标记为间断点，输出为Y=(1,0)T，否则输出为Y=(0,1)T。

Step5：重复2～4 步，每次重复需改变点xi的位置，网格大小h。

Step6：重复1～5 步，每次重复使用不同的函数u(x)。

需要注意的是，为了使得所生成样本更符合物理现象，在生成训练集与测试集中的间断函数数据时，默认参数a、b保持相同的正负号。

根据上述所列步骤，分别选取三角函数、一次函数、二次函数和间断函数等生成训练数据集T 与测试数据集V，如表1 所示。

表1 用于模型训练和测试的数据集Table 1 Data sets for model training and testing

2 神经网络模型

本节采用深度前馈神经网络（feedforward neural network,FNN）[12]构建激波识别器并开展神经网络模型的训练及验证，最后选取预测精度最高的模型作为最终实际使用的模型。

2.1 深度前馈神经网络模型

深度前馈神经网络也叫做前馈神经网络或者多层感知机，是典型的深度学习模型。前馈网络的目标是近似某个函数f*，通过定义一个映射y=f(x,θ)，并且学习参数 θ的值，使它能够得到最佳的函数近似。

在前馈神经网络中，各神经元分别属于不同的层，每一层的神经元可以接收前一层神经元的信号，并产生信号输出到下一层。第0 层称为输入层，最后一层为输出层，其他中间层称为隐藏层。整个网络中没有反馈，信号从输入层向输出层单向传播。其中网络的层数和各隐藏层的单元个数均为超参数，通常超参数是按照人为经验设定，并不断试错调整最终确定[12]。

本文所构建的深度前馈神经网络模型，选择7 层神经网络，输入层为7 个神经元，输出层为2 个神经元，中间有5 个隐藏层，其神经元个数分别为256、128、64、32、16。此网络结构中的参数还包括W={W1,W2,W3,W4,W5,W6} 和b={b1,b2,b3,b4,b5,b6}。其中W1∈R256×7、W2∈R128×256、W3∈R64×128、W4∈R32×64、W5∈R16×32、W6∈R2×16，b1∈R256、b2∈R128、b3∈R64、b4∈R32、b5∈R16、b6∈R2。激活函数选择LeakyReLU 函数。网络结构如图3 所示。

图3 FNN 模型网络结构示意图Fig.3 FNN model network structure diagram

2.2 训练神经网络

本文选择交叉熵损失作为模型训练的损失函数：

式（1）中，S是一次训练所使用的样本数；为所给样本第k个网格单元是否为间断点的实际标签，其值只有0 和1 之分；分别为模型预测该样本第k个网格单元为间断点或光滑点的概率，其中。通过交叉熵损失可以计算出两个概率分布与实际值之间的误差，最小化训练样本中预测值与真实值之间的误差，从而达到提高模型预测精度的目的。

模型训练的优化算法选择Adam 算法，此外为了防止模型过拟合，提高泛化能力，同时应用正则化与提前终止方法。训练的残差变化如图4 所示。

图4 模型训练残差图Fig.4 Model training residuals

在实际训练中，为了避免标准梯度下降中存在的无法充分利用计算机并行计算能力的问题，使用小批量梯度下降法减少参数更新的波动并节省计算资源，每一小批量选择200 个训练样本以实现效果更好更稳定的收敛。其中经过每一小批训练样本后的测试精度如图5 所示，选取其中测试精度最高的模型参数作为此次训练的结果。

图5 模型训练中的测试精度变化图Fig.5 Test accuracy variation in model training

模型训练中选择重复训练参数R=10，即按照上述流程重复训练10 次模型，选取其中测试数据集预测精度最高的参数作为最终模型参数。

3 模型验证

本节采用多个算例（Lax 激波问题、Sod 激波管问题及二维超声速无黏圆柱绕流算例）对FNN 模型进行验证，检验所训练模型的预测精度及其适用性。

3.1 一维算例

本小节首先将深度前馈神经网络模型（FNN）和原Cn2 格式[10,13]中的激波识别方法（此处简称为Cn2）进行对比，通过上述生成的训练和测试数据集检验两种激波判别式的识别效果。接着将FNN 模型嵌入Cn2 格式中以识别激波位置，并形成新的数值格式（此处简称为FNN 格式），将该格式与原Cn2 格式进行对比。通过求解一维的Lax 激波问题以及Sod 激波管问题[14]以验证FNN 模型在数值格式中的有效性。

原Cn2 格式中的激波识别方法[10]如下：

以密度场为例，令：

定义该网格单元Ij为光滑区域，否则标记为间断点。

3.1.1 训练和测试数据集

本文生成的训练数据集包含20 000 个光滑样本和20 000 个间断样本，测试数据集包含12 000 个光滑样本和12 000 个间断样本。应用FNN 和Cn2 激波识别方法对训练数据集和测试数据集中的光滑样本和间断样本分别进行预测，并对比两种激波识别方法的预测精度，如表2 所示。

表2 FNN 和Cn2 激波识别方法的预测精度对比Table 2 Comparison of predictive accuracy between FNN and Cn2 shock wave identification methods

一般神经网络模型在训练过程中，各类样本的权重是相等的，从而使得模型在对各类样本进行预测时保持预测精度基本一致。但针对激波识别问题，为保证数值求解的稳定性，通常要求在所有可能出现大梯度或者间断的流场处引入数值耗散来抑制数值非物理振荡；同时为保证流场模拟的精度，通常希望在光滑区域采用高分辨率低耗散格式。因此在FNN 模型训练过程中，经验性地设置光滑样本和间断样本的权重分别为[0.04,1.0]，在保证神经网络模型将所有间断样本都能识别出来的情况下，尽可能提高其对光滑样本的预测精度。

从表1 可以看到，通过在模型训练过程中给定各个类别样本不同的权重，放宽对光滑样本的预测精度，从而使得模型在预测间断样本时的预测精度更高。相较于Cn2 激波识别方法，FNN 模型提高了在光滑样本处的预测精度，并保持了正确识别间断样本的特性。

3.1.2 Lax 激波问题

计算域为[0,1]，计算网格数为100，初始条件为：

计算终止时间t=1.3。

图6 中给出了FNN 与Cn2 格式求解Lax 激波问题时的无量纲密度在计算域[0,1]的分布曲线，其中红色圆点为1 表示识别的激波位置。可以看到，FNN 格式分别将两个间断附近的7 个网格点和6 个网格点识别为激波区域，Cn2 格式分别将两个间断附近的10 个网格点和7 个网格点识别为激波区域。相较而言，FNN 格式识别出的激波位置更加精准，Cn2 格式在光滑区域会产生更多误判，这和3.1.1小节从训练及测试数据集中得到的结论是一致的。

图6 Lax 问题密度分布曲线Fig.6 Density distributions of Lax’s problem

3.1.3 Sod 激波管问题

计算域为[0,1]，计算网格数为100，初始条件为：

计算终止时间为t=0.2。

图7 中给出了FNN 与Cn2 格式求解Sod 激波管问题时的无量纲密度在计算域[0,1]的分布曲线。可以看到，FNN 格式分别将两个间断附近的4 个网格点识别为激波区域，Cn2 格式分别将两个间断附近的8 个网格点和6 个网格点识别为激波区域。两种激波识别方法均能很好地识别激波位置，但FNN 格式相较而言更加精准。

图7 Sod 问题密度分布曲线Fig.7 Density distributions of Sod’s problem

3.1.4 一维算例结果讨论

表3 为FNN 和Cn2 格式在求解上述Lax 激波问题和Sod 激波管问题时的耗时情况，选择CPU_TIME函数估计两种格式在同一计算平台上的计算时间。

表3 运用FNN 和Cn2 格式的计算时间对比Table 3 Comparison of computing time between FNN and Cn2 schemes

从表3 的对比可以看到，相对于Cn2 格式，FNN格式在求解Lax 激波问题时，消耗时间增加了22%；在求解Sod 激波管问题时，消耗时间增加了25%。由于FNN 模型在每个网格点处判别其是否为间断时，均需要重新计算498 个神经元的值才能给出判别结果，而Cn2 格式中激波识别方法为一次代数表达式，相较而言其消耗时间更少。

对于本文提到的Lax 和Sod 问题，FNN 模型所增加的时间消耗较少，完全可以忽略不计。对于二维或三维复杂流场数值模拟，由于流场网格数量激增以及可能出现激波反射、激波干扰等现象，采用FNN 方法相较于Cn2 方法会产生更大的时间代价；当时间代价达到一定量级时，可以使用多步判断以节省数值模拟时间。此外，可以考虑修改神经网络模型结构，通过适当减少网络层数及神经元个数，以减小程序运行的时间消耗。

3.2 二维超声速无黏圆柱绕流算例

本小节采用超声速无黏圆柱绕流[15-16]（不同来流马赫数、不同网格分布）来验证深度前馈神经网络模型（FNN）在二维真实流场中判别激波间断的有效性。数值计算方法采用Roe 格式加剪切黏性[17]，剪切黏性激活函数采用本文的FNN 激波识别方法。

图8 为来流马赫数Ma分别为5、10、15、20 时二维圆柱流场示意图。其中上半部为压力云图，下半部为使用FNN 模型识别的流场激波结果。由图中可以看出，训练得到的FNN 模型对于不同来流马赫数的流场都有很好的激波识别能力。其中网格采用161 ×161（周向×壁面法向）均匀分布网格。

图8 不同马赫数下超声速二维圆柱绕流流场Fig.8 Flow field around a supersonic two-dimensional cylinder at different Mach numbers

图9 为不同网格分布下超声速二维圆柱计算所得的流场图（Ma=20）。其中网格所示的第一个参数表示周向网格数，第二个参数表示壁面法向网格数，从左至右、从上至下网格长宽比分别为8、1/8、4、1/4、4、1/4、2、1/2、1。可以看到不同网格分布下，FNN 模型都能准确地识别出激波位置。其中图9（i）为161 ×161（周向×壁面法向）的均匀分布网格，其识别效果最好。

图9 不同网格分布下超声速二维圆柱绕流流场（Ma=20）Fig.9 Flow field around a supersonic two-dimensional circular cylinder with different grid distribution (Ma=20)

图9（b、d、f、h）分别对应长宽比为1/8、1/4、1/4、1/2 的网格计算结果。使用FNN 模型在中段正激波处识别激波位置较为精确，在末段斜激波处识别范围较大。这是由于壁面法向网格布点较多而周向网格布点较少，导致在壁面法向垂直方向的激波分辨率较高而在周向垂直方向激波分辨率较差，在末段斜激波附近即体现为所识别激波范围较大，因此激波识别精度依赖于网格分辨率。