物理模拟相关理论的简单介绍;
刚体动力学
只处理平移与旋转,以及碰撞;
平移:力、质量->加速度->速度->位移; 旋转:力、惯性张量->扭矩、转动惯量->角加速度->加速度->角度;
基于积分的方法
- 显示积分:两阶误差,不稳定性,误差累加;
- 隐式积分:两阶误差,稳定性好,求解隐式方程(非线性优化算法);
- 半隐式积分:两阶误差,稳定性好;
- 韦尔莱积分:四阶误差,二阶精度,On2复杂度?稳定性好;不需要单独储存速度,粒子的速度和加速度都是隐式处理,稳定且高效。
- 龙格库塔积分:分二阶龙格库塔积分与四阶龙格库塔积分,常用的为四阶龙格库塔积分,On4复杂度?,五阶误差,四阶精度;
碰撞的处理
- 使用sdf判断内外;
- 如果顶点在内,可以使用在顶点处施加反相冲量的方法来改变刚体物体的运动;
- 冲量大小的计算,需要根据目标顶点改变后的速度来求得;
- 得到冲量后,将冲量作用于顶点处即可;
- 对于有多个点发生碰撞时,需要使用点的平均变化量来求取冲量;
基于粒子的方法
常用的基于粒子的敢提动力学模拟方法为shapematching;可以很好的集成到PBD框架中;
- 先基于动力学方程进行粒子的运动模拟;
- 再进行原模型还原,将粒子恢复成原始刚体模型的形状;首先由粒子平均值求质心,然后令各点到质心的距离平方和最小,也即令平方和的导数为零来求解变换矩阵,并从矩阵中提取旋转矩阵来得到旋转的变化量;
碰撞的处理
- 使用sdf判断内外;
- 由于直接使用粒子进行模拟,那么使用逐个对点进行碰撞处理即可;
布料模拟
常用弹簧质点模型,结构网格常用分布为水平、竖直、米形斜线(及其简化分布)、弯曲弹簧(及其简化分布),非结构网格分布为默认网格及相邻网格对角线(弯曲弹簧);
基于积分的方法:
积分方法与刚体动力学一致;显式积分具有不稳定问题;
比较难的是解决隐式积分问题:分为直接求解法与迭代求解法,直接求解有lu分解法等,迭代法有牛顿迭代法、jacobi迭代法等;在实时基于物理动画模拟应用较少;
Bending Spring Issue问题
直接使用弯曲弹簧会有Bending Spring Issue问题,弯曲弹簧的替代模型有多种,这里提两种:Dihedral Angle Model以及Quadratic Bending Model;前者直接使用弯曲离为弯曲角度的变化量;后者直接替换了能量函数,更容易计算H矩阵,便于隐式积分实现;
Locking Issues问题
拉伸与弯曲产生耦合,导致自由度丢失。产生的问题为所有弹簧拉伸刚度很大时,丢失了一个自由度的弯曲,像纸这种情况就无法模拟;
基于PBD的方法;
PBD并不会计算力与加速度,而是直接使用约束来限制模型的位置(相当于内力模拟),不过在进行约束之前,仍然会对外力进行动力学的更新(隐式、显式等等);
对于单个弹簧,直接约束弹簧的长度为原始长度即可;
对于多个弹簧,有多种计算方式:
- Gauss-Seidel Approach:通过逐个遍历约束来逐个限制弹簧长度;
- Jacobi Approach:遍历完约束过后,取delt_x的平均来限制顶点位置;
其它约束方式:
- Spring Strain Limit:限制弹簧的长度为原来的一定比例范围,不需要与原长一致;
- Triangle Area Limit:限制三角形的面积为原来的一定比例范围,不需要与原面积一致;
心得:
- pbd基于约束的迭代,来满足约束,迭代的次数越多,越接近理想约束效果;
- pbd算法并不是基于物理来进行计算,delttime的大小与迭代次数的多少会直接影响计算的效果,如果运行时delttime不一致会导致计算结果的不稳定性;
- Strain Limiting 方法有助于解决 locking issue;
- Jacobi迭代法与Gauss-Seidel迭代法相比,其结果收敛偏慢,但是计算结果与顺序无关,因此更适合进行GPU并行化;
有限元法面向弹性体
物理量:格林应变张量、应变能量密度(Saint Venant-Kirchhoff模型假设)、能量(应变能量密度乘以面积)、应变能量密度对能量的求导(Second Piola-Kirchhoff tress tensor);
公式推导: 力为能量对位置求导,最终得到力为其它参数的方程;
计算流程: 先计算力,然后转换为动力学积分的计算;
有限体积法面向弹性体
有限体积法,使用散度理论来进行积分,从而得到当前点所受力的公式;与有限元法不同的是,公式中所用的应力张量为当前状态(当前Normal,当前traction state)的应力张量,而不是运动前的状态; 因此有限体积法计算的过程,需要考虑应力状态从上一状态切换到当前状态的变换,这一步会相对复杂;
有限元法与有限体积法的区别:有限元法从微分的角度来求力,有限体积法从积分的方法来求力;
浅波方程的计算
公式推导:
由高度速度方程、速度压强方程推导出浅波方程;
然后使用有限差分法来带入浅波方程,即可得到高度场随时间的变化;
引入体积不变准则,引入压强与高度的函数,来调整高度变化公式,同时考虑粘性来调整公式;最终得到我们的迭代公式;
边界条件处理
- Dirichlet boundary:假设边界的高度与内部一致,用于开放边界;
- Neumann boundary:假设边界的导数与内部一致,用于闭合边界;
物体对流体的影响
首先知道物体排出水的高度,为了得到水被排出的效果,在排出水的区域会添加一个虚拟高度,来使得模拟后水体的高度为我们预期的高度,从而让效果满足要求;虚拟高度的计算可以使用被排出的高度来进行反算,得到虚拟高度后,进行正常的模拟即可;
需要注意的是虚拟高度影响并不会引起水量的增加或减少,因为其影响的那部分计算指影响压强,并不影响总水量;
流体对物体的影响
知道高度差后,使用阿基米德定律来继续每个格子的浮力,然后使用浮力来积分计算物体的运动;
欧拉方式NS方程的计算
首先使用差分法来表示物理量,速度表示在格子的面上,其它物理量表示在格子中心上;
V的更新通过分步计算NS方程来更新,其中advection部分要额外处理,会使用semi-lagrangian方法来更新速度,压强的更新也需要额外处理,压强的计算需要使用不可压缩方程来构建压强的线性方程组,通过解方程组来计算压强,然后通过压强来更新速度;
更新完速度后,可以使用semi-lagrangian来更新其它物理量;
SPH NS方程的计算
使用平滑函数的确定,用来计算由粒子插值出的空间任一位置处的物理量;
基于粒子来进行更新,外部力(重力等)先更新,压力的更新(基于粒子的压力的计算,这里采用了经验公式来计算,先计算粒子的压力,再使用平滑函数来计算目标粒子的压力梯度),粘性力的更新(由速度的laplace来计算);
计算完受力后,使用显示或隐式等方法积分计算即可;
PBF方法进行流体计算
PBF为pbd与sph方法的结合版,pbf将内部力相关的计算转化为密度约束与表面张力约束,还有部分旋度及粘性的考虑;
加速算法
使用SPH算法最主要的问题是,属性的插值计算需要采样周围的粒子;算法的复杂度是非常高的,为了解决这个问题会使用一些空间划分的算法来加速的粒子的查找与遍历;在GPU上运行的算法一般会使用均匀网格划分的算法;