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

《CFD理论过渡到编程的傻瓜入门教程.docx》由会员分享，可在线阅读，更多相关《CFD理论过渡到编程的傻瓜入门教程.docx（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精品word 可编辑资料 - - - - - - - - - - - - -.CFD 理论过渡到编程的傻瓜入门教程注：这是一篇不知道谁写的介绍一维无粘可压缩Euler 方程，以及如何编程实现求解该方程的论文；作者从最根本的概念动身， 深化浅出的讲解了掌握方程，有限体积格式，MSUCL 方法，限制器， Roe格式等相关学问；这篇论文我觉得有利于大家学习CFD 编程的相关学问，所以举荐给大家；文章的后面附有我写的程序C 语言，用于求解一维激波管问题， 大家有爱好可以看看 程序中加了注释说明胡偶 2021借宝地写几个小短文， 介绍 CFD 的一些实际的入门学问；主要是由于这里支持 Latex，写起来

2、比拟便利；CFD ，运算流体力学，是一个挺难的学科，涉及流体力学、数值分析和运算机算法，仍有运算机图形学的一些学问；特别是有关偏微分方程数值分析的东西，不是那么简洁入门；大多数图书，片中数学原理而不重实际动手，由于作者都把读者当做已经把握根底学问的科班同学了；所以数学根底不那么好的读者往往看得很吃力，看了仍不知道怎么实现；本人当年虽说是学航天工程的，但是那时本科训练已经退步，根底的流体力学课被砍得只剩下一维气体动力学了，因此自学 CFD 的时候也是头晕眼花； 不知道怎么实现，也很难找到教学代码那时候网络仍不兴盛，只在教研室的故纸堆里搜罗到一些完全没有注释，编程风格也不好的冗长代码，硬着头皮分析

3、；后来网上淘到一些代码研读，结合书籍论文才渐渐入门；可以说-.可修编 - .第 1 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.中间没有老师教， 后来赌博士为了混学分上过CFD 特地课程， 不过那时候我已经都把握课堂上那些了；回想自己入门艰辛， 不免有一个想法写点通俗易懂的CFD 入门短文给师弟师妹们；本人不准备搞得很系统，而是期望能结合实际，说明一些最根本的概念和手段，其中一些复杂的道理只是点到为止；目前也没有具体的方案，想到哪里写到哪里，因此可能会很零散；但是我争取让初学CFD 的人能够明白一些根本的

4、东西，看过之后， 会知道一个 CFD 代码怎么炼成的这“炼字好像很流行啊；欢送大家提出看法，这样我尽可能的可 以追加一些修改和说明；言归正传，第一局部，我准备介绍一个最根本的算例，一维激波管问题；说白了就是一根两端封闭的管子，中间有个隔板，隔板左边和右边的气体状态密度、速度、压力不一样，突然把隔板抽去，管子面的气体怎么运动；这是个一维问题，被称作黎曼连续问题 ，好像是黎曼最初讨论双曲微分方程 的时候提出的一个问题，用一维无粘可压缩Euler 方程就可以描述了；这里这个方程就是描述的气体密度、动量和能量随时间的变化-.可修编 - .第 2 页，共 18 页 - - - - - - - - - -

5、精品word 可编辑资料 - - - - - - - - - - - - -.与它们各自的 流量 密度流量，动量流量，能量流量随空间变化的关系；在 CFD 常把这个方程写成矢量形式这里进一步可以写成散度形式肯定要熟识这种矢量形式以上是掌握方程 ，下面说说 求解思路 ；可压缩流淌 运算中，有限体积FVM 是最广泛使用的算法，其他算法多多少少都和FVM 有些联系或者共通的 思路；明白的FVM ，学习其他高级点的算法比方目前比拟热门的连续有 限元、谱 FVM 、谱 FDM 就好说点了；针对一个微元掌握体，把 Euler 方程在空间积分用微积分学问可以得到也就是说 掌握体气体状态平均值的变化是掌握体界

6、面上流通量的结果；因此我们要运算的演化，关键问题是运算掌握体界面上的；FVM 就是以这个积分关系式动身，把整个流场划分为很多小掌握体，每个掌握体和四周相邻的某个掌握体共享一个界面，通过运算每个界面上的通-.可修编 - .第 3 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.量来得到相邻掌握体之间的影响，一旦每个掌握体的变化得到，整个流场的变化也就知道了；所以，再强调一次，关键问题是运算掌握体界面上的；初学者会说， 这个不难， 把界面上的插值得到， 然后就可以运算；有道理！咱们画个图，有三个小掌握体i-1 到

