求解二维不可压Navier-Stokes方程的PINN算法

PINN Algorithm for Solving Two-Dimensional Incompressible Navier-Stokes Equations

摘要: 本文采用物理信息神经网络(PINN)来求解不可压缩湍流Navier-Stokes方程。本研究引入了动态权重调整策略，使得各项误差在训练过程中得到适当的平衡，从而避免了某些误差项主导整个训练过程的问题。此外，为了加速训练收敛并提高精度，本研究还对网络结构进行了优化，结合物理约束优化过程，改变了优化方法，提高了模型的训练效率。

Abstract: In this paper, physical information neural networks (PINN) are used to solve the Navier-Stokes equations of incompressible turbulence. In this study, the dynamic weighting adjustment strategy is presented to make the errors properly balanced in the training process, so as to avoid the problem that some error terms dominate the whole training process. In addition, in order to accelerate the training convergence and improve the accuracy, this study also optimized the network structure, combining with the physical constraint optimization process and changing the optimization method to improve the training efficiency of the model.

1. 引言

物理信息神经网络(PINN)是近年来由Raissi等人[1]提出的一种结合物理约束的深度学习方法。该方法的核心思想是将物理方程作为损失函数的一项加入到神经网络中，通过神经网络的训练过程求解偏微分方程(PDE)。与传统的基于数据的深度学习方法不同，PINN不依赖大量的实验数据，而是直接利用物理定律(如Navier-Stokes方程)作为先验知识，在训练过程中结合物理信息。PINN方法不仅可以有效地求解流体力学中的基础方程，还能够处理复杂的边界条件、非线性问题和多尺度现象，因此在流体力学中有广阔的应用前景。

PINNs已经历了诸多扩展和改进，以增强其适用性并克服原始框架的局限性。包括融合其他网络架构如卷积神经网络(CNNs) [2]和循环神经网络(RNNs) [3]，进一步拓宽了PINNs的应用范围。神经网络具有强大的非线性映射能力和自适应特性，能够有效捕捉复杂系统中的潜在规律。在求解复杂流体问题，处理复杂的边界条件、噪声和不适定问题时，表现出较好的鲁棒性[4]。

因此，结合物理信息神经网络(PINN)的方法，为不可压缩流体的Navier-Stokes方程提供了一种新的解法。通过结合深度学习的强大推理能力和物理方程的约束，PINN有望克服现有数值方法的瓶颈，提供更为高效、精确的解决方案。

当前，PINNs中标准的优化算法是Adam和拟牛顿方法(如BFGS)。BFGS通过引入可训练变量空间的二阶信息进行预处理，显著提升了PINNs的性能。尽管这两种优化器因其理论和数值特性不仅在PINNs中，而且在其他优化问题中也广为应用，但近期在PINNs中的实践表明，对于许多问题，其他算法可能更具优势。

2. 物理信息神经网络模型

2.1. 模型框架

简单的PINN网络主要有三个部分，输入层为时空坐标$\left(x,t\right)$ ，中间层为全连接层，输出层为速度和压力$\left(u,p\right)$ 。模型的损失函数的一般形式为：

$L_{PINN}={\lambda }_1{L}_{IC}+{\lambda }_2{L}_{BC}+{\lambda }_3{L}_f+{\lambda }_4{L}_{data}$ (1)

