📝 机器学习技术增强的CFD求解器  Enhancing CFD solver with Machine Learning techniques

摘要

本研究解决了在流体流动模拟中，由于捕捉时间尺度和长度尺度的需求，导致的计算网格的挑战。我们的方法专注于压力求解器，因为这是计算流体动力学 (CFD) 求解器中资源密集的组件。我们通过将机器学习 (ML) 代理模型与不可压缩流体流求解器集成来实现这一目标。我们创建了两种增强型 CFD 求解器的变体，它们能够在非稳态流体模拟中减少 CFD 压力求解器所需的迭代次数。因此，模拟产生了可比较的阻力系数和斯特劳哈尔数，伴随着八倍的执行时间减少。这些性能提升归因于每次时间迭代的计算工作量减少以及使用基于 ML 的代理模型时在模拟动力学行为上的早期强制力。本研究引入了一种通过集成代理模型来提高 CFD 模拟压力求解器计算效率的方法。我们提出了一个混合 CFD 求解器，即物理信息的求解器增强了数据驱动的代理模型。

1. 引言

人工智能 (AI) 应用正扩展到多个领域，科学计算也不例外。AI 驱动模型在科学计算中的应用范围可以大致分为三类： (i) 数据驱动模型，这些模型无需显式了解底层物理定律，并且往往需要大量数据集来训练；(ii) 物理信息模型，旨在遵守控制方程和边界条件；(iii) 混合方法，试图结合两者的优势。

AI 可以作为一种多功能工具，能够创建解决复杂科学问题的模型，例如增强 CFD 求解器。它还可以充当代理模型，替代或增强物理求解器，以加速计算工作流程。甚至有人假设，最近的发展可能导致用机器学习 (ML) 基于的解决方案部分替代 Navier-Stokes (N-S) 方程的 CFD 计算。

机器学习 (ML)，包括深度学习 (DL) 的子集，已经在科学计算中发挥了重要作用。这一贡献扩展到多个领域，如计算力学、生物力学、物理学、热力学及 CFD。在 CFD 应用中，一些示例包括：流体-粒子相互作用建模、改善传统流体模拟、海岸和海洋工程以及预测湍流量。

这些 ML 应用在特别是 CFD 中表现出显著的潜力，能够显著减少计算时间和成本，使其特别适用于设计目的。

流体流动代表了一个具有重大意义的复杂现象，涉及众多工程应用。在模拟这些物理过程时，通常是将数值模拟与实验研究相结合起到至关重要的作用。N-S 方程为描述流体流动提供了全面的框架，但对于大多数场景来说，达到时间和空间尺度的求解是不现实的。对于不可压缩流体流动，关键步骤是生成与质量守恒相一致的压力场，而这一现象由泊松方程支配。

泊松方程适用于包括电磁学、天文学和传热在内的多个领域的物理系统。在不可压缩流动中，压力泊松方程解决了大量的计算成本，因此有可能通过使用 ML 代理模型解决这一组件来增强不可压缩 CFD 求解器。

近期的研究已成功应用 ML 技术，包括卷积神经网络 (CNN) 和图神经网络 (GNN) 来解决泊松方程。在 CFD 领域的例子包括 Özbay 等人的工作，他们为二维笛卡尔网格计算压力场，Illarramendi 等人提出了一种方法来加速经典的多网格方法。尽管如此，仍然存在未解决的问题，例如，如何使这些代理模型推广到适应不同几何形状。我们之前引入了一种算法，用于创建用于解决泊松方程的 ML 代理模型。我们通过消除传统 CFD 方法中的神经网络与 CFD 求解器的耦合，能够在广泛的流域上应用训练的代理模型。

在本研究中，我们集成了其中一种代理模型进入 CFD 求解器，展示了其在非稳态流体模拟中的实际应用。为了实现这一目标，我们选择了 OpenFOAM 工具箱，并修改了 pisoFoam 不可压缩流体求解器以整合代理模型。我们系统地评估了该方法的性能，通过围绕二维圆柱体进行的初步研究来对其进行了现代化处理。

2. 方法 | 数学模型

2.1 流体求解器

2.1.1 控制方程

