一种干燥颗粒流实时仿真与交互方法与流程

本发明属于计算机图形学与人机交互领域，具体涉及一种基于浅沙方程的干燥颗粒流实时仿真与交互方法。

背景技术：

颗粒介质是地球上继水之后人类第二大最常接触的物质。严格地说，颗粒状物质是指包含大量不连续的固体颗粒的集合，其间隙填充有流体或气体。一方面它表现出类似流体的行为，另一方面其也具备一定固体的特性。沙子一方面可以像流体一样流动形成任意形状，另一方面则可以在静止的时候像固体形成不同高度的沙堆。尽管研究人员已经对颗粒流的物理学开展了大量研究，其运动机理依然没有被完全摸清。实验表明，颗粒流的动力行为不仅取决于颗粒质量，还取决于其他因素，如颗粒大小，颗粒形状，含水饱和度等。如何高效模拟可能包含数百万乃至数亿颗粒的颗粒流，长期以来一直被认为是数值模拟的难题。

本发明主要针对一种特殊的粒状介质——干燥颗粒流进行实时仿真。如果颗粒大小足够大(如尺寸大于100μm)且颗粒间的粘度足够小，则本发明可以忽略所有微观的力(如静电，空气动力学和毛细管力)对干燥颗粒流的动力学进行简化。此时，本发明可以忽略凝聚力，从而仅考虑动量转换和能量耗散来建模干燥颗粒流的动力学。在连续介质力学中，干粒状流近似的采用单侧不可压缩流体模型[rahulnarain,abhinavgolas,andmingclin.2010.free-flowinggranularmaterialswithtwo-waysolidcoupling.acmtransactionsongraphics(tog)29,6,2010]或基于drucker-prager的弹塑性材料模型[gergelykl′ar,theodoregast,andrepradhana,chuyuanfu,craigschroeder,chenfanfujiang,andjosephteran.2016.drucker-pragerelastoplasticityforsandanimation.acmtrans.graph.35,4,2016]进行模拟。然而这些方法只适用于三维场景。一旦颗粒流场景较大，其仿真的效率将会受到严重影响，无法满足实时应用对于仿真效率的要求。

技术实现要素：

本发明主要针对传统方法仿真干燥颗粒流效率低的问题，提供一种干燥颗粒流实时仿真与交互方法。本发明利用浅沙方程对干燥颗粒流动力学进行建模，以及通过降维技术将三维运动控制方程转化为二维运动控制方程，从而极大的降低了仿真计算量，同时通过采用gpu对算法进行并行求解，使得大规模干燥颗粒流的实时仿真成为可能。

具体来说，本发明的技术方案如下：

一种基于浅沙方程的干颗粒流实时仿真与交互方法，其具体步骤包括：

(1)针对干燥颗粒流的特性，提出浅沙方程对干燥颗粒流动力学进行建模，得到干燥颗粒流动力学模型；

(2)利用二维高度场对得到的干燥颗粒流动力学模型进行离散化，得到高度场每个网格单元存储干颗粒流的深度及水平运动速度；

(3)针对步骤(2)得到的高度场，对距离高度场表面设定深度h以内的区域进一步采用三维均匀网格进行离散化，其中每个网格单元包含0个或者1个粒子用于存储干燥颗粒的位置；

(4)针对步骤(2)得到的高度场施加压强力，并更新得到新的速度；

(5)针对步骤(4)得到的高度场施加摩擦力，并更新得到新的速度；

(6)利用步骤(5)得到的最终速度，对颗粒流沿其水平速度方向进行平流操作；

(7)利用步骤(2)中的水平运动速度，对步骤(3)中的该设定区域内离散化后各网格单元中的干燥颗粒沿颗粒流的水平速度方向进行平流操作；

(8)对步骤(7)中运动之后的干燥颗粒的垂直坐标进行变换，使得干燥颗粒沿着高度场的表面运动；

(9)如存在交互物体，利用交互物体当前的运动速度对高度场的运动状态进行更新；

(10)仿真运行过程中重复步骤(4)到(9)之间的操作并更新时间，直到达到终止时间或者其他仿真终止条件。

进一步的，步骤(1)结合干燥颗粒流的特性，提出了两个假设：一、干燥颗粒流仿真区域的垂直尺度要远小于水平尺度；二、整个高度场沿垂直方向被分成两层——稀疏层和稠密层，其中稀疏层为靠近高度场表面的部分，其运动遵从浅水方程所述的流体动力学；稠密层为靠近高度场底部的部分，其中内部压强保持恒定，且不受外力的影响。

进一步的，高度场的干燥颗粒流动力学模型可以利用如下的浅沙方程进行描述：

其中h是高度场高度，其沿垂直方向被分解为稀疏层和稠密层两部分(如图2所示)，h表示高度场的稀疏层的高度，s是高度场表面的垂直坐标，u和v分别表示干燥颗粒流沿x和y方向的运动速度，g是重力加速度，μ是摩擦系数，t表示时间，和分别代表水平运动速度归一向量的分量。

