FEM物理引擎解析

2024-02-05 21:45:22

FEM 把物体分为多个基本图元，如四面体。
我们可以用四面体的顶点集合来表示这个物体。

Loop:
通过受力更新顶点速度
碰撞检测<先导顶点位置更新>
顶点位置更新

为了达到更好的效果使用后向欧拉法，
x表示位置，t表示时间，可以导出两个公式：

下时刻位置=这时刻位置+这时刻到下时刻的速度*时间差；
这时刻到下时刻的速度=这时刻速度+物体下时刻位置对应受力/质量*时间差
（比较准确的速度表示方法是向后差分）

对应数学表达

\[x_{t+1} = x_t + v_{t+1}\Delta t \\ v_{t+1} = v_t + \frac{f(x_{t+1})}{M}\Delta t \]

消去x_{t+1}

\[v_{t+1} = v_t + \frac{f(x_t + v_{t+1}\Delta t)}{M}\Delta t \]

我们在x_t点泰勒展开一阶来近似

\[v_{t+1} = v_t + (f(x_t) + \frac{\partial f}{\partial x}(x_t) + v_{t+1}\Delta t) \frac{\Delta t}M \]

化简得

\[(1 - \frac{\frac{\partial f}{\partial x}(x_t) \Delta t^2}M)v_{t+1} = v_t + \frac{f(x_t)}M\Delta t \]

这里只要求一个自变量v_{t+1}，使用迭代法求解

一个四面体有四个顶点，可以看做一个 4 * 3 的矩阵 X
现在我们可以通过一个 3 * 3 的变形矩阵 R 以及平移向量 b 映射，得到变形后的四面体

\[X_{old} R + b = X_{new} \]

变形的内力指向形变减少的方向，可以通过偏导计算变形位移，并用胡克定律计算
但是这里并不直接计算内力，而是用一个势能场来计算势能 E，以更好地泛化问题。

\[E = \psi(\frac{\partial X_{new}}{\partial X_{old}}) = \psi(GR^T) \]

上面的 G 表示一个形状与 X 一样的全为 1 的矩阵
这里可以获得整个物体的内力势能，即所有四面体的势能之和:

\[U = \sum_i E_i \]

由于内力会作用在整个物体势能减小的方向，
因此，物体全局势能 U 对每一个顶点的负梯度就等效于该点的受力

\[f_i = - \frac{\partial U}{\partial x_i} \]

有限元模型就是在探讨如何选择势能场的问题，常用的比如:

\[F = GR^T \\ J = \frac{V_{new}}{V_{old}} = |R| \\ \psi(F) = \frac{\mu}{2}\sum_i[(F^TF)_{ii}-1] - \mu log(J) + \frac{\lambda}2 log^2(J) \]

\[\psi(F) = \mu \sum_i[(\sigma_i-1)^2] + \frac{\lambda}2(J-1)^2 \]

这里需要一个自动微分求解库来完成操作，可以使用Torch，
加法后面惩罚图元体积的变化，以接近现实情况
再根据重力、接触力等即可模拟物体运动

码农公寓