CN112163384A - 一种面向自由表面流的固体边界提取方法 - Google Patents

Abstract

本发明提供一种面向自由表面流的固体边界提取方法，包括获取自由表面流体中固体边界附近的若干流体粒子组成的流体计算域，计算固体边界的固体体积积分；将固体体积积分转化为固体边界与各流体粒子支撑域相交面上的表面积分，并使用若干三角形代表相应表面积分；生成一均匀网格结构存储若干三角形，获取每一均匀网格单元存储的三角面片；通过各流体粒子所属的均匀网格单元进行邻域查找，获取邻域内的各三角形索引，并对各三角形按照所在平面进行聚类；计算每一流体粒子的每一类三角形相应的固体边界积分，得到自由表面流的固体边界。本发明通过对所有三角形聚类，避免了在相邻三角形之间引入不连续性，能够处理复杂的实心边界。

Description

一种面向自由表面流的固体边界提取方法

技术领域

本发明属于计算机图形学领域，具体涉及基于物理的流体仿真动画应用中一种面向自由表面流的固体边界提取方法。

背景技术

平滑粒子流体动力学(SPH)非常适合模拟猛烈的现象，包括水溅，多相流以及粘滞流体。然而，在SPH中，处理坚固的墙体边界仍然是最具挑战性的部分之一。在过去的三十年中，为了在SPH中实施实心墙，已开发出各种不同的策略(例如，幽灵粒子[Nadir Akinci,Markus Ihmsen,Gizem Akinci,Barbara Solenthaler and MatthiasTeschner.2012.Versatile rigid-fluid coupling for incompressible sph.ACMTransactions on Graphics(TOG)31,4,2012]，密度图[Dan Koschier,JanBender.2017.Density maps for improved sph boundary handling.In Proceedings ofSCA,2017]，半分析边界[Martin Ferrand,Dominique Laurence,Benedict D Rogers,Damien Violeau,Christophe Kassiotis.2013.Unified semi-analytical wallboundary conditions for inviscid,laminar or turbulent flows in the meshlesssph method.International Journal for Numerical Methods in Fluids71,4,2013]等)，每个都有其优点和缺点。

根据标准SPH公式，如果粒子位于固体边界附近，只有流体域内部的粒子才有助于粒子相互作用的总和。这种单面的粒子近似不能给出正确的解，因为在固体边界内的场变量通常不能按照零计算(如在计算粒子密度时)。要开发一种实用的固体边界处理技术，应优先考虑准确性和效率。在所有固体边界处理技术中，幽灵粒子是最常用的。但是，应用幽灵粒子并不是处理较大平面区域的理想方法。此外，要完全消除颗粒不足的问题，需要使用多层幽灵粒子进行采样。不幸的是，大量的幽灵粒子对仿真性能有很大的副作用。半分析方法旨在将边界上的体积积分转换为表面积分以减少计算成本。但是现有方法不能处理复杂的边界几何或仅适合于采用梯度项的边界积分。

发明内容

本发明提出一种面向自由表面流的固体边界提取方法，可以处理复杂边界几何、避免在相邻三角形之间引入不连续性、提升计算效率。在SPH方法自由表面流体中，对于固体边界附近的流体粒子，计算固体边界时采用体积积分，然后将其转换为表面积分，以便于代表边界的计算机辅助设计(CAD)网格文件自然地集成到粒子模拟中。

具体来说，本发明的技术方案如下：

一种面向自由表面流的固体边界提取方法，其步骤包括：

1)获取自由表面流体中固体边界附近的若干流体粒子组成的流体计算域，计算固体边界的固体体积积分；

2)将固体体积积分转化为固体边界与各流体粒子支撑域相交面上的表面积分，并使用若干三角形代表相应表面积分；

3)生成一均匀网格结构存储若干三角形，获取每一均匀网格单元存储的三角面片；

4)通过各流体粒子所属的均匀网格单元进行邻域查找，获取邻域内的三角面片，并对每一流体粒子的各三角面片对应的三角形按照所在平面进行聚类；

