CN110992456A - 一种基于位置动力学的雪崩模拟方法 - Google Patents

Abstract

本发明公开了一种基于位置动力学的雪崩场景模拟方法。其步骤为：1)对所有顶点粒子进行了属性初始化，根据外力预测更新当前速度，根据时间步长及当前速度更新预测位置；2)使用符号距离场进行碰撞检测，更新位置及速度；3)迭代求解约束直到约束条件收敛，根据校正位置更新速度；本发明解决了在计算机图形学中雪崩灾害场景模拟缺失的问题，能够快速模拟雪崩发生过程；通过凝聚力、摩擦力约束方程，模拟雪流及雪块等不同形态的颗粒流运动过程；通过符号距离场实现了与石块和地面等刚体的耦合现象；使用并行计算框架提升了模拟效率。

Description

一种基于位置动力学的雪崩模拟方法

技术领域

本发明属于流体仿真领域，尤其涉及一种基于位置动力学的雪崩模拟方法。

背景技术

基于物理的自然场景仿真一直是计算机图形学领域研究的一大热点。其内容主要可分为用于游戏娱乐的实时仿真以及用于影视特效的离线仿真两大分支。对于仿真效率以及仿真效果的评判角度而言，前者更注重计算的效率，即实时性；而后者更为注重视觉效果，即准确感。但是即便如此，在影视界的准确感特效制作中，若能够提高计算绘制效率、减少艺术创作者的等待时间也是十分关键的，这有助于缩短制作周期、大大降低人工制作成本，从而提高市场竞争力。

雪崩是一种常见的自然灾害，是指斜坡上的雪块在摩擦力和凝聚力不足的情况下，发生整体移动而产生的灾难性地质现象。雪崩包含广泛的地面运动，包括干雪崩、湿雪崩等。重力是雪崩发生的主要驱动力，但影响边坡稳定性的其他因素也会导致雪层破坏从而发生雪崩。在许多现象情况下，雪崩是由某一特定事件引起的，如声音，爆破，融化等，但是确定雪崩的原因往往是十分困难的。

在计算机图形学中基于物理的动画领域的研究与寻找模拟物理现象(如刚体，可变形物体或流体流动)的新方法有关。与侧重于准确性的计算科学相反，此处的主要问题是稳定性，鲁棒性和速度，而结果应在视觉上保持合理。因此，计算科学中的现有方法不能一对一地采用。实际上，在计算机图形学中进行基于物理的仿真研究的主要理由是提出专门针对特定领域需求的方法。该方法属于这一类。

模拟动态对象的传统方法是使用力。在每个时间步开始时，都会累积内部和外部力。内部力的示例包括可变形物体中的弹力或流体中的粘度和压力。重力和碰撞力是外部力的例子。牛顿的第二运动定律与通过质量的加速度有关。根据顶点的密度或集中质量，将力转换为加速度。然后，可以使用任何时间积分方案首先根据加速度计算速度，然后根据速度计算位置。一些方法使用冲力而不是力量来控制动画。由于冲力直接改变速度，因此可以跳过一级积分。

在计算机图形学中，尤其是在计算机游戏中，经常需要直接控制对象或网格顶点的位置。用户可能希望附加平均运动学对象，或确保顶点始终保持在碰撞对象之外。本发明提出的方法在此处直接作用于使此类操作容易的位置。此外，基于位置的方法可以直接控制集成，从而避免与显式集成相关的过冲和能量获取问题。因此，基于位置的动力学的主要特征和优点是(1)基于位置的仿真可控制显式集成，并消除了典型的不稳定问题。(2)在仿真过程中，可以直接操纵顶点的位置和对象的一部分。(3)提出的公式允许处理基于位置的设置中的一般约束。(4)基于显式位置的求解器易于理解和实现。

发明内容

本发明的目的在于克服图形学基于物理的自然场景仿真与绘制领域对雪崩现象的高效仿真的方法缺失，提供一种基于位置动力学的雪崩模拟方法，其针对雪崩场景特点实现了雪流不同参数化模型，并且使用多种技术，实现雪崩模拟结果。

