0 引　言

流体模拟一直是研究的热点，广泛应用于工业过程建模、城市数字孪生系统以及军事仿真等领域。本文的流体模拟方法以拉格朗日粒子法^[1]为基础，粒子法不像欧拉网络法^[2]受网格约束，对固液边界、多相流体模拟以及流体的自由表面具有更好的模拟效果。主流的流体模拟方法有传统物理模拟方法和数据驱动模拟方法。传统物理方法，如SPH^[3−4]，本质上依赖一组可行的物理参数和边界条件进行数值计算，由于模拟过程受物理方程约束，结果较精确，但对于不同的场景需设置不同参数，以及进行高精度细节模拟时需增多粒子数量^[5−6]，导致需要高昂的计算成本和存储成本；数据驱动方法^[7]，利用神经网络取代求解复杂物理方程的耗时过程，并且训练好的神经网络具有一定自适应能力，避免为不同场景单独设置参数，但模拟过程中由于缺乏对局部稀疏流体的捕获能力，导致局部细节模拟能力较差，容易造成穿模现象，难以长时间维持模拟的稳定性。

本文结合传统物理模拟方法和数据驱动方法的优点，设计一套具有高效、高精度的流体模拟框架。仿真实验结果表明，本文所提出的方法在保持高效计算的同时，模拟精度比数据驱动方法高10%左右。

1 基于连续卷积的流-固特征提取网络 1.1 概　述

为了更好地提取粒子的运动特征，将粒子分为3种类型进行训练，分别是固体粒子、流体粒子及将固体粒子和流体粒子混合，统一为全局粒子，设置全局粒子的原因是为了更好提取全局特征。流体粒子用$p_i^{n*}$表示，固体粒子用${s_j}$表示。通过设计3个分支网络学习流体的运动特征以及流体粒子与环境（固体粒子）的交互现象，分支1为基于连续卷积^[8]的流体粒子和固体粒子的运动特征提取网络；分支2为基于MLP的全局粒子特征提取网络及加强局部少量流体模拟的局部特征提取网络。

流体运动特征提取网络如图1所示，网络的输出是位置的修正量$\Delta x$，通过位置修正量得出下一帧粒子的具体位置。流体粒子在投入网络之前做插值处理，预测出大概的速度及位置；固体粒子由于不会运动，不需要做额外处理。网络的输入是流体粒子的速度、位置以及粒子本身携带的物理量。时间步为n，粒子为i的速度及位置的插值计算式为：

$ {\boldsymbol{v}}_i^{n*} = {\boldsymbol{v}}_i^n + {{\Delta }}t{a_{{\mathrm{ext}}}}，$

(1)

$ {\boldsymbol{x}}_i^{{{n*}}} = {\boldsymbol{x}}_i^{{n}} + {{\Delta }}t\frac{{{\boldsymbol{v}}_i^n + {\boldsymbol{v}}_i^{n*}}}{2}。$

(2)

式中：矢量${{\boldsymbol{a}}_{{\mathrm{ext}}}}$为加速度，通过加速度可应用外力来控制流体或简单的应用重力。${\boldsymbol{v}}_i^{n*}$为时间步n、粒子i的预测速度。${\boldsymbol{x}}_i^{{{n*}}}$为时间步n、粒子i的预测位置。此时，流体粒子的输入为$({\boldsymbol{x}}_i^{{n}},{\boldsymbol{v}}_i^n,t)$，即插值预测的三维位置、三维速度；固体粒子的输入为$({\boldsymbol{x}}_i^{{n}},{\boldsymbol{x}}_{{\mathrm{normal}}}^{{n}})$，即粒子本身携带的三维位置以及其法向量；全局粒子的输入是$({\boldsymbol{x}}_i^{{n}},{\boldsymbol{v}}_i^n,{\boldsymbol{x}}_{{\mathrm{normal}}}^{{n}},t)$，包含了流体粒子插值后三维位置、三维速度、固体粒子本身携带的三维位置、法向量与时间戳。经过网络后，输出位置的修正量$\Delta x$，基于插值预测出的位置与速度，结合该位置修正量得出n+1帧粒子的准确预测位置，计算式为：

$ {\boldsymbol{x}}_i^{{{n}} + 1} = {\boldsymbol{x}}_i^{{{n*}}} + {{\Delta }}{{\boldsymbol{x}}_i}，$

