【游戏流体力学基础及Unity代码（十五）】线性有限元及弹性物体模拟

BlaXuan · 发表于 2021-5-9 15:21

上面就是我用有限元方法模拟果冻的结果。更多种类的果冻请看这个视频。代码在github上，或者是gitee上
数值积分

Galerkin法中需要算权函数的积分，然而手算太慢。因此可以用数值积分来算。积分就是算面积，将曲线下不规则的区域用长方形去拟合，就是数值积分。
此时我们用到的数值积分主要分为三种，第一种，这些长方形的底部宽度是相等的。比如两阶精度的梯形公式，三阶精度的simson公式等。比如三阶精度的公式如下。其中h就是这些等距节点间的距离，f0,f1,f2就是取值节点。

$\int_{x0}^{x2}f(x) = \frac{h}{3}(f_0 +4 f_1 + f_2) + O(x^3) \tag{1}$
比如f(x) = x^2，要从0积分到2。那么取三个点x=0,1,2，它们之间的距离为1。那么数值积分结果就是

$\int_{0}^{2}x^2 = \frac{1}{3}(0 +4 + 4) = \frac{8}{3} \tag{2}$
于此同时五阶精度与七阶精度的数值公式分别如下

$\int_{x0}^{x4}f(x)dx = \frac{2h}{45}(7f_0 + 32f_1 + 12f_2 + 32f_3 + 7f_4)\tag{3}$

$\int f(x)dx = h(\frac{41}{140}(f_0+f_6) + \frac{54}{35}(f_1+f_5) + \frac{27}{140}(f_2+f_4) + \frac{68}{35}f_3) \tag{4}$
这种方法的证明请看《数值方法(matlab版)》by John H.Mathews第7.2节。直接用搜索引擎搜索一些大学的讲义也可以，但是不容易搜到且质量参差不齐，最好还是看这些完整的教科书。
第一种方法各节点之间的距离都是相等的。也可根据勒让德多项式(legendre polynomials)，推导出另一种方式的积分，也就是高斯积分。直接贴结论了。同样来源于上面那本书。下图左侧的横坐标xn，就是勒让德多项式的根。
同样算f(x) = x^2，要从-1积分到1。选的两个点为-0.57735和0.57735，计算结果如下

$\int_{-1}^{1}x^2 dx \approx (-0.57735)^2 + (0.57735)^2 = 0.66666604 \tag{5}$ N点高斯积分的精度为2N-1，这意味着2点高斯积分算x^3的积分也能算很准确的。比如x^3从零积分到2，那么节点位置也得改一下

$\int_0^2x^3dx \approx (1-0.57735)^3 + (1.57735)^3 \approx 4 \tag{6}$ 第三种方法为蒙特卡洛积分(monte carlo integral)。前两种方法都是选特殊位置的节点，然后计算函数值，函数值还有独特的权重。蒙特卡洛积分就比较放飞自我，就是随便选些点，然后每个节点的函数的权重都是一样的。比如算x^2从a=0积分到b=2的值，用蒙特卡洛积分写成python代码如下
import numpy as np
def f(x):# 需要积分的函数
return x*x
a = 0 # 积分下边界
b = 2 # 积分上边界
n = 100 # 随机选择点的数量
summ = 0
r <span class="o">= np.zeros((n))
for i in range(0,n):
r = np.random.rand(1)
r = r*(b - a) + a # 随机产生a到b之间的值
summ += f(r)*(b-a)/n # 每个节点函数的权重为(b-a)/n用上面这段代码运行的时候，尽管我选了100个点，但最终结果仍然随机处于2.6~3之间，与上面两种方法比起来费时费劲还不准确。但在对面图形学中的光照渲染领域，蒙特卡洛积分经常被用在路径追踪(Path Tracing)等光照算法中。因为光照算法中，那个大名鼎鼎的渲染方程(Rendering Equation)积分起来实在太复杂了。并且我开始学光照渲染的时候，很多教程都是直接把蒙特卡洛积分扔上来，很少提到其它的数值积分方法。然而当时的我凭借着仅有的高等数学知识，理解起来也非常困难并且进展很慢。但现在我学了一大堆偏微分方程数值解法，在学到有限元的时候突然发现还有数值积分这么个玩意，再把新学到的等距数值积分和高斯数值积分，与之前一直无法理解的蒙特卡洛积分对比一下，很多东西就豁然开朗了。
多项式插值

