流体模拟《Fluid Simulation for Computer Graphics 》——从入门到放弃(三)半拉格朗日对流/插值/MacCormack方法/压力求解

从入门到放弃系列

┬─┬ ノ( ゜-゜ノ)啊啊啊啊啊啊啊啊啊啊数值分析好难啊啊啊啊啊啊啊啊啊((╯‵□′)╯︵┻━┻在线暴躁,,真学不动了嘛就QWQ

上一篇文章:流体模拟《Fluid Simulation for Computer Graphics
》——从入门到放弃(二)模拟视角/Navier-Stokes流体方程/MAC法
,有些提到的东西可以在前置文章里找到,哪里理解错误啥的请指出3Q


开头再放上前一篇博客里的Navier-Stokes方程和变量:

  • 向量u表示流体速度
  • ρ表示流体密度
  • p表示压力,1\ρ ∇p体现了压力对流体的影响
  • 向量g表示重力加速度
  • ν表示运动粘度
  • q表示通过流体传输的模拟量之一(比如温度)
  • t表示时间

5 对流

5.1 半拉格朗日对流

首先要寻找出发点。为了在某向量x处得到q的新值,我们可以找到理论上以向量x结尾的粒子然后得到对应的q值。

我们已经知道了这个向量结尾的粒子在哪,并且这个粒子移动速度为向量u,把我们关注的这个网格点记为向量x_G,那么我们要求的就是(q_G)^(n+1)(这个网格点G在n+1时刻对应的q值)。如果一个旧粒子也是x_G结尾并且它的q值为(q_P)^n,那么时间经过Δt后,(q_G)^(n+1)=(q_P)^n

就是通过前向欧拉法求解这个出发点:。但是前向欧拉法是一阶精度的,文中使用的是更高精度的Runge-Kutta方法(龙格-库塔方法)。

前向欧拉法这个可以看这个:Forward and Backward Euler Methods,,一点点都看不进去,,总之这个东西是拿来解常微分方程的就是了

5.2 插值

大多数情况下这个x_P可能不在网格点上,这样就不能获得确切的值,但是可以通过对附近点的插值来获得近似值。

影响数值方法稳定性的关键因素是Δt,半拉格朗日对流对Δt无条件稳定,无论粒子从哪到哪,都可以从q的旧值差值得到新的q值,所以q保持有界。

如果只考虑模拟的实时速率,那么Δt可以设置成帧与帧的时间间隔。在大多数情况下需要对Δt进行限制(不能太大),Δt≤3Δx/u_max其中u_max是估计的流体的最大速度。

插值的时候利用上一步的值的加权平均值,平均操作可以平滑或模糊尖锐的特征(这样会有数值耗散),使用简单的半拉格朗日方法的尝试求解不具有粘性的对流方程时,结果看起来像是在模拟具有粘性的流体(因为都被平滑掉了好像粘的流体流不动一样)。

可以用Catmull-Rom插值进行处理,这样可以将精度提高到二阶,并显著降低数值耗散。但是这样的缺点是稳定性不如线性插值。

5.3 MacCormack方法

半拉格朗日对流对动画非常有用,因为它是无条件稳定的,所以较大的时间步长也可以顺利模拟。但是它会导致不必要的数值平滑化,使水看起来很粘稠或或者烟雾失去细节。

为了获得更高阶的精度,可以使用MacCormack方法。

给定一个量φ,利用前向对流算子A和后向对流算子A^R得到两个中间量,最后进行一个误差分析,得到φ^(N+1)的结果,这个方法可以获得更高的精度。公式描述为:

和5.1的半拉格朗日方法不同,这个MacCormack方法不是无条件稳定的,所以对结果φ^(N+1)要用限制条件进行约束,确保它落在合理范围内。在实践中,也就是要找到最接近采样点的点,并将最终值限制在这些点能获得的最小值和最大值之间。

6 压力求解

6.1 投影和线性系统求解

可以分成如下两步:

  1. 减去压力梯度

  1. 确保结果满足流体的不可压缩性以及实心墙边界条件

如果边界是自由表面,此时压力为0(狄利克雷边界条件),如果边界是实心墙,那么流体速度法向分量等于固体速度法向分量,即(14)式。

三个维度差异是:),通过MAC网格的中心差可以计算(13)式,结合(12)和(13)式可以得到这样的泊松问题:

如果用MAC网格在二维(三维也同理的)上对(15)式进行数值逼近会得到:(实话实说,,这部分我一点也没看懂)

通过(16)式可以得到下面几点:

  • 如果是自由表面边界处的,直接删去p就可以
  • 如果是实心墙边界处的,删去p然后改改p_(i,j,k)系数啥的
  • 用一个涉及流体和固体速度差的项来增加在右侧测得的负散度

用这些乱七八糟的,最后为p的求解创建一个大型线性方程组:Ap=b

A代表压力系数矩阵,p就是未知的压力值,b就是每个流体单元网格中的负散度组成的向量。

为了解这个方程,使用MICCG(0)(modified incomplete Cholesky
conjugate gradient,0级),共轭梯度(CG,conjugate gradient)是一种迭代方法。预处理共轭梯度(PCG,preconditioned conjugate gradient)本质上是相同的算法,但是可以通过一点东西加速迭代过程。

加速是这样的:对于任意可逆矩阵MAp=b的解与MAp=Mb相同,如果MA的倒数,那么MA就是单位矩阵,这样通过CG就可以快速求解p

接下来需要考虑的就是怎么求M,用线性代数学的Cholesky分解(就是平方根法那个),如果A是实对称正定矩阵就可以用,那么有A=L·L^T

最后总结一下压力求解的过程:

  1. 计算在实心墙边界处的负散度
  2. 设置压力系数矩阵A
  3. 预处理矩阵M
  4. Ap=b求解
  5. 根据压力梯度计算新的粒子速度

今天先学这么点,to be continued……

打赏
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2018-2020 LeFlacon

奶茶一杯 快乐起飞

支付宝
微信