用有限体积方法求解欧拉方程

有限体积法求解二维可压缩Euler方程——计算流体力学课程大作业老师：夏健、刘学强学生：徐锡虎学号：SQ09018013018日期：2010年2月5日-1-目录一、内容摘要……………………………………………………（2）二、流动控制方程………………………………………………（2）三、有限体积法的空间离散……………………………………（2）四、人工耗散……………………………………………………（3）五、时间离散……………………………………………………（4）六、边界条件……………………………………………………（5）七、计算结果……………………………………………………（8）八、结论与展望…………………………………………………（11）参考文献…………………………………………………………（11）-2-一、内容摘要本文通过运用JAMESON有限体积法求解了二维定常和非定常可压缩Euler方程。程序实现语言为C++。其中，使用的网格是三角形非结构网格。在时间推进上使用的是四步龙—库塔推进格式。推进的时间步长取的是当地的时间步长。为了消除迭代误差、round-off等误差，本文采用了添加人工耗散项的办法。另外，本文计算了NACA0012翼型在跨音速下不同迎角的情况，并与fluent软件的计算结果进行了比较，来验证程序的准确性。二、流动控制方程守恒形式的Euler方程：0GdxFdywdtS（1）其中x和y代表笛卡儿坐标系。W是守恒变量。EVUW(2)F,G表示通量UHUVPUUF2,VHPVUVVG2(3),P,H和E表示密度，压强，单元总焓和单元总能量。U,V表示笛卡儿坐标系下的速度矢量。这些量由理想气体的单位体积的总能量和总焓相互联系。2/122）（）（VUPE(4)PEH(5)三、有限体积法的空间离散计算域被划分为互不重叠的单元。在每个单元运用守恒形式的Euler方程。由于每个单元相对于时间都是不变的，所以等式（1）可以写成：dGdxFdytWS)(（8）其中和S是单元的体积和边界。W是单元的平均值。在对上述方程进行时间离散前，先对空间进行离散，则方程（6）可以写为：kkkQdtdW（9）-3-其中k表示第k个单元的体积，kW是第k个单元的守恒变量。kQ表示第k个单元的通量。方程（7）的右边项可以写成：kedgesiikxGyFQ1)(（10）其中abiabiyyyxxx,（11）（8）式中的求和是对第k个单元的所有边进行的。守恒参数的量是单元中心值，在求通量时，第I条边的守恒参数值是用左右单元的平均来表示的：2/)(WpkiWW（12）引入变量：iiiiixVyUZ（13）则第k单元的Euler方程可以写为：kedgesiikkHZxPVZyPUZZdtdW11（14）在本文中，采用的是JAMENSON有限体积法，为了减少存储的相关信息的量，其存储的方式选择的是按边存储的方法。在存储的每条边的信息中，包含了这条边的边号，左右单元号和边的端点。在计算通量时采用按边循环的方式：doI=1,nedgek=connmatrix(I,1)a=connmatrix(I,2)b=connmatrix(I,3)p=connmatrix(I,4)flux=function(k,a,b,p)sum(k)=sum(k)+fluxsum(p)=sum(p)-fluxenddo这里给出的是FOTRAN语言的形式，我编写采用的是C，具体表现在上交的程序中。在计算时间步长、人工耗散项等也可用象这样按边循环，从此处我们可以看出求解时与单元的形状无关。四、人工耗散人工粘性模型对方法的成功应用起着关键作用，人工粘性抑制解在激波附近的振荡，又阻尼解在光滑区域内的高阶误差，对解的线性稳定和收敛于定态是很重要的。本文在方程（14）的右端加入了人工耗散项，如对于单元k，其表达式可以表示为：kkkkDQdtdW（15）-4-在有限体积法中，耗散项的公式可以表示为：kedgesiikedgesiikddD1)4(1)2(（16）其中：ikpiiiikpiiiWWdWWd22)4()4()2()2(（17）其中的I表示单元k和p的公共边，2定义为：kedgesikjk（18）上面的j表示与k相邻的单元。kedgesiikpkedgesiikpkPPPP11)()(（19）),0max(),max()2()4()4()2()2(iiikpikk（20）其中的量)4()2(,kk的范围是：0.121,3212561)2()4(kk。在计算时发现上面方法得到的人工耗散项并不太适合。其在光滑区域耗散项太大，而在大剃度区域又显得太小，为了弥补上面的不足，作下面的修改：kpkpiPPPP（21）自适应系数为：),0max()2()4()4()2()2(iiiikk（22）尺度系数为：22yxcxVyUi（23）其中的U,V表示边上的值，C表示当地声速。五、时间离散方程最后的稳定解是通过时间上的迭代得到的，可以写为：kkRdtdW（24）右边项的表达式为：-5-kkkkDQR/（25）为了加速收敛，时间迭代使用的是4步龙—库塔推进格式。格式如下：)4(1)1()0()()0(4..1.....WWtomfortR（26）其中的n表示的是当前的时间步，n+1表示的是新的时间步：)0()()(DQRmm（27）1,21,31,414321（28）为了减小计算时间，人工耗散项的计算只在第一步进行，在下面几步的迭代中保持不变。运用上面的方法计算，可以发现CFL数可以取到22，本文中使用的是2.0。使用显示格式迭代的主要缺点是由于稳定区域的限制，所以不能使用过大的时间步长。可以用近似的方法估算时间步长，对于任意形状的网格，可以使用下面的方法：kedgesiiiiiiiikkyxcxVyUCFLt122（29）六、边界条件1固面边界条件对与无粘流动，固面边界条件无穿透条件，设其法向的速度通量为零，即0iZ。由于压强项的影响，x－向和y－向的动量通量并不为零。固面的压强近似的取为其相邻的单元的单元中心压强。广泛的数值研究证明，如果贴近壁面的单元足够小，并且人工耗散项运用正确，用这种方法取得的压强对结果的精度不会产生太大的影响。2远场边界条件本文提到的远场，实际是人为的有界边界，对于流场中的扰动会传到很远的地方。因而对于远场边界条件，情况比较复杂，它不能直接给定具体的流场值，需要与流场内的值来共同确定远边场的流场值。如果边界取得过小，则通常采用环量修正。一般情况下，我们采用无反射边界条件。为了保证扰动波不会反射回流场，应用A.Jamson提出的远场边界法向一维特征分析方法，来建立无反射的远场边界条件。一维均熵流动的Euler方程可写成：0xtUU(30)02xxtaUUU(31)这里，动量方程中除去了压强项，将上式写成矩阵形式为：-6-0xwAtw(32)这里：Uw，UaUA2(33)A的特征值为：aU1，aU2因此，上式的两族特征线为：aUdtdx(34)aUdtdx(35)(34)为C特征线，(35)为C特征线。沿C，C给出了众所周知的Riemann不变量R，R：12aqRn(36)12aqRn(37)这里不变量R，沿入流特征线C是常数，可以用来流条件计算得到；R沿出流特征线是常数，可以用流场内部向外插值计算：12aqRn(38)12eneeaqR(39)上式中下标“”表示来流值，下标“e”表示以计算域内部参数外插获得的值。通过Riemann不变量R，R的加减，可获得远场的法向速度nq和音速a：RRqen21(40)RRae41(41)根据Riemann不变量，按边界附近信息传播的性质把远场边界条件分成以下四种情况：A、亚音速入流条件0nU，它有三条入流特征线，需规定三个条件：1212aqaqnn(42)1212aqaqnn(43))(ttqq(44)ss(45)其中下标t表示切向，n代表方向，代表自由流，e代表从流场内到边界的外插值，上式的右端皆为已知，可解出边界上的sqcqtn,,,值，再由s和c求出p和。-7-B、亚音速出流条件0nU，它有一条入流特征线，需规定一个条件1212aqaqnn(46)1212eennaqaq(47)ettqq)((48)ess(49)C、超音速出流条件0nU，无入流特征线，不需在边界上规定边界条件lWW（50）其中的W表示边上的守恒变量，lW表示与此边相邻元素的守恒变量值。D、超音速入流条件0nU，它有四条入流特征线，需规定全部四个条件WW（51）其中的W表示边上的守恒变量，W表示来流值。-8-七、计算结果本文计算了三个算例，一个是攻角0，马赫数0.80的情况，二是攻角25.1，马赫数0.80的情况，三是攻角2.5°，马赫数1.5的情况XY-10-50510-10-50510Frame00128Jun1987XY00.51-0.500.5Frame00128Jun2007翼型网格示意图1．80.00Ma，＝表面压强系数分布等压线等马赫线压力云图-9-马赫数云图升力系数CL=-3.27408E-005阻力系数CD=1.14998670E-0022.80.025.1Ma，上下表面压强系数分布图等压线图等马赫图-10-压力云图马赫数云图升力系数为CL＝2.741130755E-001阻力系数为CD=2.15905376E-0023.5.15.2Ma，上下表面压强系数分布压力云图马赫数云图升力系数为CL=0.1352582564阻力系数为CD=0.1037842242-11-八、结论与展望本文主要介绍了求解二维非结构网格欧拉方程定常解的问题。采用的是有限体积空间离散方法，单元形状可以是任意多边形。使用显式多步推进过程，求解非定常方程得到定常解。一些标准的加速技术可以加快收敛速度。Jameson格式是学习CFD编程的入门格式，其精度和收敛速度都不高，但通过实际编写程序可以学到很多方法和技巧。因此，学习Jameson格式为以后编写其他格式如Roe，AUSM等等以及求解三维方程和求解N-S方程都打下了良好的基础。本次编程采用的是C语言，并且存在着一个比较大的问题，CFL数不能取到大于1.2的数值，在超音速计算时候CFL数能取的大值变的更小，估计是时间步长上的问题但是目前还没有找到问题所在。忘老师指点。参考文献1、L.STOLCIS*andL.J.JOHNSTON*,SolutionoftheEulerequationsonunstructuredgridsfortwo-dimensionalcompressibleflow,Aeronautics/AerospaceDepartmentVonKarmanInstituteforFluidDynamicsRhode-St-Genese,Belgium.2、J.Blazek,ComputationalFluidDynamicsPrinciplesandApplications,ELSEVIER.