一种基于IISPH提高不可压缩水体模拟效率的实现方法与流程

本发明涉及计算机图形学领域，具体涉及一种基于iisph(implicitincompressibilitysmoothparticlehydrodynamics，隐式不可压缩光滑粒子流体动力学)提高不可压缩水体模拟效率的实现方法。

背景技术：

受限于计算机的计算性能与计算效率，水体模拟技术长期处于理论研究阶段。近年来随着计算机硬件性能以及计算效率的提升，水体模拟得以实现，同时计算机图形学的发展进一步促进了基于物理的水体模拟技术的发展，使其逐渐成为备受关注的研究方向之一。光滑粒子流体动力学(smoothparticlehydrodynamics，sph)方法是一种通过求解质点组动力学方程计算流体(包括水体)整体运动状态的无网格方法，因其存在能够保持质量守恒、易于捕捉水花和气泡等诸多优点，被广泛应用于电影、游戏、广告以及军事等各个领域。

传统sph方法使用理想气体的状态方程求解水体粒子受到的压力，使得水体具有明显的可压缩性，从而导致在模拟大规模水体以及大尺度运动时，水体表面出现明显的不自然震荡，造成明显的视觉失真，因此实现水体的不可压缩性对于提高水体模拟真实感具有重要意义。

现有的实现水体不可压缩性的思路主要分为两种：显式不可压缩方法以及隐式不可压缩方法。显式不可压缩方法通过直接求解液体状态方程或者压力泊松方程(pressurepoissonequation，ppe)求解压力，而隐式不可压缩方程多采用预测校正的方式求解压力。

1.显式不可压缩方法

显式不可压缩方法主要通过直接求解液体状态方程或者压力泊松方程，得到水体粒子受到的压力，从而实现水体的不可压缩性。

2007年，becker等人提出了微可压缩sph(weaklycompressiblesph，wcsph)方法(参考文件1：monaghan,j.j.simulatingfreesurfaceflowswithsph[j].journalofcomputationalphysics,1994,110(2):399-406)，该方法使用液体状态方程代替理想气体的状态方程，并且通过使用大的刚性系数使得水体粒子的密度波动量不超过1％，从而实现了水体的不可压缩性。wcsph保证了水体模拟的不可压缩性，显著提升了水体模拟的真实感。但是水体的刚度支配着cfl(courant-friedrichs-levy)条件，大的刚性系数只允许使用小的时间步长，其计算代价随着压缩性的降低而增加，因此wcsph方法无法高效模拟复杂水体运动。

同年，adams等人提出了不可压缩sph(incompressiblesph，isph)方法(参考文件2：x.y.hu,n.a.adams.anincompressiblemulti-phasesphmethod[j].journalofcomputationalphysics,2007,227(1):264–278.)。该方法通过共轭梯度直接求解压力泊松方程，从而得到粒子受到的压力。isph方法首先计算压力以外的其他力作为合力，其次通过该合力计算粒子的中间速度，接着对泊松方程进行离散，计算出粒子受到的压力。不同于wcsph，isph允许使用较大的时间步长，但是该方法在每个时间步长内需要对复杂公式以及方程组进行求解，计算开销巨大，无法实现不可压缩水体的高效模拟。

2015年，kang等人将无散度条件引入了isph方法(参考文件3：kangn,sagongd.incompressiblesphusingthedivergence-freecondition[j].computergraphicsforum,2015,33(7):219-228)，该方法在求解动量方程时，实施无散度条件。保证体积不可压缩的同时保持一个无散度速度场，在位置级和速度级上同时执行不可压缩性。相比于isph，该方法能够极大地加快压力求解的收敛速度，减少求解器的迭代次数，从而提高水体模拟效率。但是该方法在模拟过程中，容易出现水体自由表面处粒子密度被低估，从而产生粒子聚集等现象。

2.隐式不可压缩方法

隐式不可压缩方法主要通过间接求解压力泊松方程，得到水体粒子受到的压力，从而实现水体的不可压缩性。

