注册 登录
Simwe仿真论坛(forum.simwe.com),CAE/CAD/CAM/,FEA/FEM/有限元分析论坛---(手机验证注册) 返回首页

Fluid的个人空间 https://home.simwe.com/?467291 [收藏] [复制] [分享] [RSS]

日志

[推荐]CFD入门傻瓜版[2]

已有 747 次阅读2010-7-10 19:42 |个人分类:CFD理论

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


路过

雷人

握手

鲜花

鸡蛋

评论 (0 个评论)

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 注册

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 16:39 , Processed in 0.021967 second(s), 13 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

返回顶部