ROE-Scheme Physics-Augmented Graph Neural Networks in Solving Eulerian and Laminar Flow Incompressible NS Equations
-
摘要: 近年来，融合物理信息的深度学习方法为偏微分方程的求解提供了一个新的思路. 然而，到目前为止，大多数工作在解空间存在间断的问题上的计算精度不高，时间外推能力差. 针对以上两个问题，该文提出了使用图神经网络结合流体计算领域的ROE格式融合方程或数据信息的模型——ROE-PIGNN. 数值实验表明，该模型在求解由Euler方程控制的激波管问题时，可达到与传统ROE格式相当的计算精度，并具备一定时间范围的外推能力. 最后，对由Navier-Stokes(NS)方程控制的二维圆柱绕流问题进行了求解，实验结果表明：模型可以预测后续的周期性流动，并实现对部分关键位置流动结构的更精确的复现，相比纯数据驱动，误差降低了60%.Abstract: In recent years, the deep learning method incorporating physical information provided a new idea for solving partial differential equations. However, most of the studies so far has low computational accuracy and poor time extrapolation for problems with discontinuities in the solution space. To address the above 2 problems, the ROE-PIGNN model was proposed for fusing equations or data information with the graph neural networks and the ROE scheme in computational fluid dynamics. Numerical experiments show that, the model achieves a computational accuracy comparable to that of the ROE scheme in solving the shock tube problem controlled by the Eulerian equation, and has the ability of extrapolation over a certain time range. Finally, the 2D cylindrical bypass flow traditional problem controlled by the Navier-Stokes (NS) equations was solved. The experimental results show that, the model can predict the subsequent periodic flow and reproduce the flow structure more accurately at some key positions, with an error reduction of 60% compared to the purely data-driven approach.
-
Key words:
- deep learning /
- graph neural network /
- ROE scheme /
- Eulerian equation /
- NS equation
-
0. 引言
在自然科学和工程领域中，偏微分方程(partial differential equation, PDE)扮演着至关重要的角色，能够描述各种现象，具有较大的研究价值. 被广泛研究的PDE包括Navier-Stokes(NS)方程、Schrödinger方程、Maxwell方程等，但这些方程的求解仍然具有挑战性. 以NS方程为例，NS方程只在少数简单的流动上有理论解，大部分问题的求解仍需借助传统的数值求解. 但数值求解又面临求解效率和精度的矛盾.
深度学习(deep learning, DL)为数值计算等传统研究手段提供了强有力的补充或替代. Raissi等[1]提出了著名的物理信息神经网络(physics-informed neural network, PINN)来处理与PDE相关的正问题和逆问题，后续又涌现了大量与PINN相关的工作[2-4]. 相较于传统数值求解，神经网络求解PDE时无需了解方程本身的性质，无需选择离散算法，借助机器学习的算法可以很容易地实现并行计算和自动化调参. 另外，神经网络能够融合方程和各种数据，实现物理和数据双驱动的高精度求解.
但是，目前所用到的求解方法大多适用于较为简单的PDE. 众所周知，即使初始条件(initial condition, IC)是光滑的，Euler和NS方程的解也会出现不连续性[5-6]. 针对这类方程已有大量基于传统数值方法的工作，例如文献[7-9]和其中的文献. 其中，一些工作是将深度神经网络(deep neural network, DNN)的监督学习与传统的数值方法相结合[10-13]. Liu等[14]揭示了PINN在求解强非线性间断问题时存在的一种限制迭代收敛的悖论状态，并提出了一套新的间断问题求解框架，用以提升PINN类方法的间断求解能力. 通过“过渡点”压缩进光滑区，而挤出“间断”的思路，具有一定的参考价值. Cao等[15]将计算流体力学(computational fluid dynamics，CFD)中的结构化网格变换技术与PINN相结合，使用神经网络学习计算空间而不是物理空间中的流场，成功地解决了无黏翼型绕流的精确高效计算问题，并初步展示了多状态批量化求解的优势.
用DNN求解PDE的研究[16]表明，由自动微分定义的时空函数逼近器对于逼近非线性问题是有效的，例如Burgers方程、热传导方程以及NS方程. 同时，研究人员利用PINN通过数据驱动方法计算反问题，可以实现对具有高度不连续过渡的非线性双曲型方程组的求解，如Riemann问题. 然而，在未知物理场中何处存在间断解的情况下，使用DNN解决正问题仍然是一个挑战. 在本研究中，我们将通量差分分裂格式中的ROE数值格式结合到图神经网络(graph neural network, GNN)中，建立了损失函数逼近器以替代传统深度学习框架的自动微分，从而提高了捕获Riemann问题中间断解的精度.
对于可能包含激波和接触间断的守恒定律的数值方法的研究已有很多工作，如ROE格式[17]、Godunov格式[18]、MUSCL格式[6, 19]、ENO格式[20-21]和WENO格式[21-22]、分层重建[23-25]等. 数值技术是研究高速空气动力流动的一种重要研究方法，高速空气动力流动在飞机设计、燃烧问题和天文物理学中发挥着重要作用[26-29]. 然而，目前融合各种成熟的数值求解技术与机器学习技术的研究仍处于初级阶段，充分利用两者各自的优势有望在CFD中发挥更大的价值.
有部分工作是基于神经网络和加权基本无振荡算法. Bezgin等[13]在WENO方案中使用神经网络作为加权函数，解决了这些缺点. Bar-Sinai等[30]介绍了数据驱动离散化，这是一种基于已知基本方程的实际解学习PDE优化近似的方法. 该方法使用神经网络来估计空间导数，这些导数是端到端优化的，以最好地满足低分辨率网格上的方程. 以上两个工作是数值方法与神经网络的深度结合，但是在求解方程时依赖数值求解器. Magiera等[10]提出两种策略来生成约束感知神经网络，通过施加(强或弱)物理约束提高模拟的准确性. 在处理复杂问题时，计算效率也有很大提高，例如在多组分流、涉及动力学子模型的分散流或燃烧问题/化学反应问题等方面. 高普阳等[31]基于卷积神经网络模型，对传统有限体积法格式中的权重系数进行优化，以得到在粗粒度网格下具有较高精度的新数值格式，从而更适用于复杂问题的求解. 该网络模型可以准确、有效地求解Burgers方程和level set方程，数值结果稳定, 且具有较高数值精度. 这部分的工作完成了对存在间断问题的高精度求解，但是模型不具备时间外推能力.
我们期望借助GNN的局部信息归纳推理与数值离散技术，实现常用流体力学方程的求解与时间泛化. 在我们以前的研究中[32-33]，通过最小二乘方法和有限差分求解空间微分算子，据此计算PIGNN(physics-informed graph neural network)的损失函数来求解Burgers方程、传热方程等，但是在求解Euler和NS方程时无法得到可接受的近似解. 为了解决这个问题，我们以ROE格式代替最小二乘法的计算梯度. 模型求得的预测解与从相应数值计算中获得的解相当，本文中称为ROE-PIGNN. 相比于最小二乘法，ROE格式可以消除间断附近的振荡，实现更为精确的计算，误差降低了80%. 另外，我们研究了融合数据信息的NS方程的求解，结果表明，引入部分数据约束后，降低了模型的训练难度，并且具备更优的泛化能力，相比于只借助数据信息，可以复现更为精确和丰富的流动细节，误差降低了60%.
1. 方程与数值格式
1.1 NS方程与Euler方程
二维可压缩流的质量守恒、动量守恒和能量守恒可以用NS方程来描述，方程可以写成以下守恒形式：
$$ \frac{\partial \boldsymbol{U}}{\partial t}+\frac{\partial \boldsymbol{F}_1(\boldsymbol{U})}{\partial x}+\frac{\partial \boldsymbol{F}_2(\boldsymbol{U})}{\partial y}=\frac{\partial \boldsymbol{G}_1(\boldsymbol{U})}{\partial x}+\frac{\partial \boldsymbol{G}_2(\boldsymbol{U})}{\partial y}, $$ (1) 其中U =[ρ, ρu, ρv, E]T，ρ是密度，p是压力，u, v是二维x和y方向的速度分量，E是总能量；F1(U)=[ρu, ρu2+p, ρuv, u(E+p)]T，F2(U)=[ρv, ρuv, ρv2+p, v(E+p)]T，$\boldsymbol{G}_1(\boldsymbol{U})=\left[0, \tau_{11}, \tau_{12}, k \frac{\partial T}{\partial x}+u \tau_{11}+v \tau_{12}\right]^{\mathrm{T}}$，$\boldsymbol{G}_2(\boldsymbol{U})=\left[0, \tau_{21}, \tau_{22}, k \frac{\partial T}{\partial y}+u \tau_{21}+v \tau_{22}+w \tau_{23}\right]$. 要封闭这些方程，我们还需要一个方程，即描述压力和能量关系的状态方程，由方程(2)给出：
$$ E=\frac{p}{\gamma-1}+\frac{1}{2} \rho\left(u^2+v^2\right), $$ (2) 其中γ=1.4是绝热指数.
若黏性项全为零，则方程变为Euler方程. 守恒形式的各种数值方法研究较多，早期的尝试主要集中在激波捕获方法. 最近，已经发展了高阶方法，如间断Galerkin方法、基本无振荡方法及其变体、加权基本无振荡法、谱差分法等.
对于低速流动问题，二维不可压NS方程的形式如下：
$$ \left\{\begin{array}{l} \nabla \cdot \boldsymbol{V}=0, \\ \frac{\partial \boldsymbol{V}}{\partial t}+\boldsymbol{V} \cdot \nabla \boldsymbol{V}=-\nabla p+\frac{1}{Re} \nabla^2 \boldsymbol{V}, \end{array}\right. $$ (3) 其中$\boldsymbol{V}=\left\{\begin{array}{l} u \\ v \end{array}\right\}$，Reynolds数$R e=\frac{\rho U L}{\mu}$，U表示参考速度. 对于非定常状态通常需要内迭代且效率很低，由于不可压的假设，缺乏确定压力分布的显式方程，给求解带来一定的困难. 目前求解低速不可压问题的主要方法有涡流函数法、虚拟压缩法和压力修正法等. 使用GNN结合方程损失可以避免这一困难，利用神经网络的泛化能力可以实现二维不可压NS方程的快速求解与预测.
1.2 ROE格式
ROE格式基于Godunov方法发展而来，通过对变量进行ROE平均，将原始的无黏通量Jacobi矩阵用交界面左右两侧流场参数构成的常数矩阵来代替，在交界面求解近似Riemann问题，从而提高了计算效率. ROE格式数值耗散小，精度高，对流场间断有很高的分辨率. 因此，在CFD数值求解中应用非常广泛，也是高阶精度格式通常采用的通量计算方法. 采用ROE格式，网格交界面处的无黏通量可以表示为
$$ \boldsymbol{F}_{j+1 / 2}=\frac{1}{2}\left[\boldsymbol{F}\left(\boldsymbol{U}_{\mathrm{L}}\right)+\boldsymbol{F}\left(\boldsymbol{U}_{\mathrm{R}}\right)\right]-\frac{1}{2}|\tilde{\boldsymbol{A}}|\left(\boldsymbol{U}_{\mathrm{R}}-\boldsymbol{U}_{\mathrm{L}}\right), $$ (4) 其中第一项是交界面左侧通量和右侧通量的平均，第二项可以认为是耗散项，$|\tilde{\boldsymbol{A}}|$是通过ROE平均构造的近似无黏通量Jacobi矩阵，QL和QR分别表示交界面左右两侧的流场参数，耗散项的具体形式为
$$ \left|\tilde{\boldsymbol{A}}\left(\boldsymbol{U}_{\mathrm{R}}, \boldsymbol{U}_{\mathrm{L}}\right)\right|=\boldsymbol{S}^{-1}|\boldsymbol{\varLambda}| \boldsymbol{S} . $$ (5) 以平均值代替$\tilde{\boldsymbol{A}}$，这里的平均值不是简单的算术平均而是ROE平均，具体计算公式为
$$ \left\{\begin{array}{l} \bar{\rho}=\left[\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{L}}} / 2\right]^2, \\ \bar{u}=\left(\sqrt{\rho_{\mathrm{L}}} u_{\mathrm{L}}+\sqrt{\rho_{\mathrm{R}}} u_{\mathrm{R}}\right) /\left(\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}\right), \\ \bar{H}=\left(\sqrt{\rho_{\mathrm{L}}} H_{\mathrm{L}}+\sqrt{\rho_{\mathrm{R}}} H_{\mathrm{R}}\right) /\left(\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}\right), \end{array}\right. $$ (6) 其他量可以用这三个量计算
$$ \bar{p}=\frac{\gamma-1}{\gamma}\left(\bar{\rho} \bar{H}-\frac{1}{2} \bar{\rho} \bar{u}^2\right), $$ (7) $$ \bar{c}^2=(\gamma-1)\left(\bar{H}-\frac{\bar{u}^2}{2}\right). $$ (8) 计算步骤如下：
已知n时刻所有网格点上的物理量，对于j点：
1) 利用差分格式计算UR, UL；
2) 采用ROE平均公式计算ROE平均值；
3) 将Jacobi矩阵进行特征分解：$\boldsymbol{A}(\bar{\boldsymbol{U}})=\boldsymbol{S}^{-1} \boldsymbol{\varLambda} \boldsymbol{S}$；
4) 计算$\boldsymbol{S}^{-1}, \boldsymbol{\varLambda}, \boldsymbol{S}$；
5) 计算$\left|\tilde{\boldsymbol{A}}\left(\boldsymbol{U}_{\mathrm{R}}, \boldsymbol{U}_{\mathrm{L}}\right)\right|=\boldsymbol{S}^{-1}|\boldsymbol{\boldsymbol{\varLambda}}| \boldsymbol{S}$，其中$|\boldsymbol{\varLambda}|=\operatorname{diag}\left(\left|\lambda_1\right|, \left|\lambda_2\right|, \left|\lambda_3\right|\right)$；
6) 计算$\boldsymbol{f}_{j+1 / 2}=\hat{\boldsymbol{f}}\left(\boldsymbol{U}_{\mathrm{R}}, \boldsymbol{U}_{\mathrm{L}}\right)=\frac{1}{2}\left[\boldsymbol{f}\left(\boldsymbol{U}_{\mathrm{R}}\right)+\boldsymbol{f}\left(\boldsymbol{U}_{\mathrm{L}}\right)\right]-\frac{1}{2}\left|\tilde{\boldsymbol{A}}\left(\boldsymbol{U}_{\mathrm{R}}, \boldsymbol{U}_{\mathrm{L}}\right)\right|\left(\boldsymbol{U}_{\mathrm{R}}-\boldsymbol{U}_{\mathrm{L}}\right)$；
7) 计算空间导数$\left.\frac{\partial \boldsymbol{f}(\boldsymbol{U})}{\partial x}\right|_j=\frac{1}{\Delta x}\left(\boldsymbol{f}_{j+1 / 2}-\boldsymbol{f}_{j-1 / 2}\right)$.
2. ROE-PIGNN模型
本文提出了一种新的ROE-PIGNN方法，以求解上述PDE系统. 给定初始条件、边界条件(boundary condition, BC)和PDE方程，应用ROE-PIGNN模型求得PDE的近似解. 基于PDE方程本身的残差损失和边界损失构造物理信息损失函数，整个学习过程无监督. 如有数据可利用，可再增加数据损失，进一步指导方程求解. PIGNN原理如图 1所示. 通过在模型中引入数值仿真数据、专门设计的损失函数、GNN结构和消息传递机制，将物理先验知识嵌入深度模型，以帮助我们的模型学习空间和局部的动态随机变化和推理关系，从而使模型具有一定时间范围的外推能力.
模型的运行过程如下：首先将方程转化为离散的形式，包括求解域、初始状态和边界条件. 初始状态作为模型第一个时间步的输入，然后输出第二个时间步的物理场. 第二个时间步的物理场再输入到模型中输出第三个时间步的物理场，以此类推输出后续的物理场. 为充分利用方程本身的信息，边界条件在每次的输出中可以作为损失加入到损失函数中，也可直接将边界条件赋值给预测的物理场. 损失函数最重要的一部分是求解域中方程的残差信息，各个方程的残差信息以取绝对值并加权的形式组成损失函数.
2.1 模型架构
模型输入为当前时间步t方程的解，输出为下一时间步t+Δt方程的解，我们借鉴了Gilmer等[34]提出的一种称为消息传递神经网络(message-passing neural network, MPNN)的GNN，以建模从上一时间步到下一时间步的解的变化. 具体地，ROE-PIGNN采用与Sanchez-Gonzalez等[35]相同的编码器-处理器-解码器架构，如图 2所示.
编码器(encoder)部分，网络将储存在网格点上的数据信息转化为图信息，并设置节点和边的特征，即Glt=fENCONDER(u(X, t)). 处理器(processor)建模节点和边之间的相互作用模式，并通过L个图网络块(graph network, GN)将每一步挖掘的节点和边信息不断传递. 运行L步生成网络图序列Gt=(G1t, …, GLt)，其中每一步的网络图都是基于前一GN块的输出来更新，即Glt=fGN(Gl－1t). 在每个GN模块中，分三步更新边特征和节点特征：首先，更新每条边的特征，然后将更新后的边特征聚合到当前节点上，最后使用聚合后的边特征更新节点特征. 数学上，首先更新边特征，即$\tilde{e}_{i j, (l)}=\phi^e\left(e_{i j, (l)}, v_{i, (l)}, v_{j, (l)}\right)$, 根据更新后的边特征聚合每个节点的边缘特征，即$\bar{e}_{i, (l)}=\bigoplus\limits_{j \in N(i)} \tilde{e}_{i j, (l)}$，更新节点特征$\tilde{v}_{i, (l)}=\phi^v\left(e_{i j, (l)}, v_{i, (l)}\right)$. 节点特征包含时间t处的节点值，以及一个节点类型向量，表明它是内部节点还是边界节点. 而边特征由相邻两个节点i和j的相对位移及其范数组成. 即节点特征$v_i=\left\{u(x, t), n_{\text {type }}\right\}$，边特征$e_{i j}=\left\{x_{i j}, \left\|x_{i j}\right\|\right\}$. 其中$\phi^e$，$\phi^v$分别是边和节点更新函数，它们被映射到所有边/节点上以计算每边/每节点更新.
解码器(decoder)从处理器输出的最终图GLt的节点特征vi, (L)中提取节点的动力学信息，即hi, (L)=fDECODER(vi, (L))，本文使用MLP作为解码器函数. 基于解码器所得动力学特征，使用向前Euler积分方式计算下一时间步的解u(Xi, t+Δt)=u(Xi, t)+hi, (L).
2.2 损失函数
为更好地在模型中嵌入物理信息，除了模型架构的设计之外，损失函数的设计也尤为重要. 常用的技术之一是将物理方程显式地加入模型的损失函数中，通过模型训练引导模型一方面满足PDE的约束，另一方面能够拟合数值模拟数据. 因此，物理信息损失函数由物理方程的残差损失和数据损失组成(若有观测数据)，具体形式由下式给出：
$$ L\left(\theta ; D_p, D_d\right)=(1-\gamma) L_{\mathrm{PDE}}\left(\theta ; D_p\right)+\gamma L_{\mathrm{data}}\left(\theta ; D_d\right), $$ (9) 其中γ表示加权系数，θ表示模型的可训练参数，$D_p=\left\{x_i^p, t_k\right\}_{i=1}^{N_p}\left(k=1, 2, \cdots, N_{t_k}\right)$为一系列配置点，$D_d=\left\{y\left(x_i^d, t_p\right)\right\}_{i=1}^{N_d}$是观测数据集. 在离散时间中，从内部网格节点和时刻点中选择配置点$\left\{x_i^p, t_k\right\}$，以此计算PDE损失项$L_{\mathrm{PDE}}\left(\theta ; D_p\right)$. 具体形式为$L_{\mathrm{PDE}}\left(\theta ; D_p\right)=E\left[\left\|R\left(x^p, t_k ; \theta\right)\right\|_2^2\right]$，其中$\|\cdot\|_2$代表l2范数，PDE残差为
$$ R\left(x_i^p, t_k ; \theta\right)=\frac{\partial u\left(x_i^p, t_k ; \theta\right)}{\partial t}+F\left[u, \nabla_x u, \Delta_x u, \cdots\right], $$ (10) 这里空间导数的计算方法为第1.2小节提到的ROE格式. 数据损失为观测值$y\left(x_i^d, t_p\right)$和ROE-PIGNN预测值$u_{\mathrm{pred}}\left(x_i^d, t_p ; \theta\right)$之间的均方误差，具体形式为$L_{\mathrm{data}}\left(\theta ; D_d\right)=E\left[\left\|y\left(x^d, t_p\right)\right\|_2^2\right]$.
3. Euler方程算例
3.1 Riemann问题
Euler方程的间断初值问题称为激波管问题，也叫Riemann问题. 该问题描述的是无黏理想气体受一根管道中压力或者速度驱动在突变初始条件下的变化过程. 一维激波管问题是为数不多的存在精确解的CFD问题，常常作为检测算例来检验数值格式. 除此之外，很多数值方法都包含了Riemann精确求解器或近似求解器，以提高捕捉激波和处理光滑区域的能力. 因此，对一维激波管问题进行学习和实现是最基本的也是具有现实意义的. 本文选取典型工况作为模型的验证算例：
$$ \left\{\begin{array}{l} \boldsymbol{W}_{\mathrm{L}}=\left[\begin{array}{lll} 1 & 0 & 1 \end{array}\right], \\ \boldsymbol{W}_{\mathrm{R}}=\left[\begin{array}{lll} 0.125 & 0 & 0.1 \end{array}\right], \end{array}\right. $$ (11) 其中$\boldsymbol{W}_{\mathrm{L}, \mathrm{R}}=\left[\rho_{\mathrm{L}, \mathrm{R}}, u_{\mathrm{L}, \mathrm{R}}, p_{\mathrm{L}, \mathrm{R}}\right]$.
在此工况下，右侧有一道激波冲出，激波扫过的地方，压强、速度、密度都开始增加. 左侧第一道波为膨胀波，膨胀波扫过的地方压强、密度均会下降(膨胀波内部物理量连续、光滑，头尾物理量连续(并非突变)，但导数不连续(弱间断)). 计算使用的网格为在x方向和y方向从0到1的范围内100等分，共计104网格.
模型的输入为当前时刻的流场数据，包括ρ, u, v, p和网格节点信息，输出为下一时刻的流场数据. 残差损失中的空间梯度采用两种计算方法，一种是最小二乘法[32]，一种是ROE格式. 最小二乘法计算简单，在没有激波间断的算例中表现良好，但是如果流场出现间断，会因为数值格式的不稳定导致数值振荡.
训练过程中选取的时间步长为Δt=2×10－3，训练时间步数为20，模型迭代优化步数为25 000步，最终四个方程均已收敛，能量方程的残差略高为10－2，连续方程残差为10－3. 如图 3(a)所示，红色虚线表示的初始场，20个时间步后到绿色线表示的位置，实线表示模型的预测结果，虚线表示精确解. 前20个时间步为模型的训练过程，后续时间步为模型的时间外推结果. 可以看出：膨胀波内预测较为准确，但在激波后压强产生了剧烈的波动，随着时间的推进波动进一步加大，到第40个时间步，振荡的幅值占比到60%. 图 3(b)、3(c)为二维云图结果，横纵坐标为网格边界节点编号，每一组图的左侧为准确值，右侧图为模型预测值. 从图 3云图中也可以看出，采用该方法预测得到的值与精确值的差异较大. 使用ROE格式代替最小二乘法计算模型损失，结果如图 4(a)所示. 结果显示激波后的振荡明显消除，在训练工况之外，也是就时间推到第60个时间步，采用ROE格式的模型预测精度更高. 从图 4(b)、4(c)也可看出PIGNN的计算值与精确值较接近. 为量化误差，定义误差：
$$ L=\sqrt{\sum\limits_{k=1}^N\left(u_k-\bar{u}_k\right)} / N. $$ 
图 3 SOD激波管问题中，以最小二乘法计算PDE损失的模型获得的压强p结果注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.Figure 3. Pressure results obtained for the model with PDE losses with the least squares method in the SOD shock tube problem
图 4 SOD激波管问题中，以ROE格式计算PDE损失的模型获得的压强p结果Figure 4. Pressure results obtained for the model with PDE losses in the ROE format for the SOD shock tube problem误差随时间的变化曲线如图 5所示，采用ROE格式后的模型相比于最小二乘法的模型误差降低约80%.
3.2 非均匀压力传播问题
我们在上一算例中验证了结合ROE格式的PIGNN在有间断的算例中的计算精度. 接下来将对非均匀压力场的传播问题进行计算，并在不同的工况下进行测试，以验证PIGNN的泛化能力，网格设置与前一算例一致. 具体问题为在计算区域中心位置，流场压强高于周围压强，速度为零，后续高压区域将不断向四周扩散. 如图 6所示，实线表示模型计算结果，虚线表示理论计算结果. 矩形区域中心位置处压强高，周围压强低，在压差的作用下，流体开始向外流动. 红色线表示的是中间截面上压强的分布，随着流动的进行在80个时间步后推进到绿色线位置. 在此算例中取80个时间步作为训练数据，时间步长取Δt=2×10-4.可以看出训练的结果与理论计算值吻合良好，理论值的计算采用的是显式时间推进，空间梯度的计算统一采用ROE格式. 在80步之后的模型计算结果属于外推结果. 模型可以预测后80个时间步，即预测到160步表现良好，结果如图 6中的黑色实线. 为查看流场其他位置的计算结果，图 6(b)、6(c)展示了压强p的云图分布. 与理论计算结果相比，两者较为吻合. 继续时间外推，模型的预测结果变差明显.
图 6 训练工况和时间外推下，模型计算压强场的结果Figure 6. Results of the calculated pressure field for the model under training conditions and time extrapolation改变初始压强场的分布，主要变化为增大了初始压强的峰值，计算结果如图 7所示. 模型在第80个时间步的结果与理论值有微小差异，继续外推到第160个时间步，两者差异加大，但是基本吻合. 说明模型除具备时间外推能力，还具备初始条件的外推. 误差曲线展示于图 8中，在新工况下，模型的平均误差小于0.06%，峰值误差小于2%.
图 7 新初始场下，模型计算压强场的结果Figure 7. Results of the calculated pressure field for the model with the new initial field我们将流体力学中一种基于Riemann近似解的通量差分分裂格式应用于PIGNN模型来求解Euler方程. 该格式以其数值耗散低、接触间断分辨率高、激波捕捉能力强等优点，广泛应用于CFD的计算中. 数值实验结果表明，应用了ROE格式的PIGNN计算精度高于采用最小二乘法的计算结果，并得出以下结论：差分格式的精度对存在激波间断的算例影响较大，差分格式精度越高，对间断的处理越好，PIGNN的预测精度越高；ROE-PIGNN求解Euler方程具备一定的时间泛化能力，时间预测步数约为训练步数的三倍；针对单一工况的训练可实现在附近工况上的预测，预备一定的初始条件泛化能力.
4. 不可压NS方程算例
NS方程比Euler方程多了黏性项，二阶导数的计算采用最小二乘法. 选取二维圆柱绕流作为验证模型能力的算例. 方程中的Reynolds数表征的是流体的惯性力和黏性力的关系，随着Reynolds数的增大，流动状态将由定常流动转为非定常流动. 随着Reynolds数的进一步增加，流动将由非定常层流向非定常湍流发展. 后续算例取Reynolds数Re=60，此状态为非定常层流流动. 计算区域和网格如图 9所示，圆柱表面分布100个网格节点，尾迹区加密，总网格共计8 089个节点，15 892个网格单元.
图 9 计算采用的网格的整体视图与局部视图Figure 9. Overall and local views of the grid used for the calculation对于此算例，我们在Reynolds数Re=60下以传统数值求解生成流场. 速度入口的值为10 m/s，出口边界为压力出口，上下面为有滑移壁面. 计算使用的时间步长为0.01 s，历经多个时间步后，流场进入周期性阶段. 取进入周期后的流场作为模型的初始输入. 全流场数据输入流场后输出下一时刻的流场数据，训练阶段提供方程残差信息，预测阶段无方程信息参与.
4.1 方程驱动的ROE-PIGNN
在此算例中，只考虑NS方程的损失项，和传统算法类似，不借助任何数据只用方程本身求解，属于无监督学习. 在此工况下，设定速度入口的值为10 m/s，流动的周期为80个时间步，即0.8 s. 模型训练只用PDE的损失，可以训练20步(1/4周期)，时间步长取1×10-2 s. 如图 10所示，每一组的上图为模型的预测结果，下图为传统数值求解的计算结果，流场的主要变化区域为尾迹区.
图 10 以PDE损失训练PIGNN在不同时刻的计算结果Figure 10. Computational results of training PIGNN with PDE losses at different moments预测结果说明，前文介绍的方法计算的梯度和二阶导数项与数值求解基本一致. 在前20个时间步，即训练步，速度误差在1 m/s以内，误差随时间的外推不断增大，在第50时间步误差小于2 m/s，继续外推误差增加更为明显，流场速度值逐渐发散. 由于空间和时间离散误差，因此两者的计算结果很难做到完全一致.
4.2 数据增强的ROE-PIGNN
4.2.1 小区域范围添加数据损失
相较于传统数值求解，ROE-PIGNN除了借助方程信息外，还可集成数据信息. 如式(9)所示，以加权的形式将数据损失纳入训练中，同时优化方程的残差损失和数据损失. 选取圆柱后方7倍直径长、2倍直径宽的区域作为数据损失区，即认为此区域流场数据已知(如图 11所示)，流场的其他位置使用方程损失. 在选择数据损失的加权因子时，根据接近收敛时比较两项损失的大小，调整加权因子使得两者在同一数量级，保证两者对最终结果的影响是一致的.
图 11 圆柱尾迹区添加数据损失的计算结果Figure 11. Calculated results after adding data losses to the cylindrical wake area计算结果如图 11所示(上图为预测结果，下图为数值模拟结果). 方程计算100步，步长1×10-2 s，在尾迹区计算结果较好，但是在下游10倍直径位置方程损失较大，在出口位置损失大.
4.2.2 全空间范围低时间频率添加数据损失
添加低时间频率的数据损失. 数据量在时间维度上降低到PDE损失的1/10，添加数据的区域为全流场，即间隔1×10-1 s计算一次数据损失.
模型计算100个时间步，训练30 000步后残差基本持平，不再下降. 计算结果与数值求解结果较为接近，模型可以实现300步的时间外推，计算结果如图 12所示(上图为预测结果，下图为数值模拟结果). 到第250步时，最大绝对误差小于上一算例. 但是平均误差略大于上一算例，因为误差主要出现在圆柱附近，此位置网格密度较大.
图 12 在少量的几个时刻添加数据损失的计算结果Figure 12. Calculated results with data losses added at a small number of moments通过前面算例的计算和分析可以看出，只用PDE的残差作损失外推效果较差. 这是因为方程残差只能反应两个相邻时刻的误差，若某一时刻的流场预测出现偏差，这种偏差会随之积累到下一时刻，即相当于方程的初始值发生了改变. 添加数据损失会将累积的误差进行修正，确保不会偏离原始推进方向.
图 13展示了利用全空间范围低时间频率的数据驱动模型与数据加方程损失驱动模型的计算和外推结果在圆柱尾迹区的误差. 前100步的最大训练误差保持在1 m/s以内，后续的时间外推结果的最大误差保持在2 m/s以内. 选取区域的x和y坐标范围为[2, 10], [-3, -1]，与数值求解的计算结果对比计算误差. 其中data表示模型训练只用数据损失；data+PDE表示添加数据损失和方程损失一块训练，即给上一模型提供了更丰富的训练信息. 模型B的误差要小于模型A，即模型在训练过程中结合方程损失可以得到更为精细的流动结构，误差降低了60%.
图 13 纯数据驱动和数据加方程驱动模型计算的尾迹区平均误差趋势Figure 13. Trends of mean errors in the tailrace area calculated by purely data-driven and data-plus-equation-driven models5. 结论
本文提出了一种新的基于物理信息的机器学习框架，该框架结合了CFD中常用空间离散格式方案以及GNN在求解PDE时获得的更高的精度和时间外推能力. 在没有数据可用时，模型可只借助PDE的方程残差来求解方程. 通过解决与实际问题相对应的激波管问题压力场传播等(Euler方程)，说明了此架构在求解PDE方面的精确性.
与传统的数值算法相比，ROE-PIGNN框架在求解PDE时，具有显著的优势. 首先，PIGNN框架极大地降低了求解PDE的门槛，无需用户深入了解方程本身复杂的物理原理和数学结构. 这是因为PIGNN框架通过结合物理知识和机器学习算法，可以自动学习方程的隐式特征和模式，从而实现对PDE的有效求解. 其次，PIGNN引入了机器学习的自动化特性，使得并行计算的实现更加简单和高效. 传统的并行计算往往需要程序员编写复杂的并行算法和调度代码，而在PIGNN中，由于其内置的并行处理机制，用户可以更轻松地利用多核处理器或分布式计算资源，从而大大提高计算效率. 此外，PIGNN框架还具备数据同化的能力. 数据同化是一种将模型预测与实际观测数据融合的技术，这使得模型能够利用有限的实验数据来修正和优化其预测，即使在数据不完整或存在噪声的情况下，也能提供更准确的结果.
与PINN相比我们的模型还具备外推能力，训练一定的步数后，可以预测后续的结果. 在Euler方程的算例中，预备一定的初始条件泛化能力，模型的平均误差小于0.06%，峰值误差小于2%. 在二维不可压NS方程控制的圆柱绕流问题的时间外推上，单纯的方程损失只能考虑到相邻两步的相对误差，会存在误差累积的问题. 只用方程损失作训练，训练时间不足一个周期，不能满足长时外推的需求，而结合少量的数据即可解决这一问题. 数据作为一种强约束可以将方程的误差即时修正，达到更长时间的稳定泛化. 经对比计算分析得出PIGNN可以实现不可压层流NS方程的数值求解，并具备时间外推能力，外推步数为训练步数的两倍. 给模型提供关键区域或少量几个时刻的数据信息增强PIGNN的泛化能力，实现圆柱绕流的周期性求解. 与单纯的数据驱动对比，结合方程损失后，可以实现NS方程的更高精度的求解和时间外推，误差降低了60%.
目前来看，我们的模型ROE-PIGNN还存在一些问题，除了训练慢，收敛难也是个问题. 因为PDE的残差作为损失函数中的约束，可以理解为是一种软约束. 使用梯度下降法求得的最优通常只是各个点互相“妥协”的最优. 自动微分是不存在截断误差的，因此可以比较精确地求出方程的损失. 然而，我们的模型在空间微分计算时引入了逼近误差，经过实验比较发现逼近误差的数量级大于截断误差. 如果想改进这一点，可能需要更先进的优化理论或空间微分逼近方法. 靠方程约束并不是万能的，在定义域空间内，有无数种解是满足方程的，只有初始条件和边界条件作为强约束下，该解才是唯一的，而神经网络方法并不能做到强约束. 因此若无内部“监督点”，会经常出现方程的损失很小，但解的形态明显不正确. 后续工作考虑引入模态分解等技术，在求解过程中去噪，增强模型的收敛性和稳定性. 模型的框架设计上，引入LSTM等序列结构，增强时间外推能力.
-
图 3 SOD激波管问题中，以最小二乘法计算PDE损失的模型获得的压强p结果
注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.
Figure 3. Pressure results obtained for the model with PDE losses with the least squares method in the SOD shock tube problem
图 4 SOD激波管问题中，以ROE格式计算PDE损失的模型获得的压强p结果
Figure 4. Pressure results obtained for the model with PDE losses in the ROE format for the SOD shock tube problem
图 6 训练工况和时间外推下，模型计算压强场的结果
Figure 6. Results of the calculated pressure field for the model under training conditions and time extrapolation
图 7 新初始场下，模型计算压强场的结果
Figure 7. Results of the calculated pressure field for the model with the new initial field
图 9 计算采用的网格的整体视图与局部视图
Figure 9. Overall and local views of the grid used for the calculation
图 10 以PDE损失训练PIGNN在不同时刻的计算结果
Figure 10. Computational results of training PIGNN with PDE losses at different moments
图 11 圆柱尾迹区添加数据损失的计算结果
Figure 11. Calculated results after adding data losses to the cylindrical wake area
图 12 在少量的几个时刻添加数据损失的计算结果
Figure 12. Calculated results with data losses added at a small number of moments
-
[1] RAISSI M, PERDIKARIS P, KARNIADAKIS G E. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[J]. Journal of Computational Physics, 2019, 378 : 686-707. doi: 10.1016/j.jcp.2018.10.045 [2] KARNIADAKIS G E, KEVREKIDIS I G, LU L, et al. Physics-informed machine learning[J]. Nature Reviews Physics, 2021, 3 : 422-440. doi: 10.1038/s42254-021-00314-5 [3] 王江, 陈文. 基于组合神经网络的时间分数阶扩散方程计算方法[J]. 应用数学和力学, 2019, 40 (7): 741-750. doi: 10.21656/1000-0887.390288 WANG Jiang, CHEN Wen. A combined artificial neural network method for solving time fractional diffusion equations[J]. Applied Mathematics and Mechanics, 2019, 40 (7): 741-750. (in Chinese) doi: 10.21656/1000-0887.390288 [4] 林云云, 郑素佩, 封建湖, 等. 间断问题扩散正则化的PINN反问题求解算法[J]. 应用数学和力学, 2023, 44 (1): 112-122. doi: 10.21656/1000-0887.430010LIN Yunyun, ZHENG Supei, FENG Jianhu, et al. Diffusive regularization inverse PINN solutions to discontinuous problems[J]. Applied Mathematics and Mechanics, 2023, 44 (1): 112-122. (in Chinese) doi: 10.21656/1000-0887.430010 [5] DAFERMOS C M. Hyperbolic Conservation Laws in Continuum Physics[M]. Berlin: Springer, 2016. [6] MAO Z, JAGTAP A D, KARNIADAKIS G E. Physics-informed neural networks for high-speed flows[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 360 : 112789. doi: 10.1016/j.cma.2019.112789 [7] LEVEQUE R J. Finite-Volume Methods for Hyperbolic Problems[M]. Cambridge: Cambridge University Press, 2002. [8] GODLEWSKI E, RAVIART P A. Numerical Approximation of Hyperbolic Systems of Conservation Laws[M]. New York: Springer, 1996. [9] COCKBURN B, KARNIADAKIS G E, SHU C W. Discontinuous Galerkin Methods[M]. Berlin: Springer, 2000. [10] MAGIERA J, RAY D, HESTHAVEN J S, et al. Constraint-aware neural networks for Riemann problems[J]. Journal of Computational Physics, 2020, 409 : 109345. doi: 10.1016/j.jcp.2020.109345 [11] SCHWANDER L, RAY D, HESTHAVEN J S. Controlling oscillations in spectral methods by local artificial viscosity governed by neural networks[J]. Journal of Computational Physics, 2021, 431 : 110144. doi: 10.1016/j.jcp.2021.110144 [12] BEZGIN D A, SCHMIDT S J, ADAMS N A. A data-driven physics-informed finite-volume scheme for nonclassical undercompressive shocks[J]. Journal of Computational Physics, 2021, 437 : 110324. doi: 10.1016/j.jcp.2021.110324 [13] BEZGIN D A, SCHMIDT S J, ADAMS N A. WENO3-NN: a maximum-order three-point data-driven weighted essentially non-oscillatory scheme[J]. Journal of Computational Physics, 2022, 452 : 110920. doi: 10.1016/j.jcp.2021.110920 [14] LIU L, LIU S, XIE H, et al. Discontinuity computing using physics-informed neural networks[J]. Journal of Scientific Computing, 2024, 98 : 22. doi: 10.1007/s10915-023-02412-1 [15] CAO W, SONG J, ZHANG W. A solver for subsonic flow around airfoils based on physics-informed neural networks and mesh transformation[J]. Physics of Fluids, 2024, 36 (2): 027134. doi: 10.1063/5.0188665 [16] HUANG H, LIU Y, YANG V. Neural networks with inputs based on domain of dependence and A converging sequence for solving conservation laws, part Ⅰ: 1D Riemann problems[J/OL]. 2021[2024-07-08]. https://arxiv.org/abs/2109.09316 .[17] ROE P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43 (2): 357-372. doi: 10.1016/0021-9991(81)90128-5 [18] GODUNOV S K. A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations[J]. Mathematics of the USSR-Sbornik, 1959, 47 : 271-306. [19] VAN LEER B. Towards the Ultimate Conservative Difference Scheme I. The Quest of Monotonicity[M]. Berlin: Heidelberg, 1973. [20] VAN LEER B. Towards the ultimate conservative difference scheme[J]. Journal of Computational Physics, 1997, 135 (2): 229-248. doi: 10.1006/jcph.1997.5704 [21] HARTEN A, ENGQUIST B, OSHER S, et al. Uniformly high order accurate essentially non-oscillatory schemes, Ⅲ[J]. Journal of Computational Physics, 1997, 131 (1): 3-47. doi: 10.1006/jcph.1996.5632 [22] LIU X D, OSHER S, CHAN T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115 (1): 200-212. doi: 10.1006/jcph.1994.1187 [23] LIU Y, SHU C W, TADMOR E, et al. Central discontinuous Galerkin methods on overlapping cells with a nonoscillatory hierarchical reconstruction[J]. SIAM Journal on Numerical Analysis, 2007, 45 (6): 2442-2467. doi: 10.1137/060666974 [24] LIU Y, SHU C W, TADMOR E, et al. Non-oscillatory hierarchical reconstruction for central and finite volume schemes[J]. Communications in Computational Physics, 2007, 2 (5): 933-963. [25] XU Z, LIU Y, SHU C W. Hierarchical reconstruction for discontinuous Galerkin methods on unstructured grids with a WENO-type linear reconstruction and partial neighboring cells[J]. Journal of Computational Physics, 2009, 228 (6): 2194-2212. doi: 10.1016/j.jcp.2008.11.025 [26] YANG V. Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems[J]. Proceedings of the Combustion Institute, 2000, 28 (1): 925-942. doi: 10.1016/S0082-0784(00)80299-4 [27] WANG X, YANG V. Supercritical mixing and combustion of liquid-oxygen/kerosene bi-swirl injectors[J]. Journal of Propulsion and Power, 2016, 33 (2): 316-322. [28] UNNIKRISHNAN U, HUO H, WANG X, et al. Subgrid scale modeling considerations for large eddy simulation of supercritical turbulent mixing and combustion[J]. Physics of Fluids, 2021, 33 (7): 075112. doi: 10.1063/5.0055751 [29] TEYSSIER R, COMMERÇON B. Numerical methods for simulating star formation[J]. Frontiers in Astronomy and Space Sciences, 2019, 6 : 51. doi: 10.3389/fspas.2019.00051 [30] BAR-SINAI Y, HOYER S, HICKEY J, et al. Learning data-driven discretizations for partial differential equations[J]. Proceedings of the National Academy of Sciences of the United States of America, 2019, 116 (31): 15344-15349. [31] 高普阳, 赵子桐, 杨扬. 基于卷积神经网络模型数值求解双曲型偏微分方程的研究[J]. 应用数学和力学, 2021, 42 (9): 932-947. doi: 10.21656/1000-0887.420050GAO Puyang, ZHAO Zitong, YANG Yang. Study on numerical solutions to hyperbolic partial differential equations based on the convolutional neural network model[J]. Applied Mathematics and Mechanics, 2021, 42 (9): 932-947. (in Chinese) doi: 10.21656/1000-0887.420050 [32] JIANG L, WANG L, CHU X, et al. PhyGNNet: solving spatiotemporal PDEs with physics-informed graph neural network[C]//Proceedings of the 2023 2 nd Asia Conference on Algorithms, Computing and Machine Learning. Shanghai: ACM, 2023: 143-147. [33] ZHANG H, JIANG L, CHU X, et al. Combining physics-informed graph neural network and finite difference for solving forward and inverse spatiotemporal PDEs[J/OL]. 2024[2024-07-08]. https://arxiv.org/abs/2405.20000 .[34] GILMER J, SCHOENHOLZ S S, RILEY P F, et al. Neural message passing for quantum chemistry[J/OL]. 2017[2024-07-08]. https://arxiv.org/abs/1704.01212v2 .[35] SANCHEZ-GONZALEZ A, GODWIN J, PFAFF T, et al. Learning to simulate complex physics with graph networks[C]//Proceedings of the 37 th International Conference on Machine Learning. 2020: 8459-468. -
下载:
下载:
计量
- 文章访问数: 168
- HTML全文浏览量: 99
- PDF下载量: 34
- 被引次数: 0
渝公网安备50010802005915号