上一篇算Galerkin法的权函数的时候，用的是拉格朗日多项式(Lagrangian polynomials)。如果要拟合二次函数，那么就要选三个权函数，它们的形式如下

$L_1(x) = \frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}\\ L_2(x) = \frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}\\ L_3(x) = \frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}\tag{7}$
取x1 = -1,x2 = 0,x3 = 1得到

$L_1(x) = \frac{1}{2}x^2 -\frac{1}{2} \qquad L_2(x) = 1-x^2 \qquad L_3(x) = \frac{1}{2}x^2 + \frac{1}{2} \tag{8}$
这是一组正交多项式，也就是分别用非零常数乘上其中某些多项式，然后将它们加起来永远也无法得到零。
之前说到的勒让德多项式，我最初在球谐光照(Spherical harmonic lighting)看到，看得也是一头雾水。现在再看看其多项式前三项

$P_1(x) = 1 \qquad P_2(x) = x \qquad P_3(x)=\frac{3}{2}x^2 - \frac{1}{2} \tag{9}$ 而之前高斯积分的节点，比如0.57735，就是P3(x)的根。因为它是正交的，因此可以作为基底来复合别的多项式，比如

$f(x) = x^2 = \frac{1}{3}P_1(x) + \frac{2}{3}P_3(x) \tag{10}$
流体力学中经常把流体的速度，压力等当成多项式。而用了渲染方程的光照中，则是把光照密度当成多项式，这就是为什么球谐函数中要用到连带勒让德多项式(Associated Legendre polynomials)了。
常用的还有切比雪夫多项式(chebyshev polynomials)和Hermite多项式。这些多项式各有特点，但它们都可用在Galerkin法中去算积分。在学Galerkin法前，我还不知道用这些拟合多项式有什么意义。如果我一直埋头学光照渲染，而不去看看数值偏微分数值解法的话，我可能还是不知道这些多项式意义何在。
线性扩散方程编程

接下来继续写代码。代码文件为Diffusion1d.py。一维线性扩散方程如下

$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} \tag{11}$ 把上面等式左边用Galerkin法展开

$\iint \begin{bmatrix} \phi_1 \\ \phi_2 \\ ..\end{bmatrix} \begin{bmatrix} \phi_{1t} & \phi_{2t} & ..\end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\..\end{bmatrix}dtdx \tag{12}$
我一共用了9个权函数
def pf(n,t,x):
if n == 0:
      return (t*t/2 - t/2)*(x*x/2 - x/2)
elif n == 1:
      return (t*t/2 - t/2)*(1 - x*x)
elif n == 2:
      return (t*t/2 - t/2)*(x*x/2 + x/2)
elif n == 3:
      return (1 - t*t)*(x*x/2 - x/2)
elif n == 4:
      return (1 - t*t)*(1 - x*x)
elif n == 5:
      return (1 - t*t)*(x*x/2 + x/2)
elif n == 6:
      return (t*t/2 + t/2)*(x*x/2 - x/2)
elif n == 7:
      return (t*t/2 + t/2)*(1 - x*x)
elif n == 8:
      return (t*t/2 + t/2)*(x*x/2 + x/2)我的数值积分是这么算的。
def bool5(point0,point1,point2,point3,point4,h):#五阶精度数值积分
return 2*h/45*(7*point0 + 32*point1 +
               12*point2 + 32*point3 + 7*point4)
for i in range(0,8):
for j in range(0,8):
      for xr in range(0,5):
         x0 = xr/2 - 1
         for tr in range(0,5):
            t0 = tr/2 - 1
            # pf(i,t0,x0) 是算phi
            # (pf(j,t0+1,x0) - pf(j,t0-1,x0))/2 是用中心差分算phi对t求导
            Inte0[tr] = pf(i,t0,x0)*(pf(j,t0+1,x0) - pf(j,t0-1,x0))/2
         Inte1[xr] = bool5(Inte0[0],Inte0[1],Inte0[2],Inte0[3],Inte0[4],0.5)
      Tmat[i,j] = bool5(Inte1[0],Inte1[1],Inte1[2],Inte1[3],Inte1[4],0.5)首先，(phi)乘(phi对t求导)最后算出的矩阵是8x8的。因此在上面的代码中Tmat也是8x8的。
