CN101329772B - 一种基于sph的运动物体与水交互的仿真建模方法 - Google Patents

Abstract

本发明涉及一种基于SPH的运动物体与水交互的仿真建模方法，属于计算机图形学与虚拟现实领域。本发明的方法包括以下步骤：一、对仿真进行初始化，包括粒子属性的初始配置和建立粒子的相邻粒子列表；二、对水粒子的相邻粒子列表进行更新；三、计算水粒子在下一步时间步的位置和运动状态；四、计算物体粒子在下一步时间步的位置和运动状态；五、判断粒子与障碍物是否发生碰撞，如果发生碰撞，则改变粒子的运动速率和方向，否则直接执行第六步；六、基于marching cubes算法构建水面拓朴结构；七、采用图形绘制语言，进行水面和物体的绘制和显示；八、判断仿真是否停止运行。使用本发明所述方法可以模拟运动物体与水的交互现象，获得比较真实的仿真结果。

Description

一种基于SPH的运动物体与水交互的仿真建模方法

技术领域

本发明涉及一种仿真建模方法，特别涉及一种基于SPH的运动物体与水交互的仿真建模方法，属于计算机图形学与虚拟现实领域。

背景技术

水环境仿真长期以来一直是虚拟现实领域的重要研究方向，它有助于提高虚拟环境的表现能力，在计算机游戏、影视动画、军事仿真等领域发挥着重要作用。物体与水的交互作用是自然界中的常见现象，同时也是水环境仿真的研究热点和难点。它是一个复杂的流体运动过程，需要采用物理建模方法才能真实地表现物体和水的运动。目前，计算机图形学中关于物体与水交互的物理建模方法主要有两种，一种是欧拉方法，另一种是拉格朗日方法。

欧拉方法是一种基于网格的方法，在流体仿真中应用广泛。但欧拉方法有一个固有的缺欠，就是在对流场进行网格化求解时存在数值耗散，导致流体在质量上不能守恒，亦即随着时间的推进流体质量不断减少。这对于烟雾等气态流体的仿真来说并无大碍，因为烟雾本身就具有飘散后逐渐消失的现象，但对于水流仿真来说，数值耗散将导致水流细节被平滑掉，在视觉上产生失真。拉格朗日方法是一种基于粒子的方法。与欧拉方法相比，拉格朗日方法研究的是流体粒子或微团的运动，因此，流体的对流项可利用流体质点的运动来直接模拟，不需要对整个流场进行网格化求解，因而可避免欧拉方法中的数值耗散问题。

光滑粒子流体动力学(Smoothed Particle Hydrodynamics，SPH)是流体动力学领域近20年来发展起来的一种新的纯拉格朗日方法，最初主要用在天体物理的研究中，1994年，Monaghan将其用于模拟流体的自由表面流动，引起很多学者的关注，随后，流体动力学领域纷纷展开利用SPH方法模拟自由表面流动的研究。1995年，Stam首次将SPH方法引入到计算机图形学领域，模拟了火和其他气态现象。2003年，Muller基于SPH方法模拟了向杯中倒水的过程，他从拉格朗日流体控制方程出发，将水粒子所受到的力分解为压力、重力、粘滞力，并根据SPH理论和公式求解这些力，水粒子在这些力的作用下遵循牛顿第二定律运动。同时，Muller还提到水粒子与障碍物的碰撞问题和基于marching cubes算法重构水体表面，但没有给出具体的解决方案。Muller的方法为基于SPH的水流模拟研究奠定了基础，但是，现有的基于SPH的水流模拟方法只考虑了物体处于静止状态这一种情况，将物体看作盛水的容器或是限制水的运动范围的障碍物，没有让物体运动起来并参与到水的运动过程当中，因此，不能表现运动物体与水发生交互的效果。

发明内容

为了模拟运动物体与水的交互现象，本发明提出一种基于SPH的运动物体与水交互的仿真建模方法。该方法对物体和水进行粒子化，建立物体粒子与水粒子的受力分析模型和状态更新方法，以及粒子与障碍物的碰撞检测方法，并基于Marching Cubes算法从水粒子的数据中抽取出水的自由运动表面。