计算流体动力学 (CFD) 是通过基于数值方法的解算来数学建模流体流动的过程。这是一个广为人知的框架，允许实践者和研究人员分析复杂的物理系统。由于稀缺的解析解，CFD 在工程和研究领域中成为了关键工具。根据 Schlichting 和 Gersten 的工作，牛顿流体在不可压缩流动中的动量守恒方程可表达为：

其中， 是流体速度矢量，每个分量描述沿每个空间坐标的速度， 对应于运动压力， 表示偏应力张量， 表示作用在流体上的所有体力的净加速度。偏应力张量  与应变率张量  的关系为：

其中， 是单位张量， 定义为 。不可压缩流动的特点是流动过程中密度变化可忽略不计，可以假定适用于广泛的应用场景。这不仅限于液体，而是适用于特定流动条件下的任何流体，突显了其广泛适用性。考虑到不可压缩性，速度场变为无旋，质量守恒方程给出为：

除了公式 (1) 和 (3)，还必须指定初始和空间边界条件，才能完全求解系统。

2.1.2 离散化

公式 (1) 和 (3) 的解析解稀缺，因此数值方法通常是必不可少的。与将流体建模为连续介质不同，这些方法通过空间和时间值的网格离散化流场。这一过程涉及将几何形状分解为离散元素，并随后对相关的微分方程进行离散化。在 OpenFOAM 中，这是通过有限体积法 (FVM) 实现的，它利用守恒形式的控制方程，并通过应用高斯-奥斯特罗格拉德斯基定理在空间体积上积分这些方程。

这种积分确保了每个控制体积内物理量的守恒。关于在 OpenFOAM 框架下使用 FVM 的更多详细信息，读者可以参考 Versteeg 和 Malalasekera 的研究。FVM 可以应用于任意几何形状，使用结构化和非结构化网格。离散化方程中仍然存在的梯度随后由有限差分法 (FDM) 近似，例如应力张量项。

离散化后，方程被重新安排为一组线性方程，并且可以通过数值求解器求得解。当这些技术在求解简单的对流-扩散输运方程时表现良好时，处理 N-S 方程中的复杂数学性质时会遇到挑战。这些性质包括对流项的非线性、质量与动量方程之间的耦合以及在能量分布和可压缩性可以忽略不计的应用中缺乏压力演化的显式方程。这些挑战促使了许多算法的发展，能够处理这些性质。

2.1.3 压力-动量耦合

由 Issa 引入的压力-隐式分裂操作 (PISO) 方法是解决压力-动量耦合问题的广泛应用方法。PISO 通过确保速度场符合 N-S 方程和连续性方程来解决压力和速度之间的耦合。它采用了一种与逐次求解所有耦合方程不同的方法，而是将操作符分解为隐式动量预测器和多个显式校正步骤。该方法通常用于处理时间相关流体流动方程的隐式离散化。实验证明，仅需少量校正步骤即可实现较高的精度。

一旦动量方程被离散化后，可以在动量预测步骤后得到解，通常称为速度预测器，如图 1 所示。然而，需要注意的是，此解受连续性方程约束，因为 N-S 方程的速度场必须同时满足连续性方程。离散动量方程表示为：

其中  是由 FVM-N-S 方程离散化产生的系数矩阵， 是表示体加速度的矢量。根据 PISO 方法，在每次新迭代的开始，公式 (4) 可用于根据上一次迭代的压力场  解算速度场 。因此，预测的速度场  在此阶段不保证质量守恒。定义矩阵 ，其中包含公式 (4) 中  的对角成分，公式 (4) 可重新排列为：

其中  已被定义。通过结合公式 (4) 和 (5)，并通过乘以 ，对两个项取散度，并考虑到连续性方程，最终可以得到压力的方程：

通过求解此方程，得到的压力场可用于修正速度场（确保其满足连续性方程）。此修正速度场如图 1 中所示，记作 ，其中：

然而，速度场修正后，压力方程不再满足，因为  取决于  的更新。为了满足支配方程，要求在图 1 所示的流程中执行循环。