其中，$L_{data}$ 是数据误差项，$L_{IC}$ 和$L_{BC}$ 分别是初始条件和边界条件的误差项，$L_f$ 是物理方程残差项，其中包含问题所满足的物理条件。每个残差项的系数(如${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4$ )用于平衡损失函数中各项的影响，确保训练过程中的各项条件得到适当的满足。图1展示了一个二维物理信息神经网络的结构示意图。

Figure 1. Structure of physics-information neural network

图1. 物理信息神经网络结构

2.2. 初始条件与边界条件残差的处理

在PINN中，初始条件的残差处理体现在损失函数的设计中。设定初始条件时，网络需要确保在时刻$t_0$ 时，速度场和压力场满足初始条件。为了实现这一目标，在损失函数中加入初始条件的残差项，具体形式为：

$L_{IC}=\frac{1}{N_i}{\displaystyle \sum _{n=1}^{N_i}{\Vert u^n-u_i^n\Vert }^2,}$

其中，$N_i$ 表示初始条件的训练点个数，$u_i^n$ 是初始时刻$t_0$ 第n个数据点的给定速度，$u^n$ 是PINN模型第n个数据点在初始时刻预测的速度，这里$\Vert ·\Vert$ 表示$L_2$ 模。

边界条件的残差处理方式与初始条件类似，同样通过损失函数中的残差项来确保网络输出符合给定的边界条件。损失函数中与边界条件相关的残差项通常表示为：

$L_{BC}=\frac{1}{N_b}{\displaystyle \sum _{n=1}^{N_b}{\Vert u^n-u_b^n\Vert }^2,}$

其中，$N_b$ 表示边界条件的训练点个数，$u_b^n$ 是边界$\partial \Omega$ 上第n个数据点的给定速度，$u^n$ 是PINN模型第n个数据点在该边界点预测的速度。

2.3. 方程残差与物理约束结合

首先计算区域内部点的数据误差$L_{data}$ ，计算过程如下：

$L_{data}=\frac{1}{N_d}{\displaystyle \sum _{n=1}^{N_i}{\Vert u^n-u_d^n\Vert }^2,}$

其中，$N_d$ 表示数据点训练点个数，$u_d^n$ 是第n个数据点的给定速度，$u^n$ 是PINN模型在第n个数据点预测的速度。

在PINN中，为了确保网络的输出能够满足这些物理方程，将方程的残差引入到损失函数中。具体来说，方程残差的计算过程如下：

若问题满足N个物理方程$f_n\left(x,t\right)\left(n=1,\cdots ,N\right),x\in \Omega ,t\in \left[0,T\right]$ 。可以根据神经网络的输出计算每个方程项的残差，计算过程如下：

$L_f=\frac{1}{N_f}{\displaystyle \sum _{i=1}^{N_f}{\displaystyle \sum _{n=1}^N{\left|f_n\left(x^i,t^i\right)\right|}^2,}}$

其中，$\left(x^i,t^i\right)$ 表示第i个内部点的时空坐标。

通过计算这些残差，能够确定网络的输出与物理方程之间的偏差。

3. 改进后的物理信息神经网络模型

3.1. 动态权重调整策略

在PINN的训练过程中，由于初始条件、边界条件和方程残差项的数值范围和变化速率差异较大，若固定权重，某一项可能会过度主导训练，导致网络忽视其他约束条件，最终影响求解精度。因此，采用动态权重调整策略，实时调整损失函数中各项的权重，从而使得每个约束条件在网络训练的各个阶段发挥合适的作用。

残差项的大小反映了网络在某个方面的偏差程度，较大的残差表明该项在当前训练阶段对解的影响较大，需要更多的关注，较小的残差则表明该项对解的贡献较小，可以适当减小权重。一种常见的动态权重调整策略[5]是基于残差归一化的方式来更新每个项的权重。公式如下：

${\lambda }_i=\frac{L_i}{{\displaystyle \sum _{i=1}^K{L}_i+\epsilon }}$ (2)

其中，${\lambda }_i$ 是第i个损失项(如初始条件、边界条件或方程残差)的动态权重，$L_i$ 是对应项的残差，$\epsilon$ 是一个小常数，避免分母为零，在本文中$\epsilon$ 取10⁻⁸，避免对计算结果产生显著影响。本文中动态权重的调整频率为每训练1000轮次调整一次权重，既能避免每一轮调整权重给网络带来计算量和时间的增加，又可以确保了在训练过程中，较大的残差对应较大的权重，较小的残差对应较小的权重，从而实现平衡各个约束条件的作用。动态权重调整策略的算法如算法1所示。

算法1 动态权重调整策略

输入： 每个数据点对应的空间坐标$x$ 和时间$t$ 。 输出： 每个时空点对应的速度$u$ 和压力$p$ 。 1: 初始化:${\lambda }_1={\lambda }_2={\lambda }_3={\lambda }_4=1$

2: for $epochs=1,\cdots ,10000$ do 3: $\left(u,p\right)=PINN\left(x,t\right),L_{IC}=\frac{1}{N_i}{\displaystyle \sum _{n=1}^{N_i}{\Vert u^n-u_i^n\Vert }^2},L_{BC}=\frac{1}{N_b}{\displaystyle \sum _{n=1}^{N_b}{\Vert u^n-u_b^n\Vert }^2,}$ $L_{data}=\frac{1}{N_d}{\displaystyle \sum _{n=1}^{N_i}{\Vert u^n-u_d^n\Vert }^2},L_f=\frac{1}{N_f}{\displaystyle \sum _{i=1}^{N_f}{\displaystyle \sum _{n=1}^N{\left|f_n\left(x^i,t^i\right)\right|}^2,}}$ $L_{PINN}={\lambda }_1{L}_{IC}+{\lambda }_2{L}_{BC}+{\lambda }_3{L}_f+{\lambda }_4{L}_{data}$

4: if $L_{PINN}>{10}^{-5}$ or $epochs\mathrm{mod}100=0$ then

5: ${\lambda }_i=\frac{L_i}{{\displaystyle \sum _{i=1}^4{L}_i+\epsilon }}$

6: end if

7: end for

8: return $\left(u,p\right)$

3.2. Adam优化方法

神经网络中的优化问题主要寻找一组参数(权重和偏置) ${\theta }^{*}$ ，使得损失函数$L\left(\theta \right)$ 最小化：

${\theta }^{*}=\arg \underset{\theta }{\min }L\left(\theta \right)$

Adam [6] (Adaptive Moment Estimation)优化算法是一种自适应的梯度优化方法。Adam通过计算每个参数的梯度的一阶矩(均值)和二阶矩(方差)来动态调整学习率，从而在训练过程中实现快速收敛和较小的波动。Adam算法在处理具有稀疏梯度的问题时，具有非常高的效率和稳定性。

Adam算法的更新规则如下：

${\widehat{m}}_t=\frac{m_t}{1-{\beta }_1^t}$

${\widehat{v}}_t=\frac{v_t}{1-{\beta }_2^t}$

${\theta }_t={\theta }_{t-1}-\frac{\eta {\widehat{m}}_t}{\sqrt{{\widehat{v}}_t}+\epsilon }$

其中，${\widehat{m}}_t$ 和${\widehat{v}}_t$ 分别表示梯度的偏差校正一阶矩和二阶矩，${\beta }_1$ 和${\beta }_2$ 是衰减率，$\eta$ 是学习率，$\epsilon$ 是一个小常数，用来避免分母为零。在PINN的训练过程中，Adam算法通过自适应调整每个参数的学习率，使得网络在训练过程中既能保持较快的收敛速度，又能避免梯度消失或爆炸问题。

3.3. SSBroyden优化方法

在优化问题中，需要求解的基本问题是$\underset{x\in R^n}{\min }f\left(x\right)$ 。对于这一问题，一个基本的挑战是确定从当前点

过渡到改进点的最佳方向和步长。牛顿法是机器学习中用得比较多的一种优化算法。牛顿法的基本思想是利用迭代点$x_k$ 处的一阶导数(梯度)和二阶导数(Hessian矩阵)对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小值。

拟牛顿法也是一类用于优化问题的迭代算法，旨在解决牛顿法计算复杂度高的问题，同时保留其快速收敛的特性。拟牛顿法通过构造$x_k$ 处$f$ 的Hessian矩阵的近似矩阵$B_k$ ，满足以下条件：

$B_{k+1}{s}_k=y_k,$

其中，$s_k=x_{k+1}-x_k,y_k=\nabla f\left(x_{k+1}\right)-\nabla f\left(x_k\right)$ 。通过迭代公式更新$B_k$ ，避免直接计算Hessian矩阵，特别适用于高维优化问题。

Wolkowicz等人[7] [8]提出的BFGS和其他Broyden家族更新的缩放方法家族进一步支持了缩放技术在某些优化问题中的有效性。直接和逆Hessian估计的自动缩放更新可以轻松扩展到一般Broyden家族，以下称为SSBroyden更新。这些更新由以下公式给出：

$B_{k+1}={\tau }_k\left[B_k-\frac{B_k{s}_k{s}_k^{\text{T}}{B}_k}{s_k^{\text{T}}{B}_k{s}_k}+{\theta }_k\left(s_k^{\text{T}}{B}_k{s}_k\right)w_k{w}_k^{\text{T}}\right]+\frac{y_k{y}_k^{\text{T}}}{y_k{s}_k^{\text{T}}},$

$H_{k+1}=\frac{1}{{\tau }_k}\left[H_k-\frac{H_k{y}_k{y}_k^{\text{T}}{H}_k}{y_k^{\text{T}}{H}_k{y}_k}+{\Phi }_k\left(y_k^{\text{T}}{H}_k{y}_k\right)w_k{w}_k^{\text{T}}\right]+\frac{s_k{s}_k^{\text{T}}}{y_k^{\text{T}}{s}_k},$

其中，${\theta }_k$ 是一个标量参数，并且

${\tau }_k=\min \left\{1,\frac{1}{b_k}\right\},$

$w_k=\frac{y_k}{y_k^{\text{T}}{s}_k}-\frac{B_k{s}_k}{s_k^{\text{T}}{B}_k{s}_k},$

$v_k=\frac{s_k}{y_k^{\text{T}}{s}_k}-\frac{H_k{y}_k}{y_k^{\text{T}}{H}_k{y}_k},$

${\Phi }_k=\frac{1-{\theta }_k}{1+\left(h_k{b}_k-1\right){\theta }_k},$

$b_k=\frac{s_k^{\text{T}}{B}_k{s}_k}{y_k^{\text{T}}{s}_k},$

$h_k=\frac{y_k^{\text{T}}{H}_k{y}_k}{y_k^{\text{T}}{s}_k}.$

在本文的数值实验中，我们将SSBroyden优化与自适应动态权重调整策略进行结合，观察权重策略和优化器同时改变之后损失函数的变化情况。

4. 数值实验

4.1. 数据集

本文研究的问题是二维圆柱绕流问题：该问题采用二维不可压缩的Navier-Stokes方程进行求解，方程形式为：

$\nabla u=0,$ (3)

$\frac{\partial u}{\partial t}+u\nabla u=-\frac{1}{\rho }\nabla p+\mu {\nabla }^{2}u,$ (4)

$u\left(x,y,0\right)=u_0,v\left(x,y,0\right)=0,\left(x,y\right)\in \Omega ,$

$\frac{\partial u}{\partial n}=0,\frac{\partial v}{\partial n}=0,在{\Gamma }_1上,$

$u=0,v=0,在{\Gamma }_0上,$

其中，$u=\left[u,v\right]$ 是流速场，$p$ 是压力场，$\mu$ 是运动粘度，实验中取0.01。设${\Omega }_1$ 为方形计算区域，${\Omega }_2$ 为${\Omega }_1$ 中圆的内部，记$\Omega ={\Omega }_1\backslash {\Omega }_2$ ，圆的边界记为${\Gamma }_0$ ，如图2所示。

Figure 2. Computational area of two-dimensional cylinder wake

图2. 二维圆柱绕流问题的计算区域

时空区域：计算域为二维矩形区域${\Omega }_1=\left[-15,25\right]\times \left[-8,8\right]$ ，时间域为[0, 7]。边界条件为固定边界速度和无滑移条件。

由Navier-Stokes方程(3)、(4)可得，其质量守恒和动量守恒形式如下：

$f_1\left(x,y,t\right):=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,\left(x,y\right)\in \Omega$

$f_2\left(x,y,t\right):=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\frac{\partial p}{\partial x}-0.01\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=0,\left(x,y\right)\in \Omega$

$f_3\left(x,y,t\right):=\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\frac{\partial p}{\partial x}-0.01\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}=0,\left(x,y\right)\in \Omega$