本发明提出的基于SPH的运动物体与水交互的仿真建模方法，包括如下几个步骤：

一、对仿真进行初始化，包括粒子属性的初始配置和建立粒子的相邻粒子列表。首先将物体和水进行粒子化，即将物体和水用粒子的形式来表达，建立包含物体粒子和水粒子的粒子集合，同时，对集合中的粒子属性进行初始配置：

1)确定物体粒子和水粒子的数量；

2)建立粒子集合的位置数组和速度数组；

3)初始时按照物体和水的形状等间隔分布粒子，并将每个物体粒子和水粒子的初始位置保存到位置数组中，同时将每个物体粒子和水粒子的初始速度设置为零；

4)设置每个粒子的质量、密度、压强等常数属性，同时为所有粒子设定统一的粒子支持域半径，即光滑核半径；

5)仿真的初始时间设置为t＝0，并定义时间步常量TIME_STEP。

对粒子集合进行初始配置后，为每个粒子建立粒子的相邻粒子列表，即为每个粒子寻找位于其支持域内的所有粒子，列表中包含位于粒子支持域内的所有粒子的索引值。

二、仿真时间步向前推进一步，对水粒子的相邻粒子列表进行更新。

三、计算水粒子在下一步时间步的位置和运动状态。

四、计算物体粒子在下一步时间步的位置和运动状态，具体步骤为：

1)计算物体粒子j所受浮力f_j ^float

2)计算物体粒子j的重力f_j ^gravity

3)计算物体粒子所受合力F_j

4)计算F_j作用下的平移加速度a_j

5)计算物体粒子的平移速度vel_j+1

6)计算物体粒子的旋转速度vel′_j+1

7)计算物体粒子在下一刻的速度v_j+1＝vel′_j+1+vel_j+1

8)计算物体粒子在下一刻的位置x_j+1＝x_j+v_j+1Δt

五、根据计算得到的粒子位置判断粒子与障碍物是否发生碰撞，如果发生碰撞，则改变粒子的运动速率和方向，否则直接执行第六步。

六、基于marching cubes算法构建水面拓朴结构，具体步骤为：

1)由非规则分布的粒子数据产生出规则分布的水流密度场数据，即生成规则分布的标量(密度)数据场，具体做法是：将网格节点作为流场中的规则采样点，节点的采样值通过SPH插值来计算；

2)确定水面的密度值，作为等值面的阈值；

3)基于Marching Cubes算法进行水面抽取，即从水体密度场中抽取出水体表面的密度等值面，进而构建出水体表面的拓朴结构，完成水面绘制。

七、采用图形绘制语言，如OPENGL或DIRECTX，将第六步所生成的三角形网格传输到显卡进行水面的绘制和显示。同时，以物体的重心位置为基准，以三角面片的形式绘制出物体。

八、判断是否有人为关闭仿真进程的行为或仿真结束时间到，如果有，则仿真停止运行；否则仿真继续，重复执行二～七步。

为提高效率，建立粒子的相邻粒子列表可以通过对粒子所在空间进行网格化处理来加速粒子搜索。因此，建立相邻粒子列表涉及粒子空间网格化和基于网格的支持域粒子搜索两个方面，具体步骤为：

1)粒子空间网格化

首先对粒子集中的所有粒子进行一次逐一检查，找到这些粒子在X，Y，Z三个方向上的位置坐标的最大值和最小值，通过最大值和最小值之差求出粒子集在三个方向上的分布跨度；然后对每个方向上的分布跨度以光滑核半径的长度为间距进行等间隔划分，最终建立一个以光滑核半径为网格单元长度的三维空间网格。

2)基于网格的支持域粒子搜索

对每个粒子建立粒子的相邻粒子列表，即为每个粒子寻找位于其支持域内的所有粒子，列表中包含位于粒子支持域内的所有粒子的索引值，搜索方法是分别对当前粒子所在网格单元及与该单元上、下、左、右、前、后相邻的单元(共27个)内的粒子进行搜索，并判断搜索到的粒子与当前粒子的距离是否小于光滑半径，从而决定是否将其加入相邻粒子列表。