如果时间步长足够小，通常不需要执行多次压力-速度循环（或者在 OpenFOAM 术语中称为外循环）。实际上，只需要执行几次内部循环即可获得部分收敛。读者可参考 Issa 的工作了解更多详细内容。

2.2 代理模型

本研究依赖于 Sousa 等人开发的代理模型之一，详细描述了 ML 框架的工作原理。因此，鼓励读者参考该文献以深入了解模型及其组件。代理模型的架构如图 2 所示。代理模型的输入数据是映射到一组输入场的压力场 。这些输入包括非保守速度流场，由速度预测器计算得出（如  在图 1 中），以及表示流域几何形状的附加参数。

该代理模型由三部分组成：(i) 块提取；(ii) ML 模型，用于解决压力场；(iii) 装配算法。在构建基于 ML 的代理模型时，我们的目标是创建一种能够适应任何计算网格的方法。这是通过从代理模型中训练 NN 并通过 NN 预测固定形状块的压力场实现的，并实施了装配算法来确保在邻近块中保持流体域的计算一致性。

块提取将流体区域分解为输入到 NN 模型的小块。每个块代表更大的流体区域的一小部分，并可独立处理。每个块的边界条件基于几何场，几何场以基于签距函数 (SDF) 作为基础。SDF 用于将几何信息表示为每个空间坐标  的标量值。SDF 定义为：

其中， 表示流域存在的计算域， 表示包含障碍物的域。  表示 和  的边界的距离。

该装配算法系统地重新构建压力场，确保外边界条件的输出并最大限度地确保邻近块的一致性。装配算法的步骤如图 3 所示，通过考虑块的重叠区域，确保块之间的压力一致性。

NN 框架选择了一个由多层感知器 (MLP) 组成的代理模型，并具有作为网络编码器和解码器层的主成分分析 (PCA) 变换。此架构在图 2 中示意。

3. 求解器开发

为了实现一个可以利用开发的代理模型的 OpenFOAM CFD 求解器，需要在 Python 和 OpenFOAM (C++) 之间进行高效的通信，以尽量减少计算开销。主要目标是尽可能高效地在每个时间步将 C++ 数组输入到代理模型中并将输出返回到 C++。

编程细节可参考 Sousa 的源码，基础为 PyFoam 和 PythonFOAM。本文讨论的所有模拟均使用单个 CPU 内核执行，因此避免了并行化相关的复杂性。

3.1 求解器概述

减少传统 OpenFOAM 压力求解器的计算负担，并用开发的代理模型取而代之似乎是合理的。然而，最终得出的结论是，这种方法可能会丢弃常规数值求解器提供的可靠性，并可能严重损害结果的准确性。所遵循的策略包括使用 OpenFOAM 压力求解器来纠正代理模型的非物理行为，因为压力代理模型缺乏物理意识。显著减少压力修正迭代次数的可能性取决于代理模型的稳定性。

图4 展示了所遵循算法的概述。我们之前工作的提出解决方案(33)被确定为算法2。该方案涉及在动量预测步骤后集成替代模型，用作压力方程求解过程的增强，期望提供一个更好的初始猜测。另一种解决方案，称为算法1，在测试阶段表现出有希望的结果，并在求解压力方程的迭代次数受限时展示了较高的稳定性。在此方法中，替代模型部署在算法的不同阶段，用于计算压力场，，并应用于速度预测步骤。这涉及从上一个时间步中获取修正的速度场，用于替代模型的输入，与速度预测步骤中计算的非保守速度场形成对比，因为它是按原训练用途使用的。关于算法1，最佳训练场景起初似乎包括使用而非，因此需要训练一个新的替代模型。然而，这种调整将替代模型与特定的时间步绑定。因此，决定避免重新训练，而是使用Sousa等人在文献(33)中开发的替代模型。

根据算法1和算法(图4)，生成了两个不同的求解器，即DLpisoFoam-alg1和DLpisoFoam-alg2求解器。有关实现细节，读者可参考附录A。

3.2 验证案例

为验证求解器的有效性，定义了一个验证案例，涉及二维圆柱体流动的模拟。流域的可视化表示见图5，流域边界条件见表1。速度场根据描述在平行板间完全发展的层流的动量方程的解析解得出：