5)计算每一流体粒子的每一类三角形相应的固体边界积分，得到自由表面流的固体边界。

进一步地，通过以下步骤计算流体计算域内的固体体积积分：

1)计算域上的平滑函数f(x)：

其中，Ω表示积分域，ω为形函数的连续形式，

表示流体计算域

内的加权积分，

表示固体边界

的加权积分，h为支撑半径，x′是x的支撑域内的采样点；

2)计算固定体积积分g(x,x′)≡f(x′)ω(x′-x)。

进一步地，表面积分

其中

表示相对于x′的散度算符。

进一步地，每一均匀网格单元存储一相交计数器，用以指示与该网格单元相交的三角形数量；任意p个均匀网格单元，生成一计数器数组，p≥2。

进一步地，对于静态固体边界，通过以下步骤进行邻域查找，获取邻域内的三角形索引：

1)获取与每一均匀网格单元相交的三角形数量；

2)对于每个三角形，获取与该三角形相交的均匀网格单元，将对应的相交计数器加一；

3)选择该三角形附件的均匀网格单元，依据对应的相交计数器生成一计数器数组，并在该计数器数组取一个前缀和；

4)将所有三角形与所有网格的相交总数返回为I，并分配一个大小为I的附加数组。

5)获取各三角形的索引，并将各索引存储在所述附加数组中；

6)根据每一流体粒子的所属均匀网格单元，遍历与所属均匀网格单元相邻的均匀网格单元，获取邻域内的三角形索引。

进一步地，对于动态固体边界，通过多次领域查找，得到邻域内的三角形索引。

进一步地，通过以下步骤得到自由表面流的固体边界：

1)计算每一流体粒子的每一类三角形在相应流体粒子的支撑域内面积；

2)计算每一三角形s_i相对于相应流体粒子中心位置的立体角

其中1≤i≤M，M为三角网络结构中三角形的数量，n_s为三角形s_i的法线向量，d_n为最近点x_n到相应流体粒子位置x_i的归一化矢量，A_s为聚类后的该类三角形与相应流体粒子的支撑域之间的总相交区域的面积，A₀为三角形s所在平面与相应流体粒子的支撑域之间的相交区域，

l_n为从相应流体粒子位置x_i到三角形所在平面的距离，h为相应流体粒子支撑域半径，Ω₀为A₀的立体角，

3)依据每一流体粒子相应的支撑域内面积与立体角，得到自由表面流的固体边界B的积分

其中，G(r)满足G′(r)＝g(r)r²，r＝||x′-x||的单变量函数，且

G'(r)是G对r的导数，

进一步地，通过以下步骤计算每一类三角形在相应流体粒子的支撑域内面积：

1)计算每一个三角形在相应流体粒子支撑域内的顶点数量；

2)计算每一个三角形与相应流体粒子支撑域相交的边数量；

3)通过顶点数量与边数量，判断每一个三角形所在平面与相应流体粒子支撑域的相交形态，计算每一个三角形在相应流体粒子支撑域内的面积；

4)对所述面积求和，获取每一类三角形在相应流体粒子的支撑域内面积。

一种存储介质，所述存储介质中存储有计算机程序，其中，所述计算机程序被设置为运行时执行上述所述的方法。

与现有技术相比，本发明具有如下的优点：

(1)本发明提出在固体边界内采用体积积分，在每个粒子的局部球形积分域下重新构造积分，然后将其转换为表面积分，使得本发明能够处理复杂的实心边界，例如尖角和薄壳；

(2)本发明将属于同一平面的所有三角形聚类在一起，并对每个聚类应用相同的积分方法，避免了在相邻三角形之间引入不连续性；

(3)为了加快对粒子相邻三角形的查询，使用统一的网格来存储三角形索；在三角形和网格单元之间进行相交测试，然后在发生相交时将三角形索引存储到相应的网格单元中；在仿真过程中，可以在现代GPU上并行执行相邻三角形的查询；