使用本发明所述方法可以模拟运动物体与水的交互现象，获得比较真实的仿真结果，且当粒子数量低于6000个时，可达到实时仿真的效果。该方法可应用于计算机游戏、军事仿真、影视动画等领域。

附图说明

图1是基于SPH的运动物体与水交互的仿真总体流程图。

图2是粒子集合的数据结构示意图。

图3是空间网格和粒子分配示意图。

图4是网格、粒子、相邻粒子列表的数据结构及关系示意图。

图5是相邻粒子搜索流程图。

图6是粒子与障碍物发生碰撞反弹示意图。

具体实施方式

本发明提出的基于SPH的运动物体与水交互的仿真方法，其总体流程如图1所示。

一、对仿真进行初始化，包括建立粒子集合并对粒子属性进行初始配置、建立每个粒子的相邻粒子列表。

首先，是建立粒子集合并对粒子属性的初始配置。由于SPH是一种基于粒子的方法，因此，在基于SPH进行水与物体的交互仿真时必须首先将物体和水进行粒子化，即将物体和水用粒子的形式来表达，建立包含物体粒子和水粒子的粒子集合，该集合用数组结构来表达，同时，对集合中的粒子属性进行初始配置。具体涉及如下几方面内容：

1)确定物体粒子和水粒子的数量，对于规则的物体或水体来说，是确定X、Y、Z三个方向上每个方向的粒子分布数，然后将三个方向的数量作乘积。粒子分布数可以任意选择，数量增加可以使效果逼真，但会导致耗费资源增加，数量减小可以加快处理速度，但会导致仿真效果不够真实，具体数量可以根据使用资源与期望获得的效果确定。例如，假设物体和水的形状是立方体，三个方向上的粒子数量分别为sx、sy、sz，则物体粒子的数量为sx×sy×sz。对于不规则的物体或水体来说，只能用手工方法确定每个粒子的位置。将物体粒子和水粒子的数量设置为常量，分别记为SOLID-NUM和WATER-NUM，同时，确定粒子集合的总数量，记为N-PARTICLES，且N-PARTICLES＝SOLID-NUM+WATER-NUM。

2)定义两个长度为N_PARTICLES的数组，称为粒子集合的位置数组和速度数组。数组0～SOLID-NUM-1元素用于保存物体粒子的信息，SOLID-NUM～N-PARTICLES元素用于保存水粒子的信息。同时，这些粒子按坐标排列以线性顺序存放到数组中。图2是粒子集合的数据结构示意图。

3)初始时的物体和水处于静止状态，可按照物体和水的形状等间隔分布粒子，粒子间的间隔可通过粒子分布长度除以该方向上的粒子分布数量来确定，并将每个物体粒子和水粒子的初始位置保存到位置数组中。同时，将每个物体粒子和水粒子的初始速度设置为零。

4)每个粒子都具有质量、密度等常数属性，同时所有粒子都具有统一的光滑核半径，这些属性也需要在仿真开始时给定。每个粒子的质量设为m＝0.0002克，这是领域内的相关文献指出的效果较好的经验值；水粒子的密度设为ρ_水＝1000克/立方米；木头粒子的密度设为ρ_木头＝300克/立方米；粒子的光滑核半径决定着影响粒子状态的粒子支持域的范围，在本实施例中设为h＝0.005米。

5)仿真的初始时间设置为t＝0，并定义时间步常量为TIME_STEP＝0.005秒。

其次，是建立粒子的相邻粒子列表。由于每个粒子的属性都是通过粒子支持域内的其它粒子的插值来获得的，为便于管理支持域内的粒子，我们为每个粒子建立一个相邻粒子列表，列表中包含位于粒子支持域内的所有粒子的索引值。例如，如果粒子1的支持域内包含粒子3、粒子6和粒子8，则粒子1的相邻粒子列表中就会含有3、6、8三个粒子索引值。