7、 i+1 ，中间的“| 表示界面，掌握体i 右边的界面用表示，左边的就是；| i-1 | i | i+1 |好下个问题：每个小掌握体长度都是如何插值运算界面上的？最自然的想法就是：取两边的平均值呗，但是很不幸，这是不行的；那么换个方法？直接平均得到？仍是很不行，这样也不行；我靠，这是为什么？这明明是符合微积分里面的学问啊？-.可修编 - .第 4 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.这个道理有点复杂， 说开了去可以讲一本书，可以说从 50 岁月到 70 岁月，CFD 科学家就在琢磨这个问题；这里，

8、初学者只需要记住这个结论：对于流淌问题 ，不行以这样简洁取平均值来插值或者差分；假如你非要想知道这个终究， 我现在也不想给你讲清晰，由于我眼下的目的是让你快速上手，而且该不追根究底的时候就不要追根究底，这也是初学阶段一种重要的学习方法；好了，既然目的只是为了求，我在这里，只告知你一种运算方法，也是特别重要、特别流行的一种方法；简洁的说，就是假设流淌状态在界面是不连续的，先运算出界面两边的值，和，再由它们用某种方法运算出；上述方法是特别重要的，是由一个联人Godunov 在 50 岁月首创的，后来被开展成为通用Godunov 方法，闻名的ENO/WENO就是其中的一种；好了，现在的问题是：1 怎

9、么确定和2 怎么运算对于第一个问题，Godunov 在他的论文中，是假设每个掌握体中是匀称分布的，因此其次个问题， Godunov 采纳了 精确黎曼解 来运算；什么是“精确黎曼解，就是运算这个激波管问题的精确解；既然有精确解，那仍费功夫搞这些 FVM 算法干什么？由于只有这种简洁一维问题有精确解，略微复杂-.可修编 - .第 5 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.一点就不行了；精确解也比拟麻烦，要分析5 种情形，用 牛顿法迭代求解牛顿法是什么？看数值运算的书去，哦， 算了， 现在临时可以不必看；

10、这是最初 Godunov 的方法，后来在这个思想的根底上，各种变体都出来了；也不过是在这两个问题上做文章，怎么确定，怎么运算；Godunov 假设的是每个小掌握体是匀称分布，也就是所谓分段常数(piecewise constant，)所以后来有分段线性(picewise linear或)者分段二次分布(picewise parabolic，)到后来 ENO/WENO出来，那这个假设的多项式次数就连续往上走了； 都是用 多项式近似 的，这是数值运算中的一个强大工具，你可以在很多地方看到这种近似；Godunov 用的是精确黎曼解，太复杂太慢，也不必要，所以后来就有各种近似黎曼解 ，最出名的是Roe

11、 求解器、 HLL求解器和 Osher 求解器 ，都是对精确黎曼解做的简化；这个多项式的阶数是和运算精度亲密相关的，阶数越高，误差就越小；不过一般来说，分段线性就能得到不错的结果了，所以工程中都是用这个，Fluent、Fastran以及 NASA 的 CFL3D 、OverFlow 都是用这个；而 黎曼求解器对精度的影响不是那么大，但是对整个算法的物理适用性有影响，也就是说某种近似黎曼求解器可能对某些流淌问题不相宜，比方单纯的Roe 对于钝头体的脱体激波会算得乱七八糟，后来加了熵修正才算搞定；上次 gezhi.org/node/399 说到了求解可压缩流淌的一个重要算法，通用Godunov 方

12、法；其两个主要步骤就是-.可修编 - .第 6 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.1 怎么确定和2 怎么运算这里我们给出第一点一个具体的实现方法，就是基于原始变量的MUSCL格式 以下简称MUSCL；它是一种很简洁的格式，而且具有足够的精 度，NASA 闻名的 CFL3D 软件就是使用了这个格式，大家可以去它的主页cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html 上看手册，里面空间离散那一章清晰的写着；MUSCL 假设掌握体原始变量就是的分布是一次或者二次多项式，

13、假如得到了这个多项式，就可以求出掌握体左右两个界面的一侧的值和；我们以压力p 为例来说明怎么构造这个多项式；这里我只针对二次多项式来讲解， 你看完之后确定能自己推导出一次多项式的结果假如你搞不定，那我对你的智商表示疑心；OK ，开场假设，这个假设不影响最终结论，由于你总可以把一个区间线性的变换到长度为1 的区间；假设压力 p 在掌握体 i 部的分布是一个二次多项式，掌握体 i 的中心处于处，左右两个界面就是和；这里先强调一个问题，在FVM 中，每个掌握体的求解出来的变量实际上是这个掌握体的平均值；-.可修编 - .第 7 页，共 18 页 - - - - - - - - - -精品word 可