其次我选了五阶精度的bool数值积分，是因为算的是(phi)乘(phi对t求导)，但(phi)乘(phi对t求导)的x都是二次，一共四次，所以精度需要五阶。这种数值积分需要5个点，由于是从-1积分到1，因此这五个点分别是-1，-0.5，0，0.5，1。它们之间的间隔是0.5。但是xr的五个点为0,1,2,3,4，所以x0 = xr/2-1才是真正的积分点。
而且这是二重积分，所以要算出dtdx，8x8矩阵中的每个元素，都要由25个值计算得到。这算起来很麻烦。不过矩阵中有些值是相同的，可以继续化简。
其次，算(phi对t求导)时，我直接用了中心差分，因为权函数中t都是二次的，二阶精度的中心差分足够了。当然也可以之间把权函数求导之后当成函数写出来
def pfdt(n,t,x):
if n == 0:
      return (t-1/2)*(x*x/2 - x/2)
elif n == 1:
      return (t-1/2)*(1 - x*x)
elif n == 2:
      return (t-1/2)*(x*x/2 + x/2)
...
...
Inte0[tr] = pf(i,t0,x0)*pfdt(i,t0,x0)最后部分结果如下
好了，手算一下最左上角的验算一下

$\iint \phi_0 \phi_{0t} = \int_{-1}^1dx\int_{-1}^1(\frac{1}{2}t^2-\frac{1}{2}t)(\frac{1}{2}x^2-\frac{1}{2}x)(t-\frac{1}{2})(\frac{1}{2}x^2-\frac{1}{2}x)dt = \frac{4}{15}*\frac{-1}{2} \tag{13}$
然后算d^2 u/dx^2的时候，
Inte0[xr] = -(pf(i,t0,x0+1) - pf(i,t0,x0-1))/2*(pf(j,t0,x0+1) - pf(j,t0,x0-1))/2
...
Inte1[tr] = Inte1[tr] + pf(i,t0,1)*(pf(j,t0,2) - pf(j,t0,0))/2
               - pf(i,t0,-1)*(pf(j,t0,0) - pf(j,t0,-2))/2 对应的公式如下

$\int_{-1}^{1}((\begin{bmatrix} \phi_{1} \\ \phi_{2} \\ .. \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} & .. \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ .. \end{bmatrix})_{-1}^1- \int_{-1}^1 \begin{bmatrix} \phi_{1x} \\ \phi_{2x} \\ .. \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} & .. \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ .. \end{bmatrix}dx)dt \tag{14}$
然后组装矩阵的时候
for k in range(0,total):
# 组装矩阵
colu = int(k/mmax) # t轴
row = int(k%mmax) # x轴
for t0 in range(0,3):
      for x0 in range(0,3):
         idx = (colu + t0)*mmax + row + x0
         idxx = int(t0*3 + x0)
         idxk[idxx] = idx我是这么写的。比如要求解x轴4个t轴4个一共16个值，如下
但权函数只有九个，一次只能求解九个相邻顶点，因此需要组装。第一次搞定0,1,2,4,5,6,8,9,10这九个顶点。第三个1,2,3,5,6,7,9,10,11这九个，第三次4,5,6,8,9,10,12,13,14这九个，依次类推。然后刚度矩阵具体计算方式如下
for i in range(0,num):
      for j in range(0,num):
         i0 = int(idxk)
         j0 = int(idxk[j])#u1的索引
         kmat[i0,j0] = kmat[i0,j0] + Tmat[i,j] - Dx2mat[i,j] 然后设置好初始条件，然后简单地求逆就行了。最终结果如下。横着的是x轴，竖着的是t轴。我让x轴上左右两个端点分别一直为0和8。然后经过11个时间步
