CFD理论过渡到编程的最简单教程附源程序

CFD!论过渡到编程的傻瓜入门教程（注：这是一篇不知道谁写的介绍一维无粘可压缩Euler方程，以及如何编程实现求解该方程的论文。作者从最基本的概念出发，深入浅出的讲解了控制方程，有限体积格式，MSUCL法，限制器，Roe格式等相关知识。这篇论文我觉得有利于大家学习CF两程的相关知识，所以推荐给大家。文章的后面附有我写的程序（C语言），用于求解一维激波管问题，大家有兴趣可以看看（程序中加了注释说明）胡偶2011）借宝地写几个小短文，介绍CFD的一些实际的入门知识。主要是因为这里支持Latex,写起来比较方便。CFD,计算流体力学，是一个挺难的学科，涉及流体力学、数值分析和计算机算法，还有计算机图形学的一些知识。尤其是有关偏微分方程数值分析的东西，不是那么容易入门。大多数图书，偏重数学原理而不重实际动手，因为作者都把读者当做已经掌握基础知识的科班学生了。所以数学基础不那么好的读者往往看得很吃力，看了还不知道怎么实现。本人当年虽说是学航天工程的，但是那时本科教育已经退步，基础的流体力学课被砍得只剩下一维气体动力学了，因此自学CFD的时候也是头晕眼花。不知道怎么实现，也很难找到教学代码一一那时候网络还不发达，只在教研室的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着头皮分析。后来网上淘到一些代码研读，结合书籍论文才慢慢入门。可以说中间没有老师教，后来读博士为了混学分上过CFD专门课程，不过那时候我已经都掌握课堂上那些了。回想自己入门艰辛，不免有一个想法一一写点通俗易懂的CFD入门短文给师弟师妹们。本人不打算搞得很系统，而是希望能结合实际，阐明一些最基本的概念和手段，其中一些复杂的道理只是点到为止。目前也没有具体的计划，想到哪里写到哪里，因此可能会很零散。但是我争取让初学CFD的人能够了解一些基本的东西，看过之后，会知道一个CFD代码怎么炼成的（这炼”字好像很流行啊）。欢迎大家提出意见，这样我尽可能的可以追加一些修改和解释。言归正传，第一部分，我打算介绍一个最基本的算例，一维激波管问题。说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体状态（密度、速度、压力）不一样，突然把隔板抽去，管子内的气体怎么运动。这是个一维问题，被称作黎曼间断问题，好像是黎曼最初研究双曲微分方程的时候提出的一个问题，用一维无粘可压缩Euler方程就可以描述了。[整+需=0JI讥茂Z十P）_°}dt-.dz.-Ididz这里E=+eP=（7-这个方程就是描述的气体密度P、动量PU和能量pE随时间的变化（岛）与它们各自的流量（密度流量PU,动量流量P^+P,能量流量pEu+pU）随空间变化（吴）的关系。在CFD中通常把这个方程写成矢量形式这里PQ二pu_PE_puF—puu+ppEu-Fpu进一步可以写成散度形式需+V，F=O一定要熟悉这种矢量形式以上是控制方程，下面说说求解思路。可压缩流动计算中，有限体积(FVM)是最广泛使用的算法，其他算法多多少少都和FVM有些联系或者共同的思路。了解了FVM,学习其他高级点的算法(比如目前比较热门的间断有限元、谱FVM、谱FDM)就好说点了。针对一个微元控制体把Euler方程在空间积分C卷dH+rOFdf=。用微积分知识可以得到1乳―Fa—A每■十Az-⑪也就是说控制体内气体状态平均值的变化是控制体界面上流通量的结果因此我们要计算Q的演化，关键问题是计算控制体界面上的F。FVM就是以这个积分关系式出发，把整个流场划分为许多小控制体,每个控制体和周围相邻的某个控制体共享一个界面，通过计算每个界面上的通量来得到相邻控制体之间的影响，一旦每个控制体的变化得到，整个流场的变化也就知道了。所以，再强调一次，关键问题是计算控制体界面上的F初学者会说，这个不难，把界面『上的Q/插值得到，然后就可以计算F/。有道理！咱们画个图，有三个小控制体i-1至W+1,中间的“|装示界面，控制体i右边的界面用i+3表示，左边的就是"3。Ii-1IiIi+1I好，下个问题：每个小控制体长度都是Aio如何插值计算界面1+2上的Qm?最自然的想法就是：取两边的平均值呗，□.-Q—Q炉31+/一2但是很不幸，这是不行的。那么换个方法？直接平均得到Fi+另F计一还是很不行，这样也不行。我靠，这是为什么？这明明是符合微积分里面的知识啊？这个道理有点复杂，说开了去可以讲一本书，可以说从50年代到70年代，CFD科学家就在琢磨这个问题。这里，初学者只需要记住这个结论：对于流动问题，不可以这样简单取平均值来插值或者差分。如果你非要想知道这个究竟，我现在也不想给你讲清楚，因为我眼下的目的是让你快速上手，而且该不刨根问底的时候就不要刨根问底，这也是初学阶段一种重要的学习方法。好了，既然目的只是为了求Fi+上我在这里，只告诉你一种计算方法,也是非常重要、非常流行的一种方法。简单的说，就是假设流动状态在界面/是不连续的，先计算出界面a+/两边Q的值，Q叶当和^^工,再由它们用某种方法计算出F叶上上述方法是非常重要的，是由一个苏联人Godunov在50年代首创的，后来被发展成为通用Godunov方法，著名的ENO/WENO就是其中的一种。好了，现在的问题是：i怎么确定计去和2怎么计算F计*对于第一个问题，Godunov在他的论文中，是假设每个控制体中Q是Q".=Q.均匀分布的，因此.•：•Qj-i=Q计1第二个问题，Godunov采用了精确黎曼解来计算卫计上什么是精确黎曼解"，就是计算这个激波管问题的精确解。既然有精确解，那还费功夫搞这些FVM算法干什么？因为只有这种简单一维问题有精确解，稍微复杂一点就不行了。精确解也比较麻烦，要分析5种情况，用牛顿法迭代求解（牛顿法是什么？看数值计算的书去，哦，算了，现在暂时可以不必这是最初Godunov的方法，后来在这个思想的基础上，各种变体都Q-+p出来了。也不过是在这两个问题上做文章，怎么确定叶,，怎么计算计厘。Godunov假设的是每个小控制体内是均匀分布，也就是所谓分段常数(piecewiseconstant),所以后来有分段线性(picewiselinear)或者分段二次分布(picewiseparabolic),到后来ENO/WENO出来，那这个假设的多项式次数就继续往上走了。都是用多项式近似的，这是数值计算中的一个强大工具，你可以在很多地方看到这种近似。Godunov用的是精确黎曼解，太复杂太慢，也不必要，所以后来就有各种近似黎曼解,最有名的是Roe求解器、HLL求解器和Osher求解器，都是对精确黎曼解做的简化。这个多项式的阶数是和计算精度密切相关的，阶数越高，误差就越小。不过一般来说，分段线性就能得到不错的结果了，所以工程中都是用这个，Fluent、Fastran以及NASA的CFL3D、OverFlow都是用这个。而黎曼求解器对精度的影响不是那么大，但是对整个算法的物理适用性有影响，也就是说某种近似黎曼求解器可能对某些流动问题不合适，比如单纯的Roe对于钝头体的脱体激波会算得乱七八糟，后来加了嫡修正才算搞定。上次(/node/399)说到了求解可压缩流动的一个重要算法，通用Godunov方法。其两个主要步骤就是Q-工QJ1怎么确定叶蔓和叶^2怎么计算这里我们给出第一点一个具体的实现方法，就是基于原始变量的MUSCL格式(以下简称MUSCL)。它是一种很简单的格式，而且具有足够的精度，NASA著名的CFL3D软件就是使用了这个格式，大家可以去它的主页（/Cfl3dv6/cfl3dv6.html）上看手册，里面空间离散那一章清楚的写着。MUSCL假设控制体内原始变量（就是Q）的分布是一次或者二次多项式，如果得到了这个多项式，就可以求出控制体i左右两个界面的一侧的值Q计当和Qt。我们以压力p为例来说明怎么构造这个多项式。这里我只针对二次多项式来讲解，你看完之后肯定能自己推导出一次多项式的结果（如果你搞不定，那我对你的智商表示怀疑）。OK,开始假设Ar=1,这个假设不影响最终结论，因为你总可以把一个区间线性的变换到长度为1的区间。假设压力p在控制体i内部的分布是一个二次多项式/（P）=&/+bE+c,控制体i的中心处于工二0处，左右两个界面就是-1.xl2和十2。这里先强调一个问题，在FVM中，每个控制体内的求解出来的变量Q实际上是这个控制体内的平均值--JTy。所以，Pi=（at2+bi+c）dz=ar3/3+bxa/2+ext==+c。我们知道力T,玄和瓦十1,等距网格情况下1-十之处的导数可以近似表示为言|工_产△「二论一佚7罪|计产纯二讪+1-氏那么fO（2d+&）|H=T=b-a=△一（这里错了，应该是2ax+b）『（+；）=（2a/+圳h=+产匕+a=△+由上述三个有关a,b和c的方程，我们可以得到a=^±rb=^--0=M-24十1-