建立粒子的相邻粒子列表的过程就是为每个粒子寻找位于其支持域内的所有粒子的过程。在不进行特殊处理的情况下，建立一个粒子的相邻粒子列表需要对粒子集中的所有粒子进行一次从头到尾的搜索，而建立n个粒子的相邻粒子列表就需要进行n次粒子集搜索，这无疑会耗费大量机时。为此，我们提出对粒子所在空间进行网格化处理，以此来加速粒子搜索。因此，建立相邻粒子列表涉及粒子空间网格化和基于网格的支持域粒子搜索两个方面。

1)粒子空间网格化：

首先对粒子集中的所有粒子进行一次逐一检查，找到这些粒子在X，Y，Z三个方向上的位置坐标的最大值和最小值，通过最大值和最小值之差求出粒子集在三个方向上的分布跨度；然后对每个方向上的分布跨度以光滑核半径的长度为间距进行等间隔划分，最终可建立一个以光滑核半径为网格单元长度的三维空间网格，如图3所示。网格单元按坐标顺序线性排列，每个粒子的X、Y、Z坐标分别与网格最左下角点的X、Y、Z坐标相减并除以网格单元长度，得到粒子在的网格中的单元位置。同时，每个粒子都有一个指向其相邻粒子列表的指针。图4是网格、粒子、及相邻粒子列表的数据结构与关系示意图。

2)支持域粒子搜索：

由于在对粒子运动空间进行网格划分时，网格单元的长度被设置为核函数的光滑半径长度，这使得与粒子相隔较远的网格单元因超出粒子的支持域范围而不必被搜索。因此，我们对粒子支持域内的粒子搜索限定在与粒子相邻的网格单元中进行即可。搜索过程针对粒子集中的所有粒子，搜索方法是分别对当前粒子所在网格单元及与该单元上、下、左、右、前、后相邻的单元(共27个)内的粒子进行搜索，这样做可减少参与搜索的范围和粒子数量，提高搜索效率。具体步骤如下：

1)从粒子集中取出粒子，记为i。注意，如果只为物体粒子建立相邻粒子列表，则只从粒子集中的物体粒子部分取粒子，同理，如果只为水粒子建立相邻粒子列表，则只从粒子集中的水粒子部分取粒子。

2)确定当前粒子所在的网格单元

3)确定与该单元相邻的网格单元(共27个)

4)取出一个相邻网格单元

5)取出该网格单元格内的粒子，记为j

6)计算粒子i与粒子j之间的距离，如果距离小于光滑半径，则执行第7步，否则执行第9步

7)将粒子j记入粒子i的相邻粒子列表中

8)粒子i的支持域内的粒子计数加“1”

9)判断网格单元内的粒子是否全部访问完毕，如果全部访问完毕，则执行第10步，否则执行第5步

10)判断待搜索的网格单元是否全部访问完毕，如果全部访问完毕，则执行第11步，否则执行第4步

11)判断粒子集中的粒子是否全部访问完毕，如果全部访问完毕，则执行第12步，否则执行第1步

12)搜索结束

整个搜索流程如图5所示。

需要说明的是，仿真初始阶段只建立物体粒子的相邻粒子列表。这是由于物体是一种固体，而固体粒子间的相对位置始终保持不变，因此，物体粒子的相邻粒子列表只在仿真初始阶段求取一次即可。而水粒子间的相对位置会由于水的运动而不断发生变化，因此，水粒子的相邻粒子列表需要在每一个时间步的开始时刷新，因而，不必在仿真初始阶段建立，可放在仿真运行阶段完成。

二、空间网格、水粒子相邻粒子列表的更新

此处开始进入仿真运行阶段，仿真时间步向前推进一步，即t＝t+TIME_STEP。由于水粒子间的相对位置会随着粒子水的运动而不断变化，因此，每个水粒子支持域内的粒子是不断变化的，水粒子的相邻粒子列表需要在每个时间步都刷新一次，即每个时间步都要对所有水粒子的相邻粒子列表重新建立一次。建立方法与步骤(1)中的相同，即包括粒子空间网格化和基于网格的水粒子支持域搜索两个方面。需要说明的是，虽然在初始化阶段已建立空间网格，但由于网格空间在X、Y、Z三个方向上的范围刚好等于粒子在X、Y、Z三个方向上的跨度，而粒子的不断运动导致粒子在空间中的分布跨度随着粒子的运动不断变化，因此，要对粒子空间网格进行重新定位和分割，即对网格进行刷新。空间网格的重新划分方法与具体实施方案“步骤一”中的网格划分方法相同。网格被刷新后，可根据新计算出来的粒子空间位置定位每个粒子所处的网格单元。

