The Material Point Method for Simulating Continuum Materials读书笔记

3 简介

MPM是结合PIC和FLIP而发展的一种算法。MPM对Lagrangian mesh连通性没有要求。

和PIC/FLIP类似,MPM算法在背景Eularian网格的辅助下,隐式处理自碰撞和破碎。与传统的Lagrangian方法(例如FEM算固体)Eularian方法(流体)相比,MPM好处如下:

  • 和FEM一样,MPM可以由弱动量守恒形式导出,因此求解精确。
  • 可以很容易在网格和粒子上施加边界条件、固壁碰撞和外力。
  • 自动解决自碰撞/接触。因为粒子的运动,是由网格上未被扭曲的节点运动插值而得。
  • 自动对物质进行分割/合并,因为物质是由粒子表示的。这对流体和颗粒状材料很有用。
  • 只要给粒子赋不同的性质或成分模型,就能自动解决多材质和多相耦合问题。
  • MPM也能用于模拟基于网格的Lagrangian力,而不失其优点(见12.2节)。因此可以同时耦合用点表示的物体,和用mesh表示的物体,只需要在网格上求解一次,并能自动处理碰撞。

严格来说MPM是Lagrangian方法,但有一个Eularian网格用来算导数。坏处是拓扑结构变化、碰撞、超弹性材料效果不佳。好处是不用单独处理碰撞和拓扑结构变化。

4 MPM在工业界的应用

《冰雪奇缘》里的雪是用MPM做的。《超能陆战队》和《疯狂动物城》用到了MPM。

5 运动学理论

每个MPM粒子实际代表一个体积微元,而不是什么质点,原子之类的东西。也可以把它看作FEM里面四面体的顶点(见第9章)。

5.1 连续物质的运动

记\(\mathbf{X}\)为“初始位置”,\(\mathbf{x}\)为“当前位置”。初始位置集合记为\(\Omega^0\),当前位置集合记为\(\Omega^t\),它们都是\(R^d\)的子集。\(\phi\)是二者之间的映射关系:

物质点 \(\mathbf{X}\) 的速度为

加速度为

这里的速度和加速度都是以Lagrangian方式定义的,也就是说在物质点上定义的。

5.2 形变

上面说的\(\phi\)又叫做形变映射(deformation map)。它的Jacobian可以描述一些性质,比如说弹性。这个Jacobian一般记作F:

F一般简称为形变梯度(deformation gradient),也就是一个2*2或者3*3矩阵。只有一个特殊情况:对于三维布料/薄壳,F是3*2矩阵,因为材料是二维的。可以把F看作一个映射,对于某个初始点\(\mathbf{X}\),它的\(F(\mathbf{X},t)\)是对应这块物质微元在时间t的形变梯度。

写开就是:

图1:形变梯度

F的行列式的含义是微元的体积变化,一般把这个行列式记作J。如果J>1意味着体积增大,J<1意味着体积减小。

J=0意味着体积变为零(退化)。J<0意味着“反转”,比如一个三角形,一个顶点越过对边。第6章会讨论“invertible elasticity”模型。

5.3 前推和拉回

物质的初始状态\(\Omega^0\)和当前状态 \(\Omega^t\)在\(\phi\)下同胚/微分同胚,因为 \(\phi\) 是光滑的。 \(\phi\) 存在一个逆映射。所以如果一个函数定义在\(\Omega^0\)上,就可以把它给“转移”到 \(\Omega^t\)(这叫‘前推/push forward’),当然反之也可(‘拉回/pull back’).例如:有一个函数\(G: \Omega^0\rightarrow R\),可以把它前推成\(g(\mathbf{x},t)=G(\phi^{-1}(\mathbf{x},t))\),而g可以拉回为\(G(\mathbf{X})=g(\phi(\mathbf{X},t),t)\).

有时把前推函数称为Eularian(定义在\(\mathbf{x}\)上),拉回函数称为Lagrangian(定义在\(\mathbf{x}\)上)。

速度和加速度都是Lagrangian的:

把它前推成Eularian的:

当然可以反过来写出拉回的结果:

使用a和v的定义可以写出:

单个分量就是:

最后一项是用了爱因斯坦求和约定,代表对j=1,2,3求和。

把(8)和(10)联立,用\(\phi\)表示Eularian速度:

联立(11)和(15):

可以发现

注:这里求Lagrangian速度用的是偏导数,但是求Lagrangian加速度用了全导数,所以加速度多出来一项。也就是说,Lagrangian速度和Eularian速度相同,但加速度不同。Lagrangian加速度是速度的物质导数。

5.4 物质导数

物质导数:

Lagrangian加速度是速度的物质导数:

对于任意Eularian函数,物质导数的定义是:

注意它是偏导数\(\frac{\partial}{\partial t}F\)的前推。这个F是Lagrangian的,是f的拉回。这一点很重要。

形变梯度一般认为是Lagrangian的。不过也可以把它前推为Eularian的:

这里重复下标k代表求和。推导如下:

最后一步是由对(12)求导而得。

(23)式对后面求每个粒子F的离散化更新很重要(9.4节)。

5.5 体积和面积变化

初始体积微元dV在一段时间后变成了dv=JdV,其中J=det(F)。可以用三维长度推一推得到这件事情。

所以在初始物质上积分和在一段时间后物质上积分满足关系:

为了便于理解写成:\(dx=JdX\).

二维有向面积:

大写S是初始状态,小写s是中间状态。

考虑一个线段微元\(dL\)(中间状态\(dl\)),和dS(或ds)确定一个体积:

根据上面说的\(dv=JdV\),所以

注意F是形变梯度。在微元内dl=FdL是对的。这对任何dL都成立。:

为什么是逆转置呢?把向量点积写成向量和向量转置相乘的形式就知道了。

如此一来,在二维表面上的积分形式如下:

其中H是h的拉回。

6 超弹性

Stress是应力,strain是应变。

超弹性材料的第一类Piola-Kirchoff应力张量P可以由应变能量密度函数(strain energy density function)\(\Psi(F)\)表示:

分量就是:

别忘了\(\Psi(F)\)是弹性势能密度函数(是一个标量函数)。P和它大小相同。

在工程领域经常用Cauchy应力张量\(\sigma\),它和P的关系是:

固体材料的表现就由形变映射\(\phi\)、Cauchy应力张量 \(\sigma\) (当然也可以用P)这两个量确定。

对于刚体而言,形变映射就是旋转加平移:

超弹性材料可以变形,但是会抵抗变形,这由\(\Psi\)描述。注意刚体的F=R,这是一个正交矩阵。所以对于超弹性材料,当F非正交的时候,会产生应力。正交矩阵,比如R的特性是\(R^TR=I\),所以我们可以写成:

这样看起来较为清楚,因为正交是\(F^TF=I\)的特殊情况,此时\(\tilde{\Psi}\)达到最小值.这个\(F^TF\)经常记为C(右Cauchy-Green张量)。如果材料是各向同性的,那么可以把能量写成\(F^TF\)特征多项式的三个系数(或者叫做各向同性不变量/isotropic invariants):

注:各向同性的意思,是不受坐标系旋转影响。也就是说乘以一个旋转矩阵结果一样。

这个结果可以简写为:

这里面\(F=U\Sigma V^T\)是F的“Polar SVD”.这里面U和V是旋转矩阵。而F的极分解F=RS可以通过\(R=UV^T\)和\(S=V\Sigma V^T\)得到 。这个R是距离F最近的旋转矩阵。S是对称阵。

各向同性不变量通常可以用奇异值表示。

6.2 新胡克模型

弹性势能密度函数:

d是维度,\(\mu\)和\(\lambda\)由杨氏模量E、泊松比\(\nu\)确定:

记得J是形变梯度F的行列式。如果F是旋转矩阵,\(\Psi=0\).如果没有反转(J>0),那么\(\Psi(F)\geq 0\).

为了计算力,需要给出P的形式:

为了隐式积分,需要给出\(\frac{\partial P}{\partial F}\),这个在6.4节讲。

6.3 改进Corotated本构模型(Constitutive Model)