不过我觉得这不能叫迭代，因为有限元方法是一次性把所有时间所有地方的所有变量全部算好的。到了第80个时间步的时候，就稳定下来了。它确实是线性的，1.14乘上7就等于8。
好了，一维扩散方程不带劲，来算一次二维的线性扩散方程吧。如下

$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \tag{15}$
完整代码文件为Diffusion2d。
三角形网格

上面用的都是正方形网格。不过这样怎么能突出有限元的优势呢？所以接下来用三角形网格解决问题。三角形有三个顶点，如果要用它来拟合一个二维一次函数，那么权函数如下

$\phi_1 = 1 - \xi - \eta \qquad \phi_2 = \xi \qquad \phi_3 = \eta \qquad \xi + \eta < 1 \tag{16}$
比如对于如下的三角形，从顶点1到顶点2，xi从0变化到1。从顶点1到顶点3，eta从0变化到1。
假设顶点1的坐标为(0,0)，值为0，顶点2的坐标为(1,0)，值为2，顶点3的坐标为(0,1)，值为1。那么顶点2和顶点3中间的那点坐标为(0.5,0.5)，xi = eta = 0.5，因此值等于

$(1-0.5-0.5)*0 + 0.5*2 + 0.5*1 = 1.5 \tag{17}$
也刚好符合值为2x + y的形式。但是三角形的顶点不可能都正好在网格上，所以要先将其转到xy坐标上

$\begin{bmatrix} d\phi/d\xi \\d\phi/d\eta \end{bmatrix} = \begin{bmatrix} dx/d\xi & dy/d\xi\\dx/d\eta & dx/d\eta \end{bmatrix}\begin{bmatrix} d\phi/dx \\d\phi/dy \end{bmatrix} \tag{18}$
计算雅可比矩阵，其中Area就是三角形的面积

$J = \frac{\partial(x,y)}{\partial (\xi,\eta)} = \begin{bmatrix} x2 - x1 & y2 - y1 \\x3 - x1 & y3 - y1\end{bmatrix} \qquad |J| = 2|Area| \tag{19}$
那么

$\frac{\partial \xi}{\partial x} = \frac{1}{|J|}\frac{\partial y}{\partial \eta} \qquad \frac{\partial \eta}{\partial x} = -\frac{1}{|J|}\frac{\partial y}{\partial \xi}\\ \frac{\partial \xi}{\partial y} = -\frac{1}{|J|}\frac{\partial x}{\partial \eta} \qquad \frac{\partial \eta}{\partial y} = \frac{1}{|J|}\frac{\partial x}{\partial \xi} \tag{20}$
最后

$b_1 = y_2 - y_3 \qquad b_2 = y_3 - y_1 \qquad b_3 = y_1 - y_2 \\ c_1 = x_3 - y_2 \qquad c_2 = x_1 - y_3 \qquad c_3 = x_2 - x_1 \tag{21}$
那么最终变换到xy坐标轴的权函数如下

$\begin{bmatrix} \partial \phi/\partial x \\ \partial \phi/\partial y\end{bmatrix} = \frac{1}{|J|}\begin{bmatrix} b_2 & b_3 \\ c_2 & c_3\end{bmatrix}\begin{bmatrix} \partial \phi/\partial \xi \\ \partial \phi/\partial \eta\end{bmatrix}\tag{22}$
如果用它来解下面这种偏微分方程

$\frac{\partial U}{\partial x} + a\frac{\partial U}{\partial y} + b = 0 \tag{23}$
那么最后手算Galerkin法结果如下

$\frac{b_2 + ac_2}{2A}\begin{bmatrix} -1/6 & 1/6 & 0 \\-1/6 & 1/6 & 0 \\-1/6 & 1/6 & 0\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3\end{bmatrix} + \frac{b_3 + ac_3}{2A}\begin{bmatrix} -1/6 & 0 & 1/6 \\-1/6 & 0 & 1/6 \\-1/6 & 0 & 1/6\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ u_3\end{bmatrix}+b\begin{bmatrix} 1/6 \\ 1/6 \\ 1/6\end{bmatrix} = 0 \tag{24}$
无下标的a和b是偏微分方程中的系数，而带下标的b和c则和三角形三个顶点坐标位置有关。其计算方法在之前已经写出。A则是三角形的面积。将其上面这堆东西写出代码就成了triangle.py。为了简便起见，这份代码仅仅是将之前的正方形网格拆成了左下和右上两个三角形。
网格划分