(3)

$ {\boldsymbol{v}}_i^{n + 1} = \frac{{{\boldsymbol{x}}_i^{n + 1} - {\boldsymbol{x}}_i^n}}{{\Delta t}} 。$

(4)

1.2 网络模型搭建

基于连续卷积的流、固粒子运动特征提取网络包含特征聚合器和特征提取网络的搭建。特征聚合器的作用是聚合流体粒子、固体粒子以及全局粒子的特征，与直接将特定粒子的特征输入到后续网络相比，聚合后的特征更适合学习流体粒子的运动特征以及流、固粒子之间的交互。

特征聚合器的作用在于聚合3种类型粒子的特征，分别使用单层连续卷积神经网络处理流体与固体粒子，使用MLP处理全局粒子，其网络框架图如图2所示。其中，流体粒子的输入为$({\boldsymbol{x}}_i^{{n}},{\boldsymbol{v}}_i^n,t)$，即流体粒子的三维位置、三维速度；全局粒子的输入是$({\boldsymbol{x}}_i^n,{\boldsymbol{v}}_i^n,{\boldsymbol{x}}_{{\mathrm{normal}}}^n,t)$，即流体粒子的三维位置、三维速度、固体粒子的三维位置与时间戳；固体粒子的输入为$({\boldsymbol{x}}_i^{{n}},{\boldsymbol{x}}_{{\mathrm{normal}}}^{{n}})$，即粒子的三维位置以及其法向量。流体粒子、固体粒子经过连续卷积神经网络后，均输出32维的特征；全局粒子经过MLP神经网络后，输出32维的特征；最后将流体粒子的32维特征、固体粒子的32维特征与全局粒子的32维特征concat，形成96维的综合特征。这个96维的综合特征向量包含了流体粒子、固体粒子和全局粒子的运动特征。神经网络基于综合特征向量进行学习，将比单独基于每个特征进行学习的效果更好。

针对流体粒子和固体粒子交互现象的模拟，采用连续卷积神经网络学习流固之间的交互特征，也即学习流固边界处理。为了处理流体粒子与场景的碰撞，定义固态粒子为${s_j}$，法向量${\boldsymbol{x}}_{{\mathrm{normal}}}^{}$和三维位置${{\boldsymbol{x}}_i}$作为特征向量，即${s_j} = ({\boldsymbol{x}}_i^n,{\boldsymbol{x}}_{{\mathrm{normal}}}^{})$，网络实现以下功能，如式(5)所示，网络框架图如图3所示。

$ \left[ {{{\Delta }}{{\boldsymbol{x}}_1}, \ldots ,{{\Delta }}{{\boldsymbol{x}}_N}} \right] = CConvNet\left( {\left\{ {p_1^{n*}, \ldots ,p_N^{n*}} \right\} , \left\{ {{s_1}, \ldots ,{s_M}} \right\}} \right)。$

(5)

该网络的输入是经过特征聚合器后的特征，其中输入的数据格式为(B,N,96)，96代表聚合后的粒子特征，网络一共分为3层，每层都是由连续卷积神经网络与激活函数ReLU构成。为了使网络在训练过程中更好的关注上下文特征，在第二层的输出中加入第一层的输出结果，然后一并送入第三层神经网络进行训练。该网络的优点在于不仅能够学习流体粒子的运动特征，而且对于流体与环境之间的交互特征也能够学习。其原因在于经过特征聚合器后的96维特征中，包含单独的流体粒子运动特征、固体粒子的特征以及将全局的流体运动看作一个整体（这其中包含了流体粒子与环境之间的交互）后的运动特征。然后，将这些特征送入连续卷积神经网络中训练，强化了对流体运动与交互特征能力学习。同时，在训练的过程中能够及时补充上下文信息，因此具有较好的模拟效果。

2 全局粒子运动特征提取网络 2.1 全局特征提取网络