这个模型是基于SVD定义的。“改进”的意思是,它是在计算机图形学常用的corotated线性弹性模型制上进行改进。做polar SVD分解\(F=U\Sigma V^T\),那么改进corotated模型为:

这里的\(\sigma_i\)是SVD分解里面\(\Sigma\)矩阵的对角元。J等于各个 \(\sigma_i\) 之积(行列式)。把右端\(\mu\)项展开:

可以证明:

其中F=RS是极分解。给定SVD分解\(F=U\Sigma V^T\)以后,\(R=UV^T\),\(S=V\Sigma V^T\).

因此可以计算出P:

二阶导数稍微麻烦一点。

后面的太麻烦,先不管。6.4节其余部分(包括6.4.1和6.4.2)暂略。

6.5 雪的塑性(暂略)

7 物理方程

Lagrangian方程:

还记得J是形变梯度矩阵F的行列式,代表体积变化。R是Lagrangian意义下的密度。它和Eularian意义下的密度\(\rho\)关系是:\(R(\mathbf{X},t)=\rho(\phi(\mathbf{X},t),t)\).这个\(\phi\)就是形变映射。注意质量守恒可以写作\(\frac{\partial}{\partial t}(R(\mathbf{X},t)J(\mathbf{X},t))=0\).

Eularian方程:

这里v是Eularian速度,D是物质导数:\(\frac{D}{Dt}=\frac{\partial}{\partial t}+v\cdot\nabla^{\mathbf{x}}\).

在Lagrangian方程中使用第一Piola-Kirchoff张量P更方便,在Eularian方程中则是Cauchy张量\(\sigma\)更方便。

7.1 质量守恒

\(\rho(\mathbf{x},t)\)是Eularian密度,R是Lagrangian密度(也就是\(\rho\)的拉回)。二者满足关系:

Eularian密度的定义:

这个\(B_{\epsilon}^t\)是t时刻,以\(\mathbf{x}\)为秋心,半径为\(\epsilon\)的球。注意这里面的体积微元实际上都是在初始状态意义上讲的。取极限就是:

注意\(J(\mathbf{X},0)=1\)(F是单位阵)。所以这个式子可以写成:

为了推导Eularian方程,首先把它展开:

继续展开:

其中使用了行列式对矩阵求导的公式:

其中这个\(F^{-T}\)的含义是\((F^{-1})^T\).

(99)第一个等号,脚标的意思是求和。第二个等号,是考虑到F是形变梯度\(\phi\)的Jacobian矩阵,也就是说\(F_{ij}=\partial\phi_i/\partial X_j\).而Lagrangian速度\(V=\frac{\partial\phi}{\partial t}\).所以这个\(\partial F_{ij}/\partial t\)就是\(\partial\phi_i/(\partial t\partial X_j)=\partial V_i/\partial X_j\).

第三个等号是基于5.2节的(10)式的,使用链式求导法则:\(\frac{\partial V_i}{\partial X_j}=\Sigma_k\frac{\partial v_i}{\partial x_k}\frac{\partial x_k}{\partial X_j}= \Sigma_k\frac{\partial v_i}{\partial x_k}F_{kj}\).当然(99)中用爱因斯坦求和约定省去了求和号。

第四个等号,我们考虑其中固定i,k,对j求和的三项,会发现实际上是\(F^{-1}\)的一列乘以F的一行,所以就是\(\delta_{ik}\),也就是i=k时为1,否则为0.

联立(97),(98),(99):

这个J实际可以约掉。然后对两项都做前推。第一项的前推就是\(\rho\)的物质导数,而R的前推就是\(\rho\),结果就是:

注意后面那项看起来很像物质导数的对流项,但其实不是。写开之后,实际上有三项。

7.2 动量守恒

作用在物质上的力可以分为体力(如重力)和面力(基于应力,stress-based)。

面力意味着有一个应力场(柯西应力张量):

t是面力。

动量守恒公式:

这个很好理解:体积微元内总动量的变化率等于总受力。第二个等号,用到了上面的质量守恒:\(R(X,t)J(X,t)=R(X,0)\).

