天津大学硕士学位论文.pdf

分类号 密级 天津大学硕士学位论文天津大学硕士学位论文 专 业 流 体 力 学 研 究 生 王 毅 指导教师 周恒 教授 天津大学研究生院 天津大学硕士学位论文天津大学硕士学位论文 二维超音速平板边界层扰动波二维超音速平板边界层扰动波 传播的数值模拟传播的数值模拟 申请硕士学位 专 业 流 体 力 学 研 究 生 王 毅 指导教师 周恒 教授 天津大学机械工程学院力学系 2001 年 11 月 Numerical Simulation on the Evolution of T S Wave in 2 D Supersonic Plate Boundary Layer A Thesis Submitted for the Degree of Master Major Fluid Mechanics Graduate Student Wang Yi Supervisor Prof Zhou Heng Department of Mechanics Tianjin University November 2001 天津大学硕士学位论文 摘要 摘要 摘要 近来可压缩流动的稳定性和转捩问题正在成为新的研究热点 迄今为止 大 多数的研究还限于线性理论 由于超音速流的稳定性 转捩的实验对实验设备和 实验技术要求非常之高 细致的实验结果也很缺乏 这在很大程度上也阻碍了非 线性理论的建立 目前 一些学者试图将不可压缩流动稳定性的非线性理论推广 应用于可压缩流动中 以研究有限幅值扰动的非线性演化 然而 在超音速尤其 是高超音速范围内 扰动波的幅值达到一定程度时 流场中可能产生小激波 如 果有小激波 则扰动解将变得不连续 这样 建立在光滑流场基础上的不可压缩 流动的非线性理论就不适用于带有间断的可压缩流场 因此探讨超音速边界层中 当扰动达到一定的幅值是否有小激波的出现 为建立可压缩流稳定性非线性理论 提供一定的依据是有必要的 本文分别用三种格式 对于不同马赫数下超音速边界层的扰动的演化进行了 直接数值模拟 得到以下结论 1 NND 格式尽管作为目前捕捉激波较为有效的格式 但由于其只有二阶精 度且数值耗散较大 并不适合于计算小扰动演化等稳定性问题的数值模 拟 2 在马赫数为 1 4 6 0 的范围内加入 T S 波 并无小激波产生 因此 在该 马赫数范围内 将不可压缩非线性理论推广到可压缩流的研究中是可行 的 关键词关键词 小激波 T S 波 直接数值模拟 超音速边界层 天津大学硕士学位论文 ABSTRACT ABSTRACT Recently the problem of transition and stability in supersonic boundary layer has attracted more and more attention due to its technical importance Up to now most analyses were limited to linear theory In addition because the requirements on techniques and equipments for the experiment of supersonic flow are too severe detailed experiments are very rare thus the establishment of a non linear theory lacks a solid foundation Nowadays some people have tried to extend the nonlinear theory of hydrodynamic stability for incompressible boundary layer to compressible boundary layer However shocklet may appear when disturbance is large enough in the supersonic flow If shocklet does exist the disturbance field is not continuous any more In this case the incompressible nonlinear stability theory which assumes the continuity of flow cannot be simply extended to the compressible flow Therefore in order to form the base of compressible stability theory it is necessary to analyze if shocklets exist or not in a supersonic boundary layer In this thesis three kinds of schemes have been used to simulate the evolution of T S wave in supersonic boundary layer under the condition of different Mach numbers Two conclusions can be draw 1 Although NND scheme is one of the most powerful shock capturing schemes it is not suitable for the simulation of small disturbance due to its low accuracy and large numerical dissipation 2 No shocklet appears when T S wave was introduced with the Mach number ranging from1 4 to 6 0 Therefore in the above range it is reasonable to extend the nonlinear theory of hydrodynamic stability for incompressible boundary layer to compressible boundary layer Key Words Shocklet T S wave Direct Simulation Supersonic Boundary Layer 天津大学硕士学位论文 目录 目目 录录 摘要 第一章 绪论 1 1 1 稳定性与转捩 1 1 2 可压缩流动线性稳定性理论的发展 3 1 3 非线性稳定性理论 8 1 4 本文的目的和意义 9 第二章 控制方程及数值方法 12 2 1 控制方程 12 2 2 差分格式 15 2 3 边界条件 16 第三章 小扰动演化的直接数值模拟 27 3 1 相似性解和基本流 27 3 2 特征值与特征函数 29 3 3 小扰动的演化 32 第四章 结果分析与讨论 33 参考文献 41 致谢 44 天津大学硕士学位论文 第一章 绪论 1 第一章第一章 绪论绪论 1 1 稳定性与转捩 流动稳定性是指某种性态的流体运动受到扰动后恢复原来形态的能力 若流 动能恢复原来的形态 则是稳定的 反之是不稳定的 层流流动失去稳定性是湍 流发生的第一个阶段 这个概念是 Reynolds 在 1883 年提出的 我们知道 Reynolds 在作圆管流动实验时 除了发现流体具有层流和湍流两种不同的流态外 还提出 了一个无量纲数 Reua 作为从层流到湍流的临界判据 其中 u 为圆管中心的 流速 a 为圆管半径 为运动粘性系数 这个无量纲数被称之为雷诺数 他的 实验装置可以使层流状态一直保持到 Re 13000 后来的实验工作者在入口处引 进有限扰动或采用粗糙管壁发现临界 Re 可降低到 2000 而用匀称的流动和光滑 的管壁则可使临界 Re 达到 100000 甚至更高 当 Re 低于 2000 时很难发生转捩 在高 Re 时层流状态不可能维持 在低 Re 时湍流不易发生 表明流体运动存在 着稳定性问题 从数学上看 流体运动的任何一种形态都应是 Navier Stokes 方程在一定 初值和边值条件下的解 在高 Re 下虽然也照样可得到层流解 但这种解在物理 上是不能观察到的 所以 并不是运动方程的每一个解 即使是精确解 都能 够实际上在自然界中出现 在自然界中出现的流动不仅要服从流体动力学方程 而且应该是稳定的 Landu L D Lifshitz E M 1959 为了研究流体的流 动状态是否稳定 人们从 Navier Stokes 方程在一定初 边值条件下得到的层流 解上 迭加上一个理论上无限小的扰动 然后通过分析扰动能量或扰动幅值的衰 减或增长来判断某种流动状态的稳定性 这就是目前应用较广泛的线性稳定性理 论的基本做法 本文只讨论平行或近似平行剪切流的线性稳定性问题 这类流动 包括了管道内的流动 平行平板间的流动以及边界层流动 Reynolds 的文章发表以后 转捩的概念被很好的建立起来 Rayleigh 1880 在不可压流的条件下 得到了两个关于无粘流稳定性的重要定理 天津大学硕士学位论文 第一章 绪论 2 1 流体的速度剖面 yU有一拐点 0 U 是扰动能增长的必要条件 后 来Tol1mien又证明 存在拐点也是扰动能增长的充分条件 2 若把扰动波取为如下形式 ctxi ey 其中 是实数 c为复数 并0 U 对于0 i c即中性扰动情况 至少 有一点 0 yy 这里cU 即在流体内部的某点 相速度等于平均流速 对应于了cU 这点称为临界层 然而 不可压缩流的无粘理论显然还无法解释如平面Poiseuille流和平板边界层 等没有速度拐点的流动失稳并转捩的机理 二十世纪初 Orr 1907 和Sommerfeld 1908 分别给出了粘性不可压缩 的平行剪切流稳定性的重要方程 即Orr Sommerfeld方程 以后许多有关平行 剪切流的稳定性的实质性工作都是在此基础上做的 Hesienberg 1924 在其博 士论文中克服了数学上的难题 特别是与临界层有关的困难 证明了平面 Poiseuille流动在大雷诺数下是不稳定的 这是平行剪切流粘性不稳定性的第一 个重大突破 Tollmien 1929 和Schlichting 1933 进一步发展了Heisenberg 的渐近方法 具体计算了小扰动的解 预测了在转捩过程中首先会有扰动波出现 以后这类解就称为Tollmien Schlichting波 但当时未能为实验证明 Burgers 1924 和Van der Hegge Zignen 1924 的实验得出了平板流的转捩雷诺数约 为35000 但并没有发现任何T S型的放大扰动 Dryden 1934 采用当时新发 明的热线风速仪进行了仔细的测量 发现了湍流脉动 但并未探测到理论所预测 的任何选择性的放大 Taylor 1939 甚至断言稳定性理论与边界层转捩很少或 根本没有关系 这种情况一直持续到40年代初期 美国国家标准局在Dryden的 领导下 建造了一个来流湍流度非常低 0 02 的新风洞 Schubauer和Skramstad 1947 在此风洞中进行的实验 才证实了Tollmien和Schlichting的理论预测 即T S波的存在 并且证实了中性曲线与理论值符合得很好 也找到了早期实验 不能探测到T S波的原因 即当时风洞的背景湍流度非常高 将扰动淹没了 以 天津大学硕士学位论文 第一章 绪论 3 后 Liepmann 1943 的类似风洞试验和Wartmann 1955 的水流试验均表明 了 线 性 理 论 作 为 层 流 剪 切 层 不 稳 定 性 的 机 理 是 正 确 的 对 于 求 解 Orr Sommerfeld方程 在数学上遇到了很大的困难 经过Heiseberg Tollmien Schlichting 及林家翘等人的相继努力 到1945年才使用渐近法得到该特征问题 解析解的方法有了坚实的基础 他们以及后人在此方面的贡献推动了应用数学的 发展 五十年代以后 随着计算机的发展和应用 数值计算开始应用于求解 Orr Sommerfeld方程的准确解 Thomas 1953 用有限差分法 Davey 1973 用打靶法 Orszag 1971 用正交函数法 Chebyshev函数 分别求解了平面 Poiseuille流问题 其中Orszag的结果被认为精度最高 然而 应当指出的是 线性稳定性理论也有其固有的缺陷 首先 它无法反 应背景湍流度对转捩的影响 其次 它无法反应扰动强弱对稳定性的影响 并且 它还无法说明平面Couette和圆管Poiseuille流动失稳的机理 至今为止 对这两 种流动的线性稳定性分析都只能得出稳定的结果 线性稳定性理论只预计了层 流在某一Re下失稳 但这并不等于湍流的出现 虽然基于线性稳定性理论特定 放大率的经验性转捩估算 在对二维和轴对称流动上取得令人满意的结果 但完 整的转捩理论并不存在 转捩是层流流动经过一系列复杂的非线性的时 空的演 化才最终达到非定常 无序的湍流状态 而线性理论显然无法刻画非线性影响起 主导作用的阶段 线性稳定性理论作为人们试图解释层流流动失去稳定机理的一种尝试 虽然 有其自身的局限性 但不容忽视的是该理论在绝大多数情况下 探讨流动失稳的 机理 起着十分重要的作用 它也是后来进一步发展非线性理论 如弱非线性理 论 分叉理论等的基础 并已广泛的应用于层流控制 机翼设计等实际工程中 所以 线性稳定性理论仍然是人们用来研究层流到湍流过渡过程的主要手段 1 2 可压缩流动线性稳定性理论的发展 早在40年代 Lees和Lin就利用小扰动理论为可压缩流动的稳定性理论奠 天津大学硕士学位论文 第一章 绪论 4 定了一个初步的基础 但当时以及随后的几十年里 这个理论并未象不可压缩流 动稳定性那样引起人们的足够重视 近些年来这个领域能够重新活跃起来 则完 全是由于超音速民用机研制和航天技术进一步发展的需要 而人们对这个速度段 发生转捩的物理图画 发生机理知之甚微 直到今天 这种状况也只是得到了一 定程度上的改观 作为本文工作的一个引导 我们有必要回顾一下可压缩流动稳 定性的发展历程以及比较其与不可压缩流动稳定性理论的差异 线性可压缩平行流的稳定性分析的基本方程 是由Lees和Lin 1946 最早 运用小扰动理论得到的 以后 Dunn和Lin 1955 又给出了该方程的一个简化 形式 他们的主要贡献是无粘的稳定性理论的研究 这个理论后来被证明对可压 缩流动稳定性理论的发展起了重要作用 虽然 Lees和Lin最初建立这个理论时 只考虑了边界层 但是 Lees和Gold 1964 却成功地把这个理论推广应用于 自由剪切流 早期的这些稳定性分析工作主要是采用渐近的方法 第一个直接从 粘性稳定性方程中用数值方法计算出简正模态特征值的是Brown 1962 更为 精细的数值计算工作则是由Mack 1965 来完成的 他考虑了二维和三维扰动 并得到了有粘和无粘稳定性分析的结果 这一结果也反映了当代在不可压缩和可 压缩平行流动稳定性分析的发展水平 Gropengiesser 1969 用无粘的空间模式 的稳定性理论得到更广泛的可压缩混合层的数值计算结果 许多在从前的分析工 作中没有注意到的与不可压缩流动稳定性理论的重要差别 在数值计算中被揭示 出来 比如 可压缩理论中第二模式的发现 Mack l969 就是很好的证明 在 Nagfeh和El Hady 1980 的工作中他们考虑了非平行性对稳定性的影响 Mack 1979 和Lekoudis 1979 分别用稳定性理论对无限翼展后掠翼的三维边界层 进行了稳定性分析 Malik和Orszag 1981 Malik 1982a 描述了一种由他们 提出的计算时间模式的特征值的有限差分方法 这是一种很有效的方法 这无疑 又向工程应用迈进了一步 几乎同时 Malik 1982b 又提出了更高精度的二点 四阶紧致格式 通过在无限翼展后掠翼上的三维可压缩边界层的稳定性分析的数 值试验表明 特征值计算的精确性 可靠性以及收敛速度等都有较大的提高 虽 然 此时不可压缩流动稳定性的研究已经进行得十分广泛 但可压缩流动稳定性 的研究仍然局限于二维或三维平板 Gasperas 1985 首先用Runge Kutta方法 直接求解稳定性方程 对尖锥外形进行了稳定性分析 随后 Gasperas 1987 天津大学硕士学位论文 第一章 绪论 5 又对该方法进行了进一步改进 此外Mack 1987 Ma1ik 1987 都对有尖锥 外形的轴对称边界层进行了稳定性分析 并与实验进行了比较 Balakumar和 Reed 1991 完成了旋转锥的稳定性分析工作 Mack 1978 与Malik 1987 分别讨论了压力梯度对平板边界层的第一模式与第二模式的影响 可压缩稳定性的实验研究 相对起步较晚 真正对可压缩流动稳定性进行实 质研究的实验 是始于Laufer和Vrebalovich 1960 的工作 他们在对来流马 赫数为1 6到2 2时的平板边界层做稳定性实验时 发现并研究了自激励震荡现 象 Demetriades 1960 完成了马赫数为5 8的平板边界层的稳定性实验 他发 现虽然高速时的放大率比低速时小得多 但在一定的频率段内 高速边界层稳定 性的变化趋势与低速的稳定性曲线在形式上非常相似 Kendall 1966 先后研 究了马赫数在4 5和2 2时平板边界层的人工扰动现象 他的实验数据在马赫数 为4 5时与Mack 1969 的理论分析结果吻合较好 但在马赫数为2 2时 符合 情况要差一些 Lysenko和Maslov 1984 研究了平板边界层在马赫数为2到4 时 冷壁对其稳定性的影响 他发现第一模向稳定方向发展 第二模向不稳定方 向发展 并与第一模的不稳区分离开来 这些结论与可压缩流动的稳定性理论定 性上是一致的 Kendall 1975 Demetriades 1977 Stetson Thompson Donaldson 和Siler 1983 等人先后对尖锥在零攻角情况进行了稳定性实验 Kendall的实 验用的是半锥角为4度的尖锥 来流马赫数为8 5 壁面边界层外缘的马赫数为 7 7 Demetriades的实验用相同的锥 来流马赫数为8 0 壁面边界层外缘的马赫 数为7 0 在他们实验中所测的的放大率有些与Mack的理论分析吻合很好 但 有些则符合较差 Stetson等的实验 是用半锥角为7度的尖锥 来流马赫数为 8 0以及边界层外缘的马赫数为6 8 同样 实验的结果无论是与Mack对平板的 稳定性分析 再通过Mangler转换后的结果 还是 Mack 1987 对尖锥所作的 分析结果相比较 都存在某些吻合情况较好 而另一些则不能令人满意的现象 一般说来 吻合情况较差的是放大率 对此 Mack提出了两点解释 首先 实 验中边界层感受到的并不只限于计算上单一的不稳定波 其次 组成不稳定波的 响应也不只限于计算上的单一模态 还有可能激励起展向的波 此外 Stetson 等 1983 1984 还对有攻角的尖锥及零攻角的钝锥进行了稳定性分析 至今未 见到与理论分析的比较的报告 天津大学硕士学位论文 第一章 绪论 6 我们知道可压缩流动稳定性理论分析的基本方程是连续方程 动量方程 能 量方程及状态方程 而在不可压缩流动的稳定性理论中则无需考虑能量方程与状 态方程 方程形式上的差异 必然使得可压缩流动稳定性理论在某些方面有别于 不可压缩流动的稳定性理论 其中最重要的区别是可压缩流动稳定性理论中存在 着附加解 为了说明这个现象 我们引人相对马赫数的概念 2121 22 T MWU M 其中 M U W T分别是来流中的马赫数 流向速度 展向速度和温度 分别是扰动波流向及展向的波数分量 是扰动波频率 当1M 其中性波的相速度c 用来流速度无量化后 可能介于1和 M 11 之间 这里 cos MM 是波角 这样扰动波的相速度总是小于最 多等于来流速度的结论就不成立了 更重要的是 对于给定的相速度 波数 的 唯一性也不存在了 对可压缩二维和轴对称物体的层流边界层 已进行了较仔细的研究 得到了 以下几点主要的结论 1 当5 3 M时 无粘不稳定性起主要作用 在无粘理论中 可压缩流的 dy dU dy d 与 不 可 压 缩 流 的 2 2 dy Ud 项 具 有 相 同 的 作 用 满 足 0 s yy dy dU dy d 的 s y叫作广义拐点 存在广义拐点是存在不稳定亚音 速扰动的充分条件 这里的亚音速是指相对于相速度的流速是亚音速 即相对马赫数1M时总是成立的 即1 2 M是扰动失稳的 充分条件 3 考虑三维扰动时 第一模的放大率随三维扰动波与流向夹角的增加而增 加 而高模式的放大率则随之而减小 1 3 非线性稳定性理论 尽管线性稳定性理论可以说明那种剖面不稳定 分辨出放大最快的扰动频 率 但当扰动幅值增长起来以后 将出现不能用线性理论解释的现象 因此 建 立非线性理论是十分必要的 Landau在1944年提出 当扰动增长后 非线性的影响会起作用 因此由线 性理论得到的幅值a的演化规律a dt da i 将变为 53 3 aoaAa dt da i 这就是著名的Landau方程 它只包含a的奇次项 系数 i A为Landau系数 Stuart 1960 提出了弱非线性理论 给出了求解Landau系数的方法 这一理论成了 后来许多工作的出发点 Stuart的理论在应用时受到一些限制 它只能求解存在 中性曲线 并且1 i 的情况 无法处理远离中性曲线和无中性曲线的情况 天津大学硕士学位论文 第一章 绪论 8 周恒等 1982 将非线性振动理论中的K B方法应用于流动稳定性的研究中 其实质和Stuart的方法一样 但其形式更带有普遍性 更便于应用 同年 周恒 又提出了 人为中性 的方法 使弱非线性理论可用于非中性的情况 Herbert 1980 1983 提出了一种新的确定Landau系数的方法 去掉了Stuart理论中 1 i 的限制 周恒 1989 1991 发现 虽然弱非线性理论已被广泛应用于 流动稳定性研究 然而对于转捩和自由剪切流中扰动波演化这类问题 原来的弱 非线性理论不能给出与实验相符的结果 针对这一问题 周恒 1992 1993 又 作了深入的研究 指出了原弱非线性理论带根本性的一些缺陷 并提出了改进的 办法 其主要区别在于 对每一个谐波 它的幅值均按其自身规律演化 区 别于原弱非线性理论中各次谐波幅值均与基本波幅值的某次方严格成比例 除基本波形状可由线性理论给出外 数值模拟结果证实 其它各次谐波 包括 平均流修正 形状是变化的 只能从它满足的发展方程求解 而原弱非线性理论 却认为所有各次谐波形状都不变 这些修正使得弱非线性理论可用于扰动演化问 题的研究 修正后的弱非线性理论用于平面Poiseuille流时 得到的结果与数值 模拟结果符合的相当令人满意 以上的结论都是对二维扰动分析得到的 Klebanoff 1962 Saric 1984 和Kachanov 1984 等人的实验表明 在平板 层流边界层中 虽然扰动一般先以二维形式开始 但在沿流向演化时 不可避免 的要出现三维扰动 且三维扰动的增长率大于二维T S波的增长 从明显的出现 三维扰动开始到转为完全湍流 不需要很长的过程 通常在几个T S波长范围内 就完成 根据T S波的幅值不同 有三种不同的三维扰动 分别称为C Craik 型 H Herbert 型和K Klebanoff 型 Craik 1971 提出的共振三波 理论和Herbert 1983 1984 提出的二次失稳理论分别解释了C 型和H 型 失稳 共振三波理论是在弱非线性理论的基础上 加入共振条件 满足此条件的 一个二维波和一对斜波的流向相速度相等 其能量交换关系保持不变 若幅值在 初始时增长则在较长时间内的演化不会有反复 将一直增长下去 二次失稳理论 将二次失稳的概念应用于研究三维扰动失稳问题 它认为当存在定常或准定常有 限幅值二维T S波时 三维次谐波通过参数共振而失稳 导致三维扰动的增长 赵耕夫和周恒 1988 通过一个实际的例子指出了二次失稳理论与共振三波理论 间存在某种内在的联系 Kachanov 1987 在重新分析Klebanoff等人的实验资 天津大学硕士学位论文 第一章 绪论 9 料后提出了一般共振理论 这些理论尽管还不能完全解决转捩问题 但已经形成 了不可压缩流动的非线性稳定性的较为完整的理论体系 基于与不可压缩流相似的原因 可压缩流动的非线性理论也急待建立和发 展 与不可压缩流动稳定性理论中一些较成熟的非线性理论方法向比较 可压缩 流动的非线性理论还只是刚刚起步 可压缩非线性稳定性分析在无粘理论中做了 一些工作 但也只限于象尾流这种单一模态的情形 Er1cbacher和Hussaini 1987 用直接数值模拟的方法 观察到了一个有限幅值的二维波与三维扰动相互作用触 发二次失稳的现象 以此为依据El Hady 1988 导出可压缩流动二次稳定性分 析方程 其方法同Herbert 1983 1984 在不可压流中的做法相同 但其结果 缺乏足够的实验依据 我们现在知道在不可压缩流动的稳定性理论中 被Herbert 称作Squire模的简正模态 已证明是构成亚谐失稳的主要模态 而可压缩理论的 一个重要特征就是存在许多的附加模态 而且有可能构成可压缩流动亚谐失稳的 主要模态是在附加模态中 因此 缺乏可供参考的实验资料 很难提出令人信服 的非线性理论 现有一些可压缩的非线性理论 包括二次失稳理论 Duck 1990 渐近分析理论等 作为可压缩流动稳定性非线性理论的一种探索和尝试有其重要 的价值 1 4 本文的目的和意义 近来由于航空 航天技术的发展 可压缩流动的湍流及转捩问题正在成为新 的研究热点 可压缩流动稳定性问题是其中的一个重要部分 特别是可压缩非线 性稳定性问题尚未建立起系统的理论 正急需发展 由于这类问题的复杂性 从 数学上求解析解是不可能的 在实验研究方面 可压缩流特别是超音速流的稳定 性 转捩的实验 对实验设备和实验技术要求非常之高 以致国内还不具备开展 实验研究的条件 国际上也只有 Stetson 等几个小组能在特定条件下进行一些实 验 但无法提供演化细节 由于缺少全面系统的实验研究 理论上难以构造令人 天津大学硕士学位论文 第一章 绪论 10 信服的可压缩流动稳定性的非线性理论框架 目前 国内 外都有一些学者试图 将不可压缩流动稳定性的非线性理论推广应用于可压缩流动中 以研究有限幅值 扰动的非线性演化 然而 在超音速尤其是高超音速范围内 扰动波的幅值达到 一定程度时 很可能产生小激波的运动 这一点 已经在超音速混合层的直接数 值模拟中得到了证实 傅德熏 马延文 1996 在对二维超音速混合层的直接数 值模拟时 首先发现激波的存在 曹伟 周恒 2001 系统的分析了二维超音速 混合层中小激波的存在 并用流体质点的熵的明显增加来证明小激波的存在 并 分析了其对流场结构的影响 小激波的出现将使流场发生较大变化 扰动解将变 得不连续 这样 建立在光滑流场基础上的不可压缩流动的非线性理论就不适用 于带有阶跃的可压缩流场 既然在超音速混合层中存在激波 那么很自然的我们 会提出超音速边界层中是否存在小激波这个问题 然而 国际上尚无人做过这方 面的实验或计算研究 在研究可压缩流稳定性理论时 也没有人考虑到激波的问 题 因此 验证在一定的扰动幅值下 可压缩边界层特别是超音速边界层中是否 有小激波 并研究其运动规律及对稳定性的影响是十分有意义的 袁湘江在 1999 年的博士后工作期间 曾用 NND 格式对马赫数为 4 5 时小扰动在超音速边界层的 演化做过直接数值模拟 当扰动幅值为 0 01 时 他在等马赫线图中发现第一个 周期处有等值线密集的地方 并根据密度在临界层处有一个相位的变化判断该处 有激波出现 然而这种经验性的判断并没有绝对的说服力 此外他计算出的扰动 的增长率与线性理论的结果差别较大 且缺少不同流动参数及不同扰动形式的计 天津大学硕士学位论文 第一章 绪论 11 算结果的比较 因此这样的结论还不能令人信服 本文分别用不同的格式 对于不同马赫数下超音速边界层的扰动的演化进行 了直接数值模拟 以探讨超音速边界层中当扰动达到一定的幅值是否有小激波的 出现 为建立可压缩流稳定性非线性理论提供一定的依据 天津大学硕士论文 第二章 控制方程及数值方法 12 第二章 控制方程及数值方法第二章 控制方程及数值方法 2 1 控制方程 1 二维守恒型 N S 方程 在直角坐标系下 可压缩流动 N S 方程可写为守恒型 具体形式如下 y F x E y F x E t U 2 1 表达式中 U为守恒型通量 x E 和 y F 为对流项 其中包括压力项 x E 和 y F 为粘性项 其中包括热传导项 具体表达式为 vpE pv uv v F upE uv pu u E E v u U nnn 2 2 2 2 yyyyx yy yx xxyxx xy xx qvu F qvu E 0 0 2 3 y v x u xx 2 3 2 y v x u yy 2 3 2 x v y u yxxy 1 2 1 22 p vuEn x T qx y T qy yx qqTpvu 分别为流体的流向速度 法向速度 压力 温度 密度 天津大学硕士论文 第二章 控制方程及数值方法 13 粘性系数 热传导系数 流向和法向热通量 2 方程的无量纲化 用特征长度 边界层中 取给定处的边界层动量厚度 来流速度 U 无穷 远处的密度 温度 T和粘性系数 分别对相应的参数进行无量纲化 以下 标 表示无穷远处的参数 上标 表示无量纲的量 即 PrRe1 2 2 M k U t t U v v y y x x U p p T T T U u u 粘性系数 满足 Sutherland s 公式 CT C T 2 3 1 2 4 其中 C 110 4 k T Pr 为 Prandtl 数定义为 Pr Cp k Cp 是定压热容 Pr 在计 算中取常数 0 72 为比热比 取为 1 4 将无量纲的量带入 N S 方程 得到的无量纲的 N S 方程 为方便起见略去上 标 则形式与 2 1 一样 UE F x q y q的形式不变 E F的中各个分 量的形式也基本不变 但粘性应力项要除以雷诺数 U Re 形式如下 yeyyyx eyy eyx xexyxx exy exx qRvu R R F qRvu R R E 0 0 2 5 状态方程变为 2 MTp 3 坐标变换 计算中为了提高边界层内的精度 对 y 方向靠近壁面处加密网格 我们进行 天津大学硕士论文 第二章 控制方程及数值方法 14 坐标变换 计算坐标系 与直角坐标系 x y 之间的对应关系为 yx yx 令 J 为坐标变换的雅可比行列式 即 yy xx yx J det 则坐标变换后的 N S 方程如下 vv FEFE t U 11111 2 6 JUU 1 yx FEJE 1 yvxvv FEJE 1 yx FEJF 1 yvxvv FEJF 1 其中 E F与 v E v F分别如 2 2 2 5 示 4 通量分裂 为了使用迎风差分格式进行计算 对流项应根据 Jacobian 矩阵特征值的正负 进行通量分裂 以 方向的 1 E 为例 设 Jacobian 矩阵 1 1 U E A 其特征值为 i i 1 4 则可找到S使 1 SSA 则通量可作以下分裂 1 1 1 1 111 1 11 USSUSSAUU A U U EE 111 1 1 1 EEUSSUSS A 有四个特征值 其中有两个重根 令 yx vuU 22 yx cC 这里 c 是 天津大学硕士论文 第二章 控制方程及数值方法 15 当 地 音 速 则 四 个 特 征 值 为 CU 1 U 32 CU 4 4321 diag 令 2 ii i 4321 diag 则 1 E可以 分裂为 1 E 2 2 1 2 2 1 2 2 1 2 2 1 3 22 21 3 22 21 3 22 21 31 1 nyxt yxy yxx EU P E vv uu JE 2 7 其中 421 1 2 41 2 c 23 1 2 1 22 p vuEn 1 F与 1 E具有类似形式 只须将上式中的 x y相应换为 x y即可 而CU 也相应变为 22 yxyx cCvuV 分裂后 原控制方程 2 6 变为 vv FEFFEE t U 1111111 2 8 2 2 差分格式 本文的计算中 对于对流项先后采用了NND格式 三阶弱迎风紧致格式和 五阶精度的迎风紧致格式 对粘性项 采用了中心型紧致格式 其具体形式如下 天津大学硕士论文 第二章 控制方程及数值方法 16 1 二阶精度NND格式 1 2 1 2 1 ii HH xx E 2 9 式中 2 1 2 1 2 1 iii EEH 2 1 2 1 2 1 modmin 2 1 ii i i EEEE 2 3 2 11 2 1 modmin 2 1 ii i i EEEE 其中E的上标 代表Stager分裂下正负通量 1 2 1jj j EEE 式中min mod a b 定义为 min 2 mod minba bsignasign ba 2 三阶精度的弱迎风紧致格式 ii iii ii iii EE xx E x E x E EE xx E x E x E 1 11 1 11 12 58 12 85 2 10 在靠近边界的网点上 则为三阶精度弱迎风紧致格式如下 11 1 11 1 45 2 1 2 54 2 1 2 iii ii iii ii EEE xx E x E EEE xx E x E 2 11 3 五阶精度的弱迎风紧致格式 天津大学硕士论文 第二章 控制方程及数值方法 17 2112 1 2112 1 3443612 12 1 23 3443612 12 1 32 iiiii ii iiiii ii EEEEE xx E x E EEEEE xx E x E 2 12 4 六阶精度的中心型紧致差分格式 x ffff x f x f x f iiiiiii 28 123612 2211 1 1 2 13 2 3 边界条件 由特征线理论 在边界上有两组特征波 一组为传入计算域的波 一组为传 出计算域的波 对传出计算域的波可用计算域内的信息得出 即用单边差分格式 计算导数 而对传入计算域的波则需要通过适当的边界条件计算导数 这里我们 采 用Poinsot and Lele提 出 的 局 部 一 维 无 粘 关 系 式LODI The Local One Dimensional Inviscid Relations 这个关系式是在局部忽略粘性项和平行于边 界的对流项而得到的 下面我们介绍一下LODI关系式 及本文中的具体处理方 法 1 局部一维无粘关系式LODI 略去 N S 方程 2 1 的粘性项和 y 方向的对流项 设边界平行于 y 坐标轴 得 0 x E t U 2 14 为了分析边界上的特征波的行为 我们把守恒型方程写成非守恒形式 即 将原 天津大学硕士论文 第二章 控制方程及数值方法 18 来基于守恒变量 t n EvuU 的方程化为基于基本变量 t p pvuU 的方程 0 x U U U U E t U U U p p p p 2 15 令 U E A p U U P 有 0 1 x U APP t U pp 2 16 令APPAp 1 p A可以分解成 1 ppp SSA 其中 4321 diag cu 1 u 32 cu 4 所以有 0 11 x U S t U S p p p p 2 17 其展开形式为 0 0 0 0 4 3 2 2 2 1 x u c x p t u c t p x v t v x p x c t p t c x u c x p t u c t p 2 18 令 天津大学硕士论文 第二章 控制方程及数值方法 19 x u c x p x v x p x c x u c x p 44 33 2 22 11 2 19 这里的 i 分别代表以 i 为速度传播的特征波 i 的正负代表波的传播方向 把 i 的定义式 2 19 带入式 2 18 用物理量对时间的导数表示 1 2 3 4 得 t u c t p t v t p t c t u c t p 4 3 2 2 1 2 20 这就是忽略了粘性项和平行于边界的对流项后得到的局部一维无粘关系式 LODI 由特征线理论 对于穿过边界流出计算域的波 其值完全由计算域内的信息 确定 我们用单边差分格式求得边界上基本变量的空间导数 代入 2 19 得到 流出计算域的 i 而对于穿过边界流入计算域的波 其值由计算域外的信息确定 则不能用计算域内的点差分求得 而必须给定计算边界条件 我们用LODI 2 20 结合具体的物理边界条件计算 i 本文中入口和出口处的具体做法 将在下文中 详细介绍 求得边界上所有的 i 后 代入用 i 表示对流项的N S方程中 即可 求出边界上通量的空间导数和基本变量在下一个时间层的函数值 我们还可以从另一个角度理解以上处理方法 对于双曲系统 我们一般对不 天津大学硕士论文 第二章 控制方程及数值方法 20 同方向的波分别采用迎风格式以保证计算格式的稳定性 对于流出计算域的波 若在边界上用单边差分 恰好满足迎风的要求 而对于流入计算域的波 如果还 用单边差分 则不满足迎风要求 所以需给定边界条件 2 用 i 表示对流项的N S方程 经过上节的推导知 x U S p p t 1 4321 p p pp pp ppp PS x U SPS x U APPP x U U U U E U U U U x E 11 1 31441421 2 22 3421 2 14241 2 421 2 212 1 2 4 2 2 2 2 1 2 2 1 v c u c vu c v cuu c c PS x E p 2 21 对于 y 方向 经相似的推导可得 y F 的表达式 31441421 2 22 14241 2 3421 2 421 2 212 1 2 4 2 2 1 2 2 2 2 1 u c v c vu cvv c c u c y F 2 22 天津大学硕士论文 第二章 控制方程及数值方法 21 其中 1 2 3 4 及特征函数 变为 y v c y p y u y p y c y v c y p 4 3 2 2 1 4 3 2 1 cv v v cv 4 3 2 1 3 计算坐标下的各种关系式 计算坐标下 处理边界条件所需的各种关系式的推导与直角坐标系中类似 下面我们不经推导给出本文中用到的各种表达式 1 特征值 22 22 4 3 2 1 yxyx yx yx yxyx cvu vu vu cvu 2 23 2 i 的定义式 yx yx yx yx yx yx vucJJp uv JpJ c vucJJp 22 4 22 3 2 2 22 1 4 3 2 1 1 2 24 天津大学硕士论文 第二章 控制方程及数值方法 22 3 LODI关系式 yx yx yx yx yx yx t v t ucJ t Jp t u t v t Jp t J c t v t ucJ t Jp 22 4 22 3 2 2 22 1 1 2 25 4 用 i 表示的对流项 1 E 34114 22 421 2 22 314 22 421 2 314 22 421 2 421 2 1 12 1 2 2 4 2 2 2 2 2 2 2 2 1 yx yx yx x yx y y yx x uvJ c vu c vu J c c v J c c u c E 2 26 以上各式变为 方向后形式基本不变 只需将式中的 yx 换为 yx 即 可 我们还可以看到 若将x y 代入各式 则可以变换回到直角坐标系 下的表达式 天津大学硕士论文 第二章 控制方程及数值方法 23 4 本文中边界条件的处理 边界条件分为物理边界条件和计算边界条件 物理边界条件是为了满足适定 性的要求 物理上需要的对边界上的变量的描述和限定 而对于物理边界条件没 有确定的物理量 则需要计算边界条件补充 如果对流项用紧致差分格式 则计 算边界条件还需给出边界上通量的导数值 二维欧拉方程和N S方程所需的物理 边界条件的个数见表2 1 我们根据Lele等的建议 用欧拉方程的无粘边界条件 附加上粘性边界条件来计算N S方程 二维欧拉方程 二维N S方程 超音速入口 4 4 亚音速入口 3 4 超音速出口 0 3 亚音速出口 1 3 表 2 1 二维流动满足适定性所需要的物理边界条件的个数 边界层中既有亚音速区 又有超音速区 由于亚音速和超音速的特征值的符 号不同 其特征波传播的方向不同 所需的物理边界条件的个数也不同 因此需 要分别处理 我们先介绍本文中用到的无粘边界条件 1 入口的超音速区 物理边界条件为vup 都给定 因为 1 2 3 4 均大于零 四个特征波都是 天津大学硕士论文 第二章 控制方程及数值方法 24 流入计算域的 计算对流项时 i 不能用定义式 2 24 而需用LODI关系 式 2 25 对于基本流 vup 对时间的导数都为零 代入 2 25 得 1 2 3 4 都为零 对于加扰动后的非定常流 边界上速度 温度 密度 压力随时间的变化等于扰动值 pvu对时间的导数 代入 2 25 可以求 得 1 2 3 4 2 入口的亚音速区 物理边界条件只有三个 tvu 给定 与超音速不同 1 小于零 因此 1 应该 用 2 24 计算 其中的空间导数用单边迎风格式 2 3 4 仍需用LODI 关系式 2 25 对于基本流 tvu 对时间的导数为零 从 2 25 可以推出 以下关系 12 1 0 3 14 对于加扰动后的非定常流 用与超音速区相同的处理方法求得 2 3 4 3 出口的超音速区 超音速出口不需要物理边界条件 所有边界的值由上游状态确定 因为 1 2 3 4 均大于零 四个特征波都是流出计算域的 i 用定义式 2 24 计 算 4 出口的亚音速区 亚音速出口需要一个物理边界条件 本文给定无穷远处的压力 p 2 3 4 大 天津大学硕士论文 第二章 控制方程及数值方法 25 于零 2 3 4 用定义式 2 24 计算 1 用下式计算 2 1 1 M pkppk 2 27 这是一种无反射边界条件 若k等于零 为完全无反射 若k不等于零 则 为部分无反射 这种边界条件将无穷远处压力的信息带回亚音速边界 反映 了其对计算域的影响 5 下