一种基于变光滑长度的SPH流体模拟方法与流程

本发明属于计算机图形学领域，具体涉及一种基于变光滑长度的sph流体模拟方法。

背景技术：

流体遇到冲击会产生很大的变形，针对模拟大变形流体问题，sph方法不仅可以处理大变形问题而且能够对粒子进行效果很好的追踪，非常适合对大变形流体的模拟。在sph方法中，使用粒子近似法主要依赖于一个以光滑半径为长度构造出的粒子支持域，计算效率和精度取决于邻近粒子的数量。光滑长度的选取决定局部核函数近似的精度，直接影响流体模拟的计算速度和真实性。若光滑长度过小，则在支持域中没有足够的粒子对给定粒子提供插值信息和施加作用力，导致模拟结果的精度降低。若光滑长度过大，则用于更新给定粒子的邻近粒子部分信息或局部特性将会丢失，也会导致模拟结果的精度降低。

传统sph方法从光滑长度与粒子密度关系入手，将光滑长度看成受恒定临近粒子数约束条件的独立坐标变量，提出全守恒sph方程组。这种方法存在方程组复杂、数值误差较大以及适用范围窄的问题，不具有通用性和实用性。基于对前人的研究，本发明提出一种从对称核近似角度出发建立修正后的变光滑长度sph方程组，模拟遇到冲击后的变形流体。

技术实现要素：

为了从根本上消除变光滑长度效应在求解流体受到冲击产生大变形问题时造成的偏差问题，本发明从对称核近似角度出发修正后的变光滑长度sph方程组。同时在gpu上采用一种并行邻域粒子搜索，实现对粒子搜索的加速。发明具体技术方案如下：

一种基于变光滑长度的sph流体模拟方法包括如下步骤：

步骤一，从对称核近似角度出发修正sph密度演化方程、动量方程和能量方程，建立一组修正后的变光滑长度sph方程组；

步骤二，由建立的sph方程组计算出粒子所受的压力、重力、黏性力，再根据牛顿第二定律计算出合力得到加速度；

步骤三，更新下一步粒子的速度和位置；

步骤四，基于marchingcubes算法实现流体表面提取，利用gpu纹理缓存的位置、速度和颜色数据进行三维可视化渲染；

marchingcubes流体表面提取算法涉及网格节点处的密度计算、交点坐标和法矢量计算。采用cuda并行计算网格节点处的密度值和粒子网格节点密度值，计算思想由影响域改为支持域，即当前网格节点的密度值由邻域范围内所有粒子计算所得，而不是计算一个粒子点对邻域范围内产生的影响。

所述步骤一具体修正过程如下：

(1)密度方程

传统sph法中，某点ri处密度ρ(ri)的核函数近似式为

其中，w(ri-rj,h)称为核函数，h为光滑长度描述核近似的程度，n为i粒子邻近粒子总数。

光滑长度h的选取，应尽量使得粒子的邻近数目在计算中保持不变，以获得全场一致的核近似精度，通常设定粒子i的光滑长度hi半径区域内包含恒定质量msph，即：

其中，d为维数，σ对应一维到三维分别取2、π、公式(2)的引入使得光滑长度在空间分布不一致，并随同粒子的运动变化，即光滑长度同为空间和时间的函数hi＝h(ri,t)。

对于变光滑长度核函数近似，采用对称核近似的方法可以得到更高精度的近似值，本发明采用对称光滑长度的方法来对对称函数近似，此时密度核函数近似式为：

其中，

公式(3)两边分别对时间求导数得到：

公式(4)即为用对称核近似的变光滑长度sph密度演化方程，其中νij＝νi-νj，在恒定光滑长度时该方程退化为传统sph密度演化方程。

(2)动量方程

从全守恒sph方程组思想出发，改用对称的核函数近似来推导动量方程。无耗散流体动力学方程的拉格朗日函数描述如下：

其中，ρ为密度，ν为速度，u为单位质量内能，s为熵。对绝热等熵过程，由热力学第一定律可知：