14、编辑资料 - - - - - - - - - - - - -.所以，；我们知道，和，等距网格情形下和处的导数可以近似表示为那么这里错了，应当是2ax+b由上述三个有关a， b 和 c 的方程，我们可以得到这样就可以得到f(x)的表达式了，由此可以算出和通常 MUSCL 格式写成如下形式对应我们的推导结果二次多项式假设；但是这不是最终形式；假如直接用这个公式，就会导致流场在激波 连续邻近的振荡 ；由于直接用二次多项式去靠近一个连续，会导致这样的成效；所以科学家们提出要对连续邻近的斜率有所限制，因此引入了一个特别重要的修改 斜率限制器 ；参与斜率限制器后，上述公式就有点变化；-.可修编 - .第

15、8 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.这里是 Van Albada 限制器是一个小数，以防止分母为0；密度和速度通过同样的方法来搞定；密度、速度和压力被称作原始变量，所以上述方法是基于原始变量的MUSCL；此外仍有基于特点变量的MUSCL，要复杂一点，但是被认为适合更高精度的格式；然而一般运算中，基于原始变量的MUSCL 由于具有足够的精度、简洁的形式和较低的代价而被广泛使用；OK ，搞定了；下面进入其次点，怎么求；关于这一点，我不准备做具体介绍了， 直接使用现有的近似黎曼解就可以了，都是通过和

16、运算得到；比方 Roe 由于形式简洁，而特别流行；在CFL3D 软件主页 cfl3d.larc.nasa.gov/Cfl3dv6/cfl3dv6.html 上看手册，附录C 的 C.1.3 ；想了一下，仍是把Roe 求解器略微说说吧，力求比拟完整；但是不要盼望我把 Roe 求解器说明清晰，由于这个不是很简洁三言两语说清的；Roe 求解器的数学形式是这样的明显这个公式的 第一项为哪一项一个中心差分形式，从前说过简洁的中心差分不行行，缘由是耗散缺乏导致振荡，说得通俗点就像一个弹簧，假如缺乏耗散阻尼它就会始终振荡；“耗散这个术语在激波捕获格式中是最常见的； 其次项的作用就是供应足够的耗散了；-.可修

17、编 - .第 9 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.这里了；只有和已经用 MUSCL 求得了，是仍没说过的；的定义在第一讲中已经介绍这个矩阵可以写成特点矩阵和特点向量矩阵的形式而，和的具体表达式在很多书上都有，而且这里的矩阵表达有问题，所以就不写了；是由、和代入运算得到；而、和采纳所谓 Roe 平均值这才是 Roe 求解器关键的地方！总结一下，就是用Roe 平均运算界面上的气体状态，然后运算得到，这样就可以得到了； 假如有时间，我后面会找一个代码逐句分析一下；总之，运算仍是很不直接的；构造近似黎

18、曼解是挺有学问的，需要对 气体动力学的物理和数学方面有较深的懂得；通常，假如不是做根底讨论，你只需要知道它们的特点，会用它们就可以了，而不必深究它们怎么推导出来的；-.可修编 - .第 10 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -.附录程序：新建一个 .c 类型文件，将下面的程序复制粘贴到里面， 就可以运行了/*本程序用于求解一维无粘可压缩欧拉方程激波管问题运用 Dummy Cell 处理边界条件;通量运算方式: AUSM Scheme;重构方法： MUSCL 方法 限制器： Van Albada 限

19、制器时间离散：四步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

20、 dt=2.5e-5;/ 时间步长double timesum;/ 总的运算时间/ 函数声明void output();void SolveWtoU(double W3,double U3);void SolveUtoW(double W3,double U3);/根本变量和守恒变量之间的转换函数-.可修编 - .第 11 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -./ 前后各留两个网格单元作为虚拟网格单元运算网格单元从2 到NOc-3struct cellint flag;/ 网格点的类型double

21、xc;double W3,Wp3;/conservation varaible double 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)elsecelli.U0=0.445

22、; celli.U1=0.698;celli.U2=3.528;celli.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+)-.可修编 - .第 12 页，共 18 页 - - - - - - - - -

23、 -精品word 可编辑资料 - - - - - - - - - - - - -.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 限制器void SolveReconstruction(double Um1,double UI,double Up1,double

24、 Up2,double UL,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)/(a

25、Li*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+(WR

26、2+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/2SolveReconst

27、ruction(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;-.可修编 - .第 15 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -./ 流场值更新void SolveNextstep(double ar,int ir)i

28、nt 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(); getcha

29、r();printf(Please enter any key to continue.);break;-.可修编 - .第 16 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -./ 结果输出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,VARIABLE

30、S = 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(co

31、nnot 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

32、=%fn,celli.xc,celli.U0,celli.U1,celli.U2);getchar();fclose(fpd); fclose(fpp); fclose(fpu);-.可修编 - .第 17 页，共 18 页 - - - - - - - - - -精品word 可编辑资料 - - - - - - - - - - - - -./ 网格类型判定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();-.可修编 - .第 18 页，共 18 页 - - - - - - - - - -