现在我们要针对不同的形状开始画网格了。比如求一块非常平的方形平板的传热问题，那么直接用正方形网格简单又快捷。但是如果平板上有个圆形的洞，或者平板边缘是弯曲的，那么无论如何摆放正方形网格，都不可能完全覆盖到平板上每一处区域。这时候三角形网格就非常有用了。
比如对于上面有三个洞的三角形平板，划分三角形网格可如下。
这里主要介绍波前推进法(Advancing Front Method)。这种的方法需要先确定一些前线，然后遍历每条前线生成三角形。比如下图深蓝色线就是所谓前线。
现在要生成第一个三角形，那么选择第一个前线，它的端点是1和2。然后计算这条线段的法向量与中点，用于生成顶点p。然后连接顶点1p和p2。然后线段1p和线段p2就是新前线，而线段12则不再是前线。
有个小问题，如何保证线段12的法向量指向内侧呢？我的方法，绘制外围顶点的时候，是逆时针。而绘制平板内的顶点，比如画那些洞的时候，用顺时针。规定好方向后，就能方便地算出法向量了。
// 平板外围的圆
point.Add(new Vector3(CircleCenter.x + CircleRadius * Mathf.Cos(2 * Mathf.PI / CircleParticle * i), 0.0f, CircleCenter.y + CircleRadius * Mathf.Sin(2 * Mathf.PI / CircleParticle * i)));
// 平板内的圆
point.Add(new Vector3(CircleCenter2.x + CircleRadius2 * Mathf.Cos(2 * Mathf.PI / CircleParticle2 * i), 0.0f, CircleCenter2.y - CircleRadius2 * Mathf.Sin(2 * Mathf.PI / CircleParticle2 * i)));
此时，我们还需要记录一个特殊的数组PointFront，它的作用是统计每个顶点所连接的前线数量。由于一开始生成的前线线段是闭合图形，所以顶点1到9的值为2。按上面那种方法新生成的p点的值也为2，而顶点1和2的值不变。
好了，波前推进法的核心部分已经完成了。接下来处理各种特殊情况。第一种，处理顶点为2和3的前线时，按照之前的方法，新形成的顶点p，与原有的某一顶点相距太近。比如下图新形成的顶点与顶点4相距太近，那么如果不处理的话，下次就会生成顶点为3p4的三角形，这个三角形太细长，虽然是在二维平面中，但却几乎只剩下一个维度的信息。所以尽量不要太细长的三角形，而是尽量生成等边三角形。因此对付这种情况的策略就是，直接让4为新三角形的顶点。
如果有很多个顶点的话，那么就选择离顶点p最近的那个就行了。此时顶点24的PointFront值不变，顶点3的值减2。
不过这紧接着又导致另一个问题，如下，推进线段12的时候，生成了顶点3，以及三角形123。而推进线段45的时候，由于上面的策略，没有新生成顶点，而是选择3作为新三角形的一个顶点。这样就多出来一个空隙。
此时线段32是前线，线段24是前线，线段43也是前线。按照之前的方法，这里会生成三个重叠的三角形，但显然一个三角形就足够了。达成这个目标有很多种方法，我的方法是生成第一个三角形的时候，判断三角形三个边是否都是前线，如果都是的话，那么生成第一个三角形，将它们的前线标志都取消。为了找到这些边，我用的是链式前向星。而且我在添加顶点的时候非常注意顺序，对于线段12，生成的时候，起点是1，终点是2。生成新三角形时，添加两条新边，第一条起点1，终点3。第二条起点3，终点2。所以下面的程序能够准确无误执行。同时三个顶点的PointFront的值都减2。
// 比如对于上面这种情况，st是3，ed是2
int st = EdgeFrom[EdgeIndex];
int ed = EdgeTo[EdgeIndex];
for (int i = Head[ed]; i != -1; i = EdgeNext)
{
if (EdgeFront == false) continue; // 如果不是前线，直接跳
int to = EdgeTo; // 寻找2为起点的线段，这里会找到顶点4，也就是to = 4
for (int j = Head[to]; j != -1; j = EdgeNext[j])
{
      if (EdgeFront[j] == false) continue;
      if (EdgeTo[j] == st)
      { // 寻找4为起点的线段，如果能够连到st = 3，那么说明可以形成一个三角形
         Triangles.Add(st);Triangles.Add(to);Triangles.Add(ed);
         // 然后将这些线段的前线标志去掉
         EdgeFront = false; EdgeFront[j] = false;EdgeFront[EdgeIndex] = false;
         EdgeIndex++;
         return;
      }
}
}
上面这种情况是三个前线刚好组成三角形的。而如果新生成的三角形的边原来有两条是前线的话，那么除了去掉那两条前线的标志外，还需要新增一条线。如下图，对于前线线段23，起点2，终点3，要生成的三角形的第三个顶点为4，那么就要新增一条前线线段，起点为2，终点为4。此时除了内侧顶点3的PointFront值减2外，其它的不用变化。
最后一种难处理的情况就是下面这种凹块了。深蓝色为已经画好的三角形。起初线段12，线段23，线段34都是前线线段。然后画第一个三角形132，那么仅剩下线段34是前线了。然后按照之前的算法，基本上都会画出三角形324。与原来的三角形交叉重叠，显然不能要。
这时候，前面说的PointFront就派上用场了。在画完三角形132后，顶点2显然不再连接任何前线线段了，因此PointFront值也变成了零。那么之后我们只要禁止有新的三角形顶点为这些PointFront值为零的顶点就行了。
不过我的代码还是会出一点bug，那就是有时三角形新生成的点会进入到别的三角形里面去，特别是在上面那样的凹块处。我能想到的解决办法，计算量都非常大。因此暂时不打算解决。
再细节处请看代码吧。接下来是欣赏时间。
代码在这里。以及一个慢速视频。我觉得如果把网格划分的步骤可视化出来其实挺好看的，像多米诺骨牌那样，令人极度舒适。然而现在的平板只有二维的圆形方向三角形，颜色也仅仅是白色。如果弄出一些不一样的形状，设计一些特殊的网格生成顺序，再加点炫酷的特效，就可以扔到ArtStation或ShaderToy或B站上拿赞去，或者说得高端一点叫拿正反馈。
二维网格划分除了推进波前法，各种各样的Delaunay三角形剖分算法也可以用在网格划分上，就是那个用来画Vornoio图的三角形算法。这种算法的原理在互联网上很容易搜到，比如https://www.habrador.com/tutorials/math/11-delaunay/ ，这是一系列很不错的在unity上处理几何问题的教程，配套视频https://www.youtube.com/watch?v=Z-1ExrWMTfA，配套代码https://github.com/Habrador/Computational-geometry。
《面向科学计算的网格划分与可视化技术》作者：王成恩
https://www.osti.gov/servlets/purl/1394106 有关网格划分一份不错的讲义
https://github.com/SpeedyF/advancing-front-method_2D 此库中的那篇pdf写得很好
有限元方法解决弹性力学问题

