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

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 对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定。上次（http:/gezhi.org/node/399）说到了求解可压缩流动的一个重要算法，通用 Godunov 方法。其两个主要步骤就是1 怎么确定2 怎么计算和这里我们给出第一点一个具体的实现方法，就是基于原始变量的 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 格式写成如下形式对应我们的推导结果（二次多项式假设）。但是这不是最终形式。如果直接用这个公式，就会导致流场在激波（间断）附近的振荡。因为直接用二次多项式去逼近一个间断，会导致这样的效果。所以科学家们提出要对间断附近的斜率有所限制，因此引入了一个非常重要的修改斜率限制器。加入斜率限制器后，上述公式就有点变化。这里是 Van Albada 限制器是一个小数（），以防止分母为 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 求解器关键的地方！总结一下，就是用总结一下，就是用RoeRoe平均计算界面上的气体状态平均计算界面上的气体状态，这样，这样句分析一下。，然后计算得到，然后计算得到就可以得到了。就可以得到了。如果有时间，我后面会找一个代码逐总之，计算还是很不直接的。构造近似黎曼解是挺有学问的，需要对气体动力学的物理和数学方面有较深的理解。通常，如果不是做基础研究，你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导出来的。附录程序：（新建一个.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 W3,double U3);void SolveUtoW(double W3,double U3);/基本变量和守恒变量之间的转换函数/前后各留两个网格单元作为虚拟网格单元 计算网格单元从 2 到 NOc-3struct cellint flag;/网格点的类型double xc;double W3,Wp3;/conservation varaibledouble U3;/jibenbianliangdouble R3;double S;/entropy;struct cell cellNc;/网格生成及流场初始化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;iNc;i+)celli.xc=(xi+h/2)+i*h;/网格中心坐标celli.flag=0;x=celli.xc;if(x0.5)celli.U0=0.445;celli.U1=0.698;celli.U2=3.528;elsecelli.U0=0.5;celli.U1=0.0;celli.U2=0.571;SolveUtoW(celli.W,celli.U);/printf(x=%f,d=%f,u=%f,p=%fn,celli.xc,celli.U0,celli.U1,celli.U2);getchar();/边界条件：Dummy Cell 方法void boundarycondition()int i;double U3;/直接赋远场值for(i=0;iNc;i+)if(celli.flag=1)celli.U0=0.445;celli.U1=0.698;celli.U2=3.528;SolveUtoW(celli.W,U);else if(celli.flag=2)celli.U0=0.5;celli.U1=0.0;celli.U2=0.571;SolveUtoW(celli.W,U);/重构：MUSCL 方法+Van Albada限制器voidSolveReconstruction(doubleUm1,doubleUI,doubleUp1,doubleUp2,doubleUL,double UR)int i;double epsilon,aR3,aL3,bL3,bR3,sL3,sR3;epsilon=1.0e-5*h;for(i=0;i3;i+)aRi=Up2i-Up1i;/Dealta+U(I+1)bRi=Up1i-UIi;/Dealta_U(I+1)aLi=Up1i-UIi;/Dealta+U(I)bLi=UIi-Um1i;/Dealta_U(I)for(i=0;i3;i+)sRi=(2.0*aRi*bRi+1.0e-6)/(aRi*aRi+bRi*bRi+1.0e-6);sLi=(2.0*aLi*bLi+1.0e-6)/(aLi*aLi+bLi*bLi+1.0e-6);for(i=0;i=1.0)Ml=ML;pl=PL;else if(fabs(ML)=1.0)Mr=0.0;pr=0.0;else if(fabs(MR)deta)fMn=fabs(Mn);else fMn=0.5*(fabs(Mn)*fabs(Mn)+deta*deta)/deta;Fc0=0.5*Mn*(WL0*cl+WR0*cr)-0.5*fMn*(WR0*cr-WL0*cl);Fc1=0.5*Mn*(WL1*cl+WR1*cr)-0.5*fMn*(WR1*cr-WL1*cl)+pn;Fc2=0.5*Mn*(WL2+PL)*cl+(WR2+PR)*cr)-0.5*fMn*(WR2+PR)*cr-(WL2+PL)*cl);/残差计算 Rvoid SolveResidual()int i,j;double UL3,UR3,Fp3,Fm3;for(i=0;iNc;i+)if(celli.flag=0)/i+1/2SolveReconstruction(celli-1.U,celli.U,celli+1.U,celli+2.U,UL,UR);/MSUCL 重构SolveAUSMFlux(UL,UR,Fp);/二阶格式/SolveAUSMFlux(celli.U,celli+1.U,Fp);/一阶格式/i-1/2SolveReconstruction(celli-2.U,celli-1.U,celli.U,celli+1.U,UL,UR);/MSUCL 重构SolveAUSMFlux(UL,UR,Fm);/二阶格式/SolveAUSMFlux(celli-1.U,celli.U,Fm);/一阶格式for(j=0;j3;j+)celli.Rj=-(Fpj-Fmj)/h;/流场值更新void SolveNextstep(double ar,int ir)int i,j;for(i=0;i=Nc;i+)if(celli.flag=0)if(ir=0)for(j=0;j3;j+)celli.Wpj=celli.Wj;for(j=0;j3;j+)celli.Wj=celli.Wpj+arir*dt*celli.Rj;SolveWtoU(celli.W,celli.U);/Runge-Kutta 方法void SolveRungeKutta()double ar4=1/4.0,1/3.0,0.5,1.0;int it,ir;timesum=0.0;for(it=0;it+)/迭代步数for(ir=0;ir=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 infilen);return;fprintf(fpd,TITLE=TestCasen);fprintf(fpd,VARIABLES=x,densityn);fprintf(fpd,ZONE T=Only Zone,I=%d,F=POINTn,Nc);if(fpp=fopen(Pressure.plt,w)=NULL)printf(connot open infilen);return;fprintf(fpp,TITLE=TestCasen);fprintf(fpp,VARIABLES=x,pressuren);fprintf(fpp,ZONE T=Only Zone,I=%d,F=POINTn,Nc);if(fpu=fopen(Velocity.plt,w)=NULL)printf(connot open infilen);return;fprintf(fpu,TITLE=TestCasen);fprintf(fpu,VARIABLES=x,velocityn);fprintf(fpu,ZONE T=Only Zone,I=%d,F=POINTn,Nc);for(i=0;iNc;i+)fprintf(fpd,%ft%fn,celli.xc,celli.U0);fprintf(fpu,%ft%fn,celli.xc,celli.U1);fprintf(fpp,%ft%fn,celli.xc,celli.U2);/printf(x=%f,d=%f,u=%f,p=%fn,celli.xc,celli.U0,celli.U1,celli.U2);getchar();fclose(fpd);fclose(fpp);fclose(fpu);/网格类型判断void SetMeshFlag()int i;for(i=0;i=XS1&celli.xcXS2)celli.flag=2;elsecelli.flag=1;/printf(cell%d.flag=%dn,i,celli.flag);getchar();/熵的计算void SovleEntropy()int i;for(i=0;iNc;i+)celli.S=celli.U2/pow(celli.U0,gama);/主函数void main(void)initialsolve();SetMeshFlag();SolveRungeKutta();output();