基于物质点法的沙⼦模拟

学学位论 基于物质点法的沙模拟 作者姓名: 夏鸣 学科专业: 凝聚态物理 导师姓名: 刘利刚 教授 完成时间: 年五七 University of Science and Technology of China A dissertation for bachelors degree Material Point Method for Sand Simulation Author: Yiming Xia Speciality: Condensed Matter Physics Supervisor: Prof. Ligang Liu Finished time: May 17, 2018 i 致谢 感谢我的爸爸妈妈，一直以来是他们在背后默默地为我鼓劲，为我加油。是 他们提供给我莫大的精神支持，让我有了不断挑战困难的勇气与动力。我也会将 这份鼓励带到以后的生活中，健康快乐地迈向前方。 衷心地感谢刘利刚老师一年多来对我的教诲。 刘老师的计算机图形学课让我 受益匪浅，老师以他独到的训练方法，巧妙地将我从一个图形学的门外汉带入到 图形学的世界中。我对能遇到刘老师这样的良师感到十分幸运。老师不光教给了 我们图形学的知识，更教给了我们他严谨，科学的科研方法。毕设期间，是老师 循序渐进地引导我，才使得我克服对难题的恐惧，将问题分而治之，一步一步完 成毕设的任务。 感谢葛子恒同学的帮助， 我们针对同一类问题， 结队学习， 过程中多次交流， 互帮互助，让我对团队合作有了新的理解。 感谢胡渊鸣师兄对我做出的结果给出的建议。 感谢每一位关心，支持，帮助过我的老师和同学！ 夏一鸣 2018 年 5 月 15 日 ii 目录 致谢 目录 摘要 Abstract 第第一一章章 绪论绪论 1 1 1.1 引言1 1.2 相关工作 3 1.3 本文结构3 第二章第二章 物理基础物理基础 4 4 2.1 运动描述4 2.1.1 拉格朗日视角 4 2.1.2 欧拉视角5 2.2 动力学描述 6 2.3 本构模型 6 2.3.1 柯西应力与能量密度函数6 2.3.2 弹塑性模型7 2.3.3 屈服条件与屈服面8 2.3.4 塑性流动法则9 2.3.5 回映算法9 2.3.6 硬化 10 第三章第三章 物质点法物质点法 1 12 2 3.1 算法流程 12 3.2 算法详解 12 3.3 模拟相关技术 19 3.3.1 CFL 条件 19 3.3.2 碰撞检测与处理 19 3.3.3 Possion 采样 20 第四章第四章 算法部分修正算法部分修正 2222 4.1 APIC 插值 22 4.2 体积修正 23 4.3 插值上的细节问题 24 第五章第五章 结果展示结果展示 2 25 5 iii 第六章第六章 总结与展望总结与展望 2 28 8 6.1 总结 28 6.2 展望 28 附录附录 2 29 9 参考文献参考文献 3131 iv 摘要 沙子是日常生活中随处可见的物质， 静止时， 它的运动行为像固体， 运动时， 它的运动行为像流体。 本文不同于以往直接使用粒子来模拟沙子的方法， 而是将沙子视为一种复杂 流体，引入了一个用于描述沙子运动行为的 Drucker-Prager 弹塑性模型，再使 用物质点法进行离散求解，最后结合计算机图形学中的显示渲染技术，模拟出沙 子的运动行为。同时，针对模拟计算中插值方法略有不足的问题，本文引入了 APIC 插值方法， 改善了算法的稳定性。 本文还针对物质点法模拟， 设计了一个模 拟框架，并将其封装成库。 关键字：模拟仿真，物质点法，弹塑性模型，沙子 v Abstract Sand is a substance that can be seen everywhere in daily life. When it is stationary, its movement behaves like a solid. When it is in motion, its movement acts like a fluid. This paper is different from the previous method of directly using particles to simulate sand. Instead, sand is regarded as a complex fluid. A Drucker-Prager elastoplastity model is introduced to describe the motion behavior of sand. Finally, combined with the display rendering technology in computer graphics, the motion of sand was simulated. At the same time, for the problem that the interpolation method is slightly insufficient in the simulation calculation, this paper introduces the APIC interpolation method, which improves the stability of the algorithm. This article also aims at material point method simulation, designed a simulation framework, and packaged it into a library called “Flame”. Key words: simulation, material point method, elastoplastity model,sand 中国科学技术大学本科毕业论文 1 第一章 绪论 1.1 引言 1.1.1 计算机模拟仿真 对一个问题的研究，往往分为理论与实验两部分。根据实验现象，在理论上 提出模型，预测行为，然后进行实验验证，根据新的实验结果修正理论模型，如 此往复迭代。但是有时候，实验验证的成本较高，危险性较大，实验周期长，因 此，人们往往会用计算机模拟的方法，在理论与实验之间搭建桥梁，将理论提出 的模型用数值方法在计算机中进行模拟实验， 从而提高研究效率， 节约研究成本。 计算机模拟常用的思路是：将物质的运动学、动力学方程进行时间、空间上 进行数值离散，从而能够在计算机中算出每一个离散时间节点，每一个离散空间 区域内的物理量。进而可以利用这些计算出的结果去实现不同的目的，比如：飞 机、火箭等航天器的设计(图 1.1)，水流中污染物质的扩散，桥梁架构设计的合理 性评估等等，可以应用到很多工业行业当中。同时，由于模拟的结果与现实结果 十分相近，因此，将模拟的结果可视化，便可以应用到计算机动画(图 1.2)、游戏 开发当中。 1.1.2 流体模拟仿真算法 模拟的对象不同，所采用的算法也不同。本文主要涉及关于流体的模拟算 法。流体的离散计算方法，目前主流的有网格法和粒子法两大类，网格法中， 有基于三角网格的有限元法（Finite Element Method FEM19） ，基于正方形网 格的有限差分(Finite Difference Method FDM)，格子玻尔兹曼方法（Lattice 图 1.1 飞机制造设计 图片来源17 图 1.2 计算机动画 图片来源4 中国科学技术大学本科毕业论文 2 Boltzmann Method LBM24）等等，它们将物理信息存放在静止不动的格点 上，着眼点放在了格点上物理信息的传播。格点粒子法有光滑粒子动力学方法 (Smoothed Particle Hydrodynamics SPH 8)，分子动力学方法（Molecular Dynamics MD 11） ，它们将物理信息存放在粒子上，粒子是随时间变化而运动 变化的。粒子法的着眼点放在了粒子的运动变化上。除此以外，还有仅用来模 拟水波的浅水波方程（Shallow Water Equation SWE 13） 。从分析流体的两种视 角来看，网格法是从欧拉视角分析，粒子法是从拉格朗日视角分析。从人们实 际的实践经验来看，网格法计算较为稳定，不太可能出现数值爆炸的情况。但 是，因为网格法中信息的载体是被固定的网格格点，其存在拓扑上的约束，而 粒子法的信息载体则没有这方面的约束，所以粒子法在模拟物质的变形上有它 自身的优势。 本文使用的物质点法是一种将两者的长处结合的方法，它最早在 1995 年被 提出，通过不断的发展完善，物质点法开始被广泛地应用到流体动画中。首先 是 2013 年在冰雪奇缘中，场景中的雪，大多都是用物质点法2模拟出来 的。之后在 2015 年，蒋陈凡夫5等人提出了 APIC 插值的方法，改进了物质点 法插值时角动量不守恒的问题，这一技术被用到了 2016 年的海洋奇缘中。 2016 年，物质点法结合 Drucker-Prager 弹塑性模型3，给出了干沙子模拟的算 法。随后，在 2017 年，物质点法框架下水和沙子耦合的解算方法4，布料与 其他物质的耦合的解算方法被提出。2018 年，胡渊鸣等人提出了在物质点法框 架下求解物体切割问题的新方法。物质点法还在不断被完善和发展。 物质点法的大致思想是在研究区域中加入背景网格，将物质的信息存储在 粒子上，而计算则放在背景网格上进行，两者之间的信息交互是通过粒子与其 周围的格点进行插值来进行的。总的来说，物质点法适用于连续介质的模拟， 它能够模拟的物质大多具有流动性，或者弱不可压缩。物质点法相对于其他方 法而言，能够很好地计算梯度，同时由于背景网格自身不会出现缠绕等现象， 所以物质点法可以轻易地处理碰撞问题。 1.1.3 沙子模拟 本文的模拟对象是沙子，沙子是日常生活中非常常见的东西，在物理学中， 它们被定义为颗粒体系， 而颗粒体系通常指的是宏观固体颗粒组成的多体耗散体 系。 不同堆积密度的颗粒体系可以表现出类似于液体和固体的行为。 这种复杂的 中国科学技术大学本科毕业论文 3 行为想要通过数值去定量描述是很困难的。 如果我们将沙子视为一个一个粒子，粒子与粒子间存在弹力，以及摩擦力， 然后在计算机中求解，那么我们需要计算百万量级数目的粒子，运算量极大。而 且很难调整参数，使得模拟的结果与实验中沙子的行为一致。 本论文采用34中采用的方法，将沙子视为连续介质，在物质点法的离散 求解框架下，结合 Drucker-Prager 弹塑性模型，以及 APIC 插值方法，模拟了沙 子的运动行为。 1.1.4 本文工作 本文的主要工作有： 1.实现了物质点法离散求解的框架，可以用来模拟雪，水，弹性体的运动 2.实现了34中的模拟沙子的算法 1.2 相关工作 Y Zhu, R Bridson20 将沙子视为流体，模拟了沙子的流动行为，模拟的结 果看起来带有粘滞力。 Narain21等人改进了 Bridson 等人的算法， 消除了模拟中 与不可压缩有关的内聚力。Lanenerts 和 Dutre22提出了用 SPH 框架下求解水和 多孔颗粒物耦合的算法。Chang23提出一个调整后的 Hooke 定律来处理颗粒间 的摩擦力，得到了更好的沙子模拟结果。 1.3 本文结构 第一章给出本文相关的简介。 第二章描述了在本文模拟沙子时所需要构建的物理模型。 第三章详细描述了物质点法的算法流程，并结合物理模型，给出沙子模拟流 程。 第四章针对算法中存在的一些问题给出相应的改进方法。 第五章展示了本文模拟的结果。 第六章对本文的工作进行总结。 附录针对文中一些问题的进行求解推导。 中国科学技术大学本科毕业论文 4 第二章 物理基础 2.1 运动描述 研究物质的运动，就是要给出它的位置与时间的关系，根据这些关系，我们 能够获得物质在几何上的一些特征，进而进行一些相应的计算。 对于连续介质而言，目前有两种观点对它进行描述，它们分别是拉格朗日观 点和欧拉观点。 2.1.1 拉格朗日观点 在拉格朗日观点下，观察点处在物体上。随着物体的运动，观察点也同步运 动。 1.位置、速度与加速度 在被研究的物质所占的区域0中选取一点。 随着时间进程的推进， 如图 2.1， 0演化为， 对应的演化到， 这个演化过程可以用一个映射关系描述， 写作： = (,).(1.1) 由(1.1)式我们也可以推出，点在时刻的速度为： (,) = (,) .(1.2) 由此也得到此时的加速度为： (,) = (,) = 2(,) 2 .(1.3) 图 2.1 映射关系 中国科学技术大学本科毕业论文 5 2.形变梯度 现在取点周围一段非常小的连续空间，由（1.1）式我们可以得到点附 近空间演化的对应关系： = (,) . (1.4) (,) 是一个二阶张量，一般称为形变梯度(,)， （1.4）也写成： = (,).(1.5) 形变梯度给出了物质局部形变（除了平移）的描述。由(1.4)的形式可以看 出， (,)相当于是一个 Jacobi 矩阵， 进而可以推出， 点处局部体积的变化为： = (,).(1.6 ) 式中(,)=(,)。 2.1.2 欧拉观点 在欧拉观点下，观察点是坐标系中固定的点，它不随时间变化而变化。 1.位置、速度与加速度 对于坐标系中固定的点，当有物质经过该点时，我们需要按图索骥，依据 该点位置以及当前时刻，回溯找到物质区域0中对应的点。这个过程可以表 示为： = 1(,).(1.7) 将在坐标系中固定点处，于时刻流经的物质的速度记为(,)，从定义不 难得出，(,)为0中1(,)处在时间时的速度，即： (,) = (1(,),) = (,) . (1.8) 由(1.8)推出此时处的加速度为： (,) = (,) = (,) (,) + (,) = (,) (,) + (,) .(1.9) (1.9)由两部分组成， 前式是由于观察点处的物质发生变化而引发的加速度， 后式是由于观察点处的速度的变化引发的加速度。 2.物质导数 将上面的研究对象从速度换成其他与位置以及时间相关的物理量，那么在 欧氏观点下，随时间的变化率为： (,) = (,) (,) + (,) = (,) (,) + (,) .(1.10) 式中 定义为物质导数。 中国科学技术大学本科毕业论文 6 3.形变梯度随时间的变化率 在拉格朗日视角下，其变化率为(,) ，带入(1.5)中给出的定义，并带入式 (1.2),得到： (,) = (,) = (,) = (,),) = (,) (,) = (). 即有： = ().(1.11) (1.11)式表明， 点处形变梯度的变化率与改点当前的形变梯度以及对应点 处的速度梯度有关。 2.2 动力学描述 沙子的模拟是基于弹性体运动方程的，所以整个运动方程为： = = + (1.12) 2.3 本构模型 有了运动学与动力学的描述后，我们便可以进行运动方程的求解。但是， (1.12)中的应力我们无法简单的通过运动学与动力学的描述得到，它与物质自 身的性质，以及当前应变相关，我们将这种关系进行数学化描述，过程中所建立 的物理模型称为本构模型。 2.3.1 柯西应力与能量密度函数 根据弹性理论16，在物质所占的区域分布着一个与形变梯度有关的形变 能量密度分布函数()，由此可以计算出某点周围的柯西应力： = 1 () () .(1.13) 研究的物质不同， 其能量密度分布函数()也不同， 在本文的沙子模拟中， 采用的分布函数为： () = ()2) + 1 2 ()2.(1.14) 中国科学技术大学本科毕业论文 7 式中是通过对形变梯度做 SVD 分解： = 得到的，将(1.13)中() 分子 上的()用(1.14)的形式替换，可以得到： () = (21 + ()1 ) . (1.15) （(1.15)的推导见附录） ，这样就可以计算出。 2.3.2 弹塑性模型 假设现在有一根金属丝，如果我们轻轻地用力弯曲它，然后释放，金属丝会 回弹到初始时的形状，在这个过程中，金属丝只发生了弹性形变。但是，如果我 们使用较大的力弯曲它，释放后，金属丝仍然会回弹，但是它回不到初始的形状 了，在这个过程中，金属丝既发生了塑性形变又发生了弹性形变。对于金属丝这 种既有弹性形变又有塑性形变的物质，我们要使用弹塑性模型进行分析。 如图 2.2，对物质施加外力，使得其从初始构型 A 变化到构型 C。如果释放 外力，物质会变化为构型 B。这是具有弹塑性物体常常发生的物理过程。我们关 心的问题是，当物质处在构型 C 时，它的受力情况如何？我们知道，在 C 构型 中， 能产生力的是弹性形变，即对应着弹性形变梯度，我们记作。但是从 A 到 C 的形变梯度并不是完全是弹性形变梯度，其中有一部分是塑性形变梯度， 即： = .(1.16) 我们的已知是，所以我们需要根据来计算，算法是这样的，首先，我 们假设完全是弹性形变梯度 ，但是，有的材料会有弹性极限，即它们可能不 会允许有这么大的弹性形变梯度发生，比如要求形变梯度 SVD 分解后(0,0)小 于某个值。如果 超出了这个弹性极限，那就意味着，中存在塑性形变。此 时， 我们要对 按照某种法则进行修正， 修正直到在弹性极限内为止， 这时， 即为，而= 1，这样就完成了对的分解。进而可以计算出应力 图 2.2 弹塑性模型 图片来源2 图 2.3 受力分析示意图 中国科学技术大学本科毕业论文 8 2.3.3 屈服条件和屈服面 上面说到，各种物质对的约束不同，对于沙子而言，它的约束关系是这样 的：如图 2.3，我们取中的一个小微元进行受力分析。将该微元分为正区域和 负区域，在交界面上取一点，该处的法向量为，切向量为。正区域对负区域 在点处的作用力为(,)，由弹性理论16，(,)与的关系为： (,) = .(1.17) 设摩擦系数为，由库伦摩擦定律，得到： + 0.(1.18) 将(1.17)以及几何关系 = (是顺时针旋转 90 度的矩阵)带入(1.18)中， 得到一个关于应力的不等式： + 0.(1.19) 由于是对称阵， 所以， 我们可以对做特征值分解， 得 = ， 带入(1.19)， 得： + 0. 将记为 ，得： + 0.(1.20) 图中仅仅给出了某一个剖面的情形，我们的要求是，对于任何剖面，即任意 一个 （ = 1） ，(1.20)式恒成立。这个问题可以转化为一个带约束的极值问 题，使用拉格朗日乘子法进行求解(这里我们对二维的进行求解，推导方法见附 录，三维类似)，得到需要满足的条件： 1+ 2 2 2+ 1 + | 1 2| 2 0.(1.21) 这里1、2来自S = (1 0 02)。 记 = 22+1，将(1.21)写成矩阵形式： () + () 2 0.(1.22) 推广到维有： () + () 0.(1.23) (1.23)即为需要满足的条件，成为屈服条件，将(1.21)式在二维坐标系中 画出， 得到如图 2.4(a)的结果。 可以看到， 1和2的可行域构成一个锥状区域 （图 中深蓝色区域） ，三维情形下，可行域是一个圆锥。这个主应力空间中圆锥面称 中国科学技术大学本科毕业论文 9 为 Drucker-Prager 屈服面。 不同，屈服面的锥角大小也不同。 2.3.4 塑性流动法则 根据 2.3.2 中给出的算法，当假设完全是弹性形变梯度 时，它所产生的 应力在主应力空间中可能有如图几种情况，(1)如图 2.4(b)中的蓝色点，位于屈 服面内。(2)如图中的紫色点，位于屈服面上。(3)如图 2.4(b)中的绿色点，位于 屈服面外， 但是主应力相加小于零。 (4)如图 2.4(b)中的橙色点， 位于屈服面外， 但是主应力相加大于零。对于第(1)、(2)种情况，屈服条件已经得到满足，但是 第(3)、(4)种情况违背了屈服条件，根据我们之前所推导的理论，我们对形变梯 度进行修正，从而计算出它们的实际主应力。修正需要满足保持不变、能量不 增、体积守恒的条件，这一修正过程被称为塑性流动。根据推导，我们得到，对 于第(3)种情况，实际的主应力为当前应力向屈服面，以垂直于1= 2这条直线 （图中虚线） 的方向投影， 而第(4)种情况， 实际主应力为0。 投影情况如图 2.4(c)。 2.3.5 回映算法 由于我们的模拟是在离散情况下进行的， 所以我们需要有一套离散情形下的 流动准则，这种离散情况下的流动准则，我们称之为回映算法。算法是这样的， 首先假设一个弹性形变，得到一个试探应力。然后，判断试探应力是否在屈服面 之内，如果是，则保持目前的形变梯度不变。否则，使用流动法则，计算出此时 满足条件的形变梯度。 图 2.4 Drucker-Prager 屈服面 图中的蓝色区域表示屈服面，图(b)给出了不同应力在主应力空间的分布情况， 图(c)给出了不同应力向屈服面的投影情况 注注： ( (a)a) ( (b)b) ( (c)c) 中国科学技术大学本科毕业论文 10 实际计算中我们将上述过程定义成一个函数( ,)， 首先， 对 进行 SVD 分解 = ，此时试探应力= ln ()，定义两个量： = () .(1.24) = + ().(1.25) 可以看出， (1.25)式即(1.23)式不等式的左半部分。 因为 = /22+ 1， 所以(1.25)式实际上计算了在主应力空间中到屈服面的距离。 而(1.24)式则是 垂直于主应力空间(1,1,1)向量的量。所以，可以根据这些量来决定是否需要进行 屈服面上的投影操作。 情况 1. 0且有() 0或者的范数为零。那么被 投影到屈服面圆锥顶点，即坐标原点，( ,)返回 ，情况 3.上面两种情 况没有覆盖的其余情况，如前面所述，沿直线投影到屈服面上： = .(1.26) 根据对(1.24)式以及(1.25)式的理解，可以看出，(1.26)式给出了在屈服 面上的投影点，这里( ,)返回 。 2.3.6 硬化 根据之前人们的研究理论3， 沙子的塑性形变会产生硬化， 这种硬化会对沙 粒间的摩擦力产生影响，我们将塑性形变量用表示，它与计算弹性形变梯度 时产生的塑性形变梯度的变化有关。具体的形式是: 对于 2.3.5 中的情形 1，没有塑性形变发生，有： = 0. 对于 2.3.5 中的情形 2，此时主应力全部投影到坐标原点，所以： = ,+1 . 对于 2.3.5 中的情形 3，此时主应力投影到屈服面上，有： = . 硬化对摩擦力的影响体现在对 影响上，设粒子的 为 ，在运动过程中， 由于塑性形变所产生的影响，通过下面的式子体现： +1 = + .(1.27) 中国科学技术大学本科毕业论文 11 = 0+ (1 +1 3)2 +1. (1.28) +1 =2 3 2 3 .(1.29) (1.27)式的 +1给出了硬化的状态，(1.28)式的 称为摩擦角，可以看到有 0，1，2，3四个参数，它们影响着沙堆静止时的角度，之后我们会展示0参 数的影响情况。 中国科学技术大学本科毕业论文 12 第三章 物质点法 3.1 算法流程 算法流程简图如下 3.2 算法详解 下面对算法流程中的步骤进行解释 3.2.1 插值 这里的插值与拉格朗日插值法的思想一样，已知某些点的函数值，使用插值 的方法近似得出其余点的函数值。说通俗一点，就是我们首先拥有一部分采样点 的信息，通过对它们取加权平均，估算出某个特定位置处的信息，这个权重与 到采样点的距离有关。在本文的模拟中，我采用了 B 样条函数作为权重函数。 B 样条函数的形式如下： () = 1 2 |3 2+ 2 3 0 | 1 1 6 |3+ 2 2| + 4 3 1 | 2 02 | (2.1) 输入：粒子坐标、速度、质量 1. 初始化 2. 粒子插值到网格 3. 基于网格求解格点受力 4. 更新网格各物理量 5. 网格插值到粒子 6. 更新粒子各物理量 7. 输出：粒子坐标 8. 回到第 2 步 图 3.1 算法流程图 中国科学技术大学本科毕业论文 13 在 我 们 的 符 号 体 系 下 ， 是 一 个 向 量 ， 所 以 实 际 计 算 中 () = (1)(2)(3)，其中1、2、3是的三个分量 现在我们来研究点的插值情况，记点的坐标为，根据定义，我们要需要 找到与距离小于 2 个单位长度的采样点。物质点法模拟中，背景网格的格点 即为采样点，而网格长度即为单位长度，所以()中的自变量此时要换为 ( )/。 B 样条有着优良的性质，能够在插值过程中保证一部分信息不丢失。对于点 而言，将其与周围的格点算出来的权重记为，则有以下性质： = 1.(2.2) = .(2.3) 联合(2.2)与(2.3)可以推出： ( ) = 0.(2.4) 这些性质在我们后面的一些推导中会用到。 3.2.2 初始化 因为我们的输入只有粒子的位置，初速度以及粒子的质量，缺少粒子的体积 信息，所以在初始化的过程中，我们要通过插值法估算出粒子的体积。算法是这 样的，首先，用初始粒子的质量 0插值计算格点的质量0，其结果为： 0 = 0. (2.5) 之后，我们来估算格点的密度，格点所占体积为，是问题的维度。所以 图 3.2 物质点与背景网格 图 3.3 B 样条函数曲线 中国科学技术大学本科毕业论文 14 = 0/，最后，用格点的密度插值回粒子，得到粒子的密度： = .(2.6) 由此可以得到粒子的体积: 0 = 0 .(2.7) 在我们的算法中，粒子的体积在演化过程中是恒定不变的，所以其值一直是 0。除此以外，粒子的质量也是保持初始值0不变的。 3.2.3 粒子插值到网格 首先是格点质量 ，形式上同(2.5)式： = 0. (2.8) 此处证明，插值前后，格点总质量等于粒子总质量： = 0 = ( ) 0 = 0 . 之后计算格点的速度，为了能保证插值前后总动量守恒，因此我们先计算格 点的动量，由格点的动量计算格点的速度。有： = 0. (2.9) 因为 = ( )/ ，所以得到格点的速度为： = 0 0 .(2.10) 下面证明按照这样的插值方法，格点总动量等于粒子总动量。 格点总动量为： = .(2.11) 粒子的总动量为： 中国科学技术大学本科毕业论文 15 = 0 .(2.12) 将(2.9)带入(2.11)得： = ( 0) = 0 ( ) = 0 = . 上面的推导中用到了 B 样条的性质(2.2)。 3.2.4 基于网格求解网格受力 由于网格是我们人为引入的，现实中并不存在，所以无法直接计算出格点的 受力，需要间接求出。此处我们利用虚功原理来进行计算。设物质总弹性能量为 ，由(1.14)中定义能量密度分布函数，结合步骤一中算出的粒子体积 0，我们 得到： = 0( ) .(2.13) 式中， 下标是指该量与粒子有关。 有了总能量后， 只需要在处加入一个扰动， 移动到 ，即可算出此时处的受力： = = 0 ( ) = 0 ( ) .(2.14) (2.14)中， ( ) 可以用(1.15)式计算，缺少的是 。我们暂且搁置，后面 会给出它的计算方法。 3.2.5 更新网格各物理量 本文的模拟用的是显式积分格式，所以形式较为简单： +1 = + + .(2.15) 式中， 表示粒子当前时刻受到的外力，一般是重力。是时间步长，它的值 需要满足 CFL 条件。关于 CFL 条件，后面会有描述。此处还需要注意的是，对格 点要进行碰撞检测，即针对格点到边界不同距离的情况，对格点速度做出不同的 处理。碰撞处理这一部分后面会进行详细的论述，本文仅对格点做碰撞处理，对 粒子做碰撞处理容易导致 与不同步。 中国科学技术大学本科毕业论文 16 3.2.6 格点插值到粒子 不同于前面的粒子插值到格点， 这里粒子的速度是直接插值计算出来的， 即： +1 = +1 .(2.16) 式中， 权重还是取 。 此处仍然证明下格点总动量等于粒子总动量， 利用(2.11) 与(2.12)式，并将(2.16)式带入(2.12)，有： = 0 ( +1 ) = +1 ( 0 ) = +1 = . 此处没有利用格点动量去计算粒子的动量，这是物质点法的算法使然。如果 用格点动量计算粒子动量， 势必要用格点质量插值计算粒子质量， 而前面也提到， 粒子质量在整个过程中是恒定不变的， 采用插值法估算粒子的质量虽然能保证总 质量不变，但是不能保证粒子质量守恒。所以，这里直接用格点速度插值粒子速 度。 3.2.7 更新粒子各物理量 此步骤中需要更新的有粒子的坐标 以及形变梯度 ，下面分别论述更新 方法。和格点的速度更新一样，粒子的速度采用显式积分更新： +1 = + +1. (2.17) 对于形变梯度 ，将(1.11)式进行离散，有： +1 = ( +1) . (2.18) 将(2.16)带入(2.17)， 注意， 此时 +1是常量， 只作用在 上， 所以， 有： +1 = ( +1 ( ) ) . 对上式移项可得： +1 = ( + +1 ( ) ) . (2.19) 此式即为 更新的算法。 由此， 我们可以近似算出(2.14)中较难计算的 。 算法是这样的：首先，用 来近似 +1，且此时左边不再是 +1， 而是与 有关的函数 +1，这样(2.19)就可以写成： 中国科学技术大学本科毕业论文 17 +1 = ( + ( ) ( ) ) . (2.20) 这时，对于下一时刻的形变梯度 +1而言，其自变量为 ，表示格点将要 到达的位置。(2.20)式相当于描述了形变梯度关于位置分布的关系，而 表示 格点移动一个小量后，形变梯度的变化率。结合两者的理解，利用(2.20)式得 到： = ( ) . 再利用矩阵求导的相关结论，有： = ( ). (2.21) 将(2.21)带入到(2.14)，我们得到： = 0 ( ) ( ) .(2.22) 这是对格点力计算的补充。 回到 的更新， 以上所描述的物质点法流程都与物质的本构模型无关， 是相 对独立的部分。 但是在这里， 我们需要将二者建立其联系。 根据式(1.16)的结果， 我们知道当前时刻 可以分为弹性形变 ,和塑性形变 ,两个部分，现在来求 下一时刻的形变梯度的分解结果。 图 3.4 形变梯度的更新 中国科学技术大学本科毕业论文 18 更新的思想如图 3.4，从初始时刻说起，初始时的构型我们记作 A，这个构 型下材料不受外力。来到初始时刻时，材料此时产生一个形变梯度 ，变成构型 C，我们的目标是将其分解成两个部分，从而求出中间构型 B。回顾 2.3.2 节中 的求法：我们先假设 全部导致了一个弹性形变，但是，实际而言，这种材料的 弹性形变不是任意的，是有个界限的，因此，我们要对其进行分解，如果 在这 个范围内，我们就认为此时的物质弹性形变梯度 ,就是 ，如果不在这个范围 内， 我们要按照该物质的特性， 由 计算出该物质的弹性形变梯度 , = ( )， 进而得到( ,)1 ，即为塑性形变梯度 ,。这样，就完成了弹塑性分解，从 而得到了中间构型 B。 现在来到下一时刻， 由更新的法则(2.18)可以得到： +1 = ( + +1) , ,，这样，我们得到了构型 E，由于 B 构型是 C 构型释放外 力后的构型，与一开始 A 构型的状态一样，所以这时，我们以 B 为此刻的初始 构型，得到一个 E 相对于 B 的形变梯度： +1, = ( + +1) ,. (2.23) 这时重复之前求解弹性形变梯度的步骤，有 +1, = ( +1,)，这样，我们 得到了从中间构型到最终构型的弹性形变梯度，从而得出此时的塑性形变梯度 +1, = ( +1,)1 +1。从这里的论述，我们可以看出，我们是将一个大的形 变过程分解为一个个小的形变过程， 这些小的形变过程又可以分解为弹性形变与 塑性形变，通过累积小的弹性形变得到当前的受力情况。 对于(2.22)式而言， 当()仅与 ,有关时， 要使用 ,代替(2.20)中的 ， 从而(2.22)式改写成： = 0 ( ,) , ( ,) .(2.24) 3.2.8 输出粒子信息 将粒子的坐标等信息输出，作为之后更进一步的应用做准备。本文主要是将 这些信息可视化，比如用 OpenCV 对二维粒子进行了绘制，用 Blender 对三维粒 子进行渲染绘制。也可用到工程计算，设计中。 中国科学技术大学本科毕业论文 19 3.3 模拟相关技术 3.3.1 CFL 条件(CourantFriedrichsLewy condition) CFL 条件是求解偏微分方程数值离散解的重要条件，在显式积分中尤为重 要。通俗的说，CFL 条件告诉我们，在模拟中的每个时间步长都需要满足一定 条件，才能使模拟的结果不出错。对这个条件有这样的理解：每一个时间步长内