主页 › CFD入门傻瓜版[2]

CFD入门傻瓜版[2]

周一, 2006-12-18 13:22 — 优雅的混蛋

上次（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。

Your rating: None
‹ CFD入门傻瓜版[1]向上

回复

此内容将保密，不会被其他人看见。