CFD理论过渡到编程的傻瓜入门教程.doc

. . . .
CFD理论过渡到编程的傻瓜入门教程
〔注：这是一篇不知道谁写的介绍一维无粘可压缩Euler方程，以与如何编程实现求解该方程的论文。作者从最根本的概念出发，深入浅出的讲解了控制方程，有限体积格式，MSUCL方法，限制器，Roe格式等相关知识。这篇论文我觉得有利于大家学习CFD编程的相关知识，所以推荐给大家。文章的后面附有我写的程序〔C语言〕，用于求解一维激波管问题，大家有兴趣可以看看〔程序中加了注释说明〕胡偶2011〕
借宝地写几个小短文，介绍CFD的一些实际的入门知识。主要是因为这里支持Latex，写起来比拟方便。
CFD，计算流体力学，是一个挺难的学科，涉与流体力学、数值分析和计算机算法，还有计算机图形学的一些知识。尤其是有关偏微分方程数值分析的东西，不是那么容易入门。大多数图书，片中数学原理而不重实际动手，因为作者都把读者当做已经掌握根底知识的科班学生了。所以数学根底不那么好的读者往往看得很吃力，看了还不知道怎么实现。本人当年虽说是学航天工程的，但是那时本科教育已经退步，根底的流体力学课被砍得只剩下一维气体动力学了，因此自学CFD的时候也是头晕眼花。不知道怎么实现，也很难找到教学代码——那时候网络还不兴旺，只在教研室的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着头皮分析。后来网上淘到一些代码研读，结合书籍论文才慢慢入门。可以说中间没有教师教，后来赌博士为了混学分上过CFD专门课程，不过那时候我已经都掌握课堂上那些了。
回想自己入门艰辛，不免有一个想法——写点通俗易懂的CFD入门短文给师弟师妹们。本人不打算搞得很系统，而是希望能结合实际，说明一些最根本的概念和手段，其中一些复杂的道理只是点到为止。目前也没有具体的计划，想到哪里写到哪里，因此可能会很零散。但是我争取让初学CFD的人能够了解一些根本的东西，看过之后，会知道一个CFD代码怎么炼成的〔这“炼〞字好似很流行啊〕。欢迎大家提出意见，这样我尽可能的可以追加一些修改和解释。
言归正传，第一局部，我打算介绍一个最根本的算例，一维激波管问题。说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体状态〔密度、速度、压力〕不一样，突然把隔板抽去，管子面的气体怎么运动。这是个一维问题，被称作黎曼连续问题，好似是黎曼最初研究双曲微分方程的时候提出的一个问题，用一维无粘可压缩Euler方程就可以描述了。

这里

这个方程就是描述的气体密度、动量和能量随时间的变化〔〕与它们各自的流量〔密度流量，动量流量，能量流量〕随空间变化〔〕的关系。
在CFD常把这个方程写成矢量形式

这里

进一步可以写成散度形式
一定要熟悉这种矢量形式
以上是控制方程，下面说说求解思路。可压缩流动计算中，有限体积〔FVM〕是最广泛使用的算法，其他算法多多少少都和FVM有些联系或者共通的思路。了解的FVM，学习其他高级点的算法〔比如目前比拟热门的连续有限元、谱FVM、谱FDM〕就好说点了。
针对一个微元控制体，把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假设的是每个小控制体是均匀分布，也就是所谓分段常数(piecewise constant)，所以后来有分段线性(picewise linear)或者分段二次分布(picewise parabolic)，到后来ENO/WENO出来，那这个假设的多项式次数就继续往上走了。都是用多项式近似的，这是数值计算中的一个强大工具，你可以在很多地方看到这种近似。
Godunov用的是准确黎曼解，太复杂太慢，也不必要，所以后来就有各种近似黎曼解，最有名的是Roe求解器、HLL求解器和Osher求解器，都是对准确黎曼解做的简化。
这个多项式的阶数是和计算精度密切相关的，阶数越高，误差就越小。不过一般来说，分段线性就能得到不错的结果了，所以工程中都是用这个，Fluent、Fastran以与NASA的CFL3D、OverFlow都是用这个。而黎曼求解器对精度的影响不是那么大，但是对整个算法的物理适用性有影响，也就是说某种近似黎曼求解器可能对某些流动问题不适宜，比如单纯的Roe对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定。
上次〔gezhi.org/node/399〕说到了求解可压缩流动的一个重要算法，通用Godunov方法。其两个主要步骤就是
1 怎么确定和
2 怎么计算
这里我们给出第一点一个具体的实现方法，就是基于原始变量的MUSCL格式〔以下简称MUSCL〕。它是一种很简单的格式，而且具有足够的精度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主页〔cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html〕上看手册，里面空间离散那一章清楚的写着。
MUSCL假设控制体原始变量〔就是〕的分布是一次或者二次多项式，如果得到了这个多项式，就可以求出控制体左右两个界面的一侧的值和。
我们以压力p为例来说明怎么构造这个多项式。这里我只针对二次多项式来讲解，你看完之后肯定能自己推导出一次多项式的结果〔如果你搞不定，那我对你的智商表示怀疑〕。
OK，开始
假设，这个假设不影响最终结论，因为你总可以把一个区间线性的变换到长度为1的区间。
假设压力p在控制体i部的分布是一个二次多项式，控制体i的中心处于处，左右两个界面就是和。
这里先强调一个问题，在FVM中，每个控制体的求解出来的变量实际上是这个控制体的平均值。
所以，
。
我们知道，和，等距网格情况下和处的导数可以近似表示为