进一步的，干燥颗粒流动力学模型所描述的干燥颗粒流可以利用二维高度场进行离散，且整个动力学过程完全利用高度场来完成。

进一步的，步骤(3)中新生成粒子的坐标可以利用如下的公式表示：

其中s、t、l表示三维网格的沿x、y、z方向的索引，δd表示网格间距，rand(0,1)表示生成范围为(0,1)的随机函数。

进一步的，步骤(4)中压强力的计算方式如下：

其中表示索引为(i,j)的高度场网格单元在t时刻的干颗粒流深度，表示索引为(i,j)的高度场网格单元所对应的稀疏层的高度，δx、δy分别表示x和y方向的网格间距，是防止稀疏层的高度超过高度场的总高度。

进一步的，步骤(5)中摩擦力遵从动能最大耗散原理：

其中n是沿着干燥颗粒流水平速度的归一化单位向量，和中的下标(i,j)表示高度场中单元格索引，δt表示时间步长；从而摩擦力的计算公式表述如下：

其中，(hu,hv)″i,j表示索引为(i,j)的网格单元施加压强力之后得到的水平平均速度向量，n″i,j表示索引为(i,j)的网格单元施加压强力之后得到的归一化水平平均速度向量，μmax表示最大摩擦系数。

进一步的，步骤(8)中粒子通过如下公式进行变换：

其中xp＝(xp,yp,zp)表示粒子的当前坐标，zp表示粒子的z轴坐标，si(xp)表示高度场表面坐标s在粒子xp处进行二维线性插值的得到的值，表示粒子变换后的坐标。

进一步的，步骤(9)当交互物体与颗粒流存在相交区域时，本发明预先对交互物体采用有向距离场进行离散化，同时假定高度场网格单元与交互物体相交部分以与交互物体相同的速度移动，则根据动量守恒定律，二维高度场上网格单元的新速度表示为

其中vr表示交互物体的速度，vi,j表示网格单元的速度，hi,j表示二维高度场网格单元的沙层深度，d表示网格单元顶部到交互物体表面的最近距离。

至此，完成了整个算法。

与现有技术相比，本发明具有如下的优点：

(1)本发明提出浅沙方程用于建模干燥颗粒流的动力学，其有效地将传统方法所基于的三维颗粒流运动控制方程降维到二维空间进行仿真，从而极大的降低了计算开销；

(2)本发明中每一个步骤都具有高并行性，因而适用于现有gpu架构，可以利用硬件加速对大规模干燥颗粒流进行实时模拟；

(3)本发明采用有向距离场对交互物体进行预离散化并利用gpu进行了加速，支持不同几何形态的交互物体与颗粒流的实时交互。

附图说明

图1为仿真算法流程图。

图2为沙床示意图，其由稀疏层和稠密层两部分构成。

图3为交互物体与沙床交互示意图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实现和附图对本发明做详细说明,但不构成对本发明的限制。

本发明方法的硬件平台采用了型号为inteli7-6700hq四核cpu，主频为2.6ghz，nvidiageforcegtx1060显卡，现存6gb。系统程序采用c++编写，其中对于并行计算部分采用了cuda语言进行加速，并借助microsoftvisualstudio2015对程序进行编译执行，开发过程中采用了opengl、freeglut等开源库。

本发明方法的流程图如图1所示，整个流程分为两部分。最底层是一个二维高度场，用于求解干燥颗粒流的动力学问题，而顶层是一个三维均匀网格，用于追踪颗粒流的表面及展示颗粒的运动。在每个仿真步开始之初，本发明首先通过计算压强力和摩擦力来更新存储在高度场上的运动状态。新的运动状态一方面用于更新颗粒流的深度，另一方面则用于驱动粒子沿着水平方向的运动。然后，本发明根据粒子当前所处的网格单元解决位置冲突。一旦发现多于一个粒子同时进入网格单元时，本发明只保留最后进入的粒子。最后，本发明为不包含任何粒子的网格单元重新生成粒子，并将所有粒子的垂直坐标进行转换以匹配二维高度场的运动。接下来本发明将详细介绍每一个步骤的具体实现方案。

1、理论基础。结合干燥颗粒流的运动学特性，本发明提出了两个假设：一、干燥颗粒流仿真区域的垂直尺度要远小于水平尺度；二、整个沙床沿垂直方向被分成两层——稀疏层和稠密层，其中稀疏层为靠近沙床表面的部分，其运动遵从浅水方程所述的流体动力学；稠密层为靠近底部的部分，其中内部压强保持恒定，且不受外力的影响。图2根据上述假设给出了一个沙床的示意图。假设沙床的密度为ρ且沿着z轴方向密度恒定，稀疏层的高度为则沙床任意点的压强可以表示为：

