流体模拟若干名词解释

Semi Lagrangian Advection

“全拉格朗日”是:就去算一个个粒子,网格啥的一概没有。

“半拉格朗日”就“半”在,它每次用粒子方式Advect之后,就又回到了欧拉网格上。

具体地,第n+1步的A(i,j)等于第n步的A(x,y),其中(x,y)是(i,j)这个点在第n步时应该在的位置。

MacCormack Scheme

考虑这个方程

\(\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0\)

MacCormack方法分成两步:Predictor Step和Corrector Step.

Predictor Step:

把方程离散化,就变成了

\(\frac{u_i^{\overline{n+1}}-u_i^n}{\Delta t}+a\frac{u_{i+1}^n-u_i^n}{\Delta x}=0\)

这个\(u_i^{\overline{n+1}}\)表示的是“预估的第n+1步的值”的意思。把它写成explicit形式:

\(u_i^{\overline{n+1}}=u_i^n-a\frac{\Delta t}{\Delta x}\left(u_{i+1}^n-u_i^n\right)\)

Corrector Step:

假设我们知道时间n+0.5的信息,那么就应该有

\(\frac{u_i^{n+1}-u_i^{n+\frac{1}{2}}}{1/2\Delta t}+a\frac{u_i^{\overline{n+1}}-u_{i-1}^{\overline{n+1}}}{\Delta x}=0\)

注意,Predictor Step里面,关于x偏导那一项是i+1减i(前向差分),但Corrector Step里面就是后向差分了。这两个也可以交换,先后向差分再前向差分。但是必须不一样。这是达到二阶精度所必须做的。

Corrector Step关于x偏导是用“预估值”\(u^{\overline{n+1}}\)做的,需要注意。

我们不知道\(u_i^{n+\frac{1}{2}}\)是什么,就用“预估值”和\(u^n\)进行插值:

\(u_i^{n+\frac{1}{2}}=\frac{1}{2}\left( u_i^n+u_i^{\overline{n+1}} \right)\)

代入进去就变成了

\(u_i^{n+1}=\frac{u_i^n+u_i^{\overline{n+1}}}{2}-a\frac{\Delta t}{2\Delta x}\left( u_i^{\overline{n+1}}-u_{i-1}^{\overline{n+1}}\right)\)

BFECC Scheme

BFECC是缩写,全拼是Back and Forth Error Compensation and Correction.

假设有场\(\phi\).

先用“正向”Advect的算子\(A\)算出一个预估值\(\hat{\phi}^{n+1}=A(\phi^n)\)

然后再用“反向”的算子\(A^R\)算出第n步的“估计”值\(\hat{\phi}^n=A^R(\hat{\phi}^{n+1})\)

然后我们估计得到一个误差\(e=(\hat{\phi}^n-\phi^n)/2.\)

然后用这个误差修正初始值:\(\bar{\phi}^n=\phi^n-e\)

并重新Advect:\(\phi^{n+1}=A(\bar{\phi}^n)\)

Projection

Advect并不保不可压性。假设Advect得到流场\(\textbf{v}\)。

根据亥姆霍兹分解,我们可以设想,求一个“压力”场\(q\),使得\(\nabla^2q=\nabla\cdot\textbf{v}\)

然后更新\(\textbf{v}_1=\textbf{v}-\nabla q\)

我们设\(\textbf{v}=(u,v)\).右边对\(\textbf{v}\)求散度,由于拉普拉斯算子对称,所以这里也可以用中心差分:

\((\nabla\cdot\textbf{v})_{i,j}=(u_{i+1,j}-u_{i,j})/2\Delta x+(v_{i,j+1}-v_{i,j-1})/2\Delta y\)

前面的拉普拉斯算子写出来就是

\((\nabla^2q)_{i,j}=(q_{i+1,j}-2q_{i,j}+q_{i-1,j})/\Delta x^2+(q_{i,j+1}-2q_{i,j}+q_{i,j-1})/\Delta y^2\)

我们要解的是\(q\).如果把它写成一个线性方程组的形式,那么对角线就是\(-2(1/\Delta x^2+1/\Delta y^2)\),(i,j)那一行里面,(i+1,j)和(i-1,j)对应项的元素是\(1/\Delta x^2\),(i,j-1)和(i,j+1)对应项都是\(1/\Delta y^2\).方程右边的常数向量则是展开的\((\nabla\cdot\textbf{v})\).

这里还要注意一个Boundary Condition的情况。假设我们的网格范围是(1,1)~(n-1,m-1),外面一圈是“假”的边界值。如果采取“复制型”的边界条件,就是例如q(0,j)=q(1,j)这样。这时候(1,j)这一行就少了一项。

发表评论

电子邮件地址不会被公开。 必填项已用*标注