2009年，solenthaler等人提出基于预测-校正的sph(predictive-correctiveincompressiblesph，pcisph)方法(参考文件4：solenthalerb,pajarolar.predictive-correctiveincompressiblesph[j].acmtransactionsongraphics,2009,28(3):1-6)，该方法分为两个阶段：预测阶段与校正阶段。在预测阶段，该方法使用压力外的其他力作为合力，并且根据该合力计算下一时间步长粒子的速度、位置等属性，即粒子的预测物理属性。然后执行迭代校正循环，当粒子密度的误差大于给定的误差值时，通过给粒子施加压力调整粒子的速度，若平均密度误差在给定的误差范围内则退出校正循环。最后利用求解出的压力以及其他力重新计算下一时间步长的粒子属性。该方法避免了直接求解压力泊松方程的计算开销，但间接求解压力需要经历多次迭代，在一定程度上仍然增加了计算开销。

2012年，he等人提出了一种求解不可压缩性水体的局部泊松sph(localpoissonsph，lpsph)方法(参考文件5：x.he,n.liu,s.li,etal.localpoissonsphforviscousincompressiblefluids[j].computergraphicsforum,2012,31(6):1948–1958)。该方法首先将泊松方程转化为积分形式；然后应用离散化方法将连续积分方程转化为离散求和，避免了直接求解全局压力泊松方程的计算开销。

2014年，ihmsen等人提出了一种压力投影公式，即隐式不可压缩sph(implicitincompressibilitysph，iisph)方法(参考文件6：ihmsenm,cornelisj,solenthalerb,etal.implicitincompressiblesph[j].ieeetransactionsonvisualizationandcomputergraphics,2014,20(3):426-435.)，该方法将连续性方程的对称sph压力与sph离散相结合，得到压力泊松方程的离散形式，并且采用松弛jacobi算法求解压强。与之前的投影方案相比，该方法提高了求解器的收敛速度，此外，该方法基于速度而不是位置来校正密度偏差，这在一定程度上提高了时间积分的稳定性。但是该方法计算得到的速度场存在散度误差，从而增加了密度校正的平均迭代次数，降低了模拟效率。

2014年，ihmsen等人提出iisph-lfip(implicitincompressibilitysph-fluidimplicitparticle，iisph-flip)方法(参考文件7：cornelisj,ihmsenm,peera,etal.iisph-flipforincompressiblefluids[j].computergraphicsforum,2014,33(2):255-262.)。该方法采用隐式不可压缩方法进行压力投影，并采用流体隐式粒子(fluidimplicitparticle，flip)求解器进行边界处理，实现不可压缩水体的模拟。这种新的组合解决了现有sph和flip求解器的两个问题，即flip中的质量守恒与sph中的效率问题。该方法提出的iisph-flip求解器能够模拟不可压缩水体，其可量化、不可感知的密度偏差小于0.1％。

显式不可压缩水体模拟方法直接求解液体状态方程或者压力泊松方程，计算复杂度高；隐式不可压缩水体模拟方法，如pcisph、iisph等方法采用预测校正的方式，但是求解过程中迭代次数过多，且每次迭代都需执行耗时巨大的邻近粒子查找步骤，导致计算量随着粒子数的增加急剧增加。综上，现有不可压缩水体模拟方法虽然实现了水体的不可压缩性，但在施加不可压缩性时开销太大，降低了模拟效率。

技术实现要素：

本发明提出基于iisph提高不可压缩水体模拟效率的实现方法，首先通过将速度无散度模型引入iisph方法，修正iisph方法中速度的散度误差，避免密度误差随着时间的增加不断增加，从而减少修正密度误差的平均迭代次数；然后对iisph的密度恒定模型进行改进，通过密度变化率求解密度变化量，从而使得改进的密度恒定模型能够使用速度无散度模型求解出的计算因子，进而避免了冗余计算，提高了不可压缩水体模拟效率。具体包括如下步骤：

步骤一：引入速度无散度模型，以修正iisph方法的散度误差

纳维-斯托克斯(navier-stokes，n-s)方程是用于描述粘性不可压缩流体动量守恒的运动方程。该方程主要用来描述液体和气体的运动，本质上依据牛顿第二定律，也可以表述为：无穷小体积流体的动量变化是其重力、粘力、压力以及其他作用力的效果总和。n-s方程描述的是自然界中的流体运动规律，其基于物理的本质使得它非常适合用于描述自然现象，例如水体、烟雾、气流等自然流体。

n-s方程的三维表现形式如式(1)和式(2)所示。

