高精度差分格式及湍流数值模拟(三)
高精度差分格式在机翼分离流计算中的应用
ቤተ መጻሕፍቲ ባይዱ
数 修 正 ， 如 何 修 正是 一 个 难 题 。用 计 算 手 段 研 究 雷 诺 数 变 化 对 机 翼 气 流 分 离 的 影 响 将 但
能为此提供 一种 方法 。在 文献 ［ ］ ［ ］ １和 ２ 中笔者 给出 了一种机翼边 界层及 其对位流干扰 计
式
为
， 详 见 文 献 ［ ］方 程 （ ） 方程 （ ） 以 以 ， ， ｃ １。 ２至 ５可 Ｅ和 面 作 为 变 量 求 解 ， 界 条 件 边
７ ＝ ０， ７
＝ 历 ＝
＝ ０；
７ ＝ ７
，
＝
＝
＝ Ｗｅ
＝
（） ６
在 方 程 （ ） （） ， ２ 和 ３ 中 ｂ为 无 因 次 的湍 流粘 度 系 数 。 系 数 ｍ 和 ｍ 表 达式 为
一
个 大 展 弦 比后 掠 翼 给 出 了攻 角 ｌ。 不 同 雷 诺 数 下 的 机 翼 上 表 面 气 流 分 离 线 的 计 算 结 ４时
果。
１ 计 算 方 法
机 翼 的 三维 边 界 层 方 程 是 在 机翼 表 面 的 非 正 交 曲线 坐 标 系 （ ， 。 中 建立 的 。 它 们 Ｙ， ）
关 键 词 ：分 离 流 ； 界 层 计 算 ； 分 格 式 ； 性 一无 粘 干 扰 ； 翼 边 差 粘 机
中 图 分 类 号 ：Ｖ ｌ ．１ ２ １４
文 献 标 识 码 ：Ａ
０ 引 言
雷 诺 数 是边 界 层 计 算 中较 敏 感 的 参 数 ， 其 在 强 逆 压 梯 度 下 ， 诺 数 大 小 对 气 流 分 离 尤 雷 区形 成 与 扩 展 过 程 有 很 大 影 响 。 机 翼 在 较 低 雷 诺 数 下 面 的 分 离 形 态 同实 际 飞行 雷 诺 数 下
湍流数值模拟及其在工程热力学中的应用
湍流数值模拟及其在工程热力学中的应用湍流是自然界和工程中广泛存在的一种流动状态，其具有不规则、不稳定、非线性等特点。
因此，湍流研究成为了流体力学中的一个重要分支。
湍流数值模拟（Large Eddy Simulation）是目前研究湍流问题的重要手段之一，广泛应用于工程热力学中。
湍流数值模拟技术的发展历程湍流数值模拟技术起源于20世纪50年代，当时主要应用于理论模拟。
20世纪80年代后，随着计算机技术的发展，数值模拟技术应用于实际工程中，并得到广泛应用。
近年来，由于计算机性能的不断提高和算法的不断改进，湍流数值模拟技术越来越成熟，其应用范围也更加广泛。
湍流数值模拟技术的基本原理湍流数值模拟技术的基本原理是将流场分为宏观湍流和微观湍流两部分，并通过不同方法对二者进行模拟。
具体而言，宏观湍流采用平均场方程进行模拟，微观湍流则通过小尺度涡结构之间的相互作用进行模拟。
在湍流数值模拟过程中，关键是要准确地描述湍流的能量转移和钝化机制，以便合理地模拟湍流特性。
目前，湍流数值模拟技术主要有两种方法：直接数值模拟和大涡模拟。
直接数值模拟（Direct Numerical Simulation，DNS）是最为精确的湍流数值模拟方法，它直接求解完整的Navier-Stokes方程，但计算量也是最大的。
而在工程应用中，一般采用次网格模型，采用模型对小尺度湍流进行近似处理，减少计算量。
其中，大涡模拟（Large Eddy Simulation，LES）是一种很有代表性的方法，它将外部湍流场分解为大尺度湍流和小尺度湍流两部分，对大尺度湍流进行直接数值模拟，对小尺度湍流采用模型进行处理。
湍流数值模拟在工程热力学中的应用湍流数值模拟技术在工程热力学中有着广泛的应用。
具体而言，湍流数值模拟可以用来模拟涡流管道的流动、火焰、燃烧室和喷气发动机等复杂流场问题。
下面，我们将从两个方面来介绍湍流数值模拟在工程热力学中的应用：（1）流体力学问题湍流数值模拟技术在流体力学问题中得到了广泛应用，例如现代汽车设计中对车身和车厢空气动力学的研究，对于气动设计、噪声控制和气密性等方面的分析有很大的帮助。
高精度频率域弹性波方程有限差分方法及波场模拟
地 球 物 理 学 报
CHINESE JOURNAL OF GEOPHYSICS
Vol. 49 , No. 2 Mar. , 2006
殷 文 ,印兴耀 ,吴国忱等 . 高精度频率域弹性波方程有限差分方法及波场模拟 . 地球物理学报 ,2006 , 49 (2) :561～568
时间域计算方法按时间片递推计算 , 每一个时 间片的舍入误差会累积到下一片中 , 如果计算的时
u 9 λ 9 u 9v 9w ρ9 2 = + + 9 x 9x 9y 9z 9t
2
+2 +2
9 μ9u 9 μ 9 u 9v + + 9x 9x 9y 9y 9x 9 μ9w 9 μ 9 u 9w + + 9y 9y 9x 9y 9x
562
地 球 物 理 学 报 ( Chinese J . Geophys. )
49 卷
间片较多 ,最终可能导致累积误差太大而波场值太
1 引 言
波动方程数值模拟实质上是根据地下介质条件 求解地震波波动方程 , 模拟的地震波场包含地震波 传播的各种信息 . 通过地震波场数值模拟 , 人们可 以探索地震波传播机理 、 验证修改解释方案 、 进行反 问题研究等 . 这对于地震采集施工设计 , 地震资料 的处理 、 解释以及地震勘探后的岩性评价都具有重 要意义 . 随着三维地震勘探的广泛应用及复杂介质 精细成像的普遍开展 , 研究高效的正演模拟数值方 法非常必要 . 波动方程正演模拟的数值方法有许多种 , 但大 都存在处理计算精度或边界条件问题的困难 . 虽然 有限元方法不受边界几何条件的限制 , 处理自由边 界条件较为简单 , 但有限元方法要求的大计算量和 高存 储 量 在 现 有 的 计 算 水 平 下 仍 然 是 难 以 接 受 [1 ] 的 . 有限差分法是进行波动方程数值模拟的一种 有效方法 ,具有计算速度快 、 占用内存小等优点 , 该 方法对近远场及复杂边界都有广泛的适用性 , 能够 准确地模拟波在各种介质及复杂结构地层中的传播 规律 . 但常规有限差分法通常具有较强的数值频 散现象 ,存在边界反射 , 且计算效率较低 . Claton 和 [3 ] Engquist 提出一种基于声波或弹性波方程的吸收 边界条件 ,该方法仅适用于入射角较小的情况 ,并且 当泊松比大于 210 时 , 会出现不稳定问题 ; Cerjan et
湍流的几种数值模拟方法
LES特点
抓大不放小 非常有利，有力的工具 是最近,可预见未来流体 力学研究和应用的热点 近来又出现了VLES, DES等在LES上发展而 来的工具
Will RANS survive LES? Hanjalic自问自答
会。Journal of Fluids Engineering -V127, 5, pp. 831-839 (Will RANS
Prandtl(1925)混合长度模型
也被称作零方程模型 还在被广泛应用 廉价，易收敛 基本在流场比较简单，或者对计算结果 精度要求不高或者流场形状比较复杂的 行业中，比如暖通空调，流体机械等。
Prandtl混合长度模型 缺点
最明显的缺点是:当速度梯度 为零的 时候， 消失， 这与事实不符
Launder and Li(1994), Craft and Launder (1995)
目前有很多学者在继续此方面的工作
Brian E. Launder
本科Imperial College, London 硕博 MIT 实验流体力学 1964-1976 Imperial College 讲师
涡流粘度
Eddy viscosity or turbulent viscosity
二维流场分子粘性力
为描述雷诺应力，Boussinesq 1887 定义了与之相对应的
RANS模型的核心在于给出 的数 学表达式，要求精度高，适用范围广
涡流粘度,
Prandtl 1925 Prandtl 1945 Bradshaw 1968 Kolmogorov, 1942 Hanjalic 1970 Rotta 1951 Chou 1945 Davidov 1961
求解一维对流方程的高精度紧致差分格式___
应用数学MATHEMATICA APPLICATA2019,32(3):635-642求解一维对流方程的高精度紧致差分格式侯波,葛永斌(宁夏大学数学统计学院,宁夏银川750021)摘要：本文提出数值求解一维对流方程的一种两层隐式紧致差分格式,采用泰勒级数展开法以及对截断误差余项中的三阶导数进行修正的方法对时间和空间导数进行离散.格式的截断误差为O(τ4+τ2h2+h4),即该格式在时间和空间上均可达到四阶精度.利用von Neumann方法分析得到该格式是无条件稳定的.通过数值实验验证了本文格式的精确性和稳定性.关键词：对流方程;高精度;紧致格式;无条件稳定;有限差分法中图分类号：O241.82AMS(2000)主题分类：65M06;65M12文献标识码：A文章编号：1001-9847(2019)03-0635-081.引言对流方程在生物数学、能源开发、空气动力学等许多领域都具有十分广泛的应用,因此求解该类方程具有非常重要的理论价值和实际意义.然而,由于实际问题通常十分复杂,往往难以求得精确解,因此研究其精确稳定的数值解法是十分必要的.针对对流方程国内外很多学者提出了很多的数值方法.如张天德和孙传灼[1]针对一维对流方程采用待定系数法,得到了两层四点格式和四阶六点格式,并且是无条件稳定的,该方法适用于在点数确定的前提下,得到精度高的差分格式;于志玲和朱少红[2]针对一维问题建立了中间层为两个节点的三层显格式,其截断误差为O(τ2+h2);曾文平[3]针对一维对流方程推导出了一种两层半显式格式，其截断误差为O(τ2+h2),该格式是无条件稳定的.姚朝辉等人[5]将二阶的迎风格式和中心差分格式进行加权得到了WSUC格式,该格式是无条件稳定的;但该格式时间方向和空间方向仅有二阶精度.汤寒松等人[6]通过立方插值拟质点方法(CIP方法),给出了一些保单调的CIP格式;Erdogan[9]针对一维的对流方程推导出了一种指数拟合的差分格式,其截断误差为O(τ2+h2);Bourchtein[10]构造了对流方程的三层五点中心型蛙跳格式,该格式的截断误差为O(τ4+h4);即该格式时间和空间均具有四阶精度,但是该格式是三层的,空间方向需要五个点,并且是条件稳定的;Kim[11]构造了多层无耗散的迎风蛙跳格式,即时间和空间分别具有二阶、四阶、六阶精度,但格式为三层甚至是四层的,并且六阶格式空间方向最多需要五个点,给靠近边界的内点的计算带来困难.综上所述,文献中已经有的数值计算方法大多为低阶精度的,而高精度方法涉及多个时间层,需要一个或多个时间启动步,或者空间方向的网格节点多于三个,这都给计算造成困难或不便.为此本文将构造一种紧致格式,这里紧致格式的定义为对时间导数项的离散采用不超过∗收稿日期：2018-08-10基金项目：国家自然科学基金(11772165,11361045),宁夏自然科学基金重点项目(2018AAC02003),宁夏自治区重点研发项目(2018BEE03007)作者简介：侯波,男,汉族,河南人,研究方向：偏微分方程数值解法.通讯作者：葛永斌.636应用数学2019三个时间层,而对空间导数项的离散采用不超过三个网格点,时间和空间即可以达到高阶精度（三阶及三阶以上）的格式.本文拟构造的格式时间方向仅用到两个时间层上的函数值,在每个时间层上仅涉及到三个空间网格点,格式时间和空间具有整体的四阶精度.该格式的优点是无须启动步的计算,并且在对靠近边界点的计算时,不会用到计算域以外的网格节点.此外该格式为无条件稳定的,可以采用比较大的时间步长进行计算.最后通过数值实验验证本文格式的精确性和稳定性.2.差分格式的建立考虑如下一维对流方程：∂u ∂t +a∂u∂x=f,b≤x≤c,t≥0,(2.1)给定初始条件为:u(x,0)=φ(x),b≤x≤c,(2.2)给定周期性边界条件为:u(b,t)=u(c,t),t≥0,(2.3)其中,u(x,t)为未知函数,f为非齐次项,a为对流项系数,φ(x)为已知函数.将求解区域[b,c]等距剖分为N个子区间：b=x0,x1,···,x N−1,x N=c,并且定义h=c−bN,时间也采用等距剖分,步长用τ表示.在本文中,我们利用u ni ,u n+1i,u n+12i分别表示u在(x i,t n),(x i,t n+1)和(x i,t n+12)点处的函数值.假设方程(2.1)在点(x i,t n+12)成立,简写表示为:(∂u ∂t )n+12i+a(∂u∂x)n+12i=f n+12i.(2.4)将u n+1i 和u ni在点(x i,t n+12)处做泰勒级数展开,可得:u n+1i=u n+12i+τ2(∂u∂t)n+12i+(τ2)22!(∂2u∂t2)n+12i+(τ2)33!(∂3u∂t3)n+12i+O(τ4),(2.5)u ni=u n+12i−τ2(∂u∂t)n+12i+(τ2)22!(∂2u∂t2)n+12i−(τ2)33!(∂3u∂t3)n+12i+O(τ4).(2.6)(2.5)-(2.6)可得:(∂u∂t)n+12i=δt u n+12i−τ224(∂3u∂t3)n+12i+O(τ4),(2.7)其中,δt u n+12i =u n+1i−u n iτ.同理可得:(∂u∂x)n+12i=δx u n+12i−h26(∂3u∂x3)n+12i+O(h4),(2.8)其中,δx u n+12i =un+12i+1−u n+12i−12h.将(2.7)和(2.8)代入(2.4)整理可得:δt u n+12i −τ224(∂3u∂t3)n+12i+aδx u n+12i−ah26(∂3u∂x3)n+12i=f n+12i+O(τ4+h4).(2.9)为了使该格式在时间方向和空间方向上均达到四阶精度,须对(2.9)式中的∂3u∂t3和∂3u∂x3项进行二阶的离散,同时为了保证本文格式的紧致性,即空间方向不超过三个网格点,我们对(2.1)式进行如下变形:∂u ∂t =−a∂u∂x+f,∂2u∂t2=a2∂2u∂x2−a∂f∂x+∂f∂t,第3期侯波等：求解一维对流方程的高精度紧致差分格式637∂3u ∂t 3=a 2∂3u ∂x 2∂t −a ∂2f ∂x∂t +∂2f ∂t 2,∂3u ∂x 3=−1a ∂3u ∂x 2∂t +1a ∂2f ∂x 2.(2.10)将上述∂3u ∂t 3和∂3u∂x 3的表达式(2.10)代入(2.9)并整理可得：δt u n +12i+aδx u n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t)n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+h 4).(2.11)如果对上式中的δx u n +12i 项采用时间方向算术平均,即δx u n +12i =δx u n +1i+u n i 2,则会导致格式时间退化为二阶精度,为此利用(2.5)+(2.6)可得：u n +12i =12(u n +1i +u n i )−τ28(∂2u ∂t2)n +12i +O (τ4).(2.12)从而可得:δx u n +12i =12δx (u n +1i +u n i )−τ28δx (∂2u ∂t2)n +12i +O (τ4).(2.13)将(2.13)代入(2.11)得:δt u n +12i +a 2δx (u n +1i +u n i )−aτ28δx (∂2u ∂t 2)n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t )n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+h 4).(2.14)由于δx (∂2u ∂t 2)n +12i =(∂3u ∂x∂t 2)n +12i+O (h 2),所以可得:δt u n +12i +a 2δx (u n +1i +u n i )−aτ28(∂3u ∂x∂t 2)n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t)n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+τ2h 2+h 4).又因为∂3u ∂x∂t 2=−a ∂3u∂x 2∂t +∂2f ∂x∂t ,所以有:δt u n +12i +a 2δx (u n +1i +u n i )−aτ28(−a ∂3u ∂x 2∂t +∂2f ∂x∂t )n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t )n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+τ2h 2+h 4),即,δt u n +12i +a 2δx (u n +1i +u n i )+(a 2τ212+h 26)(∂3u ∂x 2∂t )n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i −aτ212(∂2f ∂x∂t )n +12i =f n +12i +O (τ4+τ2h 2+h 4).由于(∂3u ∂x 2∂t )n +12i=δ2x (∂u ∂t )n +12i +O (h 2),所以有:u n +1i −u n i τ+a 4h(u n +1i +1−u n +1i −1+u ni +1−u n i −1)+(h 26+a 2τ212)δ2x u n +1i −u n i τ−τ224(f n +1i −2f n +12i +f n −1i (τ2)2)−h 212[(∂2f ∂x 2)n +1i +(∂2f ∂x 2)n −1i ]−aτ12[(∂f ∂x )n +1i −(∂f ∂x)n −1i ]=f n +12i +O (τ4+τ2h 2+h 4),其中,δ2xu i =u i +1−2u i +u i −1h 2,舍去O (τ4+τ2h 2+h 4),等式两边同时乘以τ,并令λ=τ/h ,整理可得:u n +1i +aλ4(u n +1i +1−u n +1i −1)+(16+a 2λ212)(u n +1i +1−2u n +1i +u n +1i −1)638应用数学2019=u n i−aλ4(u n i +1−u n i −1)+(16+a 2λ212)(u n i +1−2u n i +u ni −1)+τ6(f n +1i −2f n +12i +f n i )+τ12(f n +1i +1−2f n +1i +f n +1i −1+f n i +1−2f n i +f n i −1)+aτλ24(f n +1i +1−f n +1i −1−f n i +1+f ni −1)+τf n +12i,即,(23−a 2λ26)u n +1i +(16+aλ4+a 2λ212)u n +1i +1+(16−aλ4+a 2λ212)u n +1i −1=(23−a 2λ26)u n i +(16−aλ4+a 2λ212)u n i +1+(16+aλ4+a 2λ212)u n i −1+(τ12+aλτ24)f n +1i +1(τ12−aλτ24)f n +1i −1+(τ12−aλτ24)f n i +1+(τ12+aλτ24)f n i −1+2τ3f n +12i .(2.15)由推导过程可知,该格式的截断误差为O (τ4+τ2h 2+h 4),即格式(2.15)在时间和空间上均可达到四阶精度.我们注意到,格式为两层格式,并且格式每层仅用到三个网格点,形成的代数方程组系数矩阵为循环三对角矩阵,可采用追赶法进行求解[8],同时由于要求未知时间层上(第n +1层)中间点的系数不能等于0,即23−a 2λ26=0,因此aλ=2.3.稳定性分析下面采用von Neumann 方法分析本文所推导的差分格式(2.15)的稳定性.对于(2.15)式,舍掉非齐次项f ,即假设f 项精确成立,令u n i =ηn e Iσx i,其中,η为振幅,σ为波数,I =√−1为虚数单位,有(23−a 2λ26)ηn +1e Iσx i +(16+aλ4+a 2λ212)ηn +1e Iσx i +1+(16−aλ4+a 2λ212)ηn +1e Iσx i −1=(23−a 2λ26)ηn e Iσx i +(16−aλ4+a 2λ212)ηn e Iσx i +1+(16+aλ4+a 2λ212)ηn e Iσx i −1.(3.1)两边同时约掉e Iσx i ,并整理可得：(23−a 2λ26)ηn +1+(16+a 2λ212)ηn +1(e Iσh +e −Iσh )+aλ4ηn +1(e Iσh −e −Iσh )=(23−a 2λ26)ηn+(16+a 2λ212)ηn (e Iσh +e −Iσh )−aλ4ηn +1(e Iσh −e −Iσh ).(3.2)利用Euler 公式,即e Iσh =cos σh +I sin σh,e −Iσh =cos σh −I sin σh ,可得:(23−a 2λ26)ηn +1+[(13+a 2λ26)cos σh ]ηn +1+(aλI 2sin σh )ηn +1=(23−a 2λ26)ηn +[(13+a 2λ26)cos σh ]ηn −(aλI 2sin σh )ηn .(3.3)对上式进行化简整理有[(23−a 2λ26)+(13+a 2λ26)cos σh +aλI sin σh 2]ηn +1=[(23−a 2λ26)+(13+a 2λ26)cos σh −aλI sin σh 2]ηn .(3.4)从而可得格式(2.15)的误差放大因子为:G =ηn +1ηn =(23−a 2λ26)+(13+a 2λ26)cos σh −aλI sin σh2(23−a 2λ26)+(13+a 2λ26)cos σh +aλI sin σh2.(3.5)由von Numann 稳定性定理可知当|G |≤1时,格式是稳定的,由(3.5)可得|G |=1,因此,格式(2.15)是无条件稳定的.4.数值实验第3期侯波等：求解一维对流方程的高精度紧致差分格式639为了验证本文格式(2.15)的精确性和稳定性,现考虑以下三个具有精确解的初边值问题.分别采用Crank-Nicolson(C-N)格式,文[7]中格式和本文格式(2.15)进行计算;其中,最大绝对误差及收敛阶的定义为:L∞=maxi |u n i−u(x i,t n)|,Rate=log[L∞(h1)/L∞(h2)]log(h1/h2)L∞(h1)和L∞(h2)为空间网格步长分别为h1和h2时的最大绝对误差.问题1[7]:∂u ∂t +∂u∂x=0,0≤x≤2,t>0,u(x,0)=sin(πx),0≤x≤2,u(0,t)=u(2,t),t>0,该问题的精确解为:u(x,t)=sin[π(x−t)].表1问题1当λ=τ/h=0.5,t=1时刻的最大绝对误差及收敛阶h推进步数(n)C-N格式文[7]本文格式L∞误差Rate L∞误差Rate L∞误差Rate 1/510 2.217(-1) 4.865(-2) 1.993(-3)1/1020 5.752(-2) 1.95 1.263(-2) 1.95 1.208(-4) 4.041/2040 1.450(-2) 1.99 3.199(-3) 1.987.490(-6) 4.011/4080 3.631(-3) 2.008.038(-4) 1.99 4.672(-7) 4.001/801609.082(-4) 2.00 2.014(-4) 2.00 2.919(-8) 4.001/160320 2.271(-4) 2.00 5.041(-5) 2.00 1.824(-9) 4.00表2问题1当τ=λh,t=2时刻的最大绝对误差hτλC-N格式文献[7]本文格式1/160.050000000.8 5.290(-2) 1.292(-2) 1.574(-5) 0.10000000 1.69.013(-2) 5.095(-2) 3.198(-3) 0.20000000 3.2 2.307(-1) 1.941(-1) 6.055(-2) 0.40000000 6.4 6.874(-1) 6.597(-1) 1.746(-2)1/320.025000000.8 1.330(-2) 3.230(-3)9.814(-7) 0.20000000 6.4 2.041(-1) 1.950(-1) 1.575(-3) 0.4000000012.8 6.668(-1) 6.601(-1) 1.916(-2)图1问题1当N=32,τ=0.03125,t=0.2时刻的数值解与精确解640应用数学2019表1给出了针对问题1三种格式在不同空间步长h下,当λ=τ/h=0.5,t=1时的最大绝对误差和收敛阶.我们发现C-N格式在时间和空间上都为二阶精度,由于文[7]格式时间具有二阶精度,空间具有四阶精度,因此当取τ=O(h)时,格式空间仅有二阶精度,而本文格式时间和空间均为四阶精度.图1给出N=32,τ=0.03125,t=0.2数值解与精确解对比图,可以看出数值解与精确解吻合的很好.表2给出了当h=1/16和h=1/32时,τ=λh,t=2时刻对问题1采用三种格式计算的最大绝对误差.可以看出网格比λ最大取到12.8,计算仍然是稳定的,因此本文格式是无条件稳定的.并且本文格式在所有参数下,其计算结果比C-N格式和文[7]格式计算结果更加精确.问题2[7]:∂u ∂t +∂u∂x=0,0≤x≤2,t>0,u(x,0)=e cos(πx),0≤x≤2,u(0,t)=u(2,t),t>0,该问题的精确解为:u(x,t)=e cos[π(x−t)].表3问题2当λ=τ/h=0.5,t=1时刻的最大绝对误差及收敛阶h推进步数(n)C-N格式文[7]本文格式L∞误差Rate L∞误差Rate L∞误差Rate 1/510 6.754(-1) 1.428(-1) 5.567(-2)1/1020 2.310(-1) 1.55 3.099(-2) 2.20 3.041(-3) 4.191/2040 6.027(-2) 1.94 6.825(-3) 2.18 1.904(-4) 4.001/4080 1.492(-2) 2.01 1.658(-3) 2.04 1.165(-5) 4.031/80160 3.705(-3) 2.01 4.115(-4) 2.017.252(-7) 4.011/1603209.250(-4) 2.00 1.028(-4) 2.00 4.527(-8) 4.00表4问题2当τ=λh,t=2时刻的最大绝对误差hτλC-N格式文[7]本文格式1/160.050000000.8 2.171(-1) 5.372(-2) 3.897(-4) 0.10000000 1.6 3.450(-1) 2.056(-1)7.795(-3) 0.20000000 3.2 6.810(-1) 6.111(-1) 3.416(-1) 0.40000000 6.4 1.220 1.198 2.017(-1)1/320.025000000.8 5.575(-2) 1.325(-2) 2.449(-5) 0.20000000 6.4 6.302(-1) 6.109(-1) 2.350(-2) 0.4000000012.8 1.204 1.199 2.201(-1)表3和表4给出了针对问题2利用本文格式和C-N格式以及文[7]格式的计算结果.表3考察了格式的精度,表4验证了格式的稳定性.可以看出本文格式在时间和空间上均可达到四阶精度,并且是无条件稳定的.问题3∂u ∂t +a∂u∂x=f,0≤x≤2,t>0,u(x,0)=cos(πx),0≤x≤2,u(0,t)=u(2,t),t>0,f=π(1−a)sin[π(x−t)],该问题的精确解为:u(x,t)=cos[π(x−t)].第3期侯波等：求解一维对流方程的高精度紧致差分格式641表5问题3当λ=τ/h=0.5,a=0.5,t=1时刻的最大绝对误差及收敛阶h推进步数(n)C-N格式本文格式L∞误差Rate L∞误差Rate1/510 1.124(-1) 4.244(-4)1/1020 3.520(-2) 1.67 2.744(-5) 3.951/20409.957(-3) 1.82 1.739(-6) 3.981/4080 2.551(-3) 1.96 1.134(-7) 3.941/80160 6.413(-4) 1.99 1.351(-8) 3.07问题3为非齐次问题,由于文[7]的方程模型为齐次方程,不能计算非齐次问题,因此该问题我们采用本文格式和C-N进行计算和比较,表5给出了两种格式在不同空间步长h下,当t=1时的最大绝对误差和收敛阶.可以看出当λ=τ/h=0.5,a=0.5时,C-N格式在时间和空间上都为二阶精度,而本文格式时间和空间均为四阶精度.5.结论本文针对一维对流方程提出了一种两层隐式高精度紧致差分格式,时间和空间均采用泰勒级数展开法以及截断误差余项修正法进行处理,格式截断误差为O(τ4+τ2h2+h4),即该格式在时间和空间上均可达到四阶精度.并通过von Neumann方法分析得到该格式为无条件稳定的.最后通过三个数值算例验证了格式的精确性和稳定性.通过上述研究,我们可以得出如下结论:1.文献(如[10-11])中的高精度格式往往是时间多层格式,需要另外构造启动步的计算格式,如果采用低精度格式启动,必然会影响以后时间层的计算精度.而本文格式仅为两层格式,无须启动步的计算,时间即可达到四阶精度.2.文献(如[1,10-11])中的高精度格式空间方向上往往超过三个网格节点,导致靠近边界的内点计算困难,需要采用特殊处理,而本文格式仅用到三个网格节点,可以有效避免这一问题.3.尽管本文格式要求aλ=2,这是本文格式的一个缺陷,但是由于本文格式是无条件稳定的,从理论上讲可以采用任意网格比,因此可以很容易避开aλ=2的条件限制,使得这一缺陷并不太影响格式的使用.4.由于本文方法推导过程中涉及到∂2u∂t2,∂3u∂t3,∂3u∂x3的计算,需要用原方程进行多次求导并进行反复代入计算,在考虑对流项为变系数问题时,将涉及到a(x,t)关于x和t的二阶导数,由于我们考虑在时间半点处,即(x i,t n+12)处的函数值,即要用到(∂2a∂t2)n+12i,如果采用中心差分,则时间仅具有二阶精度,因此本文方法不适用于变系数问题.5.本文方法可直接推广到二维和三维问题中去,我们将另文报道.参考文献：[1]张天德,孙传灼.对流方程的差分格式[J].计算物理,1995,12(2):191-195.[2]于志玲,朱少红.关于对流方程一类三层显格式[J].南开大学学报(自然科学版),1998,31(3):27-30.[3]曾文平.解对流方程的加耗散项的差分格式[J].应用数学,2001,14(S1):154-158.[4]陆金甫,关治.偏微分方程数值解法[M].北京:北京大学出版社,1987.[5]姚朝晖,张锡文,任玉新等.一种低耗散、无伪振荡的实用差分格式[J].水动力学研究与进展(A辑),2001,16(02):195-199.[6]汤寒松,张德良,李椿萱.对流方程保单调CIP格式[J].水动力学研究与进展(A辑),1997(02):181-187.[7]赵飞,蔡志权,葛永斌.一维非定常对流扩散方程的有理型高阶紧致差分公式[J].江西师范大学学报(自然科学版),2014,38(4):413-418.642应用数学2019[8]李青,王能超.解循环三对角线性方程组的追赶法[J].小型微型计算机系统,2002(23):1393-1395.[9]ERDOGAN U.Improved upwind discretization of the advection equation[J].Numer.Meth.PartDiffer.Equ.,2014,30:773-787.[10]BOURCHTEIN A,BOURCHTEIN L.Explicitfinite schemes with extended stability for advectionequations[J]put.Appl.Math.,2012,236:3591-3604.[11]KIM C.Accurate multi-level schemes for advection[J].Int.J.Numer.Methods Fluids.,2003,41:471-494.A High-Order Compact Difference Scheme for Solving the1DConvection EquationHOU Bo,GE Yongbin(School of Mathematics and Statistics,Ningxia University,Yinchuan750021,China)Abstract:In this paper,a two-level implicit compact difference scheme for solving the one-dimensional convection equation is proposed.Taylor series expansion and correction for the third derivative in the truncation error remainder of the central difference scheme are used for the discretization of time and space.The local truncation error of the scheme is O(τ4+τ2h2+h4),i.e.,it has the fourth-order accuracy in both time and space.The unconditional stability is obtained by the von Neumann method. The accuracy and the stability of the present scheme are validated by some numerical experiments.Key words:Convection equation;High accuracy;Compact difference scheme;Unconditional sta-bility;Finite difference method。
流体流动的高精度高分辨率格式
流体流动的高精度高分辨率格式High Accuracy and High Resolution Schemes of Fluid Flows 采用高精度高分辨格式不仅可以降低对网格规模的要求，而且能够正确分辨其中复杂的流动现象。
近年研究表明，在提高数值模拟可靠性和有效性方面，高精度高分辨率要求已成为计算流体力学（CFD）技术发展中的一个决定性因素。
高精度高分辨数值方法的研究已极大地推动了CFD学科的快速发展和应用型人才的培养，它也已在很多实际应用领域例如航空航天、海洋船舶、核物理等的产品 (如飞行器)设计、制造和研究中发挥着重要作用。
本课程即是围绕当前CFD及应用研究中的核心高精度高分辨率数值方法而设计的。
课程主要讲授高分辨激波捕捉方法、高精度有限体积方法、高阶紧致差分方法、湍流计算方法、同伦分析方法及应用、并行计算等。
内容涉及线性双曲方程的差分格式和拟线性双曲型守恒律方程的守恒型差分格式，高分辨激波捕捉方法(TVD, ENO, WENO等)，间断Galer kin方法，高阶紧致格式和应用，结构网格及非结构网格的有限体积方法，高精度有限体积方法的最新进展，结构网格及非结构网格的有限体积方法，高精度有限体积方法的最新进展，湍流计算方法（包括RANS,LES及DNS）、转捩预测方法，并行计算，同伦分析方法的基本思想、广义同伦概念及其在力学、数学中应用等。
通过模型方程和精选的示例介绍基本原理、数学理论，阐释求解问题的思路及实践应用；以专题方式介绍最新发展、若干前沿，以及复杂流体流动直接数值模拟等知识。
课程兼具理论性、实用性和前沿性。
本课程面向力学、航空航天、数学及相关理工科专业的研究生、高年级本科生。
本课程负责人为田振夫教授，并邀请国内CFD 及其应用领域三位知名教授组成教学团队。
李新亮,研究员/博导 中国科学院力学研究所高温气体动力学国家重点实验室研究员，中国科学院大学岗位教授。
主要研究方向是计算流体力学，可压缩湍流与转捩，飞行器空气动力学等。
高精度差分格式及其数值解的逼近程度分析
三阶迎风偏斜差分
Fj
1 6
2u j1 3u j-6u j-1 u j-2
kr
1 6
cos2
4 cos
3, ki
1 6
sin
2
8sin
五阶迎风偏斜差分
Fj
1 60
3u j2
30u j1 20u j
60u j1 15u j2 2u j3
kr
1 cos3
30
6 cos2
11
二阶导数的紧致差分
• 传统型差分的截断误差项的再次离散 • 四阶紧致：
1 12
S
j 1
5 6
S
j
1 12
S j1
x2u
j
• 若边界点S0和SN已知，可用解三对角矩阵方程得
到所有网格点上的差分
12
二阶导数的紧致差分(cont.)
• 一般形式
l0
l
S jl
S jl 2
bl
l 0
u jl
2u j l2
27
耗散效应
u f
t x
m
2 m x 2 m 1
2m f x 2 m
m
2 m 1x 2 m
2m1 f x 2 m 1
x
~2
u x
x
a
2u x 2
耗散
色散
f cu,u uˆ(t) exp(ikx)
2m1 f
~2
m
2 m x 2 m 1
x 2 m 1 f
f cx
u
m
(四) 高精度差分格式及其 数值解的逼近程度分析
• 指大于二阶精度的格式 • 要求准确模拟小扰动量长时间、远距离传播的速度和幅值 • 用于计算噪声、DNS等 • 有传统高精度差分和紧致差分两种
电除尘器前烟道流量分配的数值模拟
电除尘器前烟道流量分配的数值模拟摘要：针对新疆某电厂电除尘器入口气流不均匀导致的除尘效率低的问题，通过在除尘器前烟道内的增加导流板和气流挡板的方法，来优化流量的分配。
本文使用数值模拟的方法对优化前后的电除尘器入口气流流量分配进行计算。
关键词除尘器前烟道，流量分配，数值模拟1引言影响电除尘器除尘效率的因素主要有：粉尘的比电阻影响、气流的均匀性、除尘器进口的含尘浓度、漏风、烟气量、极板间距等。
经过多个电厂对于本厂所用电除尘器的历史测量数据分析、竣工图的查阅等工作，可以认为粉尘的比电阻影响、漏风、集尘极的振动、烟气量、极板间距等因素不是造成除尘器效率低的主要原因[1,2]。
要想对烟道的结构进行优化目前比较常用的研究方法是借助计算流体动力学（CFD）方法。
CFD方法是一种有效的研究流体动力学的数值预测方法，它能够描述复杂几何体内部的三维流动现象并在整个设计过程中对设计者提供帮助。
在设计的初期，CFD方法可以快速地评价设计并做出修改，而不需要花费原型生产和反复测试的代价；在设计中期，CFD手段用来研究设计变化对于研究对象的影响，减少未预料到的负面影响；在设计完成后，CFD可以提供试验中无法测量的一些流动细节，帮助设计者分析流动现象产生的原因，改进设计。
CFD方法大大减少了设计的费用、时间同时降低了新设计带来的风险。
通过计算流体力学技术可优化电除尘器前烟道的结构，比如增加烟气入口导流板以及气流挡板等，以便得到更均匀的分布流场，为烟道的优化设计提供可靠依据，极大地节省设计时间及费用[3,4]。
本文应用计算流体动力学软件ANSYS Fluent，对新疆某大型火电厂电除尘器前烟道进行了数值分析，并根据分析结果得到了除尘器入口气流流量分配均匀的烟道结构。
2电除尘器前烟道流场的数值模拟本文按照新疆某电厂除尘器前烟道的实际尺寸（如图1所示），在Gambit 软件中建立了计算的几何模型，如图2所示。
图1除尘器前烟道的实际尺寸图2除尘器前烟道的三维计算模型使用三维非结构化网格对全流道进行网格划分。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
~m; (~10m)
最小尺度涡： 微米量级； 中小尺度涡： 0.1mm~ 1mm量级
直接数值模拟： 分辨出最小尺度涡； 网格量极为巨大： 计算量、存储量大
对高性能计算需求强烈
二、 高精度CFD软件 OpenCFD
OpenCFD: 作者开发的一套高精度、开放的CFD程序
1） OpenCFD-SC : 高精度差分 2） OpenCFD-EC: 有限体积 3） OpenCFD-Comb: 化学反应计算
computational region
200
z 300
outflow boundary
buffer region
400
500
示例2：有攻角小头钝锥边界层转捩的DNS
3） 计算结果验证
a.不同网格之间的比较 （网格收敛性） b. 不同扰动形式的结果比较 （模型正确性） d. 与Horvath静风洞实验比较 （与实验结果比较） c. 与 eN方法比较 （与理论结果比较） e. 与Stenson实验比较 （与实验结果比较）
面向工程计算开放的CFD代码(Open CFD code for Engineering Computing)
算法： 多块结构网格有限体积 + RANS 特点： 适用于复杂外形工程计算
差分-有限体积混合方法
3） 高精度化学反应模拟程序—— OpenCFD-Comb
面向化学反应的开放CFD程序： An OpenCFD code for Combustion
高精度差分格式及 湍流数值模拟 （三）
Part 3 可压缩湍流DNS
1. 背景 2. OpenCFD软件简介 3. 典型可压缩湍流的直接数值模拟示例 4. 湍流模型的评估及改进 5. 湍流燃烧的DNS 6. 小结
一、 背景
➢ 直接数值模拟 （DNS） 是研究湍流机理、 模型及 控制的重要手段
✓ 无需使用湍流模型： 准确度高 ✓ 分辨湍流全部细节： 精度高
➢发卡涡形成及发展的涡动力学机制
2） 拟序结构与切应力、壁面摩擦阻力
➢雷诺瞬时切应力集中在拟序结构周 围。 ➢新旧拟序结构干扰摩擦阻力剧烈增
加。
➢ 发卡涡导致局部高剪切—— 诱发
非线性Breakdown
瞬时雷诺切应力等值面
3) 发卡涡头部脱落
➢ 发卡涡发生断裂，头部成为孤立自 由涡。
➢ 自由涡脱离主体，沿边界层外层向 下游翻滚运动。
头激波、熵层、横流； 高Reynolds数、高Mach数； 无成熟转捩（湍流）理论及模型；
DNS: 研究该问题的有力手段 尚无他人DNS报道
a
头激波
边
界
横流
层 转
捩
线
示例2：有攻角小头钝锥边界层转捩的DNS 2） 计算参数及方法
Ren(头半径度量)
10000
Ma a(攻角) Tw(壁温) T 半锥角
M6 三维翼湍流LES
8192-24万
65.14%
M6 三维翼湍流LES
32768-380万 （含从核）
约60%
M6 三维翼湍流LES
➢ OpenCFD程序 异构并行移植及优化
➢ GPU系统 （天河-1号）
➢ MIC系统 （天河-2号）
内存
➢ 国产众核系统 （神威太湖之光）
2）多块结构网格有限体积程序—— OpenCFD-EC
三、 典型可压缩湍流直接数值模拟示例
Mach 0.7-8 平板边界层湍流DNS
Mach 0.7-6 槽道湍流DNS
最大 网格 16亿
网格3.2亿
Mach 6 有攻角钝锥边界层湍流DNS
RAE2822 超临界翼型三维绕流LES及DNS
压缩折角激波-边层干扰DNS
示例1： 平板边界层湍流中的拟序结构
6
1度
294K
79K 5
数值计算： 自主开发的程序(OpenCFD-SC)
直接数值模拟； 高精度差分法
对流项: 7阶精度WENO;
粘性项： 8阶精度中心;
a
时间项： 3阶精度Runge-Kutta
头激波
边界层
湍 流
r
120
100
80
inlet boundary
60
40
20
00
100
upper boundary
✓高精度差分方法， 适合于化学反应湍流高精度模拟
➢ Chemkin-like 化学反应机理接口 ； ➢ Chemkin 热力学模型、 物理模型 （粘性、热导、扩散系数） ➢ 高精度差分方法： OMP6, GVC, WENO7 ➢ 高阶滤波，保证计算稳定性； ➢ MPI并行， 可扩展到百万CPU规模；
➢ 边界层外层的间歇现象与自由涡的 脱落有关。
4) 压缩性效应对拟序结构的影响
Ma=0.7 （平板） Ma=0.7 （尖锥） 板）
Ma=2.25 （平
Ma=6 （平板） ➢ Ma=6 平板边界层 ：准流向涡而不是发卡涡
示例2：有攻角小头钝锥边界层转捩的DNS
1）背景 飞行器的常用头部外形； 转捩及湍流机理复杂；
瞬时温度场
流动机理、 评估、改善湍 流动控制 流模型
湍射流的涡量分布：DNS
RANS
➢ DNS对高性能计算需求强烈
年代
2080
2045
2010 2000 1990
1980
DNS LES
Splart 的 预期
非定常
RANS/DES
粘性NS方程 (RANS)、定常
无粘流/Euler方程
势流、面元法
大型飞机整机流场/气动力计算
1）高精度有限差分求解器OpenCFD-SC
Open CFD code for Scientific Computing; 面向科学计算开放的CFD代码
算法： 高精度差分 特点： 精度高、适用于湍流DNS及LES等
外形相对简单， 但流场复杂
➢ OpenCFD-SC 的大规模测试及计算
时间 2007.12 2010.12 2013.4 2014.12 2015.12
计算机 国产并行机 天河1A 国产并行机 天河-2 国产并行机
计算规模（核） 最大规模时并 问题描述 行效率
4096-58240
50%
Mach 6 钝锥转捩DNS，计 算网格点3.5亿
256-16384
32768-98万 （含从核）
90% 37.5%
Mach 8 平板边界层湍流 DNS， 网格点4亿；
计算模型：平板边界层
边界层湍流DNS, Ma=0.7, 2.25,6 扰动:
二维不稳定波 + 一对三维斜波 (自然转捩） 壁面吹吸扰动 (Bypass)
➢ 应用示例： 平板边界层湍流中的拟序结构 1）拟序结构的形成及演化规律
动画演示： 平板边界层拟序结构的形成及演化
流场中涡结构的瞬时演化图（Q=10的等值面）