CFD入门傻瓜版[2]

上次（http://gezhi.org/node/399）说到了求解可压缩流动的一个重要算法，通用Godunov方法。其两个主要步骤就是
1 怎么确定 $\mathbf{Q}^{-}_{i+\frac{1}{2}}$ 和 $\mathbf{Q}^{+}_{i+\frac{1}{2}}$
2 怎么计算 $\mathbf{F}_{i+\frac{1}{2}}$

这里我们给出第一点一个具体的实现方法，就是基于原始变量的MUSCL格式（以下简称MUSCL）。它是一种很简单的格式，而且具有足够的精度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主页（http://cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html）上看手册，里面空间离散那一章清楚的写着。

MUSCL假设控制体内原始变量（就是 $\mathbf{Q}$ ）的分布是一次或者二次多项式，如果得到了这个多项式，就可以求出控制体 $i$ 左右两个界面的一侧的值 $\mathbf{Q}^{-}_{i+\frac{1}{2}}$ 和 $\mathbf{Q}^{+}_{i-\frac{1}{2}}$ 。

我们以压力p为例来说明怎么构造这个多项式。这里我只针对二次多项式来讲解，你看完之后肯定能自己推导出一次多项式的结果（如果你搞不定，那我对你的智商表示怀疑）。

OK，开始
假设 $\triangle x =1$ ，这个假设不影响最终结论，因为你总可以把一个区间线性的变换到长度为1的区间。
假设压力p在控制体i内部的分布是一个二次多项式 $f(p)=ax^2+bx+c$ ，控制体i的中心处于 $x=0$ 处，左右两个界面就是 $-\frac{1}{2}$ 和 $+\frac{1}{2}$ 。

这里先强调一个问题，在FVM中，每个控制体内的求解出来的变量 $\mathbf{Q}$ 实际上是这个控制体内的平均值 $\mathbf{\bar{Q}}_i=\frac{1}{\triangle x}\int^{i+\frac{1}{2}}_{i-\frac{1}{2}}\mathbf{Q}dx$ 。
所以，
$\bar{p}_i=\int^{i+\frac{1}{2}}_{i-\frac{1}{2}}(ax^2+bx+c)dx$
$=\left[ax^3/3+bx^2/2+cx \right]^{+1/2}_{-1/2}$
$=\frac{a}{12}+c$ 。

我们知道 ${\bar{p}}_{i-1}$ ， ${\bar{p}}_i$ 和 ${\bar{p}}_{i+1}$ ，等距网格情况下 $i-\frac{1}{2}$ 和 $i+\frac{1}{2}$ 处的导数可以近似表示为
$\frac{\partial p}{\partial x}\mid_{i-\frac{1}{2}} \simeq \triangle^-_i = {\bar{p}_{i} - \bar{p}_{i-1}}$
$\frac{\partial p}{\partial x}\mid_{i+\frac{1}{2}} \simeq \triangle^+_i ={\bar{p}_{i+1} - \bar{p}_{i}}$
那么
$f'(-\frac{1}{2}) = (2ax^2+b)|_{x=-\frac{1}{2}}=b-a=\triangle_-$
$f'(+\frac{1}{2}) = (2ax^2+b)|_{x=+\frac{1}{2}}=b+a=\triangle_+$

由上述三个有关a，b和c的方程，我们可以得到
$a=\frac{\triangle^+_i - \triangle^-_i}{2}$
$b=\frac{\triangle^+_i + \triangle^-_i}{2}$
$c=\bar{p}_i - \frac{\triangle^+_i - \triangle^-_i}{24}$
这样就可以得到f(x)的表达式了，由此可以算出 ${p}^{-}_{i+\frac{1}{2}}$ 和 ${p}^{+}_{i-\frac{1}{2}}$

通常MUSCL格式写成如下形式
${p}^{-}_{i+\frac{1}{2}}=\bar{p}_i + \frac{1}{4}[(1-k)\triangle^-_i+(1+k)\triangle^+_i]$
${p}^{+}_{i-\frac{1}{2}}=\bar{p}_i - \frac{1}{4}[(1-k)\triangle^+_i+(1+k)\triangle^-_i]$
$k=1/3$ 对应我们的推导结果（二次多项式假设）。

但是这不是最终形式。如果直接用这个公式，就会导致流场在激波（间断）附近的振荡。因为直接用二次多项式去逼近一个间断，会导致这样的效果。所以科学家们提出要对间断附近的斜率有所限制，因此引入了一个非常重要的修改——斜率限制器。加入斜率限制器后，上述公式就有点变化。
${p}^{-}_{i+\frac{1}{2}}=\bar{p}_i + \frac{s_i}{4}[(1-ks_i)\triangle^-_i+(1+ks_i)\triangle^+_i]$
${p}^{+}_{i-\frac{1}{2}}=\bar{p}_i - \frac{s_i}{4}[(1-ks_i)\triangle^+_i+(1+ks_i)\triangle^-_i]$
这里 $s_i$ 是Van Albada限制器
$s_i=\frac{2\triangle^+_i\triangle^-_i+\epsilon}{(\triangle^-_i )^2+ (\triangle^-_i)^2+\epsilon}$
$\epsilon$ 是一个小数（ $\epsilon=1\times 10^{-6}$ ），以防止分母为0。

密度和速度通过同样的方法来搞定。

密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的MUSCL。此外还有基于特征变量的MUSCL，要复杂一点，但是被认为适合更高精度的格式。然而一般计算中，基于原始变量的MUSCL由于具有足够的精度、简单的形式和较低的代价而被广泛使用。

OK，搞定了。下面进入第二点，怎么求 $\mathbf{F}_{i+\frac{1}{2}}$ 。关于这一点，我不打算做详细介绍了，直接使用现有的近似黎曼解就可以了，都是通过 $\mathbf{Q}^{-}_{i+\frac{1}{2}}$ 和 $\mathbf{Q}^{+}_{i+\frac{1}{2}}$ 计算得到 $\mathbf{F}_{i+\frac{1}{2}}$ 。比如Roe因为形式简单，而非常流行。在CFL3D软件主页（http://cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html）上看手册，附录C的C.1.3。

格致

用户登录

导航