(105)的最后一个等号可以写成:

(106)左端前推:

最后一个等号用了高斯公式。然后对微元取极限就是:

当然也可以对(106)右端做拉回。首先注意到:

其中\(F^{-T}\)这些是5.5节中推出来的边界积分的结果。

因为第一类Piola Kirchoff应力张量P和Cauchy应力张量的关系是\(P=J\sigma F^{-T}\),所以有:

结合(106),就得到了Lagrangian形式的动量守恒方程:

取极限:

这里面的\(F^{ext}\)是\(f^{ext}\)的拉回。

7.3 弱形式

首先忽略外力,动量守恒方程为:

也就是:

这里\(R_0=R(X,0)\).然后\(P_{ij,j}\)的意思是\(P_{ij}\)对第j坐标轴求导,当然用爱因斯坦约定隐含了对j求和。

对于任意函数\(Q(\cdot,t):\Omega^0\rightarrow R^d\),把它和(113)式做点积并积分起来,得到弱形式:

注意,这里似乎应该是只对j求和,不对i求和。

其中\(P_{ij}N_j\)应该被作为边界条件给出。设它是T(X,t),那么对于任意函数Q,应有:

在MPM中,应力导数是在当前状态下计算的,所以可以把涉及应力的积分前推至Eularian的。令q是Q的前推,即\(Q(X,t)=q(\phi(X,t),t)\),反过来\(q(x,t)=Q(\phi^{-1}(x,t),t)\),则:

回忆起Cauchy应力张量的定义:

令t表示Eularian上边界的单位面积受力(即T的前推),则(118)变为:

我们一项一项分析。

等号左端,\(\rho(x,t)J(X,t)=R(X,0)\),而\(dx=J(X,t)dX\),合起来就是\(\rho(x,t)dx=R(X,0)dX\).

右端第一项,本质上就是边界受力求和,是面积上加和,至于面积变化均摊的部分实际上是算到t里面了。

右端第二项就是上面说的,Cauchy应力张量。这里应该是对k求和,而不对i求和。

(120)对任意的q都成立。我们把(120)称作Eularian下的弱受力平衡形式。MPM基于这个式子做离散化。

8 物质点

在物质点上储存质量、速度、位置。但是应力是在Eularian网格上算的,所以需要一个转移的方法。

8.1 Eularian插值函数

SPH是在粒子周围定义核函数。MPM不一样,是在Eularian格点上定义核函数(插值函数)。在格点i上定义插值函数\(N_i(x)\)。这里的i其实应该是一个向量(原文中写的是黑体,我懒得每次mathbf了)。对于某个粒子点\(x_p\),常简记\(N_i(x_p)=w_{ip}\).

其中p是这个粒子,h是网格边长。之后简记:\(w_{ip}=N_i(x_p)\),\(\nabla w_{ip}=\nabla N_i(x_p)\).

当然这里核函数N的选择很关键。我们喜欢用张量积样条(tensor product spline).分段线性核函数对于FLIP很简单,但在这里不适用。原因有二([Steffen et al., 2008]):第一,如果用分段线性,\(\nabla w_{ip}\)就不连续了,产生不连续的受力。第二,用分段线性,则在离格子远处,\(w_{ip}\)本身很小,但\(\nabla w_{ip}\)却很大,这不科学。MPM中使用的核函数最好具有C1连续性。二次B样条计算更快,三次B样条计算慢一些,但是效果更好。三次核函数为:

h是网格边长。二次核函数为:

如图3所示。

图3.三次核函数(蓝),二次核函数(红)。

虽然理论上不稳定,但某些特定情况下也可以用线性核:

计算梯度的方法就是每一维单独算:

8.2 Eularian/Lagrangian质量

每个物质点代表一小块物质,质量就是这一块物质的质量。

把质量和动量(由此可以计算出能量)从物质点转移到Eularian网格上的方法就是:

\(m_i\)就是坐标为i的网格的质量。为了质量守恒,核函数应当满足单位分解(unity partition),质量守恒率写出来如下:

8.3 Eularian/Lagrangian动量

