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

本发明专利技术是一种基于IISPH(Implicit Incompressibility Smooth particle hydrodynamics，隐式不可压缩光滑粒子流体动力学)方法提高不可压缩水体模拟效率的实现方法。首先将速度无散度模型引入IISPH，修正IISPH方法每个时间步长内速度场的散度误差；其次，通过改进的密度恒定模型，使得密度恒定模型与速度无散度模型共享计算因子，实现不可压缩水体的高效模拟。通过引入无散度速度模型，实现了散度误差的修正，解决了密度误差随时间不断增加的问题；通过改进密度恒定模型，实现了速度无散度模型与密度恒定模型的计算因子共享，解决了冗余计算的问题。本发明专利技术能够达到17k粒子规模下30.10帧/秒的实时模拟速度，与PCISPH(Predictive‑Corrective Incompressible SPH，基于预测‑校正的SPH)方法相比，模拟效率提高了24.17％；与IISPH方法相比，模拟效率提高了16.21％。

【技术实现步骤摘要】
一种基于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].ComputerGr本文档来自技高网...

【技术保护点】
1.一种基于IISPH(Implicit Incompressibility Smooth particle hydrodynamics，隐式不可压缩光滑粒子流体动力学)提高不可压缩水体模拟效率的实现方法，其特征在于实现步骤如下：步骤一：将速度无散度模型引入IISPH，构建不可压缩水体模拟框架首先，基于IISPH搭建不可压缩水体模拟框架。IISPH分为预测与校正两个阶段，预测阶段使用除压力外的其他力计算粒子密度以及速度，校正阶段通过施加压力实现水体的不可压缩性。其次，改进IISPH的密度恒定模型，使其与速度无散度模型共享计算因子，减少冗余计算，提高水体模拟效率。最后，将速度无散度模型引入不可压缩水体模拟框架中。针对IISPH求解出的速度存在散度误差的问题，将速度无散度模型引入IISPH方法中，修正粒子的散度误差。步骤二：引入速度无散度模型，以修正IISPH方法的散度误差IISPH修正的粒子速度存在散度误差，即不满足

【技术特征摘要】
1.一种基于IISPH(ImplicitIncompressibilitySmoothparticlehydrodynamics，隐式不可压缩光滑粒子流体动力学)提高不可压缩水体模拟效率的实现方法，其特征在于实现步骤如下：步骤一：将速度无散度模型引入IISPH，构建不可压缩水体模拟框架首先，基于IISPH搭建不可压缩水体模拟框架。IISPH分为预测与校正两个阶段，预测阶段使用除压力外的其他力计算粒子密度以及速度，校正阶段通过施加压力实现水体的不可压缩性。其次，改进IISPH的密度恒定模型，使其与速度无散度模型共享计算因子，减少冗余计算，提高水体模拟效率。最后，将速度无散度模型引入不可压缩水体模拟框架中。针对IISPH求解出的速度存在散度误差的问题，将速度无散度模型引入IISPH方法中，修正粒子的散度误差。步骤二：引入速度无散度模型，以修正IISPH方法的散度误差IISPH修正的粒子速度存在散度误差，即不满足本发明通过引入速度无散度模型，并使用该模型计算压力项，从而实现对粒子速度的散度误差进行修正，粒子i所受压力的大小由式(1)决定。其中表示粒子i所受的压力，表示水体粒子的体积，表示压力梯度。引入表示对粒子i的速度v进行修正的刚性系数，可得压力梯度的计算式如式(2)所示。其中mj表示邻近粒子j的质量，表示核函数的梯度，即为满足速度无散度条件的刚性系数。粒子在运动的过程中满足动量守恒。对粒子i进行受力分析可知，粒子i所受的压力为邻近粒子j的压力之和，并且粒子间的压力是一对相互作用力。由于相互作用力大小相等，方向相反，因此，粒子i所受的压力与粒子i对邻近粒子j的压力满足其中表示粒子i对粒子j的压力，这意味着压力作为一种粒子间的作用力，以内力的形式存在，并且所有粒子的压力之和为0。从水体的整体看来，该压力不产生水体本身动能的改变，满足动能守恒定律。粒子i对邻近粒子压力的大小与粒子间的距离有关，其求解公式如式(3)所示。其中xj表示粒子j的三维坐标，求解器必须通过改变压力值修正粒子的速度，从而使得粒子的速度满足无散度性。密度相对于时间的变化率如式(4)所示。其中Wij表示光滑核函数，要保证水体的不可压缩性，即密度不发生变化，则式(4)求得的密度相对于时间的变化率应为0。式(4)中vi-vj表示粒子间的相对速度差，该速度差由压强差导致，不同粒子受到的压力存在差别，因此导致粒子间的速度不一样。由压力引起的速度变化量如式(5)所示。其中Δvi表示由压力引起的速度变化量，将式(5)代入式(4)，得到式(6)。由式(1)可以计算出粒子所受压力的大小，结合式(4)密度变化率为0，即可求解出刚性系数；使用该刚性系数修正粒子速度，能够使得速度的散度为0，从而使得密度相对于时间的变化率为0。将式(2)代入式(1)，得再代入式(6)，得到式(7)。根据式(7)求解出刚性系数如式(8)所示。其中且αi的大小只与粒子位置有关。式(8)为速度无散度模型，该模型通过求解刚性系数并使用该刚性系数计算压力，从而保证水体在运动的过程中密度变化率不发生波动，即保证水体的不可压缩性。最后根据式(9)计算粒子受到的合力，并使用该合力重新计算粒子的速度。其中表示粒子i所受合力，表示除压力外其他力的合力，表示粒子...

【专利技术属性】
技术研发人员：艾明晶，李锋，郑爱玉，
申请(专利权)人：北京航空航天大学，
类型：发明
国别省市：北京,11