其中u代表粒子的速度，ρ代表粒子的密度，p表示粒子的压强，即施加在单位面积上的压力，f表示流体受到的外力(如重力等)。运算符·为空间向量的点乘，表示空间向量对空间维度的偏导，即又称为梯度，运算符为拉普拉斯算子。

式(1)描述的是流体运动时的动量守恒，其右式表示流体所有的力场量。等式右边的第一项代表粒子间压力差，第二项表示粒子之间的粘力，第三项为外力(如重力)。其物理意义为：水体粒子的运动是由压力、粘力及外力共同决定的。

式(2)称为不可压缩条件，该式表示水体在运动过程中的速度散度为0，即在某个时间步长内，水体单位体积内的流入以及流出量不发生变化，从而保证水体在运动过程中的不可压缩性。

iisph修正的粒子速度存在散度误差，即不满足本发明通过引入速度无散度模型，并使用该模型计算压力项，从而实现对粒子速度的散度误差进行修正，粒子i所受压力的大小由式(3)决定。

其中表示粒子i所受的压力，表示水体粒子的体积，表示压力梯度。引入表示对粒子i的速度v进行修正的刚性系数，可得压力梯度的计算式如式(4)所示。

其中mj表示邻近粒子j的质量，表示核函数的梯度，即为满足速度无散度条件的刚性系数。

粒子在运动的过程中满足动量守恒。对粒子i进行受力分析可知，粒子i所受的压力为邻近粒子j的压力之和，并且粒子间的压力是一对相互作用力。由于相互作用力大小相等，方向相反，因此，粒子i所受的压力与粒子i对邻近粒子j的压力满足

其中表示粒子i对粒子j的压力，这意味着压力作为一种粒子间的作用力，以内力的形式存在，并且所有粒子的压力之和为0。从水体的整体看来，该压力不产生水体本身动能的改变，满足动能守恒定律。粒子i对邻近粒子压力的大小与粒子间的距离有关，其求解公式如式(5)所示。

其中xj表示粒子j的三维坐标，求解器必须通过改变压力值修正粒子的速度，从而使得粒子的速度满足无散度性。密度相对于时间的变化率如式(6)所示。

其中wij表示光滑核函数，要保证水体的不可压缩性，即密度不发生变化，则式(6)求得的密度相对于时间的变化率应为0。式(6)中vi-vj表示粒子间的相对速度差，该速度差由压强差导致，不同粒子受到的压力存在差别，因此导致粒子间的速度不一样。由压力引起的速度变化量如式(7)所示。

其中δvi表示由压力引起的速度变化量，将式(7)代入式(6)，得到式(8)。

由式(3)可以计算出粒子所受压力的大小，结合式(6)密度变化率为0，即可求解出刚性系数；使用该刚性系数修正粒子速度，能够使得速度的散度为0，从而使得密度相对于时间的变化率为0。将式(4)代入式(3)，得再代入式(8)，得到式(9)。

根据式(9)求解出刚性系数如式(10)所示。

其中且αi的大小只与粒子位置有关。

式(10)为速度无散度模型，该模型求解出使得速度满足无散度性的刚性系数该刚性系数用于计算压力，保证水体在运动的过程中密度变化率不发生波动，从而实现水体的不可压缩性。

最后根据式(11)计算粒子受到的合力，并使用该合力重新计算粒子的速度。

其中表示粒子i所受合力，表示除压力外其他力的合力，表示粒子i受到的邻近粒子j的压力。

步骤二，改进密度恒定模型，使其与速度无散度模型共享计算因子，以修正iisph的密度误差。

在速度无散度模型中，使用满足速度无散度性的刚性系数修正粒子速度，从而保证密度变化率小于给定的散度误差值η^div。在进行无散度速度修正时需要计算α因子，本发明提出的密度恒定模型复用α计算因子，对水体粒子的密度进行进一步的修正。

本发明通过密度变化率以及时间步长计算密度变化量，则t+δt时刻的密度值ρ′i，如式(12)所示。

将速度无散度模型的密度变化率求解公式(8)代入式(12)，并结合密度恒定条件ρi＝ρ0，得到δt时间步长内的密度误差ρ′i-ρ0，如式(13)所示。

其中ρ0表示粒子初始密度，将式(3)、(4)代入式(13)。

由式(14)得到刚性系数的求解公式如式(15)所示。

