虚拟现实智能硬件2021年12月3日第2*引文：杨忠耀，吴茂林，刘世光。基于Helmholtz分解的SPH虚拟现实智能硬件，2021，3（2）：118-128DOI：10.1016/j.vrih.2021.01.003·文章·Helmholtz分解SPH杨忠亚，毛林武，石光李*天津大学计算机科学与技术学院，智能与计算系，天津300072通讯作者，lsg@tju.edu.cn投稿时间：2020年12月3日修订日期：2021年1月1日接受日期：2021年1月20日国家自然科学基金资助项目（62072328，61672375）.摘要背景SPH方法在水体场景模拟中得到了广泛的应用。作为偏微分方程的一种数值方法，SPH方法可以方便地处理变形的复杂边界。此外，SPH算法实现相对简单，且结果稳定，不易发散。然而，SPH方法也有其局限性。为了进一步提高SPH方法的性能，扩大其应用范围，一系列制约SPH发展的关键和难点问题需要改进。方法将Helmholtz分解的思想引入到光滑粒子流体动力学（SPH）框架中，提出一种新的三维水流速度投影方法。首先，我们将亥姆霍兹分解应用于三维速度场，并将其分解为三个正交的子空间。然后，我们的方法结合SPH空间导数的思想，得到离散的泊松速度方程。最后，共轭梯度（CG）被用来有效地解决泊松方程。 结果实验结果表明，该方案适用于各种情况，并具有更高的效率比目前的SPH投影方案。结论与以往的投影方案相比，本文的方案不需要通过压力投影间接修正粒子速度，而是直接通过速度场投影修正粒子速度。新格式可以很好地集成到现有的SPH框架中，并可以应用于水与静态和动态障碍物的相互作用，甚至粘性流体。关键词 水模拟; SPH;亥姆霍兹分解;共轭梯度1引言自1977年首次提出以来，光滑粒子流体动力学（SPH）方法已被广泛应用于各种流体的模拟。SPH方法是一种求解偏微分方程的数值方法，可以很容易地处理变形和复杂的边界。此外，SPH相对简单，并提供难以发散的稳定结果[2]。经过几十年的发展，SPH方法在流体[3，4]、弹性固体[5]、泡沫[6]、雪[7]、多相流[8- 10]和流体控制[11- 14]的模拟方面取得了显著的成就。2096-5796/©版权所有2021北京中科学报出版有限公司Elsevier B. V.代表KeAi Communization Co. Ltd.提供的出版服务。这是CC BY-NC-ND许可证下的开放获取文章（http：//creativecommons.org/licenses/by/4.0/）。www.vr-ih.comZhongyaoyanGetall：Helmholtzdecomomosition-basedSPH119标准SPH（SSPH）[15，16]使用状态方程（EOS）计算压力，适用于可压缩流体。弱可压缩SPH[17]对弱可压缩流体施加了严格的时间步长限制，这限制了其性能。预测-校正不可压缩SPH（PCISPH）[18]基于SSPH和EOS迭代预测和粒子位置校正计算相应的压力值。PCISPH可以处理更大的时间步长，提高计算效率。它是基于粒子的不可压缩密度。IISPH（隐式不可压缩SPH）方法[19]允许更大的时间步长和更高的收敛速度，这进一步提高了计算效率，但其精度和稳定性仍有待提高。在细节方面，SPH在流体动压计算上的不精确，导致SPH在模拟流体边界层和旋涡等局部现象时精度不高。在运算速度方面，对于每个计算单元，SPH需要30- 50个周围粒子的相互作用信息，这比基于网格的方法所需的信息大得多。例如，有限体积法[20]通常需要5- 13个相邻单元的信息。另外，相邻粒子的额外搜索时间使得SPH在小尺度上的流体模拟更加耗时。在数学方面，SPH缺乏一个完整的误差分析系统。SPH模拟一旦开始，每个粒子将处于无序但非随机的状态，这使得SPH的误差分析非常困难。因此，提高SPH方法的精度和效率仍然是一个难题。本文提出了一种新的SPH投影格式来模拟水体。将Helmholtz分解的思想应用于三维速度场，以消除由压力项引入的速度发散。与以往的投影方案相比，该方案不需要通过压力投影间接修正粒子速度，而是直接投影速度场。此外，由于压力场和速度场的解耦，共轭梯度法（CG）可以有效地求解泊松方程。各种实验结果表明，本文提出的新方案可以与现有的SPH框架相结合，适用于水与静、动障碍物的相互作用，甚至可以推广到粘性流体。同时，与目前最先进的SPH投影格式相比，该格式具有更高的计算效率。2背景在介绍所提出的方案之前，我们首先简要解释与本节中的方法相关的背景知识和概念，包括SPH和Helmholtz分解方法。2.1平滑 颗粒 流体动力学（SPH）1999年，Stam[21]首次将SPH方法引入流体动画领域，而SPH在水模拟领域的首次应用要追溯到2003年[16]。SPH的基本思想是将连续流体看作离散流体质点的集合，复杂的流体运动是由流体和刚性质点之间的相互作用引起的。SPH在水模拟中的主要过程包括以下几个步骤：(1) 为所有粒子建立空间索引。(2) 计算所有粒子的压力、密度、粘性力和外力，然后根据这些力的作用计算粒子的加速度。(3) 根据粒子的加速度和时间步长（时间积分）计算下一时刻的速度和位置，然后处理必要的碰撞。120虚拟现实智能硬件2021年12月3日第2Ji ∑ijρρiJρj我J我(4) 渲染所有粒子。2.1.1SPH 插值 和 光滑 核函数因为SPH是基于离散粒子的，所以为了计算空间中任意位置的物理量，需要对某个邻域内的粒子物理量进行插值：Ai=∑A jm jW ij（1）其中，Ai是任意位置xi处的物理量，Aj是该位置xj，Wij是以下形式的核函数：W ij=W（xi-xj /λ）（2）其中λ称为光滑核半径;因此，W也称为光滑核函数。2.1.2空间导数基于可微光滑核函数的SPH方法避免了网格方法中的网格缠绕和变形问题，从而可以有效地计算空间导数。SPH方法通常使用空间导数的近似：<$Ai=ρi∑m j（Ai/ρi2+Aj/ρj2）<$W ij（3）∇⋅Ai=-ρ1∑m jAij⋅∇Wij(4)2A =2m jAxijW ij(5)jρ j其中Aij=Ai-Aj，Aij=Ai-Aj，xij=xi-xj。2.1.3力 分析 颗粒xij2+ 0.01λ2作用在颗粒上的力直接影响颗粒的运动，从而影响整个流体的动力学行为;因此，对颗粒上的力进行准确的分析至关重要。作用在质点上的力一般由三部分组成，记作F=F外部 +F压力+F粘度（6）F外部 =mg，F压力 =-mp，F粘度 =mμm2v（7）其中，Fexternal、Fpressure和Fviscosity分别是外力、压力梯度力和粘性力，μ是粘性系数。每个质点的加速度计算遵循牛顿第二定律，即：例如，ai=Fi/mi，可以通过将等式（7）代入牛顿第二定律来获得粒子的加速度：一 =g-βpi +μΩ2V（八）2.2亥姆霍兹分解1858年，亥姆霍兹提出了流场速度分解定理（以下简称亥姆霍兹分解定理）[22]。亥姆霍兹分解定理我121ZhongyaoyanGetall：Helmholtzdecomomosition-basedSPH指出，在适当的渐近无穷大，向量场可以分解为两部分：无除数部分表示旋转，其可以表示为矢势函数的旋度。无卷曲部分表示平移和压缩/扩展，其可以表示为122虚拟现实智能硬件2021年12月3日第2我Jρ我JJ ρ2i（ tρ2j（ t）我我J ρ Jijx2+ 0.01λ2标量势函数的梯度。这意味着当向量场的散度和旋度已知时，向量场可以被唯一地构造。其逆过程，i。例如，将矢量场分解为非发散和非曲率分量的过程称为亥姆霍兹分解，其可以描述为：=其中，ω是域Ω上的光滑向量场，具有光滑边界ω Ω，ω是标量势函数，ω是向量势函数。根据定义，是无卷曲的（×= 0），×是无div-free的（× = 0）。方程（9）可以缩写为ε= d + r，其中d是发散分量，r是与边界ε Ω相切的旋转分量。Chorin和Denaro等人，分别给出了Helmholtz分解的边界条件并证明了分解的存在性、正交性和唯一性[23，24]。具体地，亥姆霍兹分解必须满足以下边界条件之一：nr=0惠n=n，n×d=0惠n×（×）=t（ 10）其中n是边界的外法向量，n和t是法向和切向分量边界，分别。上述推导仅考虑无限空间中的情况，即，谐波场为零的情况。对于有界区域，可能存在非零谐波分量。本案中Helmholtz分解将向量场分成三个正交部分：n=d+r+h（11）其中h是谐波分量，其散度和旋度都为零。在这种情况下，亥姆霍兹分解需要满足由等式（10）表示的边界条件。3亥姆霍兹 基于分解 投影方案本节提出了一种基于三维速度场的亥姆霍兹分解的投影方案。该格式的优点在于求解时压力场和速度场的解耦。具体而言，由于计算保证了分解的正交性，因此一个项目的计算误差不会反映在另一个项目的计算中。这意味着使用速度场的Helmholtz分解进行投影比求解Navier-Stokes方程更有效，速度 和 压力 是 耦合。 内核我们的方案的过程如图1所示。3.1颗粒 速度预测图1我们方案的流程图。在每一个时间步的开始，我们的方案计算粒子的加速度的基础上的外力，并更新粒子的速度与此加速度。这个过程称为粒子速度预测。首先，使用平滑核函数来近似颗粒上的压力和粘性力：F压力（t）=-miip我=-m∑m（pi（t）+pj（t））W （t）（12）F粘度（t）=m μm2v=2μ∑m jvxijWijIJ（十三）当ρi（t）可以是一个加权平均数ρ i（t）=∑m jW ij和pii时，计算出的加权平均数为EOS：pi（t）=我IJ我123ZhongyaoyanGetall：Helmholtzdecomomosition-basedSPH我Jκ（（ρi（t）/ρ0）γ-1），其中参数κ和γ用于按比例缩放压力，称为刚度系数，ρ0是流体的静态密度。那么，当前粒子加速度ai（t） 可以通过将计算的压力F压力和粘滞力F粘度代入方程（8）。最后，可以根据粒子加速度计算每个时间步的粒子速度：v*i（t+ Δ t）=vi（ t）+ Δ tai（ t）（14）其中Δt是时间步长的大小。3.2三维速度场的亥姆霍兹分解由于压力项的存在，前一步计算的速度场含有发散，在处理不可压缩流动时需要消除发散。在这项研究中，投影操作是基于亥姆霍兹分解的原则。首先，速度场v*被分解成方程（11）形式的三个部分：v*=vd+vr+vh。（十五）由于我们假设流体速度场是无旋的，因此存在以下限制：v ×v* = 0。 因此，方程（15）中的旋转分量vr因此方程（15）可以简化为v*=vd+vh。为了确保不可压缩性，需要保留谐波分量h以消除发散：vh= v*-vd= v-ε（十六）通过将发散算子应用于方程（16）的两侧，在变形之后，我们可以得到：∇ ⋅ v*= ∇ ⋅ vh+ ∇ 2ϕ.（十七）对于调和向量场，项v*为零。因此，方程（17）可以简化为泊松方程：解完这个泊松方程后，我们可以得到标量场，然后我们把它代入方程（16），就可以得到无发散速度场v=vh3.3Poisson方程的共轭梯度解法泊松方程可以用不同的方法求解，例如Jacobi迭代[20]，连续超松弛[25]和CG[26]。其中，CG算法能够以高效的交互方式求解具有复杂边界条件的Poisson方程，是一种理想的求解方法。对于简化的Poisson方程，我们使用光滑核函数来近似和得到离散泊松方程：2方程（18）可以被写为线性方程组A=b。为了求解这个方程组，我们首先设定初始条件为n0= 0，q0=r0=b。然后，我们执行以下迭代计算：rT r rT rα=LL ，LL=+α q，r=r-α Aq，q=q+L+ 1l+ 1qL（十九）TAqL+ 1l l ll+ 1ll l l+ 1lrT rl迭代的终止条件是rL≤ λ，其中λ是迭代终止阈值。迭代后我们可以得到v=v*-φ，无发散速度场v可以计算为v=v*-到目前为止，基于亥姆霍兹分解的投影方案已经被完全引入。L124虚拟现实智能硬件2021年12月3日第24结果在本节中，我们将讨论我们方案的实验结果，并比较不同方法的结果。所有实验均在配备有2.6GHzIntel i7- 6700HQ CPU、32 GB存储器和NVIDIA GeForce GTX 1060 GPU的个人计算机上进行。我们将我们的投影方案集成到Koschier等人提出的SPH框架中。[27]。对于边界处理，我们使用Akinci等人提出的方案。[28]，使用He等人提出的方法处理表面张力。[29]，并使用Akinci等人[30]提出的方法重建液体表面。这些图像是使用Arnold渲染的。对于以下所有场景，我们将我们的方案与IISPH[19]和DFSPH（无发散SPH）[31]进行了比较。所有场景中的液体颗粒的半径为0.025 m，最大体积误差率设置为0.01%。此外，所有场景中的液体都是在一个立方体水箱中模拟的。考虑到照明因素，我们没有渲染水箱。4.1二维场景图2显示了二维场景中的渲染结果，其中灰色液体是水，蓝色液体是粘性流体，被障碍物隔开。此场景包含1726个液体粒子。在该场景中，所提出的方法与其他方法之间的效率比较如表1所示。在给定的时间步长下，该方法在平均迭代次数和仿真时间上都优于DFSPH（当时间步长为5ms时，仅略低于DFSPH），当时间步长为2ms时，该方法的优势最为明显.图2二维场景。（a）-（f）依次表示帧8、16、24、32、40和48的渲染结果4.2双发射极图3显示了双发射器场景中的渲染结果，其中左发射器从第1帧开始排水，并在第52帧停止，然后右发射器开始排水。此场景包含1938年的液体125ZhongyaoyanGetall：Helmholtzdecomomosition-basedSPH表1二维场景下不同时间步长下各方法效率比较对于DFSPH，avg.iter。列表示恒定密度求解程序和无发散求解程序的迭代次数Δt/msIISPHDFSPH我们的方案avg. iter。comp. 时间/msavg. iter。comp. 时间/msavg.iter.comp. 时间/ms0.52.014.6862.0/1.05.4702.05.6310.752.014.9282.0/1.06.0252.05.48512.014.2722.0/1.05.5672.05.93222.114.5982.0/1.05.8352.05.51434.223.4512.6/1.36.3672.76.220图3双发射器场景。（a）-（d）依次表示帧10、30、50和70的渲染结果。粒子数分别为2298、3066、3834和5560。粒子的数量随着时间的推移而增加。在该场景中，所提出的方法与其他方法之间的效率比较如表2所示。该方法的效率高于IISPH，但与DFSPH相比，我们的方案只有在时间步长为0.5ms时才有优势。表2双辐射源场景Δt/msIISPHDFSPH我们的方案avg. iter。comp. 时间/msavg. iter。comp. 时间/msavg. iter。comp. 时间/ms0.52.09.1622.0/1.03.7272.03.6060.752.08.0522.0/1.03.5652.03.65312.08.1852.0/1.03.5872.03.90822.08.6452.0/1.03.6062.03.71934.012.7533.2/1.13.9922.84.0264.3双 大坝 打破 关于sphere图4显示了一个双溃坝场景的渲染结果，其中一个漂浮的绿色球体与水交互。此场景包含11492个液体粒子。在该场景中，所提出的方法与其他方法之间的效率比较如表3所示。我们的方法的计算开销是轻于IISPH和类似的DFSPH。125ZhongyaoyanGetall：Helmholtzdecomomosition-basedSPH图4球体双溃坝场景。（a）-（c）依次表示帧15、35和55的渲染结果。表3球体双溃坝场景下不同时间步长下各方法效率比较Δt/msIISPHDFSPH我们的方案avg. iter。comp. 时间/msavg. iter。comp. 时间/msavg. iter。comp. 时间/ms0.52.030.6182.0/1.015.9012.015.5160.752.030.0832.0/1.014.8122.014.92012.134.1502.0/1.015.0852.015.38424.342.2662.1/1.015.5272.115.123311.580.4755.3/1.518.9875.118.6444.4微调器图5显示了一个微调器场景中的渲染结果，其中一个连续顺时针旋转的红色长方体存在于场景的底部。此场景包含5054个液体粒子。如表4所示，与IISPH和DFSPH相比，我们的方案需要更少的迭代，并在某些时间步长上实现更高的效率。图5旋转镜头。（a）-（f）依次表示帧10、20、30、40、50和60的渲染结果。4.5限制实验结果表明，本章提出的投影方案不仅可以很容易地集成到现有的SPH框架中，而且在二维场景模拟中表现良好126虚拟现实智能硬件2021年12月3日第2表4在spinner场景Δt/msIISPHDFSPH我们的方案avg. iter。comp. 时间/msavg. iter。comp. 时间/msavg. iter。comp. 时间/ms0.52.015.5862.0/1.09.6562.09.5440.752.015.2232.0/1.010.9262.09.57812.015.5772.0/1.09.7882.09.75422.016.0972.0/1.09.8842.09.67934.120.6742.5/1.112.1942.310.429固定障碍物交互场景和移动障碍物交互场景。在计算开销方面，我们的计划取得了显着的改善，平均迭代次数和计算时间相比，IISPH。此外，我们的计划的效率和性能也与DFSPH。然而，我们的方法有一定的局限性。首先，我们的方法假设三维速度场是无旋的，因此在进行亥姆霍兹分解时可以忽略无发散分量。这意味着我们的格式不能处理流体旋度，因此它不适合于模拟涡量引起的液体效应。此外，根据Ihmsen et al.[19]，应用程序CG无法在迭代之间钳制负压以消除内聚效应。事实上，在CG中钳制负压会导致极不稳定的压力场。最后，在用CG求解线性方程组时，例如，A=b，则要求系数矩阵A是对称的。这意味着CG仅适用于模拟密度均匀的流体。因此，未来仍需要一种适用于非均匀密度流体的泊松方程的解。5结论 和 今后工作提出了一种基于亥姆霍兹速度场分解的投影方案来模拟水体。本文首先给出了三维速度场的Helmholtz分解形式，然后将其与SPH方法相结合，给出了速度场Poisson方程的近似离散形式。最后，本文利用CG算法对离散Poisson方程进行了有效的求解。多组实验表明，与目前最先进的SPH投影方法相比，该方法能够更有效地计算谐波速度场。该方案也适用于与各种类型的障碍物交互的场景，并可扩展到粘性流体模拟。对于未来的研究，我们希望提高我们的方法的稳定性，以允许更大的时间步长。此外，本研究假设速度场是无旋的，所以我们计划考虑旋度。最后，将我们的方案推广到更一般的现象，如多相流，流体-刚性耦合，甚至流体和刚性之间的转换。竞合利益我们声明我们没有利益冲突。引用1李文，李文，等.光滑粒子流体动力学：理论与应用.北京：科学出版社，2000。皇家天文学会月报，1977年，181（3）：375DOI：10.1093/mnras/181.3.3752张文辉，张文辉.用于强流固耦合的互连SPH压力求解器。ACM127ZhongyaoyanGetall：Helmholtzdecomomosition-basedSPH图形学报，2019，38（1）：5 DOI：10.1145/32849803包克，张华，郑丽玲，吴鄂华。流体动画的压力校正SPH。计算机动画与虚拟世界，2009，20（2/3）：311DOI：10.1002/cav.2994黄春，朱军，孙海琴，吴鄂华.为逼真的VR环境优化SPH流体模拟。计算机动画与虚拟世界，2015，26（1）：43DOI：10.1002/cav.15645张文辉，张文辉，张文辉，张文辉.不可压缩线弹性固体的隐式SPH公式。计算机图形学论坛，2018，37（6）：135DOI：10.1111/cgf.133176杨伟杰，李晓梅，李晓梅.含泡沫的湍流微极SPH流体。IEEE Transactions on Visualization and ComputerGraphics，2018，25（6）：2284DOI：10.1109/tvcg.2018.28051107张文辉，张文辉.雪模拟的隐式可压缩SPH求解器。ACM图形学报，2020，39（4）：36DOI：10.1145/3386569.33924318张文辉，张文辉.基于物理的流体和固体模拟的光滑粒子流体动力学技术。20209任宝，李春芳，严晓，林明春，胡世明.使用混合物模型的多流体SPH模拟。ACM图形学报，2014，33（5）：171DOI：10.1145/264570310严新，姜永涛，李春芳，马丁若若，胡世民.多相SPH模拟交互流体和固体。ACM图形学报，2016，35（4）：79DOI：10.1145/2897824.292589711冯国，刘世光.带有变形约束的细节保留SPH流体控制。计算机动画与虚拟世界，2018，29（1）：e1781DOI：10.1002/cav.178112张晓阳，刘世光.具有动态细节的并行SPH流体控制。计算机动画与虚拟世界，2018，29（2）：e1801DOI：10.1002/cav.180113张晓阳，刘世光. SPH触觉交互与多流体模拟。虚拟现实，2017，21（4）：165-175 DOI：10.1007/s10055-017-0308-114张晓阳，刘世光.具有自适应湍流细节的SPH流体控制。计算机动画与虚拟世界，2015，26（3/4）：357DOI：10.1002/cav.163715Monaghan J J.平滑粒子流体动力学。天文学与天体物理学年鉴，1992，30（1）：543574DOI：10.1146/annurev.aa.30.090192.00255116Müller M，Charypar D，Gross M.用于交互式应用的基于粒子的流体模拟。在：ACM SIGGRAPH/ Eurographics计算机动画研讨会，2003年，15417Becker M，Teschner M.自由表面流的弱可压缩SPH。在：ACM SIGGRAPH/Eurographics计算机动画研讨会，2007年，20918Soldier B，Pajarola R.预测校正不可压缩SPH。ACM图形学报，2009，28（3）：40 DOI：10.1145/1531326.153134619杨文，李文.隐式不可压缩SPH IEEE transactions on Visualization and Computer Graphics，2013，20（3）：42620Acharya S，Baliga B R，Karki K，Murthy J Y，Prakash C，Vanka S P.计算流体动力学中基于压力的有限体积方法。传热学报，2007，129（4）：407虚拟现实智能硬件2021年12月3日第2128DOI：10.1115/1.271641921斯塔姆·J，稳定性液体见：SIGGRAPH，1999，12122吴晓松，李晓松，李晓松.亥姆霍兹-霍奇分解综述。IEEE Transactions on Visualization and Computer Graphics，2013，19（8）：1386DOI：10.1109/tvcg.2012.31623放大图片作者：A.流体力学的数学入门纽约州纽约市，Springer US，1990 DOI：10.1007/978-1-4684-0364-024玛丽亚·德纳罗Helmholtz-Hodge分解在一般边界条件下不可压缩流投影方法中的应用。国际流体数值方法杂志，2003，43（1）：43-69 DOI：10.1002/fld.59825放大图片作者：Jiang D.抛物扩散加速并行SOR方法在3D任意几何体上求解非结构网格泊松方程。国际计算流体动力学杂志，2016，30（5）：370DOI：10.1080/10618562.2016.123404526Helfenstein R，Koko J.并行预处理共轭梯度算法在GPU上。计算与应用数学学报，2012，236（15）：3584DOI：10.1016/j.cam.2011.04.02527张文辉，张文辉.基于物理的流体和固体模拟的光滑粒子流体动力学技术。在：欧洲图形，2020年，128张文辉，张文辉.用于不可压缩SPH的多功能刚性流体耦合。ACM图形学报，2012，31（4）：1DOI：10.1145/2185520.218555829何晓文，王海明，张芳娟，王红，王国平，周凯.基于SPH的自由表面流中稀疏采样薄特征的鲁棒模拟。ACM图形学报，2014，34（1）：1-9DOI：10.1145/268263030张文辉，张文辉，张文辉.粒子流体的平行曲面重建。 计算机图形学论坛，2012，31（6）：1797DOI：10.1111/j.1467-8659.2012.02096.x31Bender J，Koschier D.无发散光滑粒子流体动力学。第14届ACM SIGGRAPH / Eurographics计算机动画研讨会论文集。洛杉矶加利福尼亚州，纽约州，美国，ACM，2015，147 -155DOI：10.1145/2786784.2786796
