Homework2: MPM实现沙子和水的耦合

由于本人还是苦逼本科生+疫情把final推迟到了暑假末,所以没有在ddl前赶上hw2(其实那个时候也才刚啃完MPM)。最近趁复习摸鱼时间重新把hw2做了,复现了蒋老师的两篇和沙子有关的论文:
Multi-species simulation of porous sand and water mixtures 和 Drucker-Prager Elastoplasticity for Sand Animation,用MPM实现了沙子和水的耦合。


干沙子 (10000个沙子粒子,128x128网格)
video
湿沙子+水 (10000个沙子粒子,每次迭代+50个水粒子,128x128网格,约26fps)
video
湿沙子+水 显示出project操作的不同状态(和论文中使用相同的着色)

干沙子+水,根据我的理解,论文在耦合时默认沙子初始是湿的并没有考虑毛细作用,如果我代码没有写错的话,感觉这个结果并不是非常正确。

复现的过程可谓是非常痛苦,踩了许多的坑,据不完全统计有以下这些:

  1. hardening操作的参数 h0,h1,h2,h3 文本中提到的是弧度制而参考表中给的是角度制,虽然下面有一行说明但是一开始没看到。。
  2. 一开始分别实现了沙子和水,耦合的时候忽略修正沙子和水的密度,导致数量级差距过大动量项交换异常。
  3. 应力求错了,不得不说在simulator里调试NAN相关的问题非常痛苦。
  4. 如果沙子处在屈服面的范围内,project返回的时候不考虑论文里修正体积的项
  5. 在处理 cohesion 的时候论文里把不等式右端项直接换成了 cC 但是没有给出明确的推导。。本人数理基础较差,最后是参考 这里 的方法,类似于体积修正项在 project 操作前给 e 额外加上 \frac{c_C}{d*\alpha_p} \mathbf{I} 的项,也就是:
e += (c_C0[p] * (1.0 - phi_s[p])) / (d * alpha_s[p]) * ti.Matrix.identity(float, 2) # effects of cohesion

此外还有几点对17年论文的疑惑,希望能够得到解答:

  1. section 3.1.1的应力公式貌似写错了,应该是 \frac{\lambda}{2} \mathbf{tr^2(\epsilon)} ,原文是 \frac{\lambda}{2} \mathbf{tr(\epsilon)}
  2. 在实现中,感觉体积修正项的正负号和论文里的正好相反,按论文里写貌似不对?
  3. section 4.3.4 末尾关于更新 \mathbf{F}^{s E, n+1} 的公式貌似是错的?

完整demo

代码托管在:

鸣谢:
非常感谢 archibate 学长热心帮忙解决基础问题:https://github.com/taichi-dev/taichi/issues/1777
蒋老师在 GAMESWebinar 上关于本论文的报告:https://v.qq.com/x/page/d0539skispz.html

总结:
Taichi 复现论文真的是很好用的工具,希望以后有时间多写写hw2-extra!

10 Likes

有点厉害!我自己当年都没实现过这个paper :joy: 大家战斗力太强了…

太帅了!
1.good catch
2. good catch
3. good catch

1 Like