其中为满足水体密度恒定的刚性系数，α为速度无散度模型中的计算因子，在同一个时间步长内，该计算因子仅计算一次，减少重复计算导致的时间开销。

式(15)为密度恒定模型，该模型求解出使得密度恒定的刚性系数该刚性系数用于计算粒子受到的压力，保证水体在运动的过程中密度变化量不发生波动，从而实现水体的不可压缩性。

由式(11)可知，密度恒定模型修正过后粒子所受的合力如式(16)所示。

使用矫正过后的粒子所受的合力对粒子的速度进行修正，结合牛顿第二定律可得修正后t+δt时刻粒子速度v′i的求解公式，如(17)所示。

根据上述计算得出的以及式(17)对粒子的速度vi进行修正，得到修正后的粒子速度v′i。

本发明将速度无散度模型式(11)与改进的密度恒定模型式(15)相结合，实现了对速度散度误差以及密度误差的修正，从而在减少密度修正迭代次数的同时，减少了冗余计算，能够有效提高不可压缩水体的模拟效率。测试结果表明，本发明提出的方法在模拟不可压缩水体时，模拟效率优于iisph以及pcisph方法。

附图说明

图1是基于iisph提高不可压缩水体模拟效率的实现方法流程图。

图2是改进的密度恒定模型算法图。

图3是速度无散度模型算法图。

图4是不可压缩水体系统模块设计图。

图5是pcisph、iisph以及本发明方法模拟效率对比图。

图6为本发明方法多线程效率对比图。

图7为三维场景水滴坠落模拟及效果图。

具体实施方式

下面将结合附图和实施例对本发明的技术方案及结果作进一步的详细说明。

本发明提出一种速度无散度的不可压缩水体模拟方法，首先通过将速度无散度模型引入iisph方法，修正iisph方法中速度的散度误差，避免密度误差随着时间的增加不断增加，从而减少修正密度误差的平均迭代次数；然后对iisph的密度恒定模型进行改进，通过密度变化率求解密度变化量，从而使得改进的密度恒定模型能够使用速度无散度模型求解出的计算因子，因此避免了冗余计算，可以显著提高不可压缩水体模拟效率。

步骤一：将速度无散度模型引入iisph，改进iisph密度恒定模型，构建不可压缩水体模拟框架，如图1所示。

步骤1.1，执行邻近粒子搜索，得到每个水体粒子在当前时间步长内的邻近粒子集合ni(t)。

步骤1.2，根据每个粒子的邻近粒子集合ni(t)，计算粒子受到的压力以外的合力，并根据该合力计算粒子的预测速度以及预测密度。

步骤1.3，结合密度恒定模型修正粒子速度，首先计算共享计算因子，并根据该计算因子计算刚性系数，该刚性系数能够使得密度误差小于给定的密度误差值；接着使用该刚性系数对粒子的速度进行校正，使得密度的变化量|ρ′i(t+δt)-ρ0|小于给定的密度误差值η^ρ。

步骤1.4，结合速度无散度模型修正步骤(1.3)得到的粒子速度，使用步骤(1.3)计算出的共享计算因子计算刚性系数，该刚性系数能够使得密度变化率误差小于给定的散度误差值η^div，然后使用该刚性系数校正粒子的速度，使得密度变化率小于给定的散度误差值η^div。

步骤1.5，遍历所有的粒子，根据上述计算得到粒子速度，将该速度更新为下一时间步长的粒子速度，并更新粒子的其他属性。

步骤二：结合改进的密度恒定模型，修正密度误差，如图2所示。

步骤2.1，遍历所有粒子，结合其邻近粒子集合ni(t)，计算平均密度变化量。

步骤2.2，若平均密度变化量小于给定的密度误差值(ρavg-ρ0<η^ρ)或者迭代次数iter不小于2，则退出密度恒定模型，否则执行步骤2.3。

步骤2.3，对于每个粒子，遍历其邻近粒子集合ni(t)，计算对应的刚性系数并根据该刚性系数计算该粒子以及邻近粒子所受的压力。

步骤2.4，根据步骤2.3计算得到的压力修正粒子速度，使用修正后的粒子速度，循环执行步骤2.1-2.4，直到密度变化量小于给定的密度误差(ρavg-ρ0<η^ρ)。

步骤三：结合速度无散度模型，修正散度误差，如图3所示。