为了实现上述目的，本发明采用如下技术方案：

一种基于位置动力学的雪崩场景模拟方法，包括以下步骤：

1)将雪山山体雪层离散化为粒子，对所有粒子进行属性初始化；

2)根据粒子在当前t时刻的准确位置和受力情况，引入速度阻尼，得到t+Δt时刻所有粒子的预测速度和预测位置；

3)对所有粒子进行碰撞检测，若判断为粒子与刚体发生碰撞且粒子克服摩擦力发生运动，则更新t+Δt时刻粒子的预测位置和预测速度，否则保留步骤2)得到的t+Δt时刻粒子的预测位置和预测速度；

4)以最终更新后的t+Δt时刻粒子的预测位置和预测速度为参数，遍历所有约束条件，对预测位置进行修正，得到修正后的t+Δt时刻的所有粒子的准确位置和准确速度；

5)输出雪山山体雪层在t+Δt时刻的准确位置处的图像；

6)以t+Δt时刻所有粒子的准确位置作为起始位置，重复步骤2)-步骤4)，依次输出雪山山体雪层的每一帧图像，得到雪山山体雪层在雪崩场景下的模拟动画。

进一步的，所述的步骤1)具体为：

将雪山山体雪层离散化为粒子，设置所有粒子初始化的位置动力学参数，公式如下：

其中，x_i和v_i表示第i个粒子的准确位置和准确速度，

和

表示第i个粒子的初始化位置和速度，N表示粒子的个数，m_i表示第i个粒子的质量，w_i表示第i个粒子的权重参数。

进一步的，所述的步骤2)具体为：

2.1)定义已知的当前t时刻的粒子的准确位置

和速度

2.2)根据外力作用得到粒子在t+Δt时刻的预测速度：

其中，Δt表示时间步长，f_ext(·)表示所有外力的和的向量；

2.3)引入速度阻尼，进一步更新得到t+Δt时刻所有粒子的预测速度

根据

得到t+Δt时刻所有粒子的预测位置

计算公式如下：

进一步的，所述的步骤3)具体为：

3.1)对所有粒子进行碰撞检测，通过f_SDF函数计算粒子到刚体的符号距离函数的最近距离f，若最近距离为负值，则说明预测位置在刚体内部，粒子与刚体发生碰撞，否则未发生碰撞；最近距离f的计算公式如下：

3.2)若判断为粒子与刚体未发生碰撞，则保留步骤2)得到的t+Δt时刻该粒子的预测速度和预测位置；

若判断为粒子与刚体发生碰撞，首先计算粒子相对于碰撞物体的相对速度：

v_rel＝v_pOri-v_coll

其中，v_rel为粒子相对于碰撞物体的速度，v_pOri为粒子的速度，v_coll为碰撞物体的速度；

根据相对速度计算法向速度v_n，计算公式为：

v_n＝v_rel·n

其中n为碰撞物体的法向方向，由水平集φ的梯度方向来确定:

计算粒子切向速度v_t，计算公式为:

v_t＝v_rel-nv_n

F_f＝μ_fv_n

其中，μ_f为刚体表面的摩擦系数，F_f为静摩擦力；若‖v_t‖≤F_f，该粒子不发生位移，保留步骤2)得到的t+Δt时刻该粒子的预测速度和预测位置；否则，该粒子克服摩擦力运动，更新 t+Δt时刻该粒子的预测速度和预测位置：

其中，‖v_t‖表示v_t的值，

和

分别表示更新后的t+Δt时刻该粒子的预测速度和预测位置；

根据动量守恒，同时计算出粒子对于刚体的平移速度以及角速度：

v_ref＝v_line+ω_angle×(x_i-x_mass)

其中v_line为粒子对刚体的平移速度，ω_angle是粒子对刚体作用的角速度，x_i表示第i个粒子的世界坐标向量，x_mass表示由N个粒子组成的整个物体的质心坐标。

进一步的，所述的步骤4)具体为：