则方程的残差$L_f$ 变为：

$L_f=\frac{1}{N_f}{\displaystyle \sum _{i=1}^{N_f}\left({\left|f_1\left(x^i,y^i,t^i\right)\right|}^2+{\left|f_2\left(x^i,y^i,t^i\right)\right|}^2+{\left|f_3\left(x^i,y^i,t^i\right)\right|}^2\right)},$

其中，$\left(x^i,y^i,t^i\right)$ 是计算区域内部第i个内部点的时空坐标。

4.2. 实验环境与参数

本文实验环境为Windows 10家庭版(64位)，处理器11th Gen Intel(R) Core(TM) i5-1135G7 (2.42 GHz)，内存16 GB，Python版本为3.10.8。数据集为Raissi等人[1]的二维层流圆柱绕流数值计算结果，来源于网站https://github.com/maziarraissi/PINNs。实验研究选择了全连接神经网络(FCNN)作为PINN的基础架构。网络由输入层、8个隐藏层(每层16个神经元)、输出层组成，激活函数使用tanh。输出层包括两个速度分量和一个压力值。PINN的损失函数使用公式(1)，其中包括数据误差、物理方程残差、初始条件和边界条件的误差。实验的数据点随机生成，其中内部点6000个，边界点500个，初始点400个，总采样点数量为6900个。为了确保每个部分对损失函数的贡献合理，本文使用公式(2)引入了动态权重调整策略，以平衡各部分误差在训练过程中的影响。