却和1一十1-

却和1一2-叶1-2户口=曷+膜一目午+。+能）4］?Jli=R-i［（l-k）At+（l+Mk二1/3对应我们的推导结果（二次多项式假设）但是这不是最终形式。如果直接用这个公式，就会导致流场在激波（间断）附近的振荡。因为直接用二次多项式去逼近一个间断，会导致这样的效果。所以科学家们提出要对间断附近的斜率有所限制，因此引入了一个非常重要的修改一一斜率限制器。加入斜率限制器后，上述公式就有点变化。P"1=pt+却"网+（1+AM叱工=pt-^［（1-kSi）^+（1+左5*心］这里$工是VanAlbada限制器G一+2‘一（7户转G*£是一个小数（E=1x10一外，以防止分母为0。密度和速度通过同样的方法来搞定。密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的MUSCLo此外还有基于特征变量的MUSCL,要复杂一点，但是被认为适合更高精度的格式。然而一般计算中，基于原始变量的MUSCL由于具有足够的精度、简单的形式和较低的代价而被广泛使用。OK,搞定了。下面进入第二点，怎么求F叶丸关于这一点，我不打Q-1算做详细介绍了，直接使用现有的近似黎曼解就可以了，都是通过叶多和Q'*计算得到*+寺。比如Roe因为形式简单，而非常流行。在CFL3D软件主页（/Cfl3dv6/cfl3dv6.html）上看手册，附录C的C.1.3o想了一下，还是把Roe求解器稍微说说吧，力求比较完整。但是不要指望我把Roe求解器解释清楚，因为这个不是很容易三言两语说清的。Roe求解器的数学形式是这样的F3=外F（QQ）+F（Q\J]一扑*式口3一Qq）显然这个公式的第一项是一个中心差分形式，先前说过简单的中心差分不可行，原因是耗散不足导致振荡，说得通俗点就像一个弹簧，如果缺乏耗散（阻尼）它就会一直振荡。耗散”这个术语在激波捕捉格式中是最常见的。第二项的作用就是提供足够的耗散了。这里Q叶斗和Q叶4已经用MUSCL求得了，F的定义在第一讲中已经介绍了。只有।闵是还没说过的。