其中，p为压强。

公式(5)的sph离散形式为

其中，q＝(r1,...,rn,h1,...,hn),包含坐标和光滑长度，这一近似过程可理解为把无限维连续介质相空间离散到2n维的有限粒子相空间。

对q的后n个坐标分量—光滑长度引入公式(2)，使所有粒子的光滑区域内包含恒定质量物质。其约束条件描述如下：

采用带约束的拉格朗日方程：

令由公式(9)前n个方程以及梯度核函数可得

公式(10)即为采用对称核近似的变光滑长度sph动量方程，fi为修正系数。

(3)能量方程

由热力学第一定律可知

其中，t为温度，s为熵，u为单位质量内能，p为压强，v为比体积。若没有热源产生，方程可以简化为：

公式(12)的sph离散形式为采用对称核近似的变光滑长度sph能量方程。

所述步骤二为有效计算压力，将不可压缩流体视为若弱可压缩流体。引入弱可压状态方程：

其中，参数p0为参考压强，γ为常数。p0和γ共同用于控制计算中流体密度，以保证模拟流场的不可压缩性。

所述步骤三采用蛙跳方式进行时间步长推进。

所述步骤四计算思想由影响域改为支持域后对应粒子搜索方式为并行邻域粒子搜索，实现对粒子搜索的加速。

附图说明

图1为基于变光滑长度的sph流体模拟总体流程图；

图2为并行粒子搜索示意图；

图3为实验参数描述图；

图4为楔形体入水效果图；

图5为角加速度示意图。

具体实施方式

下面将结合附图和实施例对本发明进一步说明。

参照图1，一种基于变光滑长度的sph流体模拟方法，具体具体包括以下步骤：

步骤一，修正sph密度演化方程、动量方程和能量方程，建立适合模拟流体大变形的修正后变光滑长度sph方程组；

1.1采用对称光滑长度的方法来对对称函数近似，得到密度核函数近似式。

1.2核函数两边分别对时间求导数得到密度演化方程。

1.3对全守恒流体动力学的拉格朗日方程离散化，离散过程近似理解为把无限维连续介质相空间离散到2n维的有限粒子相空间。

1.4对n个坐标分量引入光滑长度，使所有粒子的光滑区域内包含恒定质量物质。

1.5由前n个方程以及梯度核函数可得到采用对称核近似的变光滑长度sph动量方程。

1.6简化由热力学第一定律描述的内能公式，离散成为对称核近似的变光滑长度能量方程。

步骤二，由于光滑长度变化率和密度相关，而密度方程又相关于其所有邻近粒子光滑长度，这使得方程组难以进行显式求解。因此采用迭代方法求解新建立的sph方程组，并由此计算出粒子所受的压力、黏性力得到粒子加速度；

2.1考虑到光滑长度h为空间和时间的函数，对光滑长度的约束公式求时间导数得到：

2.2设定求解密度公式，得到的密度变化率

2.3求解得到的密度变化率代入约束公式求解再代入求出的到密度公式求出修正的密度变化率

2.4重复上面的迭代过程，直到不再变化后求解动量方程和能量方程。

2.5压力和黏性力计算公式如下：

其中，vi和vj分别为粒子i和j的体积。

2.6粒子加速度由求得。

步骤三，在第一个时间步长结束时，速度和内能由初始状态向前推进半个时间步长，而粒子的位置向前推进一个时间步长。在随后每一个时间步的开始，粒子的速度和内能再向前推进半个时间步长，获得整数时间步上的值，以便与粒子的位置时刻保持一致。速度和位置分别按以下公式计算得到：

步骤四，流体渲染利用gpu纹理缓存的位置、速度和颜色数据进行三维可视化渲染。数据存储在gpu全局存储器中，调用cudagrapuiccsmapresources和cudamappedpointer函数映射到vbo缓冲区，获得该缓冲区的gpu指针，然后直接进行流体渲染。