4.3. 实验评价指标

在进行PINN与传统数值方法(如有限元法和有限差分法)的比较时，本文选择了一个典型的二维流动问题作为实验对象，使用有限元法求解其数值解，并将结果与PINN求解结果进行比较。

4.4. 结果对比

本实验选择了一个二维圆柱绕流问题作为测试案例，计算域为一个矩形区域，流体的初始条件和边界条件已经定义。流动问题的求解依赖于PINN模型通过同时最小化数据误差、物理方程的残差和边界条件的残差来得到最优解。本文先根据Raissi等人[1]的文章对实验进行复现，原文采用固定权重 + Adam优化的方法，本文采用动态权重 + Adam优化的方法，对两者实验结果进行比较，观察动态权重调整策略对实验结果的影响，实验结果如表1所示。

从表1中可以看出，在训练10,000轮次之后，动态权重调整策略的总损失达到了0.00083，比静态权重策略的总损失0.0051精度提高了一个数量级，并且训练时间只增加了350秒。表明在相同优化器的情况下，动态权重调整策略相比于固定权重策略，能够更有效地降低总损失，更快地提高精度。

表1可以证明动态权重调整策略相对于固定权重策略的优越性，本文接下来研究优化器的改变对于PINN实验结果的影响，本文在都采用动态权重调整策略的情况下，分别使用Adam优化方法和SSBroyden优化方法，观察不同优化方法对实验结果的影响，实验结果如表2和图3所示。