这个矩阵可以写成特征矩阵和特征向量矩阵的形式A=R-AR+而A|=R-|A|R+A,R,R一和A的具体表达式在许多书上都有，而且这里的矩阵表达有问题，所以就不写了。⑶卜+吉是由P计；、"叶多和分制代入A计算得到。而凡制、“叶印和气忖采用所谓Roe平均值_1,十^十十，7二「1+>_1,十^十十，7二「1+>~行+京%=这才是Roe求解器关键的地方！总结一下，就是用Roe平均计算界面上的气体状态Q计九然后计算得到|』计号|,这样F计孑就可以得到了。如果有时间，我后面会找一个代码逐句分析一下。总之，计算F[+孑还是很不直接的。构造近似黎曼解是挺有学问的，需要对气体动力学的物理和数学方面有较深的理解。通常，如果不是做基础研究，你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导出来的附录程序：（新建一个.c类型文件，将下面的程序复制粘贴到里面，就可以运行了）/****************************************************本程序用于求解一维无粘可压缩欧拉方程（激波管问题）运用DummyCell处理边界条件；通量计算方式：AUSMScheme;重构方法：MUSCL方法限制器：VanAlbada限制器时间离散：四步Runge-Kutta方法****************************************************/#include"stdio.h"#include"conio.h"#include"malloc.h"#include"stdlib.h"#include"math.h"#include"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-3structcell{intflag;//网格点的类型doublexc;doubleW[3],Wp[3];//conservationvaraibledoubleU[3];//jibenbianliangdoubleR[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()(inti;doubleU[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);}elseif(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方法+VanAlbada限制器Up2[],doublevoidSolveReconstruction(doubleUm1[],doubleUI[],doubleUp1[],doubleUL[],doubleUR[])Up2[],double{inti;doubleepsilon,aR[3],aL[3],bL[3],bR[3],sL[3],sR[3];epsilon=1.0e巧*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]);))〃守恒变量转化为基本变量voidSolveWtoU(doubleW[3],doubleU[3]){U[0]=W[0];//dU[1]=W[1]/W[0];//uU[2]=(gama-1.0)*(W[2]-0.5*U[0]*U[1]*U[1]);)〃基本变量转化为守恒变量voidSolveUtoW(doubleW[3],doubleU[3]){W[0]=U[0];//dW[1]=U[1]*U[0];//uW[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>=1.0){Mr=0.0;pr=0.0;}elseif(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);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);}//残差af算RvoidSolveResidual(){inti,j;doubleUL[3],UR[3],Fp[3],Fm[3];for(i=0;i<Nc;i++){if(cell[i].flag==0){//i+1/2SolveReconstruction(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/2SolveReconstruction(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;}}}}〃流场值更新voidSolveNextstep(doublear[],intir){inti,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方法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=%f\n",it,timesum);if(timesum>=0.16){output();getchar();printf("Pleaseenteranykeytocontinue...");break;}}}〃结果输出voidoutput(){inti;FILE*fpd,*fpu,*fpp;if((fpd=fopen("Density.plt","w"))==NULL){printf("connotopeninfile\n");return;)fprintf(fpd,"TITLE=\"TestCase\"\n");fprintf(fpd,"VARIABLES=\"x\",\"density\"\n");fprintf(fpd,"ZONET=\"OnlyZone\",I=%d,F=POINT\n",Nc);if((fpp=fopen("Pressure.plt","w")