弹簧质点系统与时间积分

弹簧质点系统

一个模拟变形物体最简单的方法就是将其表示为弹簧质点系统（Mass Spring Systems）。一个弹簧质点包含了一系列由多个弹簧连接起来的质点，这样的系统的物理属性非常直接，模拟程序也很容易编写。

虽然模型简单，但是也带来了一些问题：

1.物体的行为依赖于弹簧系统的设置方法；

2.很难通过调整弹簧系数来得到想要的结果；

3.弹簧质点网格不能直接获取体效果。

在很多的应用中这些缺点可以忽略，在这种场合下，弹簧质点网格是最好的选择，因为够快够简单。弹簧质点系统可用于模拟绳索、布料、头发等弹性物体。

力是改变物体运动状态的原因，在这个系统中，主要有两种力，一是弹簧的弹力和阻尼力。

对于连接两个质点的一个弹簧，弹力是：

阻尼力可以这样算：

假设系统中有N个质点，质量为$m_i$，位置为$x_i$，速度为$v_i$ , $1 < i < N$.

这些质点由一组弹簧S连接，弹簧参数为（$i$, $j$, $l_0$, $k_s$, $k_d$）。$i$,$j$为连接的弹簧质点，$l_0$为弹簧完全放松时的长度，$k_s$为弹簧弹性系数，$k_d$为阻尼系数，由胡科定律知弹簧施加在两顶点上的力可以表示为：

$$ \overrightarrow{f}_i = \overrightarrow{f}^{s}(x_i,x_j)=k_s \frac{x_j - x_i}{|x_j - x _i|}(|x_j-x_i|-l_0) $$

$$ \overrightarrow{f}_j = \overrightarrow{f}^{s}(x_j,x_i)= -\overrightarrow{f}^{s}(x_i,x_j) = - \overrightarrow{f_i} $$

由受力守恒知$f_i+f_j = 0$. $f_i$和$f_j$的大小和弹簧伸长成正比关系。

对于阻尼的计算，除了与位移有关，还与质点速度有关：
$$ \overrightarrow{f_i} = \overrightarrow{f}^{d}(x_i,v_i,x_j,v_j)=k_d(v_j - v_i)\frac{x_j - x_i}{|x_j - x _i|} $$

$$ \overrightarrow{f_j} = \overrightarrow{f}^{s}(x_j,v_j,x_i,v_i) = -\overrightarrow{f_i} $$

大小和速度成正比，并且符合力守恒，则对于一个质点，其受合力方程为：

$$ \overrightarrow{f}(x_i,v_i)=\sum \overrightarrow{f}^{s}(x_i,x_j) + \sum \overrightarrow{f}^{d}(x_i,v_i,x_j,v_j) $$

这里 $j$ 为所有与质点 $i$ 存在弹簧连接的质点。后面讨论的运算都是矢量运算，为了方便就省略不写了。

在计算机模拟中，牛顿第二定律 $f = ma$ 是关键。在已知质量和外力的情况下，通过 $a=f/m $可以得到加速度，将二阶常微分方程写成两个一阶方程：
$$\frac{\mathrm{d}v}{dt} = \frac{f(x,v)}{m} $$

$$\frac{\mathrm{d}x}{dt} = v$$
可以得到解析解：

$$v(t) = v_0 + \int_{t_0}^{t}\frac{f(x,v)}{m} \mathrm{d}t$$

$$x(t) = x_0 + \int_{t_0}^{t}v(t) \mathrm{d}t$$

初始状态为 $v(t_0)$ = $v_0$, $x(t_0)=x_0$。积分将时间t内所有的变化加和，模拟的过程就是从$t_0$开始不断地计算$x(t)$和$v(t)$，然后更新质点的位置。

整个过程的伪代码如下：

// initialization
 forall particles i
         initialize xi , vi and mi
 endfor
// simulation loop
 loop
        forall particles i
               fi ← fg + fcoll + ∑ f(xi , vi , x j , v j )
        endfor
        forall particles i
             vi ← vi + ∆t fi /mi
             xi ← xi + ∆t vi
         endfor
         display the  every nth time
 endloop

时间积分

时间积分算法将这个积分的过程离散化，使用有限的步长去迭代下一时刻的状态。

欧拉法

前向欧拉方法 (显式时间积分)

$$v_{t + 1} = v_{t} + \Delta{t}\frac{f(x_t,v_t)}{m}$$

$$x_{t + 1} = x_{t} + \Delta{t}v_{t}$$

这个就是显式的欧拉解法，下一时刻的状态完全由当前状态决定。

半隐式欧拉方法（显式时间积分）

$$v_{t + 1} = v_{t} + \Delta{t}\frac{f(x_t,v_t)}{m}$$

$$x_{t + 1} = x_{t} + \Delta{t}v_{t + 1}$$

它和前向欧拉的差别很小。

后向欧拉方法（隐式时间积分）

$$x_{t + 1} = x_{t} + \Delta{t}v_{t + 1}$$

$$v_{t + 1} = v_{t} + \Delta{t}M^{-1}f(x_{t + 1},v_{t + 1})$$

这里 $M$ 为一个 $3n * 3n$ 的对角矩阵，矩阵对角线上依次为各个质点的质量, $diag(M)=(m_1,m_1,m_1,m_2,m_2,m_2,…,m_n,m_n,m_n)$。

这是一个非线性的方程，先对 $f(x_{t + 1},v_{t + 1})$ 进行一阶泰勒展开，得

$$f(x_{t + 1},v_{t + 1}) = f(x_{t} + \Delta{x},v_{t} + \Delta{v})=f(x_{t},v_{t}) + \frac{\partial f}{\partial x}(x_t,v_t) \Delta{x} + \frac{\partial f}{\partial v}(x_t,v_t)\Delta{v} $$

由$\Delta{x} = x_{t+1}-x_{t} = \Delta{t}v_{t+1}$, $\Delta{v} = v_{t+1}-v_{t}$。得
$$f(x_{t + 1},v_{t + 1}) = f(x_{t},v_{t}) + \frac{\partial f}{\partial x}(x_t,v_t) \Delta{t} v_{t + 1} + \frac{\partial f}{\partial v}(x_t,v_t)(v_{t+1}-v_{t})$$

用 $I$ 表示单位矩阵，代入上面的式子，对上式进行移项，有

$$(I - \Delta{t}M^{-1}\frac{\partial f}{\partial v}(x_t,v_t) - (\Delta{t})^{2}M^{-1}\frac{\partial f}{\partial v}(x_{t},v_{t}))v_{t+1} = v_{t}+\Delta{t}M^{-1}(f(x_{t},v_{t})-\frac{\partial f}{\partial v}(x_{t},v_{t})v_{t})$$

$\frac{\partial f}{\partial v}$ 与 $\frac{\partial f}{\partial x}$ 的求法见这里。

上面的式子是一个形如，$AX = B$的方程组，$A$维度比较大，而且不一定可逆，一般使用雅克比迭代等方法来解。