下面这一堆力学公式可以在大部分有限元入门书籍中找到。因此下面的内容可以仅仅说是将大部分有限元入门内容的程序写出来且扔到unity中可视化了。
如果你写过布料模拟的代码，用过弹簧质点系统，你会发现弹性力学与其差不多。最常见的用于弹簧的胡克定律是

$F = kx \tag{25}$
其中F是力的大小，k是弹性系数。x是形变。而弹性力学中经常用广义胡克定律。

$\sigma = E \varepsilon \tag{26}$
其中E是杨氏模量。sigma是应力(Stress)。应力中的正应力就是流体力学中常说的压力。应力等于力的大小F除以截面积A。varepsilon是应变(Strain)，等于物体受到力的作用后伸缩量u除以原长度。也就是

$\sigma = \frac{F}{A} \qquad \varepsilon = \frac{u}{L} \tag{27}$
如果要测量某种材料弹性的大小，应该是算力与长度变化量的关系，也就是F/u。但是对于同种材料，如果材料长度不同，截面积不同，那么F/u曲线也会变化。所以有个更方便的参数来测弹性，也就是杨氏模量E，等于应力除以应变，如下。

$E = \frac{\sigma}{\varepsilon} = \frac{FL}{Au} \qquad \frac{F}{u}= \frac{AE}{L} \tag{28}$
或者写成矩阵