4.1)以更新后的t+Δt时刻粒子的预测位置和预测速度为参数，遍历所有约束条件，所述的约束条件包括阻尼约束、摩擦力约束和凝聚力约束；对预测位置进行修正，计算公式如下：

M＝diag(m₁，m₂，m₃)

其中

表示修正后的t+Δt时刻的第i个粒子的准确位置，

表示t+Δt时刻的第i个粒子的位置误差修正向量；M表示加权矩阵，m₁，m₂，m₃分别表示阻尼约束、摩擦力约束和凝聚力约束的加权项，

表示梯度计算函数，λ_i表示刚度系数；

和

分别表示阻尼约束、摩擦力约束和凝聚力约束的位置误差修正向量；

4.2)迭代求解阻尼约束：

所述阻尼约束表示为：

其中，Δv_i表示速度误差修正向量，k_damping为阻尼系数；

求解步骤为：

计算整体运动的位置p_cm以及整体运动的速度v_cm：

求解中间变量L以及I：

其中，r_i为3x3矩阵，

表示r_i的迹；求解后得到ω角速度：

ω＝I^-1L

更新当前速度：

进一步得到阻尼约束的位置误差修正向量

4.3)迭代求解摩擦力约束：

所述摩擦力约束表示为：

求解步骤为：

根据当前时间步长中粒子的相对切向位移来计算摩擦的位置增量：

其中，

和

表示可用于计算摩擦力约束的任意一对粒子，n表示法向量；

进一步计算得到摩擦力约束的位置误差修正向量

其中，μ_k，μ_s分别是动摩擦系数和静摩擦系数，ω_i和ω_j分别表示对第i个粒子和第j个粒子的加权系数；d表示滑动摩擦的距离；

4.4)迭代求解凝聚力约束：

所述凝聚力约束表示为：

求解步骤为：

计算作用在粒子上的附着力F_a：

F_a＝aA_c

其中，a表示附着力参数，A_c表示可用于计算凝聚力约束的任意一对粒子的接触面积，采用二次插值的接触面积，根据两个粒子之间的渗透关系得到接触面积A_c；

对于法向附着力，直接采用法向位移限制，并直接进行校正；对于切向附着力，通过牛顿运动定律可得：

进一步计算得到凝聚力约束的位置误差修正向量

其中，ω_i和ω_j分别表示对第i个粒子和第j个粒子的加权系数，t表示两粒子相对位置切线方向单位向量。

本发明具备的有益效果：

(1)使用基于位置的动力学，针对雪崩场景加入摩擦力与凝聚力模型，能够高效且稳定地模拟雪崩场景；

(2)采用水平集方法进行碰撞检测耦合，能够实现雪流与地面等物体交互，并提升碰撞检测及交互的效率；

(3)针对摩擦力约束相比其他方法，同时实现了静摩擦力和动摩擦力约束，针对凝聚力约束，使用基于物理的方式提出了切向和法向的凝聚力约束计算模型；

(4)基于并行计算框架，实现了雪崩的快速模拟。

附图说明

图1本发明的步骤展示；

图2本发明的不同摩擦力约束模拟稳定时对比图展示；

图3本发明的不同凝聚力约束模拟序列帧对比图展示；

图4本发明的雪崩不同阶段序列帧效果展示。

具体实施方式

一种基于位置动力学的雪崩场景模拟方法，图1为本发明的步骤展示：

步骤一、对所有顶点粒子进行了属性初始化，根据外力预测更新当前速度，根据时间步长及当前速度更新预测位置；

将雪山山体雪层离散化为粒子，设置所有粒子初始化的位置动力学参数，公式如下：

其中，x_i和v_i表示第i个粒子的准确位置和准确速度，

和

表示第i个粒子的初始化位置和速度，N表示粒子的个数，m_i表示第i个粒子的质量，w_i表示第i个粒子的权重参数；

定义已知的当前t时刻的粒子的准确位置

和速度

根据外力作用得到粒子在t+Δt时刻的预测速度：

其中，Δt表示时间步长，f_ext(·)表示所有外力的和的向量；