那么
〔这里错了，应该是2ax+b〕
由上述三个有关a，b和c的方程，我们可以得到

这样就可以得到f(x)的表达式了，由此可以算出和
通常MUSCL格式写成如下形式
对应我们的推导结果〔二次多项式假设〕。
但是这不是最终形式。如果直接用这个公式，就会导致流场在激波〔连续〕附近的振荡。因为直接用二次多项式去逼近一个连续，会导致这样的效果。所以科学家们提出要对连续附近的斜率有所限制，因此引入了一个非常重要的修改——斜率限制器。参加斜率限制器后，上述公式就有点变化。

这里是Van Albada限制器
是一个小数〔〕，以防止分母为0。
密度和速度通过同样的方法来搞定。
密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的MUSCL。此外还有基于特征变量的MUSCL，要复杂一点，但是被认为适合更高精度的格式。然而一般计算中，基于原始变量的MUSCL由于具有足够的精度、简单的形式和较低的代价而被广泛使用。
OK，搞定了。下面进入第二点，怎么求。关于这一点，我不打算做详细介绍了，直接使用现有的近似黎曼解就可以了，都是通过和计算得到。比如Roe因为形式简单，而非常流行。在CFL3D软件主页〔cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html〕上看手册，附录C的C.1.3。
想了一下，还是把Roe求解器稍微说说吧，力求比拟完整。但是不要指望我把Roe求解器解释清楚，因为这个不是很容易三言两语说清的。
Roe求解器的数学形式是这样的

显然这个公式的第一项为哪一项一个中心差分形式，先前说过简单的中心差分不可行，原因是耗散不足导致振荡，说得通俗点就像一个弹簧，如果缺乏耗散〔阻尼〕它就会一直振荡。“耗散〞这个术语在激波捕捉格式中是最常见的。第二项的作用就是提供足够的耗散了。
这里和已经用MUSCL求得了，的定义在第一讲中已经介绍了。只有是还没说过的。

这个矩阵可以写成特征矩阵和特征向量矩阵的形式

而
，，和的具体表达式在许多书上都有，而且这里的矩阵表达有问题，所以就不写了。
是由、和代入计算得到。而、和采用所谓Roe平均值

这才是Roe求解器关键的地方！
总结一下，就是用Roe平均计算界面上的气体状态，然后计算得到，这样就可以得到了。如果有时间，我后面会找一个代码逐句分析一下。
总之，计算还是很不直接的。构造近似黎曼解是挺有学问的，需要对气体动力学的物理和数学方面有较深的理解。通常，如果不是做根底研究，你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导出来的。

附录程序：〔新建一个.c类型文件，将下面的程序复制粘贴到里面，就可以运行了〕
/****************************************************
本程序用于求解一维无粘可压缩欧拉方程〔激波管问题〕
运用Dummy Cell处理边界条件;
通量计算方式: AUSM Scheme;
重构方法：MUSCL方法
限制器：Van Albada限制器
时间离散：四步Runge-Kutta方法
****************************************************/
#include "stdio.h"
#include "conio.h"
#include "malloc.h"
#include "stdlib.h"
#include "math.h"
#include "string.h"
#define h (1/400.0) //网格步长
#define Nc 404 //网格总数：与h之间的关系 Nc=1/h+4
#define PI 3.1415927
#define It 1000
#define gama 1.4 //气体比热比
double KAKA=0.0;//限制器控制参数
double XS1,XS2;//计算域的两个端点
double dt=2.5e-5;//时间步长
double timesum;//总的计算时间
//函数声明
void output();
void SolveWtoU(double W[3],double U[3]);
void SolveUtoW(double W[3],double U[3]);//根本变量和守恒变量之间的转换函数

//前后各留两个网格单元作为虚拟网格单元 计算网格单元从2到NOc-3
struct cell
{
int flag;//网格点的类型
double xc;
double W[3],Wp[3];//conservation varaible
double U[3];//jibenbianliang
double R[3];
double S;//entropy
};
struct cell cell[Nc];

//网格生成与流场初始化
void initialsolve()
{
int i;
double x,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();
}
}

//边界条件：Dummy Cell方法
void boundarycondition()
{
int i;
double U[3];
//直接赋远场值
for(i=0;i<Nc;i++)
{
if(cell[i].flag==1)
{
cell[i].U[0]=0.445; cell[i].U[1]=0.698;cell[i].U[2]=3.528;
SolveUtoW(cell[i].W,U);
}
else if(cell[i].flag==2)
{
cell[i].U[0]=0.5;cell[i].U[1]=0.0;cell[i].U[2]=0.571;
SolveUtoW(cell[i].W,U);
}

}
}