基于连续卷积的模拟主要学习了流体粒子之间的运动特征、以及流体粒子与固体粒子之间的交互特征。为了在网络训练过程中有充足的全局信息，所以单独设计全局粒子，然后用全局特征提取网络提取其运动特征，全局粒子特征提取网络的示意图如图4所示，使用多层共享感知机（MLP）^[9]，由多个全连接层构成，其在基于数据驱动的流体模拟领域表现较好。由于流体在运动过程中，遇到诸如拐角之类的障碍物，会造成如飞溅之类的局部流体，且飞溅之类的局部流体相对于整个场景来说占比较小，所以全局特征提取网络会忽略这部分流体。为了让网络在训练过程中能够关注这些局部流体，在全局特征提取网络中加入局部特征提取网络，以保证局部流体的学习效果。在本节中，全局粒子特征提取网络的设计策略为“全局-局部-全局”。

1）对全局粒子使用MLP神经网络和局部特征提取网络，分别提取其全局特征和局部特征，然后对二者的特征进行相加；

2）对二者相加后的特征再次使用MLP神经网络，强化对全局与局部特征的学习。

2.2 局部特征提取网络

在流体模拟过程中，如果遇到如河道拐弯处等场景，会产生局部的流体细节，如流体飞溅等，这些局部流体相对于全局流体占比较小，因此全局特征提取网络可能会忽略该部分流体。但是，局部区域的模拟细节对于提升整体模拟的真实感具有重要作用，如图5中圆圈圈住的范围所示，所以对于局部流体，需要专门的网络学习其特征。

局部特征提取网络采用3层^[10]局部特征的提取方式，分别是采样层、分组层和MLP层。首先，采样层从输入点中选择一组点，这些点定义了局部区域的质心；其次，分组层然后通过在质心周围寻找邻居点来构建局部区域，相当于把整体的输入划分为多个小的局部区域；然后，通过使用MLP将局部区域模式编码为特征向量，即提取出该局部区域的整体特征；最后，将多个局部区域的特征聚合在一起，以此组成整体的特征，如图6所示。

1） 采样层

给定输入的点集$ \{{x}_{1},{x}_{2},\dots ,{x}_{n}\} $，使用最远点采样方法（FPS）选定子集$ \{{x}_{i1},{x}_{i2},\dots ,{x}_{im}\} $，其中$ {x}_{ij} $是最远的点，以构成局部区域的质心。在给定相同数量质心的情况下，最远点采样比随机采样更能覆盖整个点集。该网络的输入是经过特征聚合器的输出（B,N,96），也即N个点。在该网络中，一共采样2次，第一次选取N/2个质心点，第二次选取N/4个质心点。然后以该质心点为圆心，伴随采样半径规划出各个分组。

2） 分组层

分组方法为球查询，即以质心为原点，r为半径，示意图如图7所示，查询在半径以内的所有点（以K为上限）。球查询的局部领域保证了固定的区域尺度，从而使局部区域特征在空间上更具普遍性，这对于局部区域来说是很重要的。在流体模拟中，局部区域的细节对于提升整体模拟的真实性具有重要意义。

3） MLP层

为方便理解，可把全局流体粒子看作无序的点云集合$\left\{ {{x_1},{x_2}, \ldots ,{x_n}} \right\}$，其中${x_i} \in {R^d}$，可以把网络抽象为映射函数$f:\chi \to R$，把无序的点云集合映射为向量，如下式：