2.3)引入速度阻尼，进一步更新得到t+Δt时刻所有粒子的预测速度

根据

得到t+Δt时刻所有粒子的预测位置

计算公式如下：

步骤二、使用符号距离场进行碰撞检测，更新位置及速度；

下面先介绍已有的碰撞检测模拟方法：

任意拉格朗日-欧拉(Arbitrary Lagrangian-Eulerian)方法任意拉格朗日-欧拉方法结合了传统的拉格朗日和欧拉方法，允许边界网格或面片顶点的任意运动，从而有效跟踪物质边界的运动，在内部网格上，根据定义参数求解网格，在独立于物质实体同时使得网格不发生严重的畸变。

在数学领域，水平集方法是拥有n个实值未知数的实值函数，其公式为：

L_c(f)＝{(x₁，…，x_n)|f(x₁，…，x_n)}

当未知数的数量为3时，水平集被称为水平面，同样也可以被称为等值面。图形学领域应用比较多的水平集数据结构就是符号距离函数SDFs(Signed distance functions又可以称为符号距离场。

作为本发明的优选实施方式，对所有粒子进行碰撞检测，通过f_SDF函数计算粒子到刚体的符号距离函数的最近距离f，若最近距离为负值，则说明预测位置在刚体内部，粒子与刚体发生碰撞，否则未发生碰撞；

由于SDF的预计算是在刚体空间进行的，需要先对SDFs做刚体坐标到世界坐标的变换，使其处于同一坐标系，再进行接下来的操作，把离散后的SDFs以f_SDF来表示，通过纹理内存的性质，可以通过硬件加速的GPU插值办法，将SDFs进行插值，可以生成更为圆滑的碰撞效果，由于SDFs存在与等值面相互垂直的梯度，因此可以把梯度等同于曲面的法线，通过取值需要采样点位置上下、前后、左右三轴上的6个点来计算出曲面的法线方向，并且得益于纹理内存的局部性原则，需取值位置的附近纹理也同时被读入缓存中，因此访问这些纹理并不需要多少计算量，最近距离f的计算公式如下：

若判断为粒子与刚体未发生碰撞，则保留步骤2)得到的t+Δt时刻该粒子的预测速度和预测位置；

若判断为粒子与刚体发生碰撞，首先计算粒子相对于碰撞物体的相对速度：

v_rel＝v_pOri-v_coll

其中，v_rel为粒子相对于碰撞物体的速度，v_pOri为粒子的速度，v_coll为碰撞物体的速度；

根据相对速度计算法向速度v_n，计算公式为：

v_n＝v_rel·n

其中n为碰撞物体的法向方向，由水平集φ的梯度方向来确定:

计算粒子切向速度v_t，计算公式为:

v_t＝v_rel-nv_n

F_f＝μ_fv_n