由于初始化阶段已经建立了物体粒子的相邻粒子列表，且列表内容在每个时间步都不发生变化，因此不必更新。

有了正确的物体粒子和水粒子的相邻粒子表后，就可基于这些列表进行物体粒子和水粒子的位置和运动状态。首先计算水粒子的位置和运动状态，然后计算物体粒子的位置和运动状态。

三、计算水粒子在下一刻的位置和运动状态。

水粒子的运动状态可采用Muller提出的方法进行计算。这里只给出概要的描述，并归纳出计算步骤。

每一个水粒子i所受的力分解可为压力、重力和粘滞力，这些力的计算公式如下：

$f_i=f_i^{\mathrm{pressure}}+f_i^{\mathrm{vis}\cos \mathrm{ity}}+f_i^{\mathrm{gravit}}$

$f_i^{\mathrm{pressure}}=-\▿p_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{m}_j\frac{(p_i+p_j)}{2{\mathrm{\ρ}}_j}\▿W(x_i-x_j,h)$

(1)

$f_i^{\mathrm{vis}\cos \mathrm{ity}}=\mathrm{\μ}{\▿}^2{u}_i=-\mathrm{\μ}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{m}_j\frac{(u_j-u_i)}{{\mathrm{\ρ}}_j}{\▿}^{2}W(x_i-x_j,h)$

$f_i^{\mathrm{gravit}}={\mathrm{\ρ}}_{i}g=g\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{m}_{j}W(x_i-x_j,h)$

其中，f_i ^pressure表示水粒子i所受到的压力，f_i ^viscosity表示水粒子i所受到的粘滞力，f_i ^gravity表示水粒子i所受到的重力，ρ_i表示密度，m_i表示质量，p_i表示压强，W(x_i-x_j，h)表示光滑核函数或核函数，h代表光滑核半径。

水粒子在力f_i的作用下不断发生位置和速度的变化。为完成粒子的受力计算，必须知道压强场的分布，即粒子i所在位置的压强。Muller采用理想气态方程的修改形式进行水粒子的压强计算：

p_i＝k(ρ_i-ρ₀)(2)

其中，k是气态常数，在水流模拟中可设置为声音在水中的传播速度，ρ_i是粒子i的当前密度，ρ₀是粒子i的初始密度。

水粒子在力的作用下完成位置和速度的更新，具体步骤如下：

1)计算水粒子的密度ρ_i

2)计算水粒子所在位置的压强p_i

3)计算水粒子的重力f_i ^gravity

4)计算水粒子的压力f_i ^pressure

5)计算水粒子的粘滞力f_i ^viscosity

6)计算水粒子所受到的合力f_i

7)计算水粒子加速度a_i

8)计算水粒子在下一刻的速度：v_i+1＝v_i+a_iΔt

9)计算水粒子在下一刻的位置：x_i+1＝x_i+v_i+1Δt

四、计算物体粒子在下一刻的位置和运动状态。

物体所受到的力包括重力和水的浮力，其计算公式如下：

$F_j=f_j^{\mathrm{float}}+f_j^{\mathrm{gravity}}---\left(3\right)$

其中，F_j表示物体粒子所受到的合力，f_j ^float表示物体粒子所受到的浮力，f_j ^gravity表示物体粒子所受到的重力。

重力的计算公式如下：

$f_j^{\mathrm{gravity}}={\mathrm{\ρ}}_{\mathrm{solid}}\·g---\left(4\right)$

其中，ρ_soild表示物体粒子的密度，g表示重力加速度。ρ_soild和g均是恒定不变的常数，因此，物体粒子所受的重力f_j ^gravity是一个常量。

物体的浮力来自于水粒子的影响。但是，并不是所有物体粒子都会因水粒子的作用而产生浮力，只有那些处于水粒子支持域内的物体粒子才会受到浮力作用，对于那些不在水粒子支持域内的物体粒子来说，它们只受到重力作用。