Table 1. Change of loss function and training time under different weighting strategies

表1. 不同权重策略下的损失函数和训练时间的变化

训练轮次	固定权重策略 损失函数(L)	固定权重策略 训练时间(s)	动态权重策略 损失函数(L)	动态权重策略 训练时间(s)
1000	0.0753	153	0.0272	182
2000	0.0589	301	0.00835	375
3000	0.0431	448	0.00527	561
4000	0.0319	607	0.00352	750
5000	0.0223	755	0.00235	942
6000	0.0141	910	0.00153	1085
7000	0.0095	1061	0.00118	1276
8000	0.0072	1218	0.00101	1410
9000	0.0063	1364	0.00090	1633
10,000	0.0051	1522	0.00083	1872

Table 2. Comparison of losses of different optimizers under dynamic weighting strategy

表2. 动态权重策略下不同优化器的损失比较

训练轮次	Adam优化 损失函数(L)	Adam优化 训练时间(s)	SSBroyden优化 损失函数(L)	SSBroyden优化 训练时间(s)
1000	0.0272	182	0.0115	215
2000	0.00835	375	0.00328	440
3000	0.00527	561	0.00102	659
4000	0.00352	750	0.00043	863
5000	0.00235	942	0.00031	1097
6000	0.00153	1085	0.00024	1336
7000	0.00118	1276	0.00018	1561
8000	0.00101	1410	0.00015	1772
9000	0.00090	1633	0.00012	1995
10,000	0.00083	1872	0.00009	2246

从表2和图3中可以看出，都采用动态权重调整策略的情况下，优化器从Adam改成SSBroyden之后，在训练10,000轮次之后的总损失从0.00083降低到了0.00009，表明随着训练的进行，模型的误差逐渐减少。动态权重调整策略和SSBroyden优化方法相结合会在每个训练轮次中表现出更小的损失函数，表明该方法能够更有效地优化模型。本文结果表明，SSBroyden方法显著提升了PINNs的收敛速度和精度，特别是在处理复杂和刚性PDEs时，收敛速度更快，且在复杂优化场景中表现出更强的鲁棒性。这些发现凸显SSBroyden方法在加速PINNs训练和增强数值稳定性方面的有效性。本研究表明，SSBroyden方法是训练PINNs的高效且可靠的优化器。其处理复杂PDEs的能力和强大的收敛特性使其在基于深度学习的科学计算中具有广泛的应用前景。未来的研究可以专注于结合领域特定的适应性方法，以进一步提升PINNs在高维和多尺度问题中的准确性、效率和泛化能力。

Figure 3. Loss reduction graph of different optimizers under dynamic weighting strategy

图3. 动态权重策略下不同优化器的损失下降图

5. 总结与展望

本文在基于PINN的神经网络的结构下，使用动态权重调整策略，调整神经网络损失函数中的权重。同时，对神经网络当中的优化器进行改变，将Adam优化方法变成了SSBroyden优化方法。在优化方法改变之后，也成功地在训练轮次不变的前提条件下，将PINN的总损失降低了一个数量级，提高了模型训练效率。尽管PINN在流体力学问题中已经表现出较高的计算效率和精度，但在大规模复杂流体问题的求解中，仍然面临着计算量和精度的挑战。PINN依赖于神经网络来逼近流体力学方程的解，这需要大量的训练数据和长时间的训练过程。特别是在三维高雷诺数湍流问题或多物理场耦合问题中，计算资源消耗仍然是一个瓶颈。因此，如何进一步提升PINN的计算效率和精度，成为未来研究的一个关键方向。

NOTES

^*通讯作者。