吴蔚

摘  要：沙粒模拟是计算机图形学中重要的研究内容。文章提出一种基于CUDA的并行化沙粒模拟方法，使用PCISPH流体模拟的算法实现沙粒流动过程，并加入摩擦力模型保持沙粒稳定堆积的效果。算法利用多线程进行数值求解，充分挖掘了GPU并行计算的性能优势。实验表明，该算法可以实现高效的沙粒模拟，在动画、电影、工程等领域具有重要的应用价值。

关键词：沙粒模拟;SPH;并行计算;CUDA

中图分类号：TP391.7       文献标志码：A              文章编号：2095-2945（2020）25-0009-02

Abstract： Sand simulation is the key research content in computer graphics. A parallel sand simulation method based on CUDA is proposed， which use the PCISPH fluid simulation method to realize the free sand flow and adopt the friction model to keep in stable sandpile. The numerical solution is processed by multithread， which helps to fully excavate the performance advantage of GPU parallel computing. Experiments suggest that the method can realize the sand simulation efficiently， and bear important application value in the field of animation， movie， engineering， etc.

Keywords： sand simulation; SPH; parallel computing; CUDA

1 概述

沙粒等固体颗粒的模拟，是计算机图形学中重要的研究内容，被广泛应用在动画、电影、工程等领域。在现实中可以观察到，不稳定的沙堆可以像流体一样流动。但与流体不同的是，沙粒内部的能量可以迅速消散，当这种能量达到一个较低值时，沙堆便会保持静止的堆积状态。之所以会表现出这样复杂的特性，正是由于沙粒颗粒之间相互作用的结果。对于沙粒所表现出的流动的效果，可以采用流体模拟的方法。而沙粒最终会趋于静止，却是由于沙粒之间的摩擦力抵消掉了重力的影响。因此沙粒模拟可以采用流体模拟的方法中加入粒子间摩擦力实现。

现今对于流体的模拟通常有两类，即基于网格的欧拉方法和基于粒子的拉格朗日方法。欧拉方法将模拟区域划分为网格，基于固定的网格点研究各时刻物理量在流场中的变化规律，其计算结果具有较高的数值精度;而拉格朗日方法则跟随质点运动，追踪它们的物理量随时间变化的规律，其更适用于需要表现表面细节或发生大变形的流体模拟应用中。SPH（Smoothed Particle Hydrodynamics）方法是一种典型的基于拉格朗日视角的流體模拟方法，近年来成为最具前景的流体模拟方法之一。该方法将连续的流体离散成粒子，每个粒子携带一定的物理量，比如密度、压强等，通过追踪每个粒子的物理量变化从而表现整个流体的运动过程。在沙粒模拟中，单纯使用SPH流体模拟方法并不能保持沙粒静止堆积的状态。因此沙粒粒子之间除了受到压强力和速度差导致的黏性力等外力，还需要加入摩擦力的影响。Alduan等人[1]提出了一种基于SPH的沙粒模拟方法，使用结合单侧不可压缩性的PCISPH方法模拟沙粒的流动，并加入摩擦力模型实现了沙粒的堆积效果。

通常来说，基于粒子的数值计算方法比基于网格的方法更加高效，然而SPH方法中模拟介质通常被离散成大量的粒子，每个粒子的物理量求解又依赖于对邻域粒子的搜索以完成插值计算，采用传统的CPU进行迭代式地求解很难达到实时性的要求。因此有研究者提出使用GPU对SPH求解进行加速的方法。比如Harada等[2]将SPH计算过程全部在GPU上实现，并使用桶纹理结构实现了高效的邻域搜索算法。Junior等[3]不仅实现了GPU上的数值求解，还使用空间哈希函数改进了邻域搜索算法，实现了无边界的流体模拟。

本文基于CUDA架构实现了并行化的沙粒模拟过程，将Alduan等人[1]的PCISPH数值求解和摩擦力模型完全基于GPU进行并行化处理，避免了CPU和GPU之间不必要的数据传输，显著提升了沙粒模拟的运算效率。

2 基于CUDA架构实现沙粒模拟

2.1 PCISPH方法

采用标准SPH方法模拟的沙粒具有可压缩性，严重影响模拟的效果。本文使用PCISPH方法实现沙粒的不可压缩模拟，每个时间步的算法步骤如下：

（1）首先针对每一个粒子搜索一定范围内的邻域粒子，用于后续各种物理量的计算;