$ f\left( {{x_1},{x_2}, \ldots ,{x_n}} \right) = \gamma (_{i = 1,..,n}^\mathrm{MAX}{\text{\{ }}h\left( {{x_i}} \right){\text{)}} 。$

(6)

式中：$ \gamma $和$ h $为常用的多层感知机网络（MLP），其在流体数据驱动领域效果较好，函数$ f $保证点集的置换不变性，且可逼近任意连续集合函数。

4） 跳跃连接

在经历过采样-分组-MLP层后，需要求得每个粒子在下一帧的位置，所以希望获得所有原始点的特征。在这里采用神经网络常用的跳跃连接方法（Skip Link Concatenation）和插值法，将之前的特征与现在的特征相结合，即获得了所有原始点的特征和上下文信息。

依照上述3层局部特征提取的概念，局部特征提取网络结构图如图8所示，经过2次采样和分组，共有3层MLP层用于提取特征。在训练过程中，流体粒子会在某些区域比较稀疏，为了保证增强该区域的特征，采用多分辨率分组策略（MRG）进行特征的聚合。

3 仿真实验与分析 3.1 评估指标

本实验的评估指标通过每3帧进行初始化，利用第1帧预测连续的后两帧，然后计算后续2帧预测粒子的位置与真实值粒子位置之间的偏差。后续2帧记为n+1帧和n+2帧，误差为各粒子与第1帧对应的位置误差，记为err_n1与err_n2。另外，计算整个预测序列中预测粒子与真实值粒子之间的平均误差，记为whole_seq_err，以此衡量长期的相似性。第n帧粒子之间的误差计算式为：

$ {d^n} = \frac{1}{N}\mathop{ \mathrm{min}}\limits_{{{\boldsymbol{x}}^n} \in {{\boldsymbol{X}}^n}} {\left\| {{\boldsymbol{\hat x}}_i^n - {\boldsymbol{x}}_i^n} \right\|_2}。$

(7)

式中：${{\boldsymbol{X}}^n}$为第n帧所有预测粒子位置的集合；${\boldsymbol{\hat x}}_i^n$为第i个粒子位置的真实值；${\boldsymbol{x}}_i^n$为预测粒子的位置；$ {N} $为粒子的总数量。

3.2 仿真实验结果

本实验以基础的连续卷积方法作为对比，分析本文方法的性能。在Um等^[8]采用的连续卷积方法中，对于全局粒子运动特征的提取仅使用3层MLP层。在本文中，对该方法进行了拓展。为了使网络在提取全局粒子运动特征的同时能关注到局部区域的流体运动，特设计局部特征提取网络，即使用“采样-分组-MLP”的方法提取局部特征，最后将该局部区域的特征融入到全局粒子运动特征提取网络中，后续2帧预测误差如图9所示，具体数值如表1所示。其中，err_n1是由第n帧预测的第n+1帧各粒子的平均位置误差，err_n2是由第n帧预测的第n+2帧各粒子的平均位置误差。由表1可看出，本文提出的特征聚合器模块和局部特征提取模块与Um等^[8]连续卷积相比均有提升，证明该2个模块的有效性，且本文提出的网络框架在第n+1帧、第n+2帧及全局粒子位置误差方面均比Um等^[8]连续卷积方法低，预测精度有所提升，证明该网络框架的有效性。

图10为立方体水块掉落实验的预测效果，只需给出初始水块位置，通过本文训练好的网络框架即可模拟掉落效果，共250帧，使用Blender进行渲染。从模拟效果来看，该网络框架实现了运动流体与固体之间的交互效果，模拟过程没有出现崩溃，且效果与真实现象相似。

3.3 泛化性及性能评估

为了体现本文方法的泛化性，采用多个未训练的、且更为复杂场景进行验证，掉落水块也采取不同形式、不同粒子数目的。场景如图11所示，每个场景的粒子数如表2所示；掉落水块如图12所示，每种水块的粒子数目如表3所示。通过使用训练好的模型进行不同场景搭配不同水块的形式进行验证泛化性，评估指标一方面是直观视觉效果上是否和现实有较大出入，另一方面通过网络推理时间与数值仿真时间对比、仿真结果与数值计算结果对比，以此说明本文方法的泛化性和优势。

图13为不同形状的水块组合在不同场景中的模拟效果，每个场景均进行连续800帧的模拟，从直观视觉效果上看，模拟符合真实交互效果，流体与流体、流体与固体交互顺畅。从评估指标上看，数值计算采用DFSPH方法，该方法计算较为精准，误差以此数值计算结果为基准；本文方法耗时采用网络推理时间，结果如表4所示，在连续800帧模拟时，本文方法的耗时远小于数值计算方法，针对越复杂的场景，本文方法的耗时优势越大，且误差在可接受范围之内。（计算平台硬件如下：CPU：Intel(R) Xeon(R) E5-2637 v2 @ 3.50GHz；内存125G；GPU:NVIDIA RTX 3090 24G）。

4 结　语

针对基于数据驱动流体模拟易忽略局部细节、易崩溃的问题，本文设计了一种基于连续卷积的流、固粒子运动特征提取网络，在网络中通过加入特征聚合器和局部特征提取网络，提升对流体局部细节的特征提取，从而提高整体模拟精度与稳定性。仿真结果表明，本文方法比目前的连续卷积方法表现更好，模拟精度更好，局部细节模拟效果更好，更接近真实流体情况，且仿真时间更长，具有较优秀的性能。未来将持续改进网络，增强网络的自适应能力与迁移能力。