//重构：MUSCL方法+Van Albada限制器
void SolveReconstruction(double Um1[],double UI[],double Up1[],double Up2[],double UL[],double UR[])
{
int i;
double epsilon,aR[3],aL[3],bL[3],bR[3],sL[3],sR[3];
epsilon=1.0e-5*h;
for(i=0;i<3;i++)
{
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)
}
for(i=0;i<3;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);
}
for(i=0;i<3;i++)
{
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]);
}
}

//守恒变量转化为根本变量
void SolveWtoU(double W[3],double U[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]);
}

//根本变量转化为守恒变量
void SolveUtoW(double W[3],double U[3])
{
W[0]=U[0];//d
W[1]=U[1]*U[0];//u
W[2]=U[2]/(gama-1.0)+0.5*U[0]*U[1]*U[1];
}

//格式：AUSM格式
void SolveAUSMFlux(double UL[],double UR[],double Fc[])//UL左状态，UR右状态，FC通量
{
double WL[3]={0.0},WR[3]={0.0};
double deta,fMn;
double ML,cl,MR,cr,Ml,Mr,Mn;
double pl,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;}
else if(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>=1.0)
{Mr=0.0;pr=0.0;}
else if(fabs(MR)<1.0)
{Mr=-0.25*(MR-1)*(MR-1);pr=0.25*PR*(MR-1)*(MR-1)*(2.0+MR);}
else
{Mr=MR;pr=PR;}

Mn=Ml+Mr;
pn=pl+pr;

if(fabs(Mn)>deta) fMn=fabs(Mn);
else fMn=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
void SolveResidual()
{
int i,j;
double UL[3],UR[3],Fp[3],Fm[3];
for(i=0;i<Nc;i++)
{
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<3;j++)
{
cell[i].R[j]=-(Fp[j]-Fm[j])/h;
}
}
}
}
//流场值更新
void SolveNextstep(double ar[],int ir)
{
int i,j;
for(i=0;i<=Nc;i++)
{
if(cell[i].flag==0)
{
if(ir==0)
{
for(j=0;j<3;j++)
cell[i].Wp[j]=cell[i].W[j];
}
for(j=0;j<3;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方法
void SolveRungeKutta()
{
double ar[4]={1/4.0,1/3.0,0.5,1.0};

int it,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=%f\n",it,timesum);
if(timesum>=0.16)
{
output();
getchar();
printf("Please enter any key to continue...");
break;
}
}
}
//结果输出
void output()
{
int i;
FILE *fpd,*fpu,*fpp;
if((fpd=fopen("Density.plt","w"))==NULL)
{
printf("connot open infile\n");
return;
}
fprintf(fpd,"TITLE = \"TestCase\"\n");
fprintf(fpd,"VARIABLES = \"x\", \"density\"\n");
fprintf(fpd,"ZONE T=\"Only Zone\", I=%d, F=POINT\n",Nc);

if((fpp=fopen("Pressure.plt","w"))==NULL)
{
printf("connot open infile\n");
return;
}
fprintf(fpp,"TITLE = \"TestCase\"\n");
fprintf(fpp,"VARIABLES = \"x\", \"pressure\"\n");
fprintf(fpp,"ZONE T=\"Only Zone\", I=%d, F=POINT\n",Nc);

if((fpu=fopen("Velocity.plt","w"))==NULL)
{
printf("connot open infile\n");
return;
}
fprintf(fpu,"TITLE = \"TestCase\"\n");
fprintf(fpu,"VARIABLES = \"x\", \"velocity\"\n");
fprintf(fpu,"ZONE T=\"Only Zone\", I=%d, F=POINT\n",Nc);

for(i=0;i<Nc;i++)
{
fprintf(fpd,"%f\t%f\n",cell[i].xc,cell[i].U[0]);
fprintf(fpu,"%f\t%f\n",cell[i].xc,cell[i].U[1]);
fprintf(fpp,"%f\t%f\n",cell[i].xc,cell[i].U[2]);
// 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();
}
fclose(fpd);
fclose(fpp);
fclose(fpu);
}

//网格类型判断
void SetMeshFlag()
{
int i;
for(i=0;i<Nc;i++)
{
if(cell[i].xc>=XS1&&cell[i].xc<=XS2)
cell[i].flag=0;
else if(cell[i].xc>XS2)
cell[i].flag=2;
else
cell[i].flag=1;
// printf("cell[%d].flag=%d\n",i,cell[i].flag);getchar();
}
}
//熵的计算
void SovleEntropy()
{
int i;
for(i=0;i<Nc;i++)
{
cell[i].S=cell[i].U[2]/pow(cell[i].U[0],gama);
}
}
//主函数
void main(void)
{
initialsolve();
SetMeshFlag();
SolveRungeKutta();
output();
}
18 / 18