$\begin{bmatrix} F_{12} \\ F_{21}\end{bmatrix} = \frac{AE}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix} \tag{29}$
缩写，K就是刚度矩阵，成了下面这样子。

$b = Ku \tag{30}$
上面的矩阵是奇异矩阵。所以一般会固定住其中一个节点u1，让矩阵变成如下样子。

$\begin{bmatrix} 0 \\ F_{21}\end{bmatrix} = \frac{AE}{L}\begin{bmatrix} 1 & 0 \\ -1 & 1\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix} \tag{31}$
说了这么多，其实这和布料模拟系统中的约束大同小异，不过就是各种参数多了一层物理意义。如下u1和u2之间是正横向的约束，而u1和u4之间则是斜的。有限元方法不过是通过矩阵将这种关系表示了出来。
好了开写代码。下面是假设节点都一开始都在一条水平线上，那么节点只会和自己左右的节点产生力。编程如下。
idx = int(i * direction + 1) # 此节点y轴编号
idxleft = int(constrain[i,0] * direction + 1) # 左侧节点y轴编号
idxright = int(constrain[i,1] * direction + 1) # 右侧节点y轴编号
if idxleft >= 0:
Kmat[idx,idx] += AE/L # 组装矩阵
Kmat[idx,idxleft] += -AE/L
# 广义胡克定律，应力等于应变乘杨氏模量
force[idx] += (L - (pos[t,idx] - pos[t,idxdown]))*E
if idxright >= 0:
Kmat[idx,idx] += AE/L
Kmat[idx,idxright] += -AE/L
force[idx] += ((pos[t,idxup] - pos[t,idx]) - L)*E对于数值计算来说，计算结果不收敛是很容易遇到的问题。其中一个比较好的解决的方法就是，就是找一些特殊的能量守恒的条件。比如先把最终稳定的结果作为初始条件写进去，如果程序没有问题的话，这个最终稳定的结果是不会随着时间改变的。借此也可用于调整参数。
比如两个球u1和u2，u1在地面固定，u2在u1上方。我希望它们之间一开始的垂直距离为x1，且弹簧放松时的长度为x2，那么u2受到弹簧的向上的力就是

$\sigma = E\varepsilon = E(x2-x1) \tag{32}$
E究竟应该是多少？随便猜一个？0.2？2？20？看别人写哪个我也写哪个？或者去查查木头的杨氏模量填上去？事实证明这些办法对写出能用的代码都没什么用。
不如设x1就是系统最终稳定后u2和u1之间的距离，且x1比x2小。根据牛顿第一定律，u2受到的弹簧力，肯定和重力g大小相同，方向相反。也就是要保证程序正确运行，就让u2静止时的弹簧长度，只要x2，x1，E，g这四个参数满足下面的关系就行了

$E = \frac{g}{x2 - x1} \tag{33}$
设好参数，算好刚度矩阵后，然后算一把逆矩阵就行了。我写了两个简单的python程序，FemSpring1写的是粒子排布在一水平条线上，仅仅受到垂直方向上重力与弹力的作用。FemSpring2则是粒子排布在一条竖直线上。它们的参数基本就是这杨确定的。
我查找资料的时候看到了很多有关能量守恒的算法，我认为它们的原理除了让系统尽量符合物理外，也是给系统尽可能多的限制条件，尽量少猜或少计算变量，让系统能够快点稳定下来。
最后，很多弹性系统都会依据如下的方程。

$M \ddot{u} + D\dot{u} + f(u)^{int} = f^{ext} \tag{33}$
其中M是质量矩阵。D是阻尼矩阵。f{int}是内力，f{ext}是外力。u是位置，头上两点代表两次求导，也就是加速度，头上一点代表一次求导，也就是速度。写成下面的样子各位肯定更加熟悉一些。