类似地,可以把粒子的动量\(m_pv_p\)转移到网格上(注:我们最后实际使用10.1节所述的APIC方法):

类似地,它同样满足守恒律:

用动量除以质量就是速度:

这里没有爱因斯坦求和约定。

把Eularian值转移到Lagrangian值的过程类似。不过不需要转移质量,因为Lagrangian质量总是不变的。速度可以通过插值得到:

动量守恒律:

满足动量守恒律的关键一步,是由于格点质量乃从Lagrangian质量插值而得。

9 离散化

9.1 时间离散化

根据之前的弱形式:

这个\(t^n\)是当前时间步。把加速度A离散化,左端前推:

这个能成立,要记得之前说过的,\(\rho dx=RdX\).另外式中\(v^n\)代表第n个时间步的速度。

为了不把代表第几维的i和代表网格坐标的i混淆,接下来我们用希腊字母\(\alpha,\beta,\gamma\)代表第几维。由此重写(134)为:

9.2 空间离散化

类似FEM,把之前的q和v都用格点插值代入:

适当省略后简写为:

则受力平衡式(135)化为:

即:

其中质量矩阵

中间积分项的含义是:用格点i在x处的核函数,乘以格点j在x处的核函数。

拉回得到:

这是一个对称半正定矩阵,它可以写成\(BB^T\),其中\(B_{ip}=\sqrt{m_p}N_{ip}\).但这个矩阵是奇异矩阵,我们使用mass lumping技术来构造一个正定对角矩阵来近似它。

(141)对所有q成立。我们取q为:

注意它是定义在格点上的。黑体i就是格点坐标。

则有:

这就是离散的受力平衡方程。这个可以作为显式Eularian更新动量的方程。

经常使用mass lumping简化。方法是把\(m_{ij}\)的每一行加起来放在对角线上。这样就是一个对角矩阵。格子i对应的行之和\(\hat{m_i}\)就是:

同样做简化,类似(143):

原因是\(m_p\approx V_p^0R(X_p,0)\).也就是说\(\hat{m_i}\)可被自然地近似为\(m_i\).

把\(m_iv_i^n\)写成动量\((mv)_i^n\),则离散方程近似为:

左边是动量的改变,右边近似是受力。假定我们已经对每个Lagrangian点\(x_p^n\)计算出了张量\(\sigma_p^n\approx\sigma(x_p^n,t^n)\),则:

其中\(V_p^n\)是\(B_{\delta x,p}^{t^n}\),也就是粒子p在时间步\(t^n\)占的体积。

用两种方法计算\(V_p^n\)。第一,注意到

可近似估计

第二种方法如下。物质点p的形变梯度矩阵近似为:\(F_p^n\approx F(X_p,t^n)\).那么其体积变化\(J_p^n=det(F_p^n)\).假如我们知道初始体积\(V_p^0\),则:

使用第一类Piola Kirchoff stress:\(\sigma=\frac{1}{J}PF^T\),则可以把(147)和(148)中的\(\sigma\)用P表示:

由于大部分本构模型都是指定P,所以这个写法很有用:

9.4 形变梯度更新

形变梯度较难计算。使用如下方程(见5.4节):

注意这个式子很神奇,把形变梯度的时间导数转换成了速度的空间导数。假设只有第n+1步的Eularian速度。 则有:

做时间离散化就是:

那么更新方程就是:

当然,把速度的插值公式写出来:

结果就是:

这就是给定\(v_i^{n+1}\)和\(F_p^n\)后更新\(F_p^{n+1}\)的方法。

9.5 作为势能梯度的力

对于超弹性材料,第一类Piola Kirchoff应力张量由弹性势能密度定义:

如第6章所述,某节点处的弹性势能为:

其中\(\hat{x}\)是所有Eularian格点\(\hat{x}_i\)拼接而成的向量。这里相当于是,临时认为Eularian网格也可以动,所以用的是Eularian格点。特别地:

那么粒子处的形变梯度更新方式为:

在\(\hat{x}_{i\alpha}\)处对势能求导:

此外,根据上面F的方程,求导有:

代入(167)得:

和9.3节的结果对比,发现网格受力为:

所以动量更新方程就是:

10 显式时间积分

10.1 APIC转移

可以用Affine Particle-In-Cell把物理量从粒子转移至网格。忽略代表时间的上标。

转移方程:

其中\(B_p\)是存在粒子处的一个矩阵。\(D_p\)为:

从网格转移至粒子的方程:

使用\(D_p\)可以轻松表示MPM常用的二次(\(D_p=\frac{1}{4}\Delta x^2I\))或三次(\(D_p=\frac{1}{3}\Delta x^2I\))插值权值。

如果使用线性样条(虽然不推荐),为了避免\(D_p\)成为奇异矩阵,需要改为:

而从网格到粒子的转移方程就是:

此时,\(C_p\)是\(x_p\)处的Eularian速度梯度。我们在更新\(F_p\)的时候已经计算过了。

10.2 形变梯度矩阵的更新

根据网格速度场\(v_i^{n+1}\),则形变梯度矩阵F的更新方程为:

这个就是9.4节的结果。

10.3 状态更新

显式时间更新:

根据(182),可以把节点位置想象成关于节点速度的函数:

那么:

即:显式的力就是假设节点没有移动时的力。使用(182),可以把(181)重写为:

这样\(F_p\)就是节点位置的函数。

10.4 力

MPM在格点上储存力。正如之前所述,对于超弹性固体,这个力既可以用动量守恒方程的弱形式推出,也可以用弹性势能梯度推出。后者比较简单。

首先假设我们基于网格定义了超弹性势能密度。总势能为

其中\(V_p^0\)是粒子p的体积。

节点上的弹性力等于势能的负梯度。由(186),形变应力为:

这是基于粒子计算节点处弹性力的方法。用辛Euler方程,可以写成:

完全取决于已有的插值权值,以及粒子的物理量。如果不用能量密度,也可以用Cauchy应力张量表示这个力:

10.5 MPM整体流程

下面描述基于辛Euler方法的MPM。假设已经初始化了所有粒子(赋予质量\(m_p\),体积\(V_p^0\),初始位置\(x_p\),初速度\(v_p\),仿射矩阵\(B_p\)和其他材料相关参数)。

  1. 物理量从粒子转移至网格(Particle to Grid, P2G)。使用10.1节所述的APIC公式,把粒子上的物理量(即:质量和动量)转移至网格上。
  2. 计算格点速度。\(v_i=\frac{m_iv_i}{m_i}\)。对于质量为0的格点,把速度设为0.一般认为会在这种地方设置浮点数阈值,用来避免除以零。不过其实直接和浮点数0.0f相比也没什么问题,因为就算质量很小,导致除法运算得到一个很大的速度,但最后向粒子转移的时候,也要被这个质量再放缩回去。实际上在这里设置浮点阈值反而有害,会导致总动量不守恒之类的问题。
  3. 确定自由度。这一步对于效率很关键。把质量非零的格点标记出来。其他节点没有自由度,保持不变,不在求解器中作为未知量出现。
  4. 使用(189)式,显式计算格点上的力\(f_i^n\)。
  5. 使用(183) 式更新网格速度。这一不需要考虑边界条件和物体碰撞。对于显式积分而言,每个格点处的速度都可以被Dirichlet边界条件,或者刚体碰撞来单独设置。处理碰撞的细节具体见12.1节。
  6. 使用(181)式更新粒子的形变梯度。注意我们并不实际把格点移到计算出来的新的位置,这只是一种假设,用来计算速度。
  7. 把物理量从网格转移至粒子。这一步使用10.1节所述方法,计算出新的粒子速度\(v_p^{n+1}\)和新的仿射矩阵\(B_p^{n+1}\)。
  8. 粒子advection。基于速度advect粒子:\(x_p^{n+1}=x_p^n+\Delta tv_p^{n+1}\)。注意只有APIC是这么advect的。在FLIP或者混合PIC/FLIP方法中,advect方程是\(x_p^{n+1}=\Sigma_ix_i^{n+1}w_{ip}\)。(对于PIC或者APIC,这两个式子相同)。

11 隐式积分