根据机械能守恒，当系统不受外力作用时，系统的合力应该为零。对于两相邻粒子i、j来说应满足：

f_i+f_j＝0   (5)

因此，求取每个物体粒子所受浮力时，可先找到所有受该物体粒子影响的水粒子，计算这些水粒子与固体粒子间的作用力，并求和取反，所得的值即为物体粒子所受浮力。即：

$f_j^{\mathrm{float}}=-\stackrel{N}{\mathrm{\Σ}}{f_{\mathrm{ij}}}^{\mathrm{pressure}}+f_{\mathrm{ij}}^{\mathrm{vis}\cos \mathrm{ity}}+f_{\mathrm{ij}}^{\mathrm{gravity}}---\left(6\right)$

其中，f_j ^float表示物体粒子i_j的浮力；f_ij ^pressure表示水粒子i受物体粒子j影响所产生的压力，它是产生于物体粒子j和水粒子i之间的一个作用力；同理，f_ij ^viscosity表示水粒子i受物体粒子j影响所产生的粘滞力；f_ij ^gravity表示水粒子i受物体粒子j影响所产生的重力；N表示受物体粒子j影响的水粒子的个数。

公式(6)的物体粒子所受浮力的累加计算可在水粒子状态更新的第4步和第5步中完成。在求得固体粒子的合力后，就可以更新固体粒子的加速度和速度，但此时的速度是平移线速度，记为v_j+1。平移线速度并不是可以更新固体粒子位置的最终速度，因为物体在受到水的作用后还会产生绕自身重心自转的运动，因此，每个物体粒子还有一个旋转线速度。物体粒子旋转线速度的求取步骤如下：

1)计算物体的重心，计算方法是将所有物体粒子X、Y、Z三个方向的坐标分别累加并求平均

2)计算每个物体粒子与物体重心的距离，记为r_j

3)将物体粒子的速度vel_j+1与r_j求叉乘，得到粒子的旋转角速度w_j

4)计算所有物体粒子的平均角速度，得到物体的转动角速度W

5)对每个物体粒子，作r_j与W的叉乘运算，得到粒子的旋转线速度vel′_j+1

6)在得到物体粒子的旋转线速度后，将粒子的旋转线速度与平移速度相加，得到物体粒子的最终速度v_j+1：

v_j+1＝vel′_j+1+vel_j+1  (7)

物体粒子在力的作用下完成位置和速度的更新，具体步骤如下：

1)计算物体粒子j所受浮力f_j ^float

2)计算物体粒子j的重力f_j ^gravity

3)计算物体粒子所受合力F_j

4)计算F_j作用下的平移加速度a_j

5)计算物体粒子的平移速度vel_j+1

6)计算物体粒子的旋转速度vel′_j+1

7)计算物体粒子在下一刻的速度v_j+1＝vel′_j+1+vel_j+1

8)计算物体粒子在下一刻的位置x_j+1＝x_j+v_j+1Δt

五、粒子与障碍物的碰撞检测。

为保证仿真的真实性，必须考虑水及物体与障碍物发生碰撞的情况。由于我们的方法是将水与物体统一化作粒子来处理，只是粒子的种类不同而已，因此，无论是水的运动还是物体的运动，最终都归结为粒子的运动，这使得水和物体与障碍物的碰撞检测最终可简化为粒子与障碍物间的碰撞检测。这里，将障碍物的边界用多边形集合来描述，每一个多边形代表障碍物上的一个面。因此，水流及物体与障碍物的碰撞检测就可转化为粒子与多边形的碰撞检测。

前面我们给出了水粒子和物体粒子的位置状态更新步骤，但该步骤没有考虑粒子与周围环境发生碰撞的情况。实际上，粒子的碰撞检测应发生在粒子位置更新之前，以水粒子为例，应发生在步骤8与步骤9之间。步骤8计算得到的速度是粒子在下一时刻的速度，通过该速度我们可以预计出粒子在下一时刻所能达到的位置，如果该位置与当前时刻粒子位置的连线穿过障碍物表面，则说明如果按照当前计算出的下一时刻的速度来位移粒子，粒子就会与障碍物发生碰撞并进入障碍物体内，这显然不是我们所期望的，因此，必须改变粒子在下一时刻的运动速度。

