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

有限体积法求解二维可压缩有限体积法求解二维可压缩 Euler 方程方程 计算流体力学课程大作业计算流体力学课程大作业 老师 老师 夏健 刘学强 学生 学生 徐锡虎 学号 学号 SQ09018013018 日期 日期 2010 年 2 月 5 日 1 目目 录录 一 内容摘要 2 二 流动控制方程 2 三 有限体积法的空间离散 2 四 人工耗散 3 五 时间离散 4 六 边界条件 5 七 计算结果 8 八 结论与展望 11 参考文献 11 2 一 内容摘要一 内容摘要 本文通过运用 JAMESON 有限体积法求解了二维定常和非定常可压缩 Euler 方程 程 序实现语言为 C 其中 使用的网格是三角形非结构网格 在时间推进上使用的是四步 龙 库塔推进格式 推进的时间步长取的是当地的时间步长 为了消除迭代误差 round off 等误差 本文采用了添加人工耗散项的办法 另外 本文计算了 NACA0012 翼型在跨 音速下不同迎角的情况 并与 fluent 软件的计算结果进行了比较 来验证程序的准确性 二 流动控制方程二 流动控制方程 守恒形式的 Euler 方程 1 0 GdxFdywd t S 其中 x 和 y 代表笛卡儿坐标系 W 是守恒变量 2 E V U W F G 表示通量 3 UH UV PU U F 2 VH PV UV V G 2 P H 和 E 表示密度 压强 单元总焓和单元总能量 U V 表示笛卡儿坐标系 下的速度矢量 这些量由理想气体的单位体积的总能量和总焓相互联系 4 2 1 22 VUPE 5 PEH 三 有限体积法的空间离散三 有限体积法的空间离散 计算域被划分为互不重叠的单元 在每个单元运用守恒形式的 Euler 方程 由于每个单 元相对于时间都是不变的 所以等式 1 可以写成 8 d GdxFdy t W S 其中和 S 是单元的体积和边界 W 是单元的平均值 在对上述方程进行时间离散前 先对空间进行离散 则方程 6 可以写为 9 k k kQ dt dW 3 其中表示第 k 个单元的体积 是第 k 个单元的守恒变量 表示第 k 个单元 k k W k Q 的通量 方程 7 的右边项可以写成 10 kedges i ik xGyFQ 1 其中 abiabi yyyxxx 11 8 式中的求和是对第 k 个单元的所有边进行的 守恒参数的量是单元中心值 在 求通量时 第 I 条边的守恒参数值是用左右单元的平均来表示的 12 2 W pki WW 引入变量 13 iiiii xVyUZ 则第 k 单元的 Euler 方程可以写为 14 kedges i i k k HZ xPVZ yPUZ Z dt dW 1 1 在本文中 采用的是 JAMENSON 有限体积法 为了减少存储的相关信息的量 其存储的 方式选择的是按边存储的方法 在存储的每条边的信息中 包含了这条边的边号 左右单 元号和边的端点 在计算通量时采用按边循环的方式 do I 1 nedge k 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 flux sum p sum p flux end do 这里给出的是 FOTRAN 语言的形式 我编写采用的是 C 具体表现在上交的程序中 在计算时间步长 人工耗散项等也可用象这样按边循环 从此处我们可以看出求解时 与单元的形状无关 四 人工耗散四 人工耗散 人工粘性模型对方法的成功应用起着关键作用 人工粘性抑制解在激波附近的振荡 又阻尼解在光滑区域内的高阶误差 对解的线性稳定和收敛于定态是很重要的 本文在方 程 14 的右端加入了人工耗散项 如对于单元 k 其表达式可以表示为 15 k kk kDQ dt dW 4 在有限体积法中 耗散项的公式可以表示为 16 kedges i i kedges i ik ddD 1 4 1 2 其中 17 i kpiii i kpiii WWd WWd 22 4 4 2 2 其中的 I 表示单元 k 和 p 的公共边 定义为 2 18 kedges i kjk WWW 1 2 上面的 j 表示与 k 相邻的单元 19 kedges i ikp kedges i ikp k PP PP 1 1 20 0max max 2 4 4 2 2 ii ikpi k k 其中的量的范围是 4 2 kk0 121 3212561 2 4 kk 在计算时发现上面方法得到的人工耗散项并不太适合 其在光滑区域耗散项太大 而在大剃度区域又显得太小 为了弥补上面的不足 作下面的修改 21 kp kp i PP PP 自适应系数为 22 0max 2 4 4 2 2 ii ii k k 尺度系数为 23 22 yxcxVyU i 其中的 U V 表示边上的值 C 表示当地声速 五 时间离散五 时间离散 方程最后的稳定解是通过时间上的迭代得到的 可以写为 24 k k R dt dW 右边项的表达式为 5 25 kkkk DQR 为了加速收敛 时间迭代使用的是 4 步龙 库塔推进格式 格式如下 26 4 1 1 0 0 4 1 WW tomfortRWW WW n m m m n 其中的 n 表示的是当前的时间步 n 1 表示的是新的时间步 27 0 DQR mm 28 1 21 31 41 4321 为了减小计算时间 人工耗散项的计算只在第一步进行 在下面几步的迭代中保持 不变 运用上面的方法计算 可以发现 CFL 数可以取到 本文中使用的是 2 0 22 使用显示格式迭代的主要缺点是由于稳定区域的限制 所以不能使用过大的时间步长 可以用近似的方法估算时间步长 对于任意形状的网格 可以使用下面的方法 29 kedges i iiiiiii k k yxcxVyU CFL t 1 22 六 边界条件六 边界条件 1 1 固面边界条件固面边界条件 对与无粘流动 固面边界条件无穿透条件 设其法向的速度通量为零 即 由0 i Z 于压强项的影响 x 向和 y 向的动量通量并不为零 固面的压强近似的取为其相邻的单 元的单元中心压强 广泛的数值研究证明 如果贴近壁面的单元足够小 并且人工耗散项 运用正确 用这种方法取得的压强对结果的精度不会产生太大的影响 2 2 远场边界条件远场边界条件 本文提到的远场 实际是人为的有界边界 对于流场中的扰动会传到很远的地方 因 而对于远场边界条件 情况比较复杂 它不能直接给定具体的流场值 需要与流场内的值 来共同确定远边场的流场值 如果边界取得过小 则通常采用环量修正 一般情况下 我 们采用无反射边界条件 为了保证扰动波不会反射回流场 应用 A Jamson 提出的远场边界法向一维特征分析 方法 来建立无反射的远场边界条件 一维均熵流动的 Euler 方程可写成 30 0 xt UU 31 0 2 xxt a UUU 这里 动量方程中除去了压强项 将上式写成矩阵形式为 6 32 0 x w A t w 这里 33 U w U a U A 2 的特征值为 AaU 1 aU 2 因此 上式的两族特征线为 34 aU dt dx 35 aU dt dx 34 为特征线 35 为特征线 C C 沿 给出了众所周知的 Riemann 不变量 C C R R 36 1 2 a qR n 37 1 2 a qR n 这里不变量 沿入流特征线是常数 可以用来流条件计算得到 沿出流特 R C R 征线是常数 可以用流场内部向外插值计算 38 1 2 a qR n 39 1 2 e nee a qR 上式中下标 表示来流值 下标 表示以计算域内部参数外插获得的值 e 通过 Riemann 不变量 的加减 可获得远场的法向速度和音速 R R n qa 40 RRq en 2 1 41 RRa e 4 1 根据 Riemann 不变量 按边界附近信息传播的性质把远场边界条件分成以下四种情 况 A 亚音速入流条件亚音速入流条件 它有三条入流特征线 需规定三个条件 它有三条入流特征线 需规定三个条件 0 nU 42 1 2 1 2 a q a q nn 43 1 2 1 2 a q a q nn 44 tt qq 45 ss 7 其中下标 t 表示切向 n 代表方向 代表自由流 代表从流场内到边界的外插 e 值 上式的右端皆为已知 可解出边界上的值 再由 和 求出和 sqcq tn scp B 亚音速出流条件亚音速出流条件 它有一条入流特征线 需规定一个条件 它有一条入流特征线 需规定一个条件 0 nU 46 1 2 1 2 a q a q nn 47 1 2 1 2 e enn a q a q 48 ett qq 49 e ss C 超音速出流条件超音速出流条件 无入流特征线 不需在边界上规定边界条件 无入流特征线 不需在边界上规定边界条件 0 nU l WW 50 其中的 W 表示边上的守恒变量 表示与此边相邻元素的守恒变量值 l W D 超音速入流条件超音速入流条件 它有四条入流特征线 需规定全部四个条件 它有四条入流特征线 需规定全部四个条件 0 nU WW 51 其中的 W 表示边上的守恒变量 表示来流值 W 8 七 计算结果七 计算结果 本文计算了三个算例 一个是攻角 马赫数 0 80 的情况 二是攻角 马赫 0 25 1 数 0 80 的情况 三是攻角 2 5 马赫数 1 5 的情况 X Y 10 50510 10 5 0 5 10 Frame 001 28 Jun 1987 X Y 00 51 0 5 0 0 5 Frame 001 28 Jun 2007 翼型网格示意图 1 1 80 0 0 Ma 表面压强系数分布 等压线 9 等马赫线 压力云图 马赫数云图 升力系数 CL 3 27408E 005 阻力系数 CD 1 14998670E 002 2 80 0 25 1 Ma 上下表面压强系数分布图上下表面压强系数分布图 10 等压线图等压线图 等马赫图等马赫图 压力云图压力云图 马赫数云图马赫数云图 升力系数为 CL 2 741130755E 001 阻力系数为 CD 2 15905376E 002 3 5 15 2 Ma 上下表面压强系数分布上下表面压强系数分布 11 压力云图 马赫数云图 升力系数为 CL 0 1352582564 阻力系数为 CD 0 1037842242 八 结论与展望八 结论与展望 本文主要介绍了求解二维非结构网格欧拉方程定常解的问题 采用的是有 限体积空间离散方法 单元形状可以是任意多边形 使用显式多步推进过程 求解非定常方程得到定常解 一些标准的加速技术可以加快收敛速度 Jameson 格式是学习 CFD 编程的入门格式 其精度和收敛速度都不高 但 通过实际编写程序可以学到很多方法和技巧 因此 学习 Jameson 格式为以后 编写其他格式如 Roe AUSM 等等以及求解三维方程和求解 N S 方程都打下了良 好的基础 本次编程采用的是 C 语言 并且存在着一个比较大的问题 CFL 数不能取 到大于 1 2 的数值 在超音速计算时候 CFL 数能取的大值变的更小 估计是时 间步长上的问题但是目前还没有找到问题所在 忘老师指点 参考文献参考文献 1 L STOLCIS and L J JOHNSTON Solution of th