（2）使用SPH插值方法计算出所有非压强力，比如粒子间粘性力、重力、所受外力等。同时将压强和压强力初始化为0;

（3）使用预测-矫正策略循环更新压强和压强力的值。迭代之初，使用粒子所受合力预测粒子临时速度和临时位置，并基于预测的位置信息预估出粒子的密度信息。该密度值用于矫正压强和压强力，并作为是否继续下一迭代的依据。若当前预测的密度值过大，说明粒子之间过于紧密，需要按照如下公式增大粒子所受压强Pi：

Pi+=δρerr

其中δ为预计算的常量，ρerr为与预测的密度相关的变量。所得压强值用于计算出新的压强力，重新进入循环更新粒子的速度和位置。经过若干次迭代，粒子间距便可满足预先设定的压缩率要求;

（4）基于粒子所有更新后的力進行时间积分，得到下一时间步的粒子速度和位置。

上面步骤中第（3）步包含若干次迭代计算，是整个算法最耗时的部分。这种迭代式的方式使得每个时间步的计算比标准SPH方法花费更多的时间。然而该方法允许设置较大的时间步长，这反而极大地提升了沙粒模拟整体的运行效率。

2.2 添加摩擦力模型

与单纯的流体模拟不同，沙粒模拟还需要添加粒子间摩擦力，从而阻止粒子间的相对滑动。摩擦力模型基于最大耗散原则实现，这又可通过最小化应变速率来表达。为此，使用反推法基于粒子间相对速度得到粒子应变速率值，然后迭代地计算出摩擦应力张量值s，从而计算出粒子所受摩擦力大小：

其中θ为预定义的静止角，用来决定沙粒稳定堆积的坡度。由于该摩擦力模型本质上也是基于预测-矫正策略，因此可以很容易地集成到前述PCISPH方法中。

2.3 基于CUDA的算法实现

SPH作为一种典型的拉格朗日方法，存在大量粒子间无关的数值计算过程，因此可以充分借助GPU的性能优势实现高度并行化。算法初始化时，所有沙粒粒子的属性数据及边界信息从CPU端内存拷贝到GPU端显存。后续的数值计算完全在GPU端进行，避免了CPU和GPU之间不必要的数据传输。

沙粒模拟的每一帧都需要计算出所有粒子更新后的属性值，比如密度、压强、应变速率等，以及各种受力，比如压强力、摩擦力等。在CUDA编程模型中，可以并行执行的代码称为核函数。因此针对沙粒模拟中的邻域搜索及各种物理量的计算，本文使用核函数创建的轻量级线程与每个粒子一一对应，实现对所有粒子的并行化处理。

2.4 实例验证

为了验证上述算法，设计了沙粒模拟实验，如图1所示。初始时由沙粒颗粒组成的长方体沙堆位于长方体盒子包围下的角落位置。相比于流体模拟，沙粒模拟由于粒子间摩擦力的作用，在流动过程中逐渐呈现出明显的堆积效果，验证了本文沙粒模型的准确性。为了对本文基于CUDA的并行算法效果进行验证，针对沙粒模拟场景分别采用CPU串行和GPU并行的方式实现，每帧平均运行时间对比如表1所示。可以看到，相比于串行方式，并行化使得算法在运行效率上达到了8倍的加速比，充分验证了本文并行化算法的有效性。

3 结束语

本文基于CUDA框架实现了并行化的沙粒模拟过程，将PCISPH流体模拟的方法和摩擦力模型完全基于GPU进行数值求解，充分挖掘了GPU并行计算的性能优势，提升了沙粒模拟的运行效率。今后的研究将在本文研究的基础上继续探索CUDA优化技术在大规模沙粒模拟场景中的效果。

参考文献：

[1]Alduán I， Otaduy M A. SPH granular flow with friction and cohesion[C]//Proceedings of the 2011 ACM SIGGRAPH/Eurographics symposium on computer animation， 2011： 25-32.

[2]Harada T ， Koshizuka S ， Kawaguchi Y . Smoothed Particle Hydrodynamics on GPUs[J]. Proc.computer Graphics Int.rio De Janeiro Brazil May Jun， 2007， 4（4）：671-691.

[3]Junior J R S， Joselli M， Zamith M， et al. An architecture for real time fluid simulation using multiple GPUs[J]. XI SBGames， Brasília， 2012：8.