Abstract
我们提出了一种改进的位置动力学方法，该方法采用修正的光滑颗粒流体力学(SPH)核来模拟可变形固体。该方法采用遵循连续介质力学的应变能约束，既保持了基于位置方法的效率和稳定性，又提高了仿真的物理合理性。我们可以很容易地模拟各向异性和塑性材料的行为，因为该方法是基于物理的。与以往基于位置的连续材料模拟不同，我们使用弱结构粒子来离散模型，以便于可变形物体的切割。在这种情况下，采用修正后的SPH核函数来测量变形梯度并计算每个粒子的应变能。我们还提出了大变形条件下颗粒间反演和侵彻的解决方案。为了执行复杂的交互场景，我们提供了一种简单的碰撞检测方法。通过模拟各种场景，包括各向异性弹性变形、塑性变形、模型切割和大规模弹性碰撞，证明了该方法的灵活性、效率和鲁棒性。
Why(为什么做这个？)
The dynamic simulation of deformable solids is a classic topic in computer graphics.
可变形固体的动态模拟是计算机图形学中的一个经典课题。
能量约束PBD
Several studies have combined the continuum mechanics with the position-based approach to simulate many physical phenomena . Instead of using geometric distance to constrain the particle positions like the traditional PBD method, they used potential energy as the constraint, which enables the simulation of the physical properties of different materials by choosing the constitutive model.
一些研究将连续介质力学与基于位置的方法相结合来模拟许多物理现象。他们不像传统的PBD方法那样使用几何距离来约束粒子位置，而是使用势能作为约束，通过选择本构模型来模拟不同材料的物理性质。
采用基于粒子的方法
以往的研究大多采用基于网格的离散化方法，如有限元法（FEM），这种方法可以通过形状函数直接计算离散元素上的张量变量。然而，物体切割和断裂模拟将面临复杂的拓扑变化。
为了解决上述问题，我们采用了基于粒子的方法。
What(做的是什么事儿？)
An energy constraint of continuum mechanics is proposed in combination with particle-based discretization, which makes the PBD method follow the physical properties of the object during deformation while at the same time dealing with the cutting problem more conveniently.
将连续介质力学的能量约束与基于颗粒的离散化相结合，使PBD方法在遵循物体变形过程中的物理性质的同时，更方便地处理切削问题。
A corrected SPH kernel function is introduced to solve the interpolation problem of the deformation
gradient on each particle, which can preserve both linear and angular momenta.
引入修正的SPH核函数来解决每个粒子上变形梯度的插值问题，可以同时保持线动量和角动量。
A technique is developed to solve the problem of particle inversion and penetration in large deformation.
提出了一种解决大变形条件下颗粒反演和侵彻问题的技术。
How(作者是怎么做的？)
算法流程
PBD算法主要包括三个步骤。
首先，通过时间积分方案预测每个顶点xn+1的位置。
然后，定义一组双边约束C(x) = 0来修改预测的顶点位置。
最后，用修改后的位置更新顶点速度vn+1。
位置约束的定义
PBD方法的核心在于位置约束的定义。这里的目的是确定满足：
计算初始速度 V0
@ti.kernel
# 计算粒子初始速度 V0
def compute_V0():
for p in range(n_particles):
p_V0[p] = 0.
p_w[p] = 0.
p_dw[p] = [0., 0., 0.]
# 本实验考虑粒子的 6 个邻近粒子，即上下左右前后
for i in ti.static(range(6)):
# 检测到邻近粒子存在
if nei[p, i] != -1:
fX = x_o[p] - x_o[nei[p, i]]
# wpoly6() 为粒子插值核函数，将邻近粒子数量插入到特定粒子 ？
w = wpoly6(fX * h_dx)
p_V0[p] += w
# p_w[p] += w
# 此处加上 1e-14 是防止插值后粒子速度为 0
p_V0[p] = 1. / (p_V0[p] + 1e-14)计算形变梯度 F
V为初始粒子体积，u = x−X为粒子位移
相比之下，基于网格的方法可以直接利用每个单元的边缘向量来计算变形梯度。同时，对于基于粒子的方法，计算每个离散粒子上的变形梯度是困难的。这里，我们采用SPH的思想，使用修正后的SPH核计算每个粒子的变形梯度，如下所示:
for p in range(n_particles):
# 变形梯度 F
F[p] = ti.Matrix.identity(float, 3)
for i in ti.static(range(6)):
# nei[p, i]为neighbor（相邻）
# 本实验考虑粒子的 上下左右前后 共计六个粒子
# 所以为 range(6)
if nei[p, i] != -1:
fX = x_o[p] - x_o[nei[p, i]]
# w = quadratic_2d(fX, wdx)
fx = x[p] - x[nei[p, i]] - fX
# F[p] += fx.outer_product(fX) * w * inv_M
# dw = d_wspiky(fX * h_dx)
dw = nei_dw[p, i]
# 变形梯度 F
F[p] -= fx.outer_product(dw) * p_V0[nei[p, i]]
# 这里的 outer_product 是taichi中的计算矩阵外积的函数能量约束
if F[p].determinant() < 0.0:
# polar_decompose() 为taichi中矩阵的奇分解，
# 返回两个矩阵，第一个矩阵为上三角矩阵，第二个矩阵下三角矩阵
R, S = ti.polar_decompose(F[p], ti.f32)
e_c = S - ti.Matrix.identity(float, 3)
# 能量密度函数 Ψ(F)
# matric_cwise(A,B) 返回 矩阵||A - B||F^2 ,即 F范数的平方
C[p] = mu_0 * matrix_cwise(e_c, e_c) + lambda_0 * 0.5 * e_c.tra
# first Piola-Kirchhoff stress 第一皮奥拉-基尔霍夫 应力 P(F)
PF[p] = R @ (2 * mu_0 * e_c + lambda_0 * e_c.trace() * ti.Matr
else:
# F_E Green strain tensor 格林应变张量
F_E = 0.5 * (F[p].transpose() @ F[p] - ti.Matrix.identity(float
C[p] = (mu_0 * matrix_cwise(F_E, F_E) + lambda_0 * 0.5 * F_E.tr
PF[p] = 2 * mu_0 * F[p] @ F_E + lambda_0 * F_E.trace() * F[p] 
计算校正和函数梯度
@ti.kernel
# 计算校正梯度 ▽Wij(X)
def correct_weight(x_in: ti.template()):
for p in range(n_particles):
La = ti.Matrix.zero(float, 3, 3)
for i in ti.static(range(6)):
if nei[p, i] != -1:
fX = x_in[p] - x_in[nei[p, i]]
# fX_x = x_in[p] - (x_in[nei[p, i]])
# fX_y = x_in[p] - (x_in[nei[p, i]])
# dw = d_wspiky(fX * h_dx)
# dw = d_wspiky(fX * h_dx)
dw = d_wpoly6(fX * h_dx)
# dw_cot 校正的核函数梯度
dw_cot = (dw - p_dw[p] / p_w[p]) / p_w[p]
# La += dw_cot.outer_product(x_in[nei[p, i]]) * p_V0[nei[p, i]]
# La += dw_cot.outer_product(x_in[nei[p, i]]) * p_V0[nei[p, i]]
# 这里的 outer_product 是taichi中的计算矩阵外积的函数
# 至于为什么是 La -=
# 原因为此代码中 fx 对应关系为 Xp - Xi 对应论文公式就是 Xi - Xj
# 与论文中的 Xj - Xi 相反，故可以直接求外积后再 -= 即可
# 注意：La 的名词为 correction martix 修正矩阵
La -= dw_cot.outer_product(fX) * p_V0[nei[p, i]]
# if nei_num[p] == 2:
# print("beta_mat:", beta_mat, "beta_dw:",beta_dw)
# La = ani_mat(La)
La = La.inverse() # 矩阵取逆
# La = ti.Matrix.identity(float, 2)
for i in ti.static(range(6)):
if nei[p, i] != -1:
fX = x_in[p] - x_in[nei[p, i]]
# dw = d_wspiky(fX * h_dx)
# dw = d_wspiky(fX * h_dx)
dw = d_wpoly6(fX * h_dx)
dw_cot = (dw - p_dw[p] / p_w[p]) / p_w[p]
# 此时的 nei_dw[] 才是真正的校正梯度
# 其中 @ 为taichi中矩阵乘法运算符
nei_dw[p, i] = La @ dw_cot计算拉格朗日乘子 Δλ 和确定粒子位置
Cip_norm = 0.
Cip_p = ti.Vector([0., 0., 0.])
for i in ti.static(range(6)):
if nei[p, i] != -1:
# 校正的核函数梯度
dw = nei_dw[p, i]
# 这里先把 压力 P(F) 乘进去
Cip = PF[p] @ dw * p_V0[nei[p, i]]
# k=i 时，选择公式上部有负号的
Cip_p -= Cip
# norm_sqr() 为内置函数，计算矩阵的所有的元素的平方
Cip_norm += Cip.norm_sqr()
# k=j 时，选择公式下部为正号的
Cip_norm += Cip_p.norm_sqr()
alpha = 1. / 1. * inv_dt * inv_dt
# 论文中计算 Δλ 的公式实现代码
C[p] += alpha * p_lamb[p]
p_delta_lamb[p] = -C[p] / (Cip_norm + alpha)
p_lamb[p] += p_delta_lamb[p]
# 确定粒子位置
for i in ti.static(range(6)):
if nei[p, i] != -1:
dw = nei_dw[p, i]
# 计算 △X 将粒子 p 和 邻近粒子分开计算
v[p] -= p_delta_lamb[p] * PF[p] * p_V0[nei[p, i]] @ dw
v[nei[p, i]] += p_delta_lamb[p] * PF[p] * p_V0[nei[p, i]] @ dw 碰撞检测
for p in range(n_particles):
if x_o[p].z < obj_offset_z - 0.9 or x_o[p].z > ob
v[p] = [0., 0., 0.]
if x[p].y < bound_wall and v[p].y < 0.:
v[p].y = 0
if x[p].x < bound_wall and v[p].x < 0.:
v[p].x = 0
if x[p].z < bound_wall and v[p].z < 0.:
v[p].z = 0
if x[p].x > bound_length - bound_wall and v[p].x
v[p].x = 0
if x[p].z > bound_length - bound_wall and v[p].z
v[p].z = 0
x[p] += v[p]
for p in range(n_particles):
v[p] = [0., 0., 0.]
for i in ti.static(range(6)):
if nei_collision[p, i] != -1:
collision_idx = nei_collision[p, i]
if collision_idx >= n_particles:
fx = x[p] - x[collision_idx]
fx_dis = fx.norm()
if fx_dis < cut_pr:
move_ori = fx / fx_dis
v[p] -= (fx_dis - cut_pr) * (move
else:
fx = x[p] - x[nei_collision[p, i]]
fx_dis = fx.norm()
if fx_dis < p_r:
v[p] -= (fx_dis - p_r) * 0.5 * fx
if x_o[p].z < obj_offset_z - 0.9 or x_o[p].z > ob
v[p] = [0., 0., 0.]
if x[p].y < bound_wall and v[p].y < 0.:
v[p].y = 0
if x[p].x < bound_wall and v[p].x < 0.:
v[p].x = 0
if x[p].z < bound_wall and v[p].z < 0.:
v[p].z = 0
if x[p].x > bound_length - bound_wall and v[p].x
v[p].x = 0
if x[p].z > bound_length - bound_wall and v[p].z
v[p].z = 0
x[p] += v[p]
v[p] = [0., 0., 0.]
# @ti.kernel
# # 暂时没用上，以后在考虑
# def move_x():
# for p in range(n_particles):
# if x[p].y > 0:
# x[p].y -= dx * 0.1
更新粒子速度和位置
@ti.kernel
# 更新粒子的速度和位置
def update_velocities():
for p in range(n_particles):
# print("x - x_b [", p, "]: ", x[p] - x_o[p])
v[p] = (x[p] - x_b[p]) * inv_dt
x_b[p] = x[p]
# print("v[p]: ", v[p])演示demo（未渲染）
参考文献
An energy constraint position-based dynamics with corrected SPH kernel.pdf
https://zhuanlan.zhihu.com/p/442450577
https://aiuai.cn/aifarm2004.html
https://blog.csdn.net/weixin_43940314/article/details/129830991
评论区