流体模拟若干名词解释

Semi Lagrangian Advection

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

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

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

MacCormack Scheme

考虑这个方程

[latex]\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0[/latex]

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

Predictor Step:

把方程离散化,就变成了

[latex]\frac{u_i^{\overline{n+1}}-u_i^n}{\Delta t}+a\frac{u_{i+1}^n-u_i^n}{\Delta x}=0[/latex]

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

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

Corrector Step:

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

[latex]\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[/latex]

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

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

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

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

代入进去就变成了

[latex]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)[/latex]

BFECC Scheme

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

假设有场[latex]\phi[/latex].

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

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

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

然后用这个误差修正初始值:[latex]\bar{\phi}^n=\phi^n-e[/latex]

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

Projection

Advect并不保不可压性。假设Advect得到流场[latex]\textbf{v}[/latex]。

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

然后更新[latex]\textbf{v}_1=\textbf{v}-\nabla q[/latex]

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

[latex](\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[/latex]

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

[latex](\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[/latex]

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

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

发表评论

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