其中，为从中心线的距离，为板间距离的一半，为平均纵向速度。

为了有效捕捉靠近壁面的陡峭梯度，网格在那些区域进行了加密，如图6所示。

表征仿真案例的其余参数包括阻塞比率，其中为双维圆柱体的直径，为平行板之间的距离；以及雷诺数，定义为：

其中为运动粘度，为中心线处的速度大小。

验证案例使用10个网格加密级别执行，以评估计算规模对求解器性能的影响。这些加密级别基于网格单元数表征，见表2。每次仿真都使用来自Intel E5-2450 CPU的单核，性能为16.8 GFLOPs。

为了评估新开发求解器的性能，对其进行了系统的基准测试，比较对象为pisoFoam求解器。该比较涉及两个关键指标：计算效率，通过时钟时间测量来量化，以及准确性，通过计算平均阻力系数和Strouhal数来量化。阻力系数定义为：

其中，为障碍物所产生的阻力，为参考面积。对应于准稳态区域中平均，其中涡流脱落是周期性的，而定义为：

其中，为涡流脱落的频率。

4. 数值结果

4.1. 求解器配置

在介绍替代模型的论文中(33)，从误差场中可以很容易识别出边界伪影，代表了预测的压力场和真实值之间的偏差。这些伪影出现在预测的压力场中，当引入离散化的Navier-Stokes方程（方程(4)）时，会显著增加CFD求解器计算的残差，因为它们扰乱了局部平衡。为了解决这一问题，进行了研究，探索应用高斯滤波器平滑预测的压力场的可行性，如图7所示。该研究旨在表征平滑滤波器对归一化均方根误差（RMSE）和计算时间的影响。RMSE定义为：

其中是计算估计值的单元数，是参考值，是预测值。这些误差进一步归一化，按除以。此外，我们优化了组装算法中重叠区域的大小，以在计算效率和精度之间达到最佳平衡。图7展示了应用平滑滤波器对精度的影响较小。然而，有必要承认在OpenFOAM框架内，从每个控制体积中的局部平衡推导出的残差，平滑操作具有比全球RMSE所能表明的更显著的效果。如预期的那样，计算成本随着重叠比率和滤波器核大小的增加而上升。

在使用开发的求解器时，有几个关键变量需要调整，这些变量大致可以分为与替代模型相关的参数和与PISO压力修正相关的其他参数。第一组参数包括重叠比、Gaussian滤波器核大小和阈值（如算法4中定义），其设置如表3所述。另一方面，第二组参数，即PISO压力修正的配置，研究了每个加密级别下的影响。运行之间的关系和使用的特定参数如表4所示。

4.2. 求解器性能

两个求解器都进行了相同的测试，比较的性能指标包括执行时间，以及计算的流动相关变量和。这些变量用于与来自OpenFOAM工具箱的pisoFoam求解器进行比较。流动相关的变量源自和的时间演变，如图8和图9所示。的计算包括在波动变为周期性时的时间段内取平均值，而的计算基于这些波动的频率。定义仿真完成的标准包括验证了准稳态流动状态，即如果在拖曳系数中的时间波动仍然波动但围绕一个恒定的平均值，则认为仿真完成。为避免每次仿真都进行独立验证，进行了初步研究。该研究揭示了pisoFoam和DLpisoFoam求解器在达到准稳态状态所需时间的显著差异。基于初步研究，维数无关的终止时间为pisoFoam仿真设置为，而DLpisoFoam仿真设置为，其中，，为特征长度，为入口处的中心线速度。从图8和图9可以看出，所实现的求解器加速了准稳态状态的开始，特征为周期性涡脱落。

图10展示了所有运行的执行时间，涵盖了不同网格大小和运行配置，如表4中所指定的。然而，测量最终加密级别的确切执行时间具有挑战性，原因是在单个CPU核心上运行pisoFoam所需的显著计算时间。为解决这一挑战，我们估算了最终加密级别的pisoFoam执行时间，使用了指数拟合（对应于类型的拟合）和三阶多项式拟合。指数拟合可以涵盖整个结果范围，但倾向于提供过于乐观地估计了较大网格的情况。因此，导致对pisoFoam性能的估计倾向于低估DLpisoFoam求解器的相对性能。相反，多项式拟合可以更好地代表最大网格的执行时间，但提供了pisoFoam性能的更为保守的估计。这两种估计都在图10中用虚线表示。有关更全面的分析，请参阅附录B中的值。

