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

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

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

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

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

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

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

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

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

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

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

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

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