其中，μ_f为刚体表面的摩擦系数，F_f为静摩擦力；根据库伦-摩尔第一定律(摩擦力与作用在摩擦面上的正应力成正比，与跟表面的接触面积无关。若‖v_t‖≤F_f，粒子的作用力不超过静摩擦力，该粒子不发生位移，保留步骤一得到的t+Δt时刻该粒子的预测速度和预测位置；否则，该粒子能够克服摩擦力运动，此时摩擦力表现为动摩擦，需更新t+Δt时刻该粒子的预测速度和预测位置：

其中，‖v_t‖表示v_t的值，

和

分别表示更新后的t+Δt时刻该粒子的预测速度和预测位置；

根据动量守恒，同时计算出粒子对于刚体的平移速度以及角速度：

v_ref＝v_line+ω_angle×(x-x_mass)

其中v_line为粒子对刚体的平移速度，ω_angle是粒子对刚体作用的角速度，x_i表示第i个粒子的世界坐标向量，x_mass表示由N个点组成的整个物体的质心坐标。

步骤三、迭代求解约束直到约束条件收敛，根据校正位置更新速度；

求解某个约束的过程中，对于一个单边或双边约束，可将其写为：

C_i(x+Δx)＝0，i＝1，...，N

C_j(x+Δx)≥0，j＝1，...，N

其中x＝[x₁，x₂，…，x_n]^T是粒子位置的向量。约束通过Gauss-Seidel迭代求解时，其中每个约束使用C_i围绕x的线性化顺序求解。

方程中，位置变化量Δx被约束为沿着约束梯度方向，并且由质量矩阵的逆M＝diag(m₁，...，m_n)加权:

结合前面两个公式，λ_i可以写为：

通常，在处理每个约束后更新位置，并在多次迭代后根据总约束增量计算速度变化

作为本发明的优选实施方式，以更新后的t+Δt时刻粒子的预测位置和预测速度为参数，遍历所有约束条件，所述的约束条件包括阻尼约束、摩擦力约束和凝聚力约束；

(1)阻尼约束

作为本发明的优选实施方式，迭代求解阻尼约束

计算整体运动的位置p_cm以及整体运动的速度v_cm：

求解中间变量L以及I：

其中，r_i为3x3矩阵，

表示r_i的迹；求解后得到ω角速度：

ω＝I^-1L

更新当前速度：

进一步得到阻尼约束的位置误差修正向量

通常情况下，k_damping的值为0-1之间，此处本发明采用k_damping＝0.1进行计算。

(2)摩擦力约束

一个物体在另一个物体表面上具有相对运动趋势时，所受到的阻碍物体相对运动趋势的力叫静摩擦力。两个相接触的物体做相对运动时发生的阻碍它们相对运动的现象，称为“动摩擦”。如果物体接触的相对位置是滑动的(相对切线速度不为零)，则摩擦力必须始终与相对切线速度相反。同样，摩擦力的大小与法向力成线性比例，由两个摩擦系数之一定标，分别为静摩擦系数和动摩擦系数μ_s以及μ_k。

摩擦力约束求解过程为首先解决穿透问题，应用投影距离约束解决：

C(x_i，x_j)＝|x_ij|-r≥0

解决投影距离约束后，计算切向位移，位置增量的计算通过时间步长确定，得到Δx_⊥。然后，根据位置增量的绝对值|Δx_⊥|，判断是否小于阈值|Δx_⊥|＜μ_sd，即可确定粒子的位置改变量Δx_i。

作为本发明的优选实施方式，迭代求解摩擦力约束

根据当前时间步长中粒子的相对切向位移来计算摩擦的位置增量：

其中，

和

表示可用于计算摩擦力约束的任意一对粒子，n表示法向量；

进一步计算得到摩擦力约束的位置误差修正向量

其中，μ_k，μ_s分别是动摩擦系数和静摩擦系数，ω_i和ω_j分别表示对第i个粒子和第j个粒子的加权系数；d表示滑动摩擦的距离，此处，动摩擦系数通常为0至1之间的值，通常取0.3，静摩擦系数根据实际需要进行改变，通常取值为0至1之间。

(3)凝聚力约束

作为本发明的优选实施方式，迭代求解凝聚力约束

粒子受到剪切应力，材质中的内聚力随作用的表面积增加而增加。从提供的附着力参数a 计算作用在粒子上的附着力F_a时，应考虑到接触面A_c：

F_a＝aA_c

其中，a表示附着力参数，A_c表示可用于计算凝聚力约束的任意一对粒子的接触面积，采用二次插值的接触面积，根据两个粒子之间的渗透关系得到接触面积A_c，通常凝聚力参数a的取值为0至1之间，不同的凝聚力系数会影响雪的凝结成块状的性质；

对于法向附着力，直接采用法向位移限制，并直接进行校正；对于切向附着力，通过牛顿运动定律可得：

进一步计算得到凝聚力约束的位置误差修正向量

其中，ω_i和ω_j分别表示对第i个粒子和第j个粒子的加权系数，t表示两粒子相对位置切线方向单位向量。

对预测位置进行修正，计算公式如下：

M＝diag(m₁，m₂，m₃)

其中

表示修正后的t+Δt时刻的第i个粒子的准确位置，

表示t+Δt时刻的第i个粒子的位置误差修正向量；M表示加权矩阵，m₁，m₂，m₃分别表示阻尼约束、摩擦力约束和凝聚力约束的加权项，

表示梯度计算函数，λ_i表示刚度系数；

和

分别表示阻尼约束、摩擦力约束和凝聚力约束的位置误差修正向量。

重复上述步骤，依次输出雪山山体雪层的每一时刻的图像，得到雪山山体雪层在雪崩场景下的模拟动画。

通过以上步骤，图1-4给出了本发明在实际使用中的过程以及产生的效果。图2为本发明的不同摩擦力约束模拟稳定时对比图，其中图2(a)至2(f)表示静摩擦力约束参数从大到小的情况下，粒子运动稳定时的结果对比图；图3为本发明的不同凝聚力约束模拟序列帧对比图，其中图3(a1)至3(a4)表示凝聚力约束参数为0的情况下的模拟序列帧，图3(b1)至(b4)表示凝聚力参数为0.2的情况下的模拟序列帧，图3(c1)至(c4)表示凝聚力参数为0.5的情况下的模拟序列帧；图4为本发明雪崩不同阶段模拟序列帧效果图，其中4(a)表示初始化状态，雪层稳定，4(b)表示受外力等因素，雪层不稳定开始发生雪崩，4(c)至4(e)为雪崩中间结果，4(f) 为雪崩最终结果与影响范围；

以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (5)

1.一种基于位置动力学的雪崩场景模拟方法，其特征在于，包括以下步骤：

1)将雪山山体雪层离散化为粒子，对所有粒子进行属性初始化；