使用图10中显示的数据，DLpisoFoam-alg1和DLpisoFoam-alg2相对于基准求解器（来自OpenFOAM的pisoFoam）的加速效果展示在图11中，涵盖了各种网格大小和运行配置，如表4所定义。与执行时间结果相似，这些加速结果也包括最大网格的乐观和悲观估计。在此背景下，悲观结果由部分透明线表示。在图11中，标出了一个水平阈值，即新开发的求解器变得具有优势的转折点，该阈值在约个网格单元的仿真中被超越。

值得注意的是，PISO循环中的压力修正迭代次数对DLpisoFoam求解器的模型结果产生了显著影响。增加PISO算法的迭代次数将结果更接近于传统的pisoFoam求解器，但以增加计算工作量为代价。关于所需的PISO压力修正次数，这两种不同的算法表现出不同的行为。使用DLpisoFoam-alg2进行的仿真在压力修正迭代次数不足（小于4次内部PISO修正循环）时表现出发散，其值在附录B的表B.1中显示。相比之下，DLpisoFoam-alg1在所有情况下表现出稳健的收敛性。鉴于这些DLpisoFoam-alg2仿真的显著比例发生发散，这些特定求解器配置已从后续结果中删除。然而，详细数据可在附录B中找到。

必须强调的是，通过集成滤波器并通过组装算法扩展交叉区，可以潜在地提高精度并减少对PISO压力求解器的过度依赖。然而，这种精度改进的代价是计算时间的显著增加。对于粗网格，计算工作量的增加可能会完全消除使用开发的机器学习增强型求解器的优势。然而，对于更精细的网格，计算效率的影响可能会变得可以忽略不计。尽管每次仿真都在不同的HPC计算节点上使用等效的硬件，性能的变化仍然变得显而易见，如图12(a)中所示。需要注意的是，这种变化可能源自节点间和节点内的来源。

节点间变化源于不同计算节点之间的性能变化，尽管使用了等效的硬件，但仍可能引入仿真结果的不确定性。这种变化可能受到节点配置、网络状况和系统负载等因素的影响。在我们的研究中，我们观察到结果中的节点间变化，这导致了性能测量中的不确定性。为减轻这种节点间的变化，我们在不同节点上执行了相同的一组仿真，并采取措施使执行时间在节点之间归一化。通过归一化，确保了跨节点的平均执行时间保持一致，从而减少了节点间变化的影响。结果性能如图12(b)所示。值得注意的是，尽管采取了这些措施来减少节点间变化，仍然存在一些节点内变化，并且无法消除。节点内变化涉及同一仿真在不同运行中的性能变化，可能由局部条件、系统负载波动和资源争用引起。

4.3. 求解器验证

为了验证DLpisoFoam求解器，使用这些求解器计算了平均阻力系数和Strouhal数，并研究了在二维圆柱体上的约束流动。这些结果的相对误差，定义为误差，用作从OpenFOAM工具箱收敛计算的结果的参考（参见附录B）。这些相对误差计算并在图13和图14中分别图示了DLpisoFoam-alg1和DLpisoFoam-alg2的情况。结果还与文献中的结果进行了额外验证(47)（参见附录B中的表B.2和表B.）。

两者和在大多数情况下的相对误差的绝对值均低于0.5％，验证了开发的求解器相对于成熟的CFD求解器的准确性。对压力修正迭代次数的更严格限制导致误差增加和运行时间缩短。研究发现，对PISO修正内循环（或OpenFOAM术语中的nCorr）的最大迭代次数施加更严格的限制比仅限制实际的内循环数量提供了更好的结果。比较“4 p corr - 2 iter/corr”和“2 p corr”（分别进行了4次和2次修正迭代）表明，前者产生了更好的结果（参见图13）且运行时间较短（参见图10），同时仍然使用相同的迭代次数来求解压力方程。此外，对每个内循环迭代次数的限制被证明是一种更有效的方法，因为它允许收敛性较差的求解器（DLpisoFoam-alg2）收敛，即使每个内修正循环仅进行一次迭代（参见“4 p corr - 1 iter/corr”结果）。