隐式和显式的主要区别是更新网格速度场。

11.1 力的导数

总的弹性势能可以由势能密度积分而得:

对应力做空间离散化,等于在格点上对这个势能密度做离散近似,然后求导。当然,我们并不真正改变格点的位置。设格点i的位置是\(x_i\),那么它下一步“应该”在的位置是\(\hat{x_i}=x_i+\Delta tv_i\),为简单起见,略去下标i,那么总势能就是:

这个\(V_p^0\)就是粒子p初始时占据的体积,\(\hat{F}_p\)为:

那么应力离散化为:

这个\(f_i(\hat{x})\)就是格点i的弹性力。写成Cauchy应力张量\(\sigma_p=\frac{1}{J_p^n}\frac{\partial\Psi}{\partial F}(\hat{F_p(\hat{x})})(F_p^n)^T\):

这个\(V_p^n=J_p^nV_p^0\),就是时间步n时粒子p的体积。

我们使用势能的Hessian矩阵来做隐式更新。对于任意小扰动\(\delta u\)的Hessian矩阵为:

其中

这里面冒号的含义是:A=C:D意味着\(A_{ij}=C_{ijkl}D_{kl}\)(隐含对kl的求和)。

也可以用格点的方式推导。注意这个力在(188)中已经被推导出来,求导有:

其中\(\frac{\partial^2\Psi}{\partial F^2}\)可以由\(\Psi(F)\)导出。

11.2 反向Euler系统

和(183)式(\(v_i^{n+1}=v_i^n+\Delta tf_i(x_i^n)/m_i\))不同,反向Euler的更新方式为:

其实就是把第n步的受力换成了第n+1步的受力。也就是说,受力会隐式地依赖于网格的运动。这允许我们把时间步变得更大。

重排后得到一个非线性的格点速度公式:

11.3 牛顿法

可以用牛顿法解(200)式。初始猜测可以取上一步的速度。然后迭代步为:

每一步都是求解一个线性系统,这可以用诸如MINRES的Krylov求解器。在牛顿迭代中,\(F_p\)需要根据速度更新。需要储存上一个时间步的\(F_p^n\),以备迭代使用。

在许多情况下,只需要做一步牛顿迭代,效果就非常好,能以比显式方法多出数量级的时间步做到相同效果。这几乎和所谓“半隐式”MPM积分效果相同。

11.4 线性力

弹塑性响应可以由网格位置的虚拟“更新”确定。把这种虚拟更新写成\(\hat{x}\)关于v的函数:\(\hat{x}=\hat{x}(v)\)。我们使用如下记号:\(f_i^n=f_i(\hat{x}(0))\),\(f_i^{n+1}=f_i(\hat{x}(v_{n+1}))\),\(\frac{\partial^2 e^n}{\partial\hat{x}_i\hat{x}_j}=-\frac{\partial f_i^n}{\partial \hat{x}_j}=-\frac{\partial f_i}{\partial \hat{x}_j}(\hat{x}(0))\).

由此写出隐式速度更新:\(v_i^{n+1}=v_i^n+\Delta tm_i^{-1}((1-\beta)f_i^n+\beta f_i^{n+1})=\approx v_i^n+\Delta tm_i^{-1}(f_i^n+\beta\Delta t\Sigma_j\frac{\partial f_i^n}{\partial \hat{x}_j}v_j^{n+1})\).这是一个求解\(v_i^{n+1}\) 的对称方程组 :

其右端就是:

参数\(\beta\)取0时就是显式,取0.5时就是梯形公式,取1就是反向欧拉。这相当于以全0场为初值,做一次牛顿迭代。因此,这样每步只需要解一个线性系统。MPM半隐式时间步见图4。

图4:隐式时间积分

11.5 基于优化的时间积分

Newton-Raphson求解器性能比显式时间步好得多,但仍然需要较小的时间步才能达到稳定。

如果力是由势能函数得到的(这在大多时候成立),那么就可以让牛顿法在和CFL条件相近的任意\(\Delta t\)内收敛,方法是把远方程写成一个优化方程。

实际上,求解:

等价于最小化:

这允许我们采取大得多的时间步,节省大量运算资源。具体情况见[Gast et al., 2015].

12 其他

12.1 物体碰撞

在使用受力更新网格速度后,我们就处理碰撞问题。对于半隐式积分,这涉及到修改线性方程组的右端,并且消除涉及到碰撞的格点的自由度。可以在更新粒子坐标之前,对粒子速度\(v_p^{n+1}\)做一次碰撞检测,用来解决插值引入的新误差。所有碰撞都是非弹性的。注意基于粒子的碰撞会让粒子上的形变梯度张量自相矛盾。因此在模拟中,只应该在必要情况下使用这种碰撞。

我们使用level set表示碰撞的物体。那么碰撞检测就很简单(\(\phi\leq 0\))。首先计算出相对速度\(v_{rel}=v-v_{co}\),这个co是碰撞物体的意思。如果两个物体互相分离,那么不加碰撞。记\(v_t=v_{rel}-nv_n\)为相对速度的切向分量。如果需要一个“sticking impulse”(\(||v_t||\leq-\mu v_n\)),那么直接设\(v’_{rel}=0\),这个撇代表已经处理过碰撞了。否则就用使用动态摩擦,更新为\(v’_{rel}=v_t+\mu v_nv_t/||v_t||\),其中\(\mu\)是摩擦系数。最后基于相对速度去更新速度:\(v’=v’_{rel}+v_{co}\).

我们有两种碰撞物体:刚体和软体。刚体就储存一个固定的level set,和刚体的位置参数,用来计算每个点的\(\phi\),法向和物体速度。对于软体,在关键帧加载level set,并使用\(\phi(x,t+\gamma\Delta t)=(1-\gamma)\phi(x-\gamma\Delta tv_{co},t)+\gamma\phi(x+(1-\gamma)\Delta tv_{co},t+\Delta t)\)做时间插值。但速度则用\(v_{co}=(1-\gamma)v(x,t)+\gamma v(x,t+\Delta t)\),忽略位置变化。

最后考虑粘性碰撞。此时Coulomb摩擦是不够的,因为法向相对速度应该是0(垂直)或者正数(悬挂在下面,然后由于重力又分开)。这里我们无条件地设置\(v’_{rel}=0\)。也就是说,这等价于Dirichlet边界条件。

12.2 Lagrangian力

在传统MPM模拟中,使用(188)式计算力。在粒子处储存形变梯度有一定好处,例如不用考虑拓扑。但是在用MPM模拟传统的软体,比如橡胶之类的东西时,使用基于mesh的FEM求解器更有道理,因为它提供了连通性信息,能更精确地计算形变梯度矩阵,并且有利渲染。在此情况下,MPM可以和基于mesh的方法轻松整合在一起。

我们可以使用任何Lagrangian受力模型(例如弹簧,有限元等等),只要能写出每个点处的总势能\(W(x_p)\)。基于这个模型计算出受力\(f_p=-\frac{\partial W}{\partial x_p}\)。设对于任何粒子处的位移\(\delta u_q\),都可以计算相应受力\(\delta f_p=\Sigma q\frac{\partial f_p}{\partial x_q}\delta u_q\).这些是纯Lagrangian的,可以用常规方法做,我们就不详述了。

然后考虑如何把它们加到网格上。我们需要根据粒子处的受力\(f_p\)计算处网格受力\(f_i\),并根据\(\delta u_i\)计算\(\delta f_i\)。考虑之前粒子-网格之间互相转移的方法,我们得到\(x_p=\Sigma i w_{ip}^nx_i\).使用链式法则:

这里有三个隐含求和,不过计算起来并不难,可以一个一个做。这些是网格上的受力,那么MPM和Lagrangian方法就可以整合起来。把每个粒子标记为MPM粒子或者是mesh粒子。我们不会使用储存在mesh粒子上的形变梯度矩阵,因为直接在mesh上计算即可。这就能把Lagrangian方法的精确表面追踪和Eularian网格的自动处理碰撞给结合在了一起。

参考文献(略)

发表评论

电子邮件地址不会被公开。 必填项已用*标注