论文引用格式：Shao X Q， Zhang X， Yang S H and Jin Y Z. 2023. SPH fluid-solid interacted real-time solid fracture animation. Journal of Image and Graphics， 28（08）：2447-2460（引用格式:邵绪强， 张鑫， 杨仕和， 金佚钟. 2023. SPH固流交互中的实时固体破碎动画模拟. 中国图象图形学报， 28（08）：2447-2460）［0　引言真实地模拟现实中的流体运动是图形学的一个重要研究方向。早期由于计算机硬件性能的制约，对流体现象的模拟主要采用易于求解的基于数学或几何的方法。随着计算硬件技术的快速发展，计算性能大幅提高，以求解复杂物理模型为基础的基于物理的模拟技术发展迅速，在计算机中真实地模拟水流涌动、浓烟滚滚以及水花飞溅等复杂流体线性成为可能。基于物理的流体模拟就是通过数值方法求解描述流体运动状态的纳维—斯托克斯（Navier-Stokes equations，N-S）方程或浅水方程（shallow water equations，SWE）（邵绪强 等，2016），得到每一时刻流体的运动状态。SPH作为经典的基于物理的流体模拟方法，由Gingold和Monaghan（1977） 以及Lucy （1977） 分别提出，早期主要用于解决天体物理问题。Müller等人（2003）将SPH方法首次引入流体的模拟领域，随后基于SPH的流体方法快速发展。Becker和Teschner（2007） 提出的微可压缩SPH （weakly compressible SPH，WCSPH）方法，通过使用Tait方程替代理想气体状态方程和使用较大的刚度系数来保证流体的不可压缩性。WCSPH算法由于相对简单，在固流交互模拟、流体与多孔介质的交互模拟、多相流以及颗粒流模拟等关注流体交互问题的研究上有大量应用。流体在运动过程中不可避免地会与各种固体进行交互，如山洪、海啸和泥石流等自然现象中就包含了流体运动、固流交互、多孔介质、多相流以及固体破坏等诸多复杂物理现象，当前对这类复杂现象进行基于物理的真实感模拟的难度仍旧非常大。现有的SPH固流交互研究大多关注固流交互中固体的运动、形变、渗透和侵蚀等现象，对流体冲击引起的固体破碎模拟研究少之又少。而关于固体破碎模拟的研究仅关注固体与固体之间的碰撞交互，同样没有涉及流体冲击。在固体破碎的研究中，基于物理的模拟方法主要通过分析固体应力应变或固体的材料特性等物理特性，进而构建物理模型求解破碎的发生和发展来模拟破碎效果。早期Terzopoulos和Fleischer （1988）通过分析物体受力瞬间的形变，对固体模型直接施加切割算法，物理模型较为简单，易于实现，但破碎效果的真实感与现实世界相差甚远。Norton等人（1991）首次引入质点—弹簧模型来研究固体破碎模拟，在受力时通过对固体的质点—弹簧系统进行力学分析，计算和模拟固体的应变、断裂，对破碎后的碎块进行同样的处理，提高了模拟的真实感。O’Brien和Hodgins（1999） 通过建立四面体有限元模型，对固体进行应力分析，计算应力传播信息来模拟固体的破碎。Smith等人（2001）改进了质点—弹簧模型，提出了适用于脆性材料的高效可控的破碎算法。O’Brien （2002）将之前的工作推广到塑性材料的破碎模拟，取得了较好的效果表现。Müller等人（2004） 通过用点来标识物体的体积和表面，采用最小二乘法求解离散位移场模拟固体破碎。Pauly等人（2005）通过物体离散点采样，采用基于高动态表面采样的方法对弹性和塑性材料进行断裂模拟，并且能够动态模拟表面裂缝。Busaryev等人（2013）利用四面体有限元模型结合材料自身特性，提出一种自适应破碎算法，在纸张、布料等薄脆材料的撕裂模拟上表现良好。Hahn和Wojtan （2015，2016）采用基于边界元的方法，对脆性材料的断裂进行模拟，优化了断裂裂纹的延展细节。He等人（2018）提出一种基于粒子表示的方法，对固体形变中的断裂细节进行较为真实的模拟，但实时性仍有待提高。总体上，基于物理的方法虽然能够得到真实的、细节丰富的破碎模拟效果，但由于其往往包含复杂的应力分析和力学模型求解过程，计算复杂度高，难以在实时性要求较高的交互场景中应用。在电子游戏等对物理真实性宽容度较高、强调实时交互性的应用场景中，一些基于几何的方法因为求解快速而得到了广泛应用。Raghavachary （2002）提出一种利用Voronoi图结构对物体表面进行多边形划分的方法，能够在模拟多边形表面上的脆性断裂。Eberle等人（2003）提出一种基于对象空间点插值的方法，基于Havok引擎实现了固体的局部破碎模拟。Oda和Subramaniam （2006）提出将固体以四面体网格形式表示，使用多阶段动态细化技术，通过动态控制四面体网格的数量，实现了低计算消耗下的多粒度破碎模拟；Liu等人 （2011）采用一种基于几何的无网格方法来模拟脆性固体的破裂，实现了交互级的固体破碎模拟。Schvartzman和Otaduy（2014）提出一种基于应变能分布和Voronoi图的破碎模拟方法，在基于几何的方法中引入物理模型，使破碎效果有了一定的物理依据，一定程度上提高了模拟的真实感。赵夫群等人（2017） 提出一种用于匹配碎块断裂面的非物理方法，实现了对固体破碎及碎块复原的双向模拟。吕长建等人（2018）提出一种基于Voronoi的实时破碎模拟算法，通过以不同方式建立Voronoi种子点集来控制破碎的行为模式，丰富了刚体破碎的效果。非物理方法并不追求极致的真实物理效果，具有较高的实时性，但模拟的真实感亟待提高。在计算复杂的流体模拟领域，由于求解流体运动原本就十分复杂，要实时模拟固流交互并且求解固体的多种复杂行为，对固体破碎算法的要求也就更高。基于几何的方法具有良好的实时性但真实感较差；基于物理的方法具有良好的真实感，但因物理模型求解复杂，难以达到实时性要求。郭琳（2019）利用预计算的思想，采用预先计算固体破碎并在实时运行中以碎块的无缝拼接代替实时破碎的方法，模拟了固流交互中固体的破碎，但预计算方法依赖破碎模板，实时运行时的破碎与实际的固体受力无关和物理特性无关，难以达到真实感要求。赵建旺等人（2019）将Schvartzman和Otaduy （2014）提出的基于应变能分布的方法与流体隐式粒子法（fluid implicit particle，FLIP）和无散度SPH法（divergence-free SPH，DFSPH） 结合，采用基于固体几何质心的快速Voronoi方法模拟固体破碎，实现了FLIP-DFSPH流体冲击下固体破碎的真实感模拟，但未将固体网格离散为拉格朗日粒子形式，应变能求解方法无法与流体模拟框架高效耦合，且模拟系统应用GPU加速技术，无法满足实时性要求。总体而言，在固体模拟领域，现有研究主要关注固体与固体之间的交互，而在流体模拟领域，研究的焦点集中于对固流交互中的动力学、形变等现象的模拟，对于固流交互引发的固体破坏研究很少。少量相关研究主要以引入固体模拟领域的方法或思想为主，与SPH固流交互求解器的耦合性较差，大量的耦合开销引起模拟系统性能降低。对于固流交互过程中固体破碎现象的模拟问题，由于流体、固流交互和固体行为的物理模型都十分复杂，要兼顾真实感与实时性地进行模拟，需要结合固体求解器自身特性与固流交互模拟系统特性提高多求解器耦合性，若直接引入固体破碎方法往往会导致模拟系统构成复杂、计算量大以及耦合性差。为真实地模拟固流交互过程中固体受流体冲击发生破碎的现象，针对相关研究较少、已有研究存在真实感与实时性难以兼顾、与固流交互模拟系统耦合性差等问题，本文基于SPH统一粒子框架，提出一种物理与几何混合的实时固体破碎模拟方法。该方法理论模型的构建基于断裂力学理论，针对固体的粒子形式与SPH方法的邻域特性，通过实时分析流体粒子与固体边界粒子之间的能量传递、固体自身的能量转化，得到满足能量限界条件的粒子，并将它们作为破碎的启发点；而后以启发点集为种子点，采用可并行构建的Voronoi图作为固体碎块生成方法，对固体体积空间进行快速子空间划分，完成固体碎块的生成。由于本文方法与SPH统一粒子框架能够高度耦合，具有可并行化特性，故可以进行并行化处理并加载至GPU作大规模并行加速以提高实时性。1　SPH固流交互模拟1.1　流体运动模拟流体力学中通过N-S方程来描述不可压缩流体运动。在拉格朗日视角下，这一非线性偏微分方程的形式为dvidt=-∇pi+gi+μ∇2vi （1）式中，vi表示流体粒子i的密度，vi表示粒子i的速度，pi表示粒子i处的压强， gi表示流体粒子i所受到的合外力，一般为重力， μ表示流体动力粘性系数，dρidt为物质导数。SPH方法是一种经典的拉格朗日粒子方法，它将整个流体域离散为粒子，每个粒子拥有流体的一小块体积空间，如图1所示，并且携带流体的一些场量，任意粒子处的场量均连续地分布在整个光滑核函数的紧支域半径h内。以光滑核函数作为权重系数来考虑核紧支域半径h内的邻域粒子贡献，可以得到任意离散粒子处的场量，具体为Ai=∑jmjρjAjW(xij,h) （2）式中，Ai表示某位置为xi的粒子i处的任一场量，j为粒子i在核半径h内的某一邻域粒子，mi为邻域粒子j的质量，ρi为邻域粒子j的处的密度，W(xij,h)为核半径为h的光滑核函数，xij为粒子i到粒子j的位移向量。10.11834/jig.220176.F001图1拉格朗日视角下流体粒子的受力情况Fig.1Forces of the Lagrangian fluid particle将式（2）用于计算位置为xi处的流体粒子i的密度，可得ρi=∑jmjWpoly6(xij,h) （3）式中，Wpoly6(xij,h)为计算密度使用的特定核函数 （Müller等，2003）。根据式（1），流体粒子运动由其所受合力决定，包括粒子之间因密度差异产生的压强力、因速度差异产生的粘性力、作用于整个流体的重力以及其他力。粒子i所受压强力Fip可以表示为Fip=-mi∑jmjpiρi2+pjρj2∇Wspiky(xij,h) （4）式中，Wspiky(xij,h)为计算压强力时使用的特定核函数（Müller等，2003）， pi和pj分别为粒子i处和它的邻域粒子j处的压强。本文使用Tait方程 （Becker和Teschner，2007）求解粒子i处的压强pi，具体为pi=kρiρ0γ-1 （5）式中， k为刚度系数， ρ0为流体静止密度，一般取ρ0=998 kg/m3，γ为Tait方程参数，一般取γ=7（Becker和Teschner，2007）。粒子i所受粘性力Fiv可以表示为Fiv=μ∑jmjρjvj-vi∇2Wviscosity(xij,h) （6）式中， μ为流体的人工粘性系数，Wviscosity(xij,h)为计算粘性力时使用的特定核函数（Müller等，2003）。粒子i所受重力Fig为Fig=mig （7）根据粒子所受合力Fit，在一个时间步长Δt内作数值积分即可更新流体粒子i的速度和位置。本文采用半隐式欧拉法作时间积分，具体为vi(t+Δt)=vi(t)+FitmiΔt （8）xi(t+Δt)=xi(t)+vi(t+Δt)Δt （9）式中，vi(t)和xi(t)分别为更新前粒子i的速度和位置，vi(t+Δt)和xi(t+Δt)为更新后粒子i的速度和位置。1.2　固流交互模拟SPH统一粒子框架的思想是将场景中所有参与交互的模拟物体都以粒子形式表示，而固体需要从三角网格形式采样成粒子形式。1.2.1　固体网格的自适应粒子化本文使用子空间轴对齐包围盒（axis align bounding box，AABB）相交策略进行粒子化采样，在预处理阶段将固体采样成单层边界粒子层形式，以得到具有均匀间距的固体边界粒子分布，自适应采样方法流程如图2所示。该方法首先计算整个网格模型的AABB包围盒B0，按照给定的采样间距r0对包围盒B0进行均匀划分，并对整个包围盒按给定分辨率进行扫描，判断模型的每一个三角形面Ti是否占据某些空间单元。若有，则将所有被占用的空间单元标记为非空。最后将这些非空单元转化为粒子，即得到最终的边界粒子集合。10.11834/jig.220176.F002图2三角网格模型的自适应粒子化采样Fig.2Adaptive particle sampling for triangular mesh model与流体粒子类似，所有边界粒子都携带一些场量，用于向边界附近流体粒子贡献密度、防止流体粒子穿透固体边界、参与固流双向耦合交互力计算、固体运动求解等。1.2.2　基于权重因子矫正的固流交互力计算SPH方法中对粒子任意场量的求解都依赖其核函数半径内的邻域粒子贡献。如图3（a）所示，在流体和固体边界交界面附近，由于流体粒子的邻域流体粒子数量不足，会导致插值得到的流体密度值过小，产生边界粘附伪影。10.11834/jig.220176.F003图3边界附近的流体粒子邻域Fig.3Neighborhoods of particle near the boundary（（a） the fluid particles near the boundary haveinsufficient neighbors； （b） the boundary particlescontribute to the fluid particle field quantity）为了消除这一失真现象，提高固流交互模拟的稳定性和真实感，需要在计算边界附近流体粒子密度时引入邻域边界粒子的贡献（Akinci等，2012），如图3（b）所示，对密度进行校正，即ρi=∑jmjWpoly6(xij, h)+γ1∑bmbWpoly6(xib, h) （10）式中， mb为邻域粒子中边界粒子b的质量， γ1为密度体积权重矫正因子。由于本文使用单层边界粒子，在某边界附近流体粒子处，处于流体粒子邻域内的局部边界粒子数可能并不足以填补缺失的邻域流体粒子，所以需要通过体积权重因子自适应调节边界粒子对流体粒子的贡献权重大小。流体粒子与边界粒子之间存在的主要交互力包括压强力和摩擦力，在计算交互力时也需要引入权重校正系数。边界附近流体粒子受到的压强力为Fibp=-γ2mi∑bmbpiρi2+piρ02 ∇Wspiky(xib, h) （11）式中，γ2为压力权重校正因子。固体边界粒子可视为参与压强力计算的“静态流体粒子”，所以可假设其密度始终为流体静止密度ρ0。并且，固体粒子受到来自流体粒子的压力与流体粒子受到的来自固体粒子的压力是一对相互作用力，根据牛顿第三定律，它们大小相等方向相反。边界附近流体粒子与固体边界粒子之间的摩擦力可以等同为粒子之间的粘性力，同样这也是一对相互作用力。流体粒子在边界处受摩擦力可以表示为Fibv=μb∑bmbρ0vb-vi ∇2Wviscosity(xib, h) （12）式中，vb为固体粒子b的当前运动速度， μb为人工施加的固流边界摩擦力因子。2　物理和几何混合的固体破碎模拟2.1　破碎过程的能量平衡固体断裂、破碎的能量平衡模型是断裂力学的一个基本理论，核心是将整个固体视为一个热力学体系时，固体的各种行为诸如运动、形变和破碎等均满足能量守恒。当外应力作用在固体上时，固体的机械能增加，外应力与应力势能关系符合材料本构关系。在外应力作用下，固体会发生应变且应力势能不断增加直到达到材料断裂临界能阈值，此时破碎发生并生成新的表面，应力势能转化为固体自由表面能，余下能量转化为内能、动能等其他形式。从微观层面来讲，固体微粒之间存在相互作用力。当微粒间距处于平衡距离时，相互作用力不表现吸引或排斥属性；当微粒间距小于平衡间距时，相互作用力表现为排斥力；当微粒间距大于平衡间距，相互作用力表现为吸引力。由此可见，如果要使具有相互作用力的微观粒子间距增大，需要施加外力将其间距拉伸，使具有结合关系的微粒远离平衡位置。当这个力所做的功足够大时，微粒间结合关系会断开，相互作用力消失，即固体在此微粒处也发生了不可逆的断裂。2.2　基于物理的固体破碎计算外部应力载荷对固体系统施加应力会引起固体能量的变化，当应力势能达到临界极限时，固体发送断裂，此时用于断裂的应力势能转化为固体的自由表面能。在本文的固流交互场景中，流体作为外部载荷向固体热力学系统传递能量，其能量转化满足热力学第一定律，即EL=Eξ+Ek+Es+Ef （13）式中，EL是外载荷应力向系统传递的能量，Eξ为固体的应力势能，Ek为固体的动能，Es为固体的自由表面能，Ef为固体的内能。在流体与固体交互时，EL表现为参与交互的流体粒子的总动能变化，即EL=12∑imiΔvi)2 （14）初始时，固体处于静止状态。冲击发生时，流体携带的动能转化为固体的应力势能和内能。当应力势能达到临界点时，破碎发生并产生新的表面。自由表面能是固体断裂、表面积增大时，破坏原有分子间结合键所需的总能量。自由表面能仅在应力势能足以破坏固体内部分子间作用力时由应力势能转化，若应力势能不足以产生新的表面，则它稍后再次转化为流体的动能。因此，整个冲击破碎过程中，固体系统的能量转化关系为Eξ=Es+∑rEkr （15）式中，Ekr表示固体碎块r的动能。本文仅考虑刚体的破碎现象，其应力势能仅作为上述能量转化的中间量，不会独立留存于固体的热力学系统。稳定状态下分子的间距为平衡间距x0，当分子间距小于x0时，分子间相互作用力表现为排斥力，否则表现为吸引力。若要引发分子结合键断裂，外应力需要克服分子间作用力做功。当分子间距变化量为Δx时，分子键结合能Eσ为Eσ=∫x0x0+Δxϕ(x)dx （16）式中，ϕ(x)为分子间作用力。要使固体发生破碎，外部能量必须足以使固体内部分子间距增加至临界点以上。当分子间距大于分子结合能作用半径κ时，分子结合键不可逆转地断裂，分子间不再有相互作用力。图4展示了外力拉伸分子间距时，分子间作用力的变化情况。10.11834/jig.220176.F004图4分子间作用力与分子间距的关系示意图Fig.4The relationship between molecular force and distance假设在Δx=0及Δx=κ时，分子间作用力为0，则分子间作用力满足ϕ(x)=ξπxκ （17）式中， ξ为成键分子对所受的外应力，则分子键结合能为Eσ=∫x0x0+Δxξπxκdx （18）在SPH统一粒子框架中，由于参与固流交互的固体以离散粒子形式表示，并且粒子与粒子在核半径h具有邻域关系，因此可以将具有相互作用的分子成键关系建模为粒子与核半径内的邻域粒子的作用关系。对于固体上每一个边界粒子，取初始采样间距作为分子平衡间距x0，以核半径h内粒子数密度作为固体粒子xri处的局部分子结合键数，则固体粒子xri处的局部结合能为Eri=Eσ∑rjW(xri-xrj, h) （19）即对于以边界粒子形式表示的固体，若某粒子处因受力而与邻域固体粒子发生体积断裂，用于断裂的能量应大于粒子与其邻域粒子的局部结合能。由于断裂本身产生两个表面，外应力势能转化得到的自由表面能和分子结合能的关系为Eri=2Es （20）因此，破碎过程应力势能转化方程可以改写为Eξ-2Es=Ek （21）当Eξ2Es时，固体仅受动力学影响而不会发生断裂；若较短时间内（如NΔt个时间步长），外应力向固体边界粒子传递极大能量，以至于Eξ2Es时，则该粒子应与其邻域粒子发生体积断裂。外应力除转化为自由表面能外，余下的能量继续传递至该粒子上，用于新产生的碎块的动力学计算。2.3　基于几何的固体碎块生成三维Voronoi图是由一组空间中相邻两点之间线段的垂直平分面构成的连续多面体。三维Voronoi图的对偶图为Delaunay四面体网。Delaunay四面体的4个顶点分别是Voronoi图的相邻种子点，Voronoi图中多面体的顶点为是Delaunay四面体的外接球球心。由于它们的这种对偶关系，通常可以采取先构建Delaunay四面体网，再转化为对偶Voronoi图的间接方法构建Voronoi图。使用增量插入法构建Voronoi图的过程如图5所示。10.11834/jig.220176.F005图5构建Voronoi图的增量插入法Fig.5Insertion method for constructing the Voronoi根据本文的固体破碎计算方法，在受到来自流体的冲击力时，固体边界粒子层会收到来自流体的能量传入，这些能量将会转化为其他形式的能量并最终达到能量平衡。当破碎发生时，固体某处边界粒子由于短时间传入巨大的能量，导致其与邻域粒子的间距被拉伸至极限间距以上，局部结合关系断开，此时该粒子应当作为破碎的发生点，加入启发点集合。由于固体模型的三角形网格拓扑和体积形状往往比较复杂，无法将其作为可参数化的几何空间对Voronoi图构建施加边界约束，因此需要基于固体模型网格的顶点建立Voronoi图边界约束条件。而Voronoi图对网格模型进行细分，实质上是以Voronoi图的多面体对固体体积空间进行细分和三角形重构。由于增量插入法构建的Voronoi图中初始种子点最终会成为多面体顶点，因此可以将网格模型的顶点作为Voronoi的边界顶点，在构建完成后将最外层的Voronoi多边形与这些边界顶点连接完成边界约束施加和网格空间闭合。三角化细分后的多个子模型即为此次破碎过程生成的碎块，基于各碎块空间进行一次快速自适应粒子化采样，并将原始边界粒子匹配至碎块，原始边界粒子的非交互场量根据质心距和邻域数密度分配至新碎块边界粒子，至此完成固体碎块的生成。3　GPU并行实现SPH方法采用粒子形式离散问题域，每个粒子处的场量求解可与其余粒子独立进行，这样的特性使得SPH统一粒子框架在GPU大规模并行加速上具有得天独厚的优势。本文的固体破碎计算方法和采用的碎块生成方法基于SPH统一粒子框架特性构建，同样具有良好的可并行性。为确保模拟系统的实时性，本文将整个模拟系统进行并行优化，并针对求解过程中数据的显存分布特性进行数据结构设计，充分利用GPU强大的并行计算能力进行加速。3.1　流体的并行模拟GPU的大规模并行计算特性使得运行在它之上的程序必须额外考虑线程间数据同步的问题。在本文采用的框架中，求解过程之间的数据依赖是GPU并行化需要考虑的关键因素，若当前帧某环节对邻域粒子场量存在依赖，则必须确保该邻域粒子的对应场量在被访问之前就已经完成更新写入。根据固流交互求解环节的数据依赖，本文将求解器过程拆分为多个子过程，每个子过程映射为一个GPU核函数，由所有GPU线程并行执行。图6为根据各环节依赖关系对整个流体模拟过程进行步骤划分的结果，包括密度与压强计算、合力计算和时间积分3个具有前后数据依赖的并行步骤。各步骤并行执行细节如下：10.11834/jig.220176.F006图6流体模拟各关键环节的数据依赖Fig.6Data dependence of each stage of fluid simulation1）密度与压强计算。线程独立访问每一个粒子，根据当前帧粒子自身的位置进行邻域粒子搜索（neighborhood particle search，NPS），访问其邻域粒子位置，计算密度。根据密度结果，求解当前帧该粒子处的压强。2）合力计算。线程独立访问每一个粒子，根据当前帧粒子自身的位置进行NPS，访问邻域粒子密度、压强和位置数据，计算压强力；访问邻域粒子速度、位置数据，计算粘性力；将重力与压强力、粘性力相加获得当前帧粒子所受合力。3）时间积分计算。线程独立访问每一个粒子，根据当前帧粒子所受合力计算下一帧速度；根据求解得到的下一帧速度计算粒子下一帧位置。3.2　固体破碎的并行计算固体系统传入的能量需要由参与交互的流体粒子损失的动能来计算，而每个流体粒子损失的动能会分散传递至邻域内所有固体粒子，因此对于固体传入能量的计算必须从流体粒子视角进行。根据数据依赖关系，固体破碎计算过程可细分为能量计算、破碎计算和动能计算环节。1）能量计算。为每一个流体粒子指派一个线程，根据流体粒子在预设定的若干个时间步长内的速度变化，计算其动能损失；根据人工摩擦力和粘性力因子，排除动能中因固液交界面摩擦、流体内部粘性导致的动能损耗，余下动能向流体粒子邻域内的固体粒子按距离权重传递。2）破碎计算。为每一个固体粒子指派一个线程，根据系统传入能量，按照局部结合能计算该处传入的能量是否足以引发结合键断裂。若能量足够大，则将该固体粒子标定为断裂启发点，用于后续的Voronoi图构建。3）动能计算。为每一个固体粒子指派一个线程，根据因断裂转化为自由表面能的能量，计算固体修正速度。另外，因为固体粒子之间的相对位置仅会因破碎发生而改变，对固体粒子的分子结合能分析可以在初始化阶段预计算。在运行时，若当前帧未发生破碎，则无需更新分子结合能；若当前帧发生破碎，则需要重新计算分子结合能。3.3　Voronoi图的并行构建Ray等人（2018）提出一种Voronoi图并行化构建方法。在正式进行增量构建Voronoi图之前，先执行一次基于空间一致性网格的邻域种子点搜索，并将所有种子点按照网格单元索引进行并行排序。这样的处理与本文NPS优化方法高度类似，因此可以将NPS中的一致性网格复用于Voronoi图的并行构建，增加碎块生成算法与模拟系统的耦合性。创建Voronoi图前，首先需要对种子点集合执行NPS，确定每个种子点的邻域种子点集。与一致性网格存储SPH粒子的方式类似，对种子点集合按照网格索引执行并行排序。由于碎块上限较低，种子点的数量往往较少，在空间分布上种子点稀疏且仅存在于固体体积空间内，因此对种子点集的NPS可以采用空间一致性网格的局部子集，通过固体粒子的空间分布，控制种子点NPS的搜索范围，提高实时构建效率。3.4　流体的真实感绘制由于SPH方法模拟的流体以粒子形式存在，且流体运动过程中粒子的位置不断改变，需要重建流体表面信息才能够进行渲染。为了尽可能降低GPU和中央处理器（central processing unit， CPU）间通信开销，本文采用基于屏幕空间的流体渲染（screen space fluid rendering， SSFR）算法（van der Laan等 2009），与GPU并行计算耦合，使渲染管线可以直接访问显存中的流体粒子数据进行渲染，避免流体表面重建和网格提取过程的巨大计算和通信开销。SSFR方法的流程如图7所示，具体步骤如下：10.11834/jig.220176.F007图7SSFR方法流程Fig.7The processes of SSFR1）渲染深度纹理。通过以点精灵的形式绘制粒子并计算球径深度偏移，利用渲染到纹理 （render to texture，RTT）技术将深度渲染到高精度纹理。2）渲染厚度纹理。水体的光学透射的光衰减与光路传播深度有关，在SSFR中，某视角下流体的厚度能够反映出水体透射深度。附加球径深度偏移的点精灵绘制的流体粒子在渲染管线中以叠加方式混合能够得到厚度信息，再利用RTT技术即可存储到高精度纹理。3）重建法线。深度信息可结合屏幕空间信息，对流体表面各处进行世界坐标重建，并根据位置偏微分求解梯度用于重构流体表面法线。流体表面法线用于光照计算，采用RTT技术可以将流体表面法线数据存储至高精度纹理。4）光照计算。真实流体在光照下涉及到的光学现象有反射、折射和透射等，对流体进行光照计算是真实感绘制必不可少的环节。结合深度、法线和厚度进行流体光照计算后，即可得到最终的流体真实感绘制场景。4　实验结果及分析本文实验依托的硬件平台为Intel Core i7-9750H CPU，Nvidia GeForce GTX 1660Ti GPU，16 GB内存。在实现过程中，使用C#、HLSL（high level shader language）、C++混合编程，集成开发环境为Visual Studio 2019，操作系统为Windows 10 × 64，图形渲染基于DirectX 11完成。4.1　稳定性分析图8为采用本文方法对流体冲击坚固墙壁场景进行模拟的结果。设定该场景的固体碎块数上限（maximum of solid chunks，MSC）为120个，流体粒子数上限 （maximum of fluid particles，MFP）为32 k，墙壁断裂抗性因子（fracture resistance coefficient，FRC）为0.97。在图8所示场景中，流体发射器向墙壁发射具有一定初速度的流体，墙壁受到涌来的流体的冲击后发生破碎。从图中结果可以看出，本文提出的固体破碎模拟方法能够稳定地模拟固体受到流体冲击而发生的破碎现象。在受冲击更严重的区域，产生的碎块更细小和密集，这与现实生活中的实际规律相符合，具有良好的真实性。10.11834/jig.220176.F008图8实验场景1：流体冲击坚固墙壁Fig.8Experimental scene 1： fluid impacts solid wall图9为左右两个观察视角下，流体冲击固体脆板场景的模拟结果。该场景设定MSC为160个，MFP为52 k，FRC为0.45。与图8中的场景类似，流体以一定初速度涌向固体脆板并对其进行冲击，固体脆板在流体作用下发生破碎。从图9可以看出，相比实验场景一，由于固体的MRC降低，固体受到较小的冲击力即会发生破碎。此外，固体受流体冲击后产生的碎块能够继续与流体进行交互，在流体作用下被冲散，符合物理规律。10.11834/jig.220176.F009图9实验场景2：流体冲击固体脆板Fig.9Experimental scene 2： fluid impacts solid fragile plate （（a） left side； （b） right side）从实验场景1和实验场景2的模拟结果可以看出，本文提出的模拟框架能够稳定地模拟流体运动、固流交互现象，能够求解不同FRC固体的破碎过程，破碎产生的碎块能够与流体进行良好的交互。4.2　真实感效果分析图10为流体冲击固体脆板场景的真实感模拟结果。经过SSFR方法的渲染，流体自身颜色、反射以及折射效果真实感良好。10.11834/jig.220176.F010图10实验场景2的真实感渲染结果Fig.10Rendering results of experimental scene 2（（a） left side； （b） right side）图11为流体冲击河谷大坝场景的真实感模拟结果，该场景设定MSC为180个、MFP为650 k，FRC为0.97。从结果可以看出，流体与地形边界交互良好，流体沿河谷倾泻后冲击河谷中的大坝，大坝底部受巨大冲击导致断裂最终倒塌，符合真实物理规律。大坝模型与地形、碎块与碎块以及碎块与流体之间的动力学交互效果良好。真实感渲染场景中水体漫反射、镜面反射和折射等渲染效果良好。10.11834/jig.220176.F011图11实验场景3：流体冲毁河谷大坝Fig.11Experimental scene 3： fluid washes out valley dam图12为本文提出的固体破碎模拟方法与赵建旺等人（2019）方法的模拟效果对比。从结果可以看出，本文方法在固体受冲击更严重的地方碎块更细小和密集，更加符合真实世界物理规律。10.11834/jig.220176.F012图12本文方法与赵建旺等人（2019）方法的效果对比Fig.12The effect comparison with Zhao et al. （2019）（（a） ours； （b） Zhao et al. （2019））图13为洪水冲毁房屋场景的真实感模拟结果。该场景设定MSC为900个、MFP为3 360 k，FRC分别为左上水泥建筑0.97、左下砖瓦仓库0.72、右上砖瓦房屋0.52、右下砖瓦小屋0.33。从结果可以看出，首先受到冲击的房屋从受冲击部位发生破碎，碎块随流体继续运动，房屋逐渐垮塌。不同建筑受流体冲击的程度、位置和方向不同，因而具有不同的破碎特征。由于预设定的不同的建筑FRC不同，各个建筑的倒塌程度有所区别，高FRC的建筑能够承受更高强度的流体冲击，垮塌发生时传入固体系统的能量更大；而低FRC的建筑则在较低强度的冲击下就能够发生垮塌，良好地反映出水泥与砖瓦结构建筑的区别。此外，场景中还未受到流体冲击的房屋保持完好无损。该实验场景中各房屋与流体交互良好，破碎过程符合物理规律，具有较好的真实感。10.11834/jig.220176.F013图13实验场景4：洪水冲毁房屋Fig.13Experimental scene 4： floods destroyed houses4.3　时间效率分析上述实验场景在2 400帧（最大粒子数）时的粒子数、碎块数以及分别在CPU和GPU端的帧率 （frame per second，FPS） 统计如表1所示。10.11834/jig.220176.T001表1各场景第2 400帧的FPSTable 1The FPS of the 2 400th frame of each scene场景MFP/kMSCFPS/（帧/s）CPUGPU13212023.992.125216013.971.936501800.221.243 3609000.000 72.1注：加粗字体表示具有更优时间效率的结果。表2为基于实验场景1—场景3进行的关于本文提出的模拟方法与赵建旺等人 （2019） 采用的方法在时间效率方面的对比结果。从表1和表2可以看出，本文方法经过GPU并行优化后，在650 k粒子规模下达到了实时级的模拟帧率，在3 360 k粒子规模下的帧率仍能够满足交互需求。10.11834/jig.220176.T002表2本文方法与赵建旺等人（2019）方法的FPS对比Table 2Comparison of the FPS between the method of Zhao et al. （2019） and ours方法粒子数32 k52 k650 k赵建旺等人（2019）9.8041.7010.041本文92.171.921.2注：加粗字体表示具有更优时间效率的结果。帧/s前述4个场景从第400帧到第2 400帧（粒子数达到MFP）的逐帧统计如图14所示，每个场景的统计数据包括逐帧GPU帧率、CPU帧率以及GPU加速比。从图14可以看出，本文提出的并行模拟框架能够在百万级粒子规模下达到可交互级别的模拟帧率，具有明显的时间性能优势。从结果可以看出，随着粒子规模的增大，GPU并行加速比也越高。经过本文计算，前述4个实验场景从第400帧到第2 400帧的平均GPU加速比分别达到了3.63、4.74、108.28和3 047.96。这主要是由于在粒子规模越大的场景中，非模拟核心计算的占比越小，并行收益也就越高。10.11834/jig.220176.F014图14各场景的第400—2 400帧的帧率Fig.14The FPS of each scene from the 400th to 2 400th frame （（a） scene 1； （b） scene 2； （c） scene 3； （d） scene 4）5　结论为模拟固流交互过程中的固体破碎现象，针对固流交互模拟研究较少关注固流交互中固体破碎现象的模拟、相关研究或真实感不足或无法满足实时模拟要求、固体破碎模拟方法与SPH固流交互模拟耦合性差等问题，本文以SPH统一粒子框架为基础，提出一种物理和几何混合的固体破碎模拟方法。该方法结合固体破碎的物理学理论与SPH统一粒子框架特性，以基于物理的模型计算固体破碎的发生情况，采用基于几何的方法生成固体碎块，解决了固体破碎模拟方法与SPH固流交互模拟框架耦合性差的问题，弥补了物理方法计算开销大、几何方法真实感不足的缺点。实验结果表明该方法能够稳定、高效地模拟固流交互引发的固体破碎现象，破碎效果符合现实世界物理规律，真实感良好。为进一步提高模拟的时间效率，本文将模拟流程进行并行处理后加载至GPU进行并行计算，并采用SSFR方法渲染流体以降低CPU与GPU间的通信开销。实验结果表明模拟系统充分利用了GPU的并行计算能力，在650k粒子规模下达到了21.2帧/s的模拟性能。虽然本文实现了兼顾真实感和实时性的模拟，但在诸多方面仍有待优化。例如可以采用更高效的流体模拟方法提高模拟效率、考虑固体的多阶段破碎模拟、引入动态粒度优化实现多粒度自适应模拟以及将NPS空间优化方法与固体碰撞检测结合等，以实现更高效率、更具真实感的模拟。