将粒子和障碍物看作理想刚体，则粒子与障碍物发生碰撞后应做反弹运动，如图6所示。设N为障碍物表面的法向量，A点表示粒子的当前位置，B点表示粒子按照步骤8计算出的下一时刻的速度v′_j+1运动一个时间步后所能到达的位置。当A、B两点间的虚线穿透障碍物表面时，若不改变粒子在下一时刻的运动速度，则粒子将与障碍物发生碰撞并进入障碍物体内。此时，让粒子在碰撞点C处(AB与障碍物表面的交点)发生反弹运动，运动方向与N之间的夹角等于AB与N之间的夹角。设E点是B点相对于障碍物表面垂直对称的点，在理想情况下，粒子被反弹后最终应运动到E点，但考虑到实际的碰撞损耗，粒子反弹后只能运动到D点，速度记为v_i+1：

v_i+1＝k_damp(v′_i+1×2N(N·v′_i+1))  (8)

其中，k_damp表示碰撞后的速度衰减系数。

粒子的加速度与粒子的受力有关，碰撞作用对粒子所产生的作用力与N同向，该作用力产生的加速度分量可表示为a_N＝αN，其中α表示从碰撞作用力到加速度的转化系数，则碰撞后粒子的加速度为a_j+1＝a_j+a_N。

六、基于marching cubes算法的水面重构。

Marching Cubes算法是一种经典的体等值面抽取算法。基于Marching Cubes算法的水面抽取实质是从水体密度场中抽取出水体表面的密度等值面，进而构建出水体表面的拓朴结构，完成水面绘制。

由于Marching Cubes算法要求体数据是在三维空间场呈规则分布的离散数据，而采用SPH方法成的粒子运动数据是非规则分布的，因此，必须首先由非规则分布的粒子数据产生出规则分布的水流密度场数据，即生成规则分布的标量(密度)数据场。具体做法是：将网格节点作为流场中的规则采样点，节点的采样值通过SPH插值来计算。首先，找到以该节点为顶点的所有体元，对于网格体内部的节点来说，每个节点对应8个体元，位于网格体边界的节点对应4个体元，位于网格体4个顶角的节点对应1个体元。其次，对每个节点对应的体元内的粒子进行SPH插值，求出节点处的密度：

${\mathrm{\ρ}}_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{m}_{j}W(x_i-x_j,h)---\left(9\right)$

其中，j表示与节点i相对应的所有网格单元内的粒子。ρ_i表示密度，m_i表示质量，W(x_i-x_j，h)表示光滑核函数或核函数，h代表光滑核半径。

为采用Marching Cubes算法从水密度场中抽取出水体表面的密度等值面，还必须确定水面的密度值，作为等值面的阈值。根据水面压强等于大气压强，并根据公式(2)可得到水面的密度如下：

其中，ρ_surface表示水面密度，p_气表示大气压强，k是气态常量。

关于基于Marching Cubes算法的标量数据场等值面抽取，可参阅相关文献(如唐圣泽出版的《科学计算可视化》一书，以及Marching Cubes算法最初的创始人William E.Lorensen所发表的文献)。

七、采用图形绘制语言，如OPENGL或DIRECTX，将第6步所生成的三角形网格传输到显卡进行水面的绘制和显示。同时，以物体的重心位置为基准，以三角面片的形式绘制出物体。

前一步已经确定了表示水面的三角形网格在空间中的位置，这一步是将该三角形网格绘制到屏幕上。绘制语言可采用OPENGL或DIRECTX。绘制方法是将每一个三角形传输到显卡上进行绘制。此外，具体实施方案“步骤四”中已经计算得到物体粒子的重心，本步骤可以物体的重心位置为基准，用三角面片的形式绘制出物体。

采用三角面片的图形绘制方法这是图形学中常用的绘制方法，这里不再赘述。由于每一帧都进行绘制，因此，最终可以看到的是物体和水面连续运动的动态效果。

Claims (2)