需要注意的是，随着网格大小的增加，计算变量的误差似乎也在增加。这是预期的，因为对迭代次数的固定限制对于较大的网格变得更加严格，较大网格通常需要更多的迭代。然而，“4 p corr”结果仍然在可接受的验证范围内，同时提供了显著的加速效果。

5. 分析

在性能方面，当处理大网格时，两种求解器变体都优于传统数值模型。转折点出现在大约个单元处，并随算法使用的PISO压力修正迭代次数变化。值得注意的是，这些求解器的优势随着网格中单元数量的增加而增加，在研究范围内速度提高了近一个数量级。

速度改进来自于更快的时间迭代和更快地达到仿真的准稳态状态（由周期性涡脱落表征）。在这一点上，DLpisoFoam-alg1在这方面优于DLpisoFoam-alg2，正如图8中所示。仔细检查图8表明，理想情况下，DLpisoFoam-alg2仿真应允许进行更多的时间迭代以进行更公平的比较。

来自替代模型的压力估计的引入对数值方法的稳定性产生了积极影响。这一改进导致达到收敛所需的压力修正迭代次数减少，从而降低了总体仿真成本。显著的是，DLpisoFoam-alg1求解器通过始终以最少的压力修正迭代收敛，展示了其在PISO算法中的稳健性，即使依赖于较少的压力修正迭代次数也能实现收敛。相比之下，DLpisoFoam-alg2求解器在选择适当的压力修正迭代次数以实现收敛时需要更多的考虑。当使用较少的PISO内循环时，它表现出发散现象（如附录B的表B.1所示），尤其是在大网格的情况下。值得注意的是，使用pisoFoam进行仿真时，被推到了极限以允许公平比较。确保进一步减少迭代次数或压力修正会导致这些仿真中的发散。

选定的感兴趣指标和验证了新求解器。这些计算的误差通常低于0.5％。这种高水平的精度是预期的，因为对于每次时间迭代，PISO算法使用选定数量的压力修正迭代，从而减少了替代模型的预测，并确保了每个控制体积内离散化控制方程的最小残差。正如预期的那样，随着对PISO修正循环最大迭代次数的更严格限制，误差随着迭代次数的增加而加剧。

6. 结论

基于Sousa等人（33）开发的机器学习（ML）替代模型的不同集成结果，生成了两种不同的求解器：DLpisoFoam-alg1和DLpisoFoam-alg2。两者都能够通过使用机器学习替代模型辅助PISO算法中的压力修正步骤，大大加速不可压缩等温二维流体的仿真。这两个求解器能够生成绝对误差通常低于0.5％的精确解，同时计算时间减少了8倍。与DLpisoFoam-alg2相比，DLpisoFoam-alg1表现出更好的性能和稳定性。它表现出更高的稳定性，更少依赖PISO压力修正，并在更早阶段引入流体动力行为。

研究结果表明，随着网格大小的增加，加速因子的增长趋势保持一致，这主要是因为替代模型的计算成本在很大程度上不受网格单元数的影响，而传统求解器则明显受到影响。因此，预计替代模型的影响在涉及更大网格，尤其是在三维流动的场景中将更加显著。

本工作提供了一个可应用于OpenFOAM中任何不可压缩求解器的框架。开发的求解器是公开可用的，基于Sousa （44）的研究，已经能够在某些特定条件下（例如大网格和接近替代模型的训练条件下的仿真）从性能改进中获益。实现了一个经过良好优化的集成，保留了替代模型的准确性和求解器的改进，使得这些CFD求解器具有实际意义。因此，未来的工作将集中于改进替代模型，使这些ML增强的CFD求解器在比参考CFD求解器更广泛的条件下更具优势。最终目标是集成能够准确覆盖特定类型仿真的替代模型，使得CFD求解器在设计仿真中极大地减少计算成本。