步骤3.1，遍历所有粒子，结合其邻近粒子集合ni(t)，计算粒子的密度变化率。

步骤3.2，若密度变化率小于等于给定的散度误差值或迭代次数不小于1，则退出速度无散度模型，否则执行步骤3.3。

步骤3.3，对于每个粒子，遍历其邻近粒子集合ni(t)，计算对应的刚性系数并根据该刚性系数计算该粒子以及邻近粒子所受的压力。

步骤3.4，根据步骤3.3计算得到的压力修正粒子速度，使用修正后的粒子速度，循环执行步骤3.1-3.3，直到速度散度误差小于给定的散度误差值

步骤四：利用步骤二、三所求解的物理模型，构建不可压缩水体运动场景，输出结果并显示。

步骤二、三中求解的不可压缩水体模型，仅仅是基于数学物理求解得出的结果，从视觉上描述，是由一系列粒子构成，并不能逼真地呈现现实中的水体场景。因此还需要对上述模型进行相应的渲染操作，从而得到具有真实感的水体，本发明中采用屏幕空间渲染技术实现水体渲染，从而展示逼真的流体效果。

本发明所有测试都是在戴尔台式机上进行的，cpu型号是英特尔corei7-4790@3.6ghz，显卡型号amdradeon^tmr5240。

图4为不可压缩水体系统模块设计图，整个系统分为4个模块：初始化模块、预测校正模块、边界处理模块以及渲染模块。初始化模块主要接收用户设定的参数，生成粒子初始信息，包括粒子位置、压力、密度、速度等属性。

预测校正模块综合了速度无散度模型、密度恒定模型，该模块在每个时间步长内执行一次，生成下一时间步长的粒子速度，主要分为两个部分：预处理部分以及速度计算部分。其中预处理部分主要计算并存储只与粒子位置相关的数据，这部分数据在该时间步长内不会发生变化，通过预计算的方式可以避免重复计算。速度计算部分使用速度无散度模型以及密度恒定模型修正粒子的预测速度，使得校正过后的速度满足速度无散度性和密度恒定。

边界处理模块用于校正经预测校正模块校正之后的速度，从而避免粒子在运动的过程中穿透固体边界。首先检测粒子是否会与边界发生碰撞，若发生碰撞则通过基于排斥力的碰撞响应模型修正粒子速度，从而保证容器边界的不可穿透性；若不发生碰撞则无须执行碰撞响应模型。

渲染模块主要根据粒子位置信息实现水体渲染，包括各向异性核表面构建、深度平滑、施加光照以及绘制天空盒等步骤，提高了水体表面的平滑度以及水体渲染的真实感。

本发明将速度无散度模型引入iisph，并对iisph的密度恒定模型进行了改进，这种基于iisph提高不可压缩水体模拟效率的实现方法，可以更加高效地模拟不可压缩水体运动场景。

分别使用pcisph算法、iisph算法以及本发明方法实现对水滴坠落模型的模拟，该模型由17074个粒子组成，测试包含12次独立实验，通过计算得到12次结果的平均帧率加速比。表1为水滴坠落模型模拟帧率对比表，图5为对应的折线图。由表1中数据计算可得，pcisph算法的平均帧率为24.24fps，iisph算法的平均帧率为25.91fps，不可压缩水体加速模拟方法的平均帧率为30.10fps。相比于pcisph算法，本论文方法的计算效率提高了24.17％，相比于iisph算法，本论文方法的计算效率提高了16.21％。

表1水滴坠落模型模拟帧率对比表

本发明使用多线程提高邻近粒子搜索、核函数及其梯度求解效率。为了验证本发明程序多线程实现的有效性与高效性，分别启动1、4、8个线程对7104个粒子进行模拟，得到图6所示的模拟帧率图。由图6可知，启动1个线程时，无法达到实时模拟要求；启动4个线程时，模拟效率明显提高，且满足实时性要求；启动8个线程时，模拟效率进一步提高，从而验证了本发明多线程方法的有效性与高效性。

图7描绘了三维场景水滴坠落模拟效果图，图7(a)表示水滴初始状态，图7(b)、图7(c)为水滴坠落并与地面碰撞后多次反弹的效果，图7(d)为水体逐渐恢复稳定的效果，由图可知，水体在运动的过程中没有出现粒子穿透边界、粒子粘附边界、粒子在边界处聚集等紊乱现象。