2)根据粒子在当前t时刻的准确位置和受力情况，引入速度阻尼，得到t+Δt时刻所有粒子的预测速度和预测位置；

3)对所有粒子进行碰撞检测，若判断为粒子与刚体发生碰撞且粒子克服摩擦力发生运动，则更新t+Δt时刻粒子的预测位置和预测速度，否则保留步骤2)得到的t+Δt时刻粒子的预测位置和预测速度；

4)以最终更新后的t+Δt时刻粒子的预测位置和预测速度为参数，遍历所有约束条件，对预测位置进行修正，得到修正后的t+Δt时刻的所有粒子的准确位置和准确速度；

5)输出雪山山体雪层在t+Δt时刻的准确位置处的图像；

6)以t+Δt时刻所有粒子的准确位置作为起始位置，重复步骤2)-步骤4)，依次输出雪山山体雪层的每一帧图像，得到雪山山体雪层在雪崩场景下的模拟动画。

2.如权利要求1所述的基于位置动力学的雪崩场景模拟方法，其特征在于，所述的步骤1)具体为：

将雪山山体雪层离散化为粒子，设置所有粒子初始化的位置动力学参数，公式如下：

其中，x_i和v_i表示第i个粒子的准确位置和准确速度，

和

表示第i个粒子的初始化位置和速度，N表示粒子的个数，m_i表示第i个粒子的质量，w_i表示第i个粒子的权重参数。

3.如权利要求1所述的基于位置动力学的雪崩场景模拟方法，其特征在于，所述的步骤2)具体为：

2.1)定义已知的当前t时刻的粒子的准确位置

和速度

2.2)根据外力作用得到粒子在t+Δt时刻的预测速度：

其中，Δt表示时间步长，f_ext(·)表示所有外力的和的向量；

2.3)引入速度阻尼，进一步更新得到t+Δt时刻所有粒子的预测速度

根据

得到t+Δt时刻所有粒子的预测位置

计算公式如下：

4.如权利要求1所述的基于位置动力学的雪崩场景模拟方法，其特征在于，所述的步骤3)具体为：

3.1)对所有粒子进行碰撞检测，通过f_SDF函数计算粒子到刚体的符号距离函数的最近距离f，若最近距离为负值，则说明预测位置在刚体内部，粒子与刚体发生碰撞，否则未发生碰撞；最近距离f的计算公式如下：

3.2)若判断为粒子与刚体未发生碰撞，则保留步骤2)得到的t+Δt时刻该粒子的预测速度和预测位置；