(4)本发明可以集成到基于位置的方法和基于投影的不可压缩SPH解算器中，实验表明，本发明和用幽灵粒子得到的结果可以达到可比，且性能更好，尤其是对于大型平滑边界。

附图说明

图1为本发明的流程图。

图2为邻域三角网格查找示意图。

图3为三角形和平面圆相交可能导致的所有9种不同情况。

图4为粒子的支撑域和固体边界之间的示意图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实现和附图对本发明做详细说明,但不构成对本发明的限制。

本发明方法的硬件平台采用了型号为Intel i7四核CPU，主频为3.46GHz，NVIDIAGeForce GTX 1060显卡，现存6GB。系统程序采用C++编写，其中对于并行计算部分采用了CUDA语言进行加速，并借助Microsoft Visual Studio 2019对程序进行编译执行，开发基于开源项目PhysIKA(https://github.com/PhysikaTeam/PhysIKA)。流体粒子和三角形网格之间加入连续碰撞检测(CCD)，防止流体粒子渗透到固体中。

请参考图1，本发明方法的流程主要分为以下几个部分。在边界处理之前，使用均匀网格存储三角网格进行邻域查找加速。接下来，对于每个流体粒子，将其邻域内的三角面片按照所在平面进行聚类。对于聚类后的每一组面片，计算每个三角面片和支撑域的面积交集求和，进而计算该类面片对应的固体角以及积分后的形函数G。最后，将支撑域内每一组面片的贡献求和，得到固体边界修正。接下来将详细介绍每一个步骤的具体实现方案。

1、理论基础。在标准SPH方法中，连续性用下标域

表示的一组粒子表示。考虑在域

上定义的平滑函数f，粒子i上f的标准SPH插值由

给出，其中j表示所有相邻粒子，V_j是j粒子的体积，f_j＝f(x_j)，且形函数W_ij＝W(‖x_i-x_j‖,h)，且支撑半径为h。如果粒子位于固体壁边界附近，则函数f的积分表示具有以下恒等式:

其中x’是x的支撑域内的采样点，Ω表示积分域，ω为形函数的连续形式，

表示流体计算域

内的加权积分，

表示固体边界B的加权积分。

的计算较为简单，

难点在于有效而准确地计算

为了方便下面的讨论，记：

g(x,x′)≡f(x′)ω(x′-x)

根据[Martin Ferrand,Dominique Laurence,Benedict D Rogers,DamienVioleau,Christophe Kassiotis.2013.Unified semi-analytical wall boundaryconditions for inviscid,laminar or turbulent flows in the meshless sphmethod.International Journal for Numerical Methods in Fluids 71,4,2013]中的讨论，如果存在函数G(x,x′)满足：

其中，

表示相对于x′的散度算符，边界

上的体积积分g(x,x′)可以转化为固体边界与流体粒子支撑域的相交面

上的表面积分G(x,x′)。不幸的是，找到这样的G(x,x′)绝不是一件容易的事。但是，如果f属于径向基函数(RBF)，即f(x′)＝f(||x′-x||)，推导过程可以大大简化。

SPH的形函数ω是径向基函数，因此可以将

内部的f积分重写为：

其中Ω是积分边界

从x视角的立体角，G(r)是一个满足G′(r)＝g(r)r²，r＝||x′-x||的单变量函数，且

是微分球面坐标的公式，原点到积分点的连线与正z轴之间的天顶角为θ，与正x轴之间的方位角为

h代表支撑域半径。

2、邻域三角网格查找。在查找邻域的三角网格时，对由N个粒子和M个三角网格组成的场景进行朴素的蛮力搜索具有O(MN)的计算复杂度，这对于大规模仿真而言太重了。因此，使用类似于最近粒子查找的策略来加速三角网格查找。对于动态边界，GPU(图形处理器)上的实现可分为以下几个阶段：

(1)分配一个统一的三维网格，该网格足够大以覆盖所有边界三角形以及相同尺寸的数组，以存储相交计数器，指示每个网格单元有多少个三角形相交；

(2)对于每个三角形，找到所有与三角形相交的网格单元，并将对应的相交计数器加一；

(3)在计数器数组上取一个前缀和，并将所有的三角形与所有网格的相交总数返回为I；

(4)如图2所示，分配一个大小为I的附加数组。对于每个三角形，相交计数器不清零，再次执行步骤(2)，并将三角形索引存储在一个附加数组中，以提供之后三角面片邻域查找使用邻域。对于每个粒子，找到其所属的网格单元，遍历存储在最近的27个(在二维情况下为9个)网格单元中的所有三角形索引，并删除重复的三角形索引。

对于静态边界，前四个步骤仅执行一次，并且可以在仿真开始时进行预先计算。

3、三角网格按照所在平面聚类。因为采样点选取的问题，位于同一平面内的三角网格边界处可能会产生不连续的现象。因此将粒子支撑域内所有三角面片按照其所在平面进行聚类，聚类后同一类三角面片选取相同采样点。下文中所有固体边界积分采样点均为使用聚类后的三角面片采样点。

4、计算三角网格在支撑域内的面积。首先在粒子i的支撑域和包含三角形s的平面之间进行相交测试。如果从粒子i到三角所在平面的距离大于或等于h，则粒子i的支撑域和三角形之间不应发生交集。否则，可以简化平面圆和三角形之间的交点问题。图3展示了三角形和平面圆相交可能导致的所有9种不同情况。计算主要步骤如下：

(1)计算网格三角位于支撑域内的顶点个数；

(2)计算网格三角与支撑域相交的边个数；

(3)如果顶点个数、交边个数均为0，特殊判断三角所在平面是否与支撑域相交以区分图3所示的①、②情况，接下来按照对应情况进行面积计算；

(4)重复步骤(1)到步骤(3)并对结果求和，直到聚类在同一类的不同三角全部计算完成。

5、计算固积分。如图4所示，考虑在支撑域内的三角面片，通过应用单点高斯正交规则，边界积分

在无穷小三角形s上为以下近似值：

其中x_s表示单点高斯积分的采样点，Ω_s表示三角形s相对于粒子位置x_i的立体角。可以注意到，

的精度取决于如何选择采样点x_s以及如何计算Ω_s。计算Ω_s的标准公式为：

其中n_s是三角形s的法线向量，A_s表示三角形面积，l_s表示从x_i到三角形s的距离，d_s是从三角形s上最近点指向x_i的归一化矢量。但是实验表明这种计算方法受三角面片分辨率影响较大。为了在面片稀梳的情况下仍然能够精确计算，描述固体三角形s相对于粒子位置x_i的立体角Ω_s时采用：

其中n_s是三角形s的法线向量，A_s表示聚类后的该类三角形与粒子的支撑域之间的总相交区域的面积，A₀是三角形s所在平面与粒子的支撑域之间的相交区域，Ω₀表示A₀的立体角，x_n是粒子位置x_i到三角形s所在平面上的最近点，d_n表示最近点x_n到粒子位置x_i的归一化矢量。计算A₀和Ω₀的公式写为：

其中l_n表示从粒子位置x_i到三角形所在平面的距离，h代表支撑域半径。接下来，固体边界表面积分的计算方法就可以表示为：

其中，G(r)满足其G′(r)＝g(r)r²，r＝||x′-x||的单变量函数，且

G'(r)是G对r的导数，计算G(r)所用的采样点x_s取三角形所在平面上的最近点x_s＝x_n。

以上实施方示例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求所述为准。

Claims (9)

1.一种面向自由表面流的固体边界提取方法，其步骤包括：

1)获取自由表面流体中固体边界附近的若干流体粒子组成的流体计算域，计算固体边界的固体体积积分；