1.一种基于SPH的运动物体与水交互的仿真建模方法，包括以下步骤：

一、对仿真进行初始化，包括粒子属性的初始配置和建立粒子的相邻粒子列表，首先将物体和水进行粒子化，即将物体和水用粒子的形式来表达，建立包含物体粒子和水粒子的粒子集合，同时对集合中的粒子属性进行初始配置：

1)确定物体粒子和水粒子的数量；

2)建立粒子集合的位置数组和速度数组；

3)初始时按照物体和水的形状等间隔分布粒子，并将每个物体粒子和水粒子的初始位置保存到位置数组中，同时将每个物体粒子和水粒子的初始速度设置为零；

4)设置每个粒子的质量、密度、压强常数属性，同时为所有粒子设定统一的粒子支持域半径，即光滑核半径；

5)仿真的初始时间设置为t＝0，并定义时间步常量TIME_STEP；

对粒子集合进行初始配置后，为每个粒子建立粒子的相邻粒子列表，即为每个粒子寻找位于其支持域内的所有粒子，列表中包含位于粒子支持域内的所有粒子的索引值；建立粒子的相邻粒子列表通过对粒子所在空间进行网格化处理来加速粒子搜索，具体步骤为：

1)粒子空间网格化

首先对粒子集中的所有粒子进行一次逐一检查，找到这些粒子在X，Y，Z三个方向上的位置坐标的最大值和最小值，通过最大值和最小值之差求出粒子集在三个方向上的分布跨度；然后对每个方向上的分布跨度以光滑核半径的长度为间距进行等间隔划分，最终建立一个以光滑核半径为网格单元长度的三维空间网格；

2)基于网格的支持域粒子搜索

对每个粒子建立粒子的相邻粒子列表，即为每个粒子寻找位于其支持域内的所有粒子，列表中包含位于粒子支持域内的所有粒子的索引值，搜索方法是分别对当前粒子所在网格单元及与该单元上、下、左、右、前、后相邻的单元内的粒子进行搜索，并判断搜索到的粒子与当前粒子的距离是否小于光滑半径，从而决定是否将其加入相邻粒子列表；

二、仿真时间步向前推进一步，对水粒子的相邻粒子列表进行更新；

三、计算水粒子在下一步时间步的位置和运动状态；

四、计算物体粒子在下一步时间步的位置和运动状态，具体步骤为：

1)计算物体粒子j所受浮力f_j ^float

2)计算物体粒子j的重力f_j ^gravity

3)计算物体粒子所受合力F_j

4)计算Fj作用下的平移加速度a_j

5)计算物体粒子的平移速度vel_j+1

6)计算物体粒子的旋转速度vel′_j+1

7)计算物体粒子在下一刻的速度v_j+1＝vel′_j+1+vel_j+1

8)计算物体粒子在下一刻的位置x_j+1＝x_j+v_j+1Δt

五、根据计算得到的粒子位置判断粒子与障碍物是否发生碰撞，如果发生碰撞，则改变粒子的运动速率和方向，否则直接执行第六步；

六、基于marching cubes算法构建水面拓朴结构，具体步骤为：

1)由非规则分布的粒子数据产生出规则分布的水流密度场数据，即生成规则分布的标量密度数据场，具体做法是：将网格节点作为流场中的规则采样点，节点的采样值通过SPH插值来计算；

2)确定水面的密度值，作为等值面的阈值；

3)基于Marching Cubes算法进行水面抽取，即从水体密度场中抽取出水体表面的密度等值面，进而构建出水体表面的拓朴结构，完成水面绘制；

七、采用图形绘制语言将第六步所生成的三角形网格传输到显卡进行水面的绘制和显示；同时，以物体的重心位置为基准，以三角面片的形式绘制出物体；

八、判断是否有人为关闭仿真进程的行为或仿真结束时间到，如果有，则仿真停止运行；否则仿真继续，重复执行二～七步。

2.根据权利要求1所述的仿真建模方法，其特征在于，在对仿真进行初始化的步骤中，只建立物体粒子的相邻粒子列表，而将水粒子的相邻粒子列表放在仿真运行阶段建立和更新。

Images