$M\frac{\partial^2 u}{\partial x^2} + D\frac{\partial u}{\partial x} = f \tag{34}$
接下来是认参数时间，之后会用的参数越来越多，与其把它们当成黑盒瞎调，不如从现在开始好好弄清每个参数的含义。
阻尼(Damping)是指任何系统在振动的时候，由于某些原因导致振动幅度逐渐下降的特性，这和流体力学中的粘性差不多。阻尼越大越难以振动。那么设置速度阻尼为很小，小，以及大的情况，并且应用到上面的代码中，如下图。
这些球相当于布料或可变形物体的节点。节点之间还有的各种弹簧没画出来。
杨氏模量(young&#39;s modulus)是指所受力的大小，与物体形变程度的比较。杨氏模量越大，物体越难以变形。如下，左侧三个点固定，且从左往右杨氏模量依次增大。
泊松比(Possion Ratio)是指物体一个方向受到压力时，另一个方向形变的程度。比如用手用力按馒头或面团这类东西，它在高度变小的同时还会变宽。如下图泊松比从左往右增大。这些粒子都仅受到重力和弹力作用。泊松比低的物体，尽量只往受重力那个方向也就是y轴变形，而泊松比高的物体，却还会往重力的垂直方向也就是x轴变形。不过这个写成代码并不复杂，之前用的一种是u1u2或u1u3之间那种正向的弹簧，一种是u1u4之间那种斜着的弹簧，泊松比大的物体，那种斜着的弹簧的刚度更大罢了。
上面这三种比较的代码在这里。好了，本章最后，用有限元方法模拟果冻或史莱姆。效果如下。我先用推进波前法生成三角形网格，然后施加重力让其撞击地面。
虽然说是地面，但我的地面其实是一堆不动的粒子组成的。类似于光滑粒子动力学里的墙壁。这些地面的粒子，也会对果冻上的粒子施加弹性影响。如果不希望物体撞击地面后形变这么大。那么就让果冻内部的刚度更大，从而不易变形。更多的请看这个视频。
这种形变效果，用布料模拟系统改一下，让它别那么柔软也可模拟出来，或者光滑粒子动力学或混合网格粒子法也可，不是必须用有限元，只是各有优缺点。而且也不是必须用二维的三角形或三维的四面体，二维的正方形和三维的体素也经常用。
只是有关弹性固体结构力学的许多知识都融到很多有限元书籍中，各种用有限元方法模拟弹性材料的资料也很多，搜&#34;Finite Element Deformable&#34;可以搜出一大堆资料与代码，可以反过来补充布料模拟的知识。
不过我的代码并不怎么好，没有优化过性能，如果把步长设置得太小，或者力设置得太大会出Bug。
根据《Gpu Pro 360 Guide to Geometry Manipulation》by Wolfgang F Engel第16章来看，想要画可以变形的物体，还有一种纯视觉的方法，那就是弄个体素化，然后调调B样条曲线之类的。
部分参考如下
根据虚幻4引擎官方发布的视频来看，他们也有个FEM有限元物理模拟系统的插件 https://www.bilibili.com/video/BV1EJ411J7WP
SIGGRAPH 2012 Course FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction.  http://www.femdefo.org/
Physically Based Deformable Models in Computer Graphics
Real-time Reduced Large-Deformation Models and Distributed Contact for Computer Graphics and Haptics
https://github.com/alecjacobson/geometry-processing-deformation
https://www.brown.edu/Departments/Engineering/Courses/En2340/Notes/notes.html
本章扯了这么一大堆，还其实仅仅是线性有限元的部分，算刚度矩阵只要算一次逆矩阵就行了。但是大自然中大部分现象都是非线性的，这时刚度矩阵就是非线性的，要用到非线性方程组的一些解法了。也就是下篇的内容。
上一篇：光影帽子：【游戏流体力学基础及Unity代码（十四）】舌尖上的有限元Galerkin法
下一篇：光影帽子：【游戏流体力学基础及Unity代码（十六）】非线性有限元及牛顿迭代法

		自动登录	找回密码
密码			立即注册

[笔记] 【游戏流体力学基础及Unity代码（十五）】线性有限元及弹性物体模拟

本帖子中包含更多资源

浏览过的版块