完整版CFD理论过渡到编程的傻瓜入门教程Word格式.docx
《完整版CFD理论过渡到编程的傻瓜入门教程Word格式.docx》由会员分享，可在线阅读，更多相关《完整版CFD理论过渡到编程的傻瓜入门教程Word格式.docx（24页珍藏版）》请在冰点文库上搜索。
针对一个微元控制体
，把Euler方程在空间积分
用微积分知识可以得到
也就是说控制体内气体状态平均值的变化是控制体界面上流通量的结果。
因此我们要计算
的演化，关键问题是计算控制体界面上的
。
FVM就是以这个积分关系式出发，把整个流场划分为许多小控制体，每个控制体和周围相邻的某个控制体共享一个界面，通过计算每个界面上的通量来得到相邻控制体之间的影响，一旦每个控制体的变化得到，整个流场的变化也就知道了。
所以，再强调一次，关键问题是计算控制体界面上的
初学者会说，这个不难，把界面
上的
插值得到，然后就可以计算
有道理！
咱们画个图，有三个小控制体i-1到i+1，中间的“|”表示界面，控制体i右边的界面用
表示，左边的就是
|i-1|i|i+1|
好下个问题：
每个小控制体长度都是
如何插值计算界面
？
最自然的想法就是：
取两边的平均值呗，
但是很不幸，这是不行的。
那么换个方法？
直接平均得到
还是很不行，这样也不行。
我靠，这是为什么？
这明明是符合微积分里面的知识啊？
这个道理有点复杂，说开了去可以讲一本书，可以说从50年代到70年代，CFD科学家就在琢磨这个问题。
这里，初学者只需要记住这个结论：
对于流动问题，不可以这样简单取平均值来插值或者差分。
如果你非要想知道这个究竟，我现在也不想给你讲清楚，因为我眼下的目的是让你快速上手，而且该不刨根问底的时候就不要刨根问底，这也是初学阶段一种重要的学习方法。
好了，既然目的只是为了求
，我在这里，只告诉你一种计算方法，也是非常重要、非常流行的一种方法。
简单的说，就是假设流动状态在界面
是不连续的，先计算出界面
两边
的值，
和
，再由它们用某种方法计算出
上述方法是非常重要的，是由一个苏联人Godunov在50年代首创的，后来被发展成为通用Godunov方法，著名的ENO/WENO就是其中的一种。
好了，现在的问题是：
1怎么确定
2怎么计算
对于第一个问题，Godunov在他的论文中，是假设每个控制体中
是均匀分布的，因此
第二个问题，Godunov采用了精确黎曼解来计算
什么是“精确黎曼解”，就是计算这个激波管问题的精确解。
既然有精确解，那还费功夫搞这些FVM算法干什么？
因为只有这种简单一维问题有精确解，稍微复杂一点就不行了。
精确解也比较麻烦，要分析5种情况，用牛顿法迭代求解（牛顿法是什么？
看数值计算的书去，哦，算了，现在暂时可以不必看）。
这是最初Godunov的方法，后来在这个思想的基础上，各种变体都出来了。
也不过是在这两个问题上做文章，怎么确定
，怎么计算
Godunov假设的是每个小控制体内是均匀分布，也就是所谓分段常数（piecewiseconstant），所以后来有分段线性（picewiselinear）或者分段二次分布（picewiseparabolic），到后来ENO/WENO出来，那这个假设的多项式次数就继续往上走了。
都是用多项式近似的，这是数值计算中的一个强大工具，你可以在很多地方看到这种近似。
Godunov用的是精确黎曼解，太复杂太慢，也不必要，所以后来就有各种近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是对精确黎曼解做的简化。
这个多项式的阶数是和计算精度密切相关的，阶数越高，误差就越小。
不过一般来说，分段线性就能得到不错的结果了，所以工程中都是用这个，Fluent、Fastran以及NASA的CFL3D、OverFlow都是用这个。
而黎曼求解器对精度的影响不是那么大，但是对整个算法的物理适用性有影响，也就是说某种近似黎曼求解器可能对某些流动问题不合适，比如单纯的Roe对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定。
上次（http:
//gezhi.org/node/399）说到了求解可压缩流动的一个重要算法，通用Godunov方法。
其两个主要步骤就是
这里我们给出第一点一个具体的实现方法，就是基于原始变量的MUSCL格式（以下简称MUSCL）。
它是一种很简单的格式，而且具有足够的精度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主页（http:
//cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html）上看手册，里面空间离散那一章清楚的写着。
MUSCL假设控制体内原始变量（就是
）的分布是一次或者二次多项式，如果得到了这个多项式，就可以求出控制体
左右两个界面的一侧的值
我们以压力p为例来说明怎么构造这个多项式。
这里我只针对二次多项式来讲解，你看完之后肯定能自己推导出一次多项式的结果（如果你搞不定，那我对你的智商表示怀疑）。
OK，开始
假设
，这个假设不影响最终结论，因为你总可以把一个区间线性的变换到长度为1的区间。
假设压力p在控制体i内部的分布是一个二次多项式
，控制体i的中心处于
处，左右两个界面就是
这里先强调一个问题，在FVM中，每个控制体内的求解出来的变量
实际上是这个控制体内的平均值
所以，
我们知道
，
，等距网格情况下
处的导数可以近似表示为
那么
（这里错了，应该是2ax+b）
由上述三个有关a，b和c的方程，我们可以得到
这样就可以得到f（x）的表达式了，由此可以算出
通常MUSCL格式写成如下形式
对应我们的推导结果（二次多项式假设）。
但是这不是最终形式。
如果直接用这个公式，就会导致流场在激波（间断）附近的振荡。
因为直接用二次多项式去逼近一个间断，会导致这样的效果。
所以科学家们提出要对间断附近的斜率有所限制，因此引入了一个非常重要的修改——斜率限制器。
加入斜率限制器后，上述公式就有点变化。
是VanAlbada限制器
是一个小数（
），以防止分母为0。
密度和速度通过同样的方法来搞定。
密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的MUSCL。
此外还有基于特征变量的MUSCL，要复杂一点，但是被认为适合更高精度的格式。
然而一般计算中，基于原始变量的MUSCL由于具有足够的精度、简单的形式和较低的代价而被广泛使用。
OK，搞定了。
下面进入第二点，怎么求
关于这一点，我不打算做详细介绍了，直接使用现有的近似黎曼解就可以了，都是通过
计算得到
比如Roe因为形式简单，而非常流行。
在CFL3D软件主页（http:
//cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html）上看手册，附录C的C.1.3。
想了一下，还是把Roe求解器稍微说说吧，力求比较完整。
但是不要指望我把Roe求解器解释清楚，因为这个不是很容易三言两语说清的。
Roe求解器的数学形式是这样的
显然这个公式的第一项是一个中心差分形式，先前说过简单的中心差分不可行，原因是耗散不足导致振荡，说得通俗点就像一个弹簧，如果缺乏耗散（阻尼）它就会一直振荡。
“耗散”这个术语在激波捕捉格式中是最常见的。
第二项的作用就是提供足够的耗散了。
已经用MUSCL求得了，
的定义在第一讲中已经介绍了。
只有
是还没说过的。
这个矩阵可以写成特征矩阵和特征向量矩阵的形式
而
的具体表达式在许多书上都有，而且这里的矩阵表达有问题，所以就不写了。
是由
、
代入
计算得到。
采用所谓Roe平均值
这才是Roe求解器关键的地方！
总结一下，就是用Roe平均计算界面上的气体状态
，然后计算得到
，这样
就可以得到了。
如果有时间，我后面会找一个代码逐句分析一下。
总之，计算
还是很不直接的。
构造近似黎曼解是挺有学问的，需要对气体动力学的物理和数学方面有较深的理解。
通常，如果不是做基础研究，你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导出来的。
附录程序：
（新建一个.c类型文件，将下面的程序复制粘贴到里面，就可以运行了）
/****************************************************
本程序用于求解一维无粘可压缩欧拉方程（激波管问题）
运用DummyCell处理边界条件;
通量计算方式:
AUSMScheme;
重构方法：
MUSCL方法
限制器：
VanAlbada限制器
时间离散：
四步Runge-Kutta方法
****************************************************/
#include"
stdio.h"
conio.h"
malloc.h"
stdlib.h"
math.h"
string.h"
#defineh（1/400.0）//网格步长
#defineNc404//网格总数：
与h之间的关系Nc=1/h+4
#definePI3.1415927
#defineIt1000
#definegama1.4//气体比热比
doubleKAKA=0.0;
//限制器控制参数
doubleXS1,XS2;
//计算域的两个端点
doubledt=2.5e-5;
//时间步长
doubletimesum;
//总的计算时间
//函数声明
voidoutput（）;
voidSolveWtoU（doubleW[3],doubleU[3]）;
voidSolveUtoW（doubleW[3],doubleU[3]）;
//基本变量和守恒变量之间的转换函数
//前后各留两个网格单元作为虚拟网格单元计算网格单元从2到NOc-3
structcell
{
intflag;
//网格点的类型
doublexc;
doubleW[3],Wp[3];
//conservationvaraible
doubleU[3];
//jibenbianliang
doubleR[3];
doubleS;
//entropy
};
structcellcell[Nc];
//网格生成及流场初始化
voidinitialsolve（）
{
inti;
doublex,xi,xe;
XS1=-0.0;
XS2=1.0;
xi=XS1-2*h;
xe=XS2+2*h;
for（i=0;
i<
Nc;
i++）
{
cell[i].xc=（xi+h/2）+i*h;
//网格中心坐标
cell[i].flag=0;
x=cell[i].xc;
if（x<
0.5）
{
cell[i].U[0]=0.445;
cell[i].U[1]=0.698;
cell[i].U[2]=3.528;
}
else
cell[i].U[0]=0.5;
cell[i].U[1]=0.0;
cell[i].U[2]=0.571;
SolveUtoW（cell[i].W,cell[i].U）;
//printf（"
x=%f,d=%f,u=%f,p=%f\n"
cell[i].xc,cell[i].U[0],cell[i].U[1],cell[i].U[2]）;
getchar（）;
}
}
//边界条件：
DummyCell方法
voidboundarycondition（）
//直接赋远场值
if（cell[i].flag==1）
SolveUtoW（cell[i].W,U）;
elseif（cell[i].flag==2）
//重构：
MUSCL方法+VanAlbada限制器
voidSolveReconstruction（doubleUm1[],doubleUI[],doubleUp1[],doubleUp2[],doubleUL[],doubleUR[]）
doubleepsilon,aR[3],aL[3],bL[3],bR[3],sL[3],sR[3];
epsilon=1.0e-5*h;
3;
{
aR[i]=Up2[i]-Up1[i];
//Dealta+U（I+1）
bR[i]=Up1[i]-UI[i];
//Dealta_U（I+1）
aL[i]=Up1[i]-UI[i];
//Dealta+U（I）
bL[i]=UI[i]-Um1[i];
//Dealta_U（I）
sR[i]=（2.0*aR[i]*bR[i]+1.0e-6）/（aR[i]*aR[i]+bR[i]*bR[i]+1.0e-6）;
sL[i]=（2.0*aL[i]*bL[i]+1.0e-6）/（aL[i]*aL[i]+bL[i]*bL[i]+1.0e-6）;
UR[i]=Up1[i]-0.25*sR[i]*（（1+KAKA*sR[i]）*bR[i]+（1-KAKA*sR[i]）*aR[i]）;
UL[i]=UI[i]+0.25*sL[i]*（（1+KAKA*sL[i]）*aL[i]+（1-KAKA*sL[i]）*bL[i]）;
//守恒变量转化为基本变量
voidSolveWtoU（doubleW[3],doubleU[3]）
{
U[0]=W[0];
//d
U[1]=W[1]/W[0];
//u
U[2]=（gama-1.0）*（W[2]-0.5*U[0]*U[1]*U[1]）;
//基本变量转化为守恒变量
voidSolveUtoW（doubleW[3],doubleU[3]）
W[0]=U[0];
W[1]=U[1]*U[0];
//u
W[2]=U[2]/（gama-1.0）+0.5*U[0]*U[1]*U[1];
//格式：
AUSM格式
voidSolveAUSMFlux（doubleUL[],doubleUR[],doubleFc[]）//UL左状态，UR右状态，FC通量
doubleWL[3]={0.0},WR[3]={0.0};
doubledeta,fMn;
doubleML,cl,MR,cr,Ml,Mr,Mn;
doublepl,pr,pn,PL,PR;
deta=0.25;
SolveUtoW（WL,UL）;
SolveUtoW（WR,UR）;
PL=UL[2];
//计算左单元的压力
cl=sqrt（gama*PL/UL[0]）;
PR=UR[2];
//计算右单元的压力
cr=sqrt（gama*PR/UR[0]）;
ML=UL[1]/cl;
MR=UR[1]/cr;
if（ML>
=1.0）
{Ml=ML;
pl=PL;
}
elseif（fabs（ML）<
1.0）
{Ml=0.25*（ML+1）*（ML+1）;
pl=0.25*PL*（ML+1）*（ML+1）*（2.0-ML）;
else
{Ml=0.0;
pl=0.0;
if（MR>
{Mr=0.0;
pr=0.0;
elseif（fabs（MR）<
{Mr=-0.25*（MR-1）*（MR-1）;
pr=0.25*PR*（MR-1）*（MR-1）*（2.0+MR）;
{Mr=MR;
pr=PR;
Mn=Ml+Mr;
pn=pl+pr;
if（fabs（Mn）>
deta）fMn=fabs（Mn）;
elsefMn=0.5*（fabs（Mn）*fabs（Mn）+deta*deta）/deta;
Fc[0]=0.5*Mn*（WL[0]*cl+WR[0]*cr）-0.5*fMn*（WR[0]*cr-WL[0]*cl）;
Fc[1]=0.5*Mn*（WL[1]*cl+WR[1]*cr）-0.5*fMn*（WR[1]*cr-WL[1]*cl）+pn;
Fc[2]=0.5*Mn*（（WL[2]+PL）*cl+（WR[2]+PR）*cr）-0.5*fMn*（（WR[2]+PR）*cr-（WL[2]+PL）*cl）;
//残差计算R
voidSolveResidual（）
inti,j;
doubleUL[3],UR[3],Fp[3],Fm[3];
if（cell[i].flag==0）
//i+1/2
SolveReconstruction（cell[i-1].U,cell[i].U,cell[i+1].U,cell[i+2].U,UL,UR）;
//MSUCL重构
SolveAUSMFlux（UL,UR,Fp）;
//二阶格式
//SolveAUSMFlux（cell[i].U,cell[i+1].U,Fp）;
//一阶格式
//i-1/2
SolveReconstruction（cell[i-2].U,cell[i-1].U,cell[i].U,cell[i+1].U,UL,UR）;
//MSUCL重构
SolveAUSMFlux（UL,UR,Fm）;
//二阶格式
//SolveAUSMFlux（cell[i-1].U,cell[i].U,Fm）;
for（j=0;
j<
j++）
{
cell[i].R[j]=-（Fp[j]-Fm[j]）/h;
}
//流场值更新
voidSolveNextstep（doublear[],intir）
=Nc;
{
if（ir==0）
{
for（j=0;
j++）
cell[i].Wp[j]=cell[i].W[j];
}
j++）
cell[i].W[j]=cell[i].Wp[j]+ar[ir]*dt*cell[i].R[j];
SolveWtoU（cell[i].W,cell[i].U）;
//Runge-Kutta方法
voidSolveRungeKutta（）
doublear[4]={1/4.0,1/3.0,0.5,1.0};
intit,ir;
timesum=0.0;
for（it=0;
;
it++）//迭代步数
for（ir=0;
ir<
4;
ir++）//四步Rouger-Kutta法
boundarycondition（）;
SolveResidual（）;
SolveNextstep（ar,ir）;
timesum+=dt;
printf（"
it=%d,timesum
升级会员 