其中b表示地面坐标，s表示沙床表面的z坐标，pa表示大气压强，g表示重力加速度。为了简便起见，以下讨论忽略x、y及z坐标参数，则压强力的计算可以通过沿z轴方向的积分得到，即

在计算摩擦力f^f的过程中，本发明假设f^f与f^p之间存在如下的关系式：

‖f^f‖＝μ‖f^p‖cot(θ)

其中μ表示摩擦系数，θ表示摩擦角。因此摩擦力f^f的计算公式可以表示为

n表示流体的运动方向，上式摩擦力的方向始终与颗粒流的相对运动速度相反。

综上，干燥颗粒流的动力学可以利用如下的浅沙方程进行统一建模：

其中h是沙床高度，表示稀疏层的高度，s是沙床表面的垂直坐标，u和v分别表示水平运动速度的两个分量，g是重力加速度，μ是摩擦系数，t表示时间。

2、离散化。这里包含对两部分数据的离散化，第一部分为动力学计算网格的离散化，采用均匀的二维高度场对浅沙方程中包括h、u、v在内的标量场进行离散化；第二部分为表面粒子的离散化，采用均匀三维网格对靠近高度场表面的区域进行离散，其中每个网格单元存储最多一个粒子，用于记录粒子当前所处的位置。三维网格与二维高度场的顶部齐平，主要用于颗粒流的可视化。

3、更新速度。整个速度更新的过程采用预测校正的方案来完成，即首先在不考虑摩擦力的前提下计算压强力，并更新流体的运动状态。通过采用中心差分的形式，压强力的离散形式可以表示为：

其中是t时刻(i+1,j)网格单元的高度场表面z坐标，是t时刻(i-1,j)网格单元的高度场表面z坐标，是t时刻(i,j+1)网格单元的高度场表面z坐标，是t时刻(i,j-1)网格单元的高度场表面z坐标，δx、δy分别表示x和y方向的网格间距，是防止稀疏层的高度超过总沙床的高度，g表示重力加速度。在施加压强力得到临时速度之后，本发明进一步考虑实际所需摩擦力的大小。

根据连续介质力学最大动能耗散原理，实际摩擦应力的大小应该能保证流体最终状态的总动能最小。通过这一原理，每个网格单元的实际摩擦力可以通过求解以下优化问题来获取，即：

其中n是沿着干燥颗粒流水平速度的归一化单位向量，δt表示时间步长；

从而摩擦力的计算公式可以表述如下：

其中μmax表示最大摩擦系数。

4、平流。平流操作主要用于更新二维高度场的沙床高度及速度，具体表达式如下：

其中表示网格单元i+1,j和网格单元i,j之间的通量。具体计算方式可以参考[alexanderkurganov,guerganapetrova.asecond-orderwell-balancedpositivitypreservingcentral-upwindschemeforthesaint-venantsystem.communicationsinmathematicalsciences5,1(2007),133–160]。

5、移动粒子。在得到网格新的速度之后，三维网格中粒子会根据二维高度场中沙床的水平速度进行运动，速度更新的公式如下

x'p＝+δt[ui(xp),vi(xp),0]^t

其中ui(xp),vi(xp)表示粒子p所处的当前位置xp所对应的沙床水平速度。

6、解决粒子冲突。在得到了新的位置之后x'p之后，本发明需要根据其位置确定其所在的网格单元，新的网格单元的索引可以根据以下公式计算得出

公式表示粒子运动只沿着水平方向进行，因而其垂直方向的下标始终保持恒定。一旦发现多于1个粒子进入同一网格单元时，本发明只保留最后进入的粒子，从而保证每个网格单元依然只存储一个粒子。这样做的目的是为了保证计算的并行性。

7、生成新的粒子。当发现某网格单元中不存在任何粒子时，本发明将根据以下公式生成一个新的粒子，新生成粒子的坐标可以利用如下的公式表示：

其中s、t、l表示三维网格的沿x、y、z方向的索引，δd表示网格间距，rand(0,1)表示生成范围为(0,1)的随机函数。

8、粒子坐标变换。由于粒子的运动只是沿着二维水平方向进行，而实际沙床的形态是三维的。为了真实表现颗粒流的三维运动，本发明对粒子的坐标进行如下的变换：

其中xp＝(xp,yp,zp)表示粒子变换之前的坐标，si(xp)表示高度场表面坐标s在xp处进行二维线性插值而得到的值。如此一来可以实现沙堆逼真的三维仿真效果。

9、交互。对于需要交互的场景，如图3所示，本发明首先将交互物体用有向距离场进行离散。当交互物体与颗粒流存在相交区域时，本发明假定网格单元与交互物体相交部分以与刚体相同的速度移动，则根据动量守恒定律，二维高度场上网格单元的新速度表示为

其中vr表示交互物体的速度，vi,j表示网格单元的速度，d表示网格单元顶部到交互物体表面的最近距离。

以上实施方示例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求所述为准。