流体表面提取基于marchingcubes算法，首先通过调用cudagraphicsglregisterbuffer指令，使用cuda中openglvbo缓冲区，同时该函数返回一个指向vbo缓冲区的句柄。具体提取本文设计四个kernel函数，第一个kernel函数根据粒子的位置将粒子插入至三维空间网格中，构成粒子—网格对。第二个kernel函数为并行计数排序函数，使粒子—网格对按照网格编号进行排序，形成网格—粒子对。第三个kernel函数将它们放到同一内核中求解，并进行粒子更新。第四个kernel函数在求取各流体粒子内力、外力后，便可进行边界处理提取出流体轮廓，这样设计多个核函数的优点是每次计算中只需要进行一次优化估值。

但是marchingcubes算法需耗费大量的计算资源，该算法最为耗时的操作在于数据处理阶段，特别是大量的空体素占据了计算资源。在计算思想由影响域改为支持域后对应粒子搜索方式为并行邻域粒子搜索，实现对粒子搜索的加速。邻域粒子搜索算法采用三维空间网格对粒子进行如图2所示的网格划分后，查找近邻粒子时仅搜索粒子所在网格邻域内的27个网格，步骤如下：

(1)对三维空间进行均匀网格划分，并根据每个网格在空间中的位置给每个网格分配一个整型变量作为网格编号，分配若干网格缓冲区，包括粒子数、偏移地址。

(2)分配若干粒子缓冲区，包括位置、速度、颜色等，并按粒子的空间位置插入至网格中，建立粒子编号与网格编号之间的关系。

(3)采用并行前缀求和算法求得每个网格中的粒子数量和偏移地址，接着采用并行计数排序算法按网格编号排列粒子缓冲区。其中，并行计数排序算法和前缀求和算法均采用cuda并行加速策略，以提高程序的性能。

(4)求取计算粒子的邻近27个网格单元，对于网格单元中的每个近邻粒子，判断其是否在支持域内，若近邻粒子与计算粒子的距离小于光滑核半径，则进行求和计算；否则，该近邻粒子不参与物理计算。

实施例

楔形体入水冲击是一种典型的流固耦合问题，带初始倾斜角的入水冲击问题耦合了大变形、多自由度等现象，其模拟难度较大。其中角加速度是需要重点关注的信息，选取的计算模型及相关参数见图3。

首先采用传统sph方法对粒子均匀分布模型进行模拟，此时粒子间距均为0.005，光滑长度均为h＝1.23δx，δx＝0.01，粒子总数为265277个，时间步长取10^-5s，状态方程中水的声速值取14(gh)^0.5，h为水深1.8m。固壁边界压力采用sph粒子近似法得到。然后在实验条件不变的情况下，采用本文变光滑长度方法模拟楔形体入水冲击。

图4模拟了楔形体入水冲击后水面产生大变形实验效果图。图左和图右展示的分别是本文变光滑长度方法和传统sph恒定光滑长度方法的模拟结果。在自由表面处，大量粒子基于本文算法动态调整光滑长度，提升了插值精度。相比图右，捕捉到了更丰富的流体细节和变形程度，增强了流体模拟的真实感。

同样，图5可得恒定光滑长度算法得到的计算结果较差，这是由于粒子非均匀分布导致角加速度不稳定性造成的。但当采用了本文变光滑长度方法后，角加速度与实际情况计吻合得都很好，在0.38s后流场紊乱程度在图中体现出来。这是因为粒子相互作用对称化思想，给每个粒子单独配置光滑长度，其支持域内的粒子数量能够保证基本一致。

最后，对楔形体入水场景分别利用cpu和gpu模拟，其中cpu和gpu版本物理计算的邻域粒子搜索算法分别为邻近网格的粒子搜索法和cuda并行邻域粒子搜索法。物理计算时间和加速比如表1所示。

表1不同粒子数的sph物理计算时间

可以看出通过并行化处理，本发明获得了最高接近57倍的加速比，提升了模拟效率。体现出并行邻域搜索算法的快速性，同时也说明采用并行搜索算法可以成为模拟变光滑长度流体的有效工具。