若判断为粒子与刚体发生碰撞，首先计算粒子相对于碰撞物体的相对速度：

v_rel＝v_pOri-v_coll

其中，v_rel为粒子相对于碰撞物体的速度，v_pOri为粒子的速度，v_coll为碰撞物体的速度；

根据相对速度计算法向速度v_n，计算公式为：

v_n＝v_rel·n

其中n为碰撞物体的法向方向，由水平集φ的梯度方向来确定:

计算粒子切向速度v_t，计算公式为:

v_t＝v_rel-nv_n

F_f＝μ_fv_n

其中，μ_f为刚体表面的摩擦系数，F_f为静摩擦力；若‖v_t‖≤F_f，该粒子不发生位移，保留步骤2)得到的t+Δt时刻该粒子的预测速度和预测位置；否则，该粒子克服摩擦力运动，更新t+Δt时刻该粒子的预测速度和预测位置：

其中，‖v_t‖表示v_t的值，

和

分别表示更新后的t+Δt时刻该粒子的预测速度和预测位置；

根据动量守恒，同时计算出粒子对于刚体的平移速度以及角速度：

v_ref＝v_line+ω_angle×(x_i-x_mass)

其中v_line为粒子对刚体的平移速度，ω_angle是粒子对刚体作用的角速度，x_i表示第i个粒子的世界坐标向量，x_mass表示由N个粒子组成的整个物体的质心坐标。

5.如权利要求1所述的基于位置动力学的雪崩场景模拟方法，其特征在于，所述的步骤4)具体为：

4.1)以更新后的t+Δt时刻粒子的预测位置和预测速度为参数，遍历所有约束条件，所述的约束条件包括阻尼约束、摩擦力约束和凝聚力约束；对预测位置进行修正，计算公式如下：

M＝diag(m₁，m₂，m₃)

其中

表示修正后的t+Δt时刻的第i个粒子的准确位置，

表示t+Δt时刻的第i个粒子的位置误差修正向量；M表示加权矩阵，m₁，m₂，m₃分别表示阻尼约束、摩擦力约束和凝聚力约束的加权项，

表示梯度计算函数，λ_i表示刚度系数；

和

分别表示阻尼约束、摩擦力约束和凝聚力约束的位置误差修正向量；

4.2)迭代求解阻尼约束：

所述阻尼约束表示为：

其中，Δv_i表示速度误差修正向量，k_damping为阻尼系数；

求解步骤为：

计算整体运动的位置p_cm以及整体运动的速度v_cm：

求解中间变量L以及I：

其中，r_i为3x3矩阵，

表示r_i的迹；求解后得到ω角速度：

ω＝I^-1L

更新当前速度：

进一步得到阻尼约束的位置误差修正向量

4.3)迭代求解摩擦力约束：

所述摩擦力约束表示为：

求解步骤为：

根据当前时间步长中粒子的相对切向位移来计算摩擦的位置增量：

其中，

和

表示可用于计算摩擦力约束的任意一对粒子，n表示法向量；

进一步计算得到摩擦力约束的位置误差修正向量

其中，μ_k，μ_s分别是动摩擦系数和静摩擦系数，ω_i和ω_j分别表示对第i个粒子和第j个粒子的加权系数；d表示滑动摩擦的距离；

4.4)迭代求解凝聚力约束：

所述凝聚力约束表示为：

求解步骤为：

计算作用在粒子上的附着力F_a：

F_a＝aA_c

其中，a表示附着力参数，A_c表示可用于计算凝聚力约束的任意一对粒子的接触面积，采用二次插值的接触面积，根据两个粒子之间的渗透关系得到接触面积A_c；

对于法向附着力，直接采用法向位移限制，并直接进行校正；对于切向附着力，通过牛顿运动定律可得：

进一步计算得到凝聚力约束的位置误差修正向量

其中，ω_i和ω_j分别表示对第i个粒子和第j个粒子的加权系数，t表示两粒子相对位置切线方向单位向量。

Images