2)将固体体积积分转化为固体边界与各流体粒子支撑域相交面上的表面积分，并使用若干三角形代表相应表面积分；

3)生成一均匀网格结构存储若干三角形，获取每一均匀网格单元存储的三角面片；

4)通过各流体粒子所属的均匀网格单元进行邻域查找，获取邻域内的各三角形索引，并对各三角形按照所在平面进行聚类；

5)计算每一流体粒子的每一类三角形相应的固体边界积分，得到自由表面流的固体边界。

2.如权利要求1所述的方法，其特征在于，通过以下步骤计算流体计算域内的固体体积积分：

1)计算域上的平滑函数f(x)：

其中，Ω表示积分域，ω为形函数的连续形式，

表示流体计算域

内的加权积分，

表示固体边界

的加权积分，h为支撑半径，x′是x的支撑域内的采样点；

2)计算固定体积积分g(x,x′)≡f(x′)ω(x′-x)。

3.如权利要求2所述的方法，其特征在于，表面积分G(x,x′)＝g(x,x′)/▽_x′，其中▽_x′表示相对于x′的散度算符。

4.如权利要求1所述的方法，其特征在于，每一均匀网格单元存储一相交计数器，用以指示与该网格单元相交的三角形数量；任意p个均匀网格单元，生成一计数器数组，p≥2。

5.如权利要求4所述的方法，其特征在于，对于静态固体边界，通过以下步骤进行邻域查找，获取邻域内的三角形索引：

1)获取与每一均匀网格单元相交的三角形数量；

2)对于每个三角形，获取与该三角形相交的均匀网格单元，将对应的相交计数器加一；

3)选择该三角形附件的均匀网格单元，依据对应的相交计数器生成一计数器数组，并在该计数器数组取一个前缀和；

4)将所有三角形与所有网格的相交总数返回为I，并分配一个大小为I的附加数组。

5)获取各三角形的索引，并将各索引存储在所述附加数组中；

6)根据每一流体粒子的所属均匀网格单元，遍历与所属均匀网格单元相邻的均匀网格单元，获取邻域内的三角形索引。

6.如权利要求5所述的方法，其特征在于，对于动态固体边界，通过多次领域查找，得到邻域内的三角形索引。

7.如权利要求1所述的方法，其特征在于，通过以下步骤得到自由表面流的固体边界：

1)计算每一流体粒子的每一类三角形在相应流体粒子的支撑域内面积；

2)计算每一三角形s_i相对于相应流体粒子中心位置的立体角

其中1≤i≤M，M为三角网络结构中三角形的数量，n_s为三角形s_i的法线向量，d_n为最近点x_n到相应流体粒子位置x_i的归一化矢量，A_s为聚类后的该类三角形与相应流体粒子的支撑域之间的总相交区域的面积，A₀为三角形s所在平面与相应流体粒子的支撑域之间的相交区域，

l_n为从相应流体粒子位置x_i到三角形所在平面的距离，h为相应流体粒子支撑域半径，Ω₀为A₀的立体角，

3)依据每一流体粒子相应的支撑域内面积与立体角，得到自由表面流的固体边界B的积分

其中，G(r)满足G′(r)＝g(r)r²，r＝||x′-x||的单变量函数，且

G'(r)是G对r的导数。

8.如权利要求7所述的方法，其特征在于，通过以下步骤计算每一类三角形在相应流体粒子的支撑域内面积：

1)计算每一个三角形在相应流体粒子支撑域内的顶点数量；

2)计算每一个三角形与相应流体粒子支撑域相交的边数量；

3)通过顶点数量与边数量，判断每一个三角形所在平面与相应流体粒子支撑域的相交形态，计算每一个三角形在相应流体粒子支撑域内的面积；

4)对所述面积求和，获取每一类三角形在相应流体粒子的支撑域内面积。

9.一种存储介质，所述存储介质中存储有计算机程序，其中，所述计算机程序被设置为运行时执行权利要求1-8中任一所述方法。

Images