《推荐飞行器流动仿真讲稿第7章拟一维喷管流动的数值模拟》由会员分享，可在线阅读，更多相关《推荐飞行器流动仿真讲稿第7章拟一维喷管流动的数值模拟（105页珍藏版）》请在金锄头文库上搜索。
1、2013年年6月月飞行器流动仿真飞行器流动仿真Computational Fluid DynamicsThe Basics with Applications1第第7 7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法第第7.5节节 含激波喷管流动的含激波喷管流动的C
2、FD解法解法2p本章以拟一维喷管流动为例，讨论本章以拟一维喷管流动为例，讨论CFD方法的具体应用；方法的具体应用； p通过实际编写拟一维喷管流动数值模拟程序，掌握程序的调通过实际编写拟一维喷管流动数值模拟程序，掌握程序的调试手段与技巧，包括错误排除等，体会试手段与技巧，包括错误排除等，体会CFD求解问题的基本求解问题的基本过程和步骤，加深对过程和步骤，加深对CFD方法的理解。为进一步学习和运用方法的理解。为进一步学习和运用计算流体力学，自行开发流体力学数值程序或使用商业软件计算流体力学，自行开发流体力学数值程序或使用商业软件奠定必要的基础；奠定必要的基础； p本章对三个拟一维喷管定常流动问题，
3、使用本章对三个拟一维喷管定常流动问题，使用MacCormack两两步显式方法计算：步显式方法计算：亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法；解法；全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法；解法；用激波捕捉法求解含激波的喷管流动。用激波捕捉法求解含激波的喷管流动。 第第7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟3第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题p拉瓦尔喷管横截面积拉瓦尔喷管横截面积A是轴向距离是轴向距离x的函数：的函数：第第7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟一、流动模型一、流动模型喷管
4、入口接驻室，其截面积足够大喷管入口接驻室，其截面积足够大(理论上无穷大理论上无穷大)，驻室，驻室内流速内流速V0驻室内压强和温度驻室内压强和温度p0、T0为总压和总温；为总压和总温；喷管出口处流体参数喷管出口处流体参数为：为：pe、Te、Ve、Mae等；等；喷管出口外环境压强喷管出口外环境压强（反压）为（反压）为pa。4p喷管内部流动完全由进口总压喷管内部流动完全由进口总压p0和出口反压和出口反压pa的差值决定。的差值决定。设喷管任意截面流动参数均匀，流动仅随设喷管任意截面流动参数均匀，流动仅随x变化变化拟一维定拟一维定常流动常流动；不考虑流体黏性；不考虑流体黏性等熵流动等熵流动。 p对于本章
5、研究的三个问题，喷管内的流动状态分别为：对于本章研究的三个问题，喷管内的流动状态分别为：亚声速亚声速超声速等熵流动超声速等熵流动pa远小于远小于p0，内部流动由进口，内部流动由进口p0和面积决定，称为和面积决定，称为欠膨胀状态：欠膨胀状态：收敛段为收敛段为Ma1膨胀流；膨胀流；反压变化不影响内部流动，也不影响反压变化不影响内部流动，也不影响Mae 1和和pe只改变出只改变出口外的膨胀；口外的膨胀； 第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题5第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题含激波喷管流动含激波喷管流动pa比比p0小一些小一些
6、，扩张段出现激波，流动由进口扩张段出现激波，流动由进口p0和出口反压和出口反压pa的差值及面积变化的差值及面积变化决定，也是决定，也是过度膨胀状态过度膨胀状态：收敛：收敛段为段为Ma1膨胀流，激波后为膨胀流，激波后为Ma1压缩流；反压变化影响激波压缩流；反压变化影响激波位置；激波前后各为熵值不同的位置；激波前后各为熵值不同的等熵流动。等熵流动。全亚声速等熵流动全亚声速等熵流动pa比比p0小得不多小得不多，喷管内部流动由，喷管内部流动由进口进口p0和出口反压和出口反压pa的差值及面积变化决定，称为的差值及面积变化决定，称为过度膨过度膨胀状态胀状态：收敛段为：收敛段为Ma1膨胀流动；喉道膨胀流动；
7、喉道Ma1；扩张段为；扩张段为Ma1压缩流动；压缩流动；6p任意两个截面流动参数的关系（控制方程组）任意两个截面流动参数的关系（控制方程组）二、拟一维定常等熵流动的解析解二、拟一维定常等熵流动的解析解第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题7p等熵流动的解析解等熵流动的解析解第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题面积比公式面积比公式(7-6)(7-7)(7-8)(7-9)(7-10a)8第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题已知总温、总压和临界截面积已知总温、总压和临界截面积A*由由(7-6)式
8、解出待求截式解出待求截面的马赫数面的马赫数Ma，再用，再用(7-7)式式(7-9a)式求出其它参数；式求出其它参数； 已知总温、总压和待求截面的压强已知总温、总压和待求截面的压强由由(7-7)式解出待求式解出待求截面的马赫数截面的马赫数Ma，再用，再用(7-7)式式 (7-9a)式求出其它参数；式求出其它参数；已知总温、总压和某截面的马赫数已知总温、总压和某截面的马赫数由由(7-6)式解出临界式解出临界截面积截面积A*，再用，再用(7-6)式解出待求截面的马赫数式解出待求截面的马赫数Ma，用，用(7-7)式式 (7-9a)式求出其它参数。式求出其它参数。p等熵流动三类问题等熵流动三类问题9p守
9、恒守恒形式控制方程组形式控制方程组三、拟一维无黏流动的控制方程组三、拟一维无黏流动的控制方程组第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题双双曲曲型型方方程程组组10p非守恒非守恒形式控制方程组形式控制方程组第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题量热状态方程量热状态方程e=cvT，cv=Const双双曲曲型型方方程程组组11p无量纲化更容易观察流动，并可把流动规律推广到同类问题。无量纲化更容易观察流动，并可把流动规律推广到同类问题。合理使用无量纲量，使参与计算的各量具有差不多的量级，合理使用无量纲量，使参与计算的各量具有差不多的量级
10、，还可避免在计算中出现大小相差悬殊的数值；还可避免在计算中出现大小相差悬殊的数值；p使用无量纲控制方程组计算完毕后，根据需要可将无量纲量使用无量纲控制方程组计算完毕后，根据需要可将无量纲量重新化为有量纲量。重新化为有量纲量。p选择如下参考量：选择如下参考量：四、控制方程组的无量纲化四、控制方程组的无量纲化第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题温度温度：滞止温度：滞止温度T0无量纲内能无量纲内能：e=e/e0无量纲温度无量纲温度：T=T/T0密度密度：滞止密度：滞止密度0无量纲密度无量纲密度：=/0内能内能：滞止内能：滞止内能e0= cvT0长度长度：喷管全长：
11、喷管全长L无量纲长度无量纲长度：x=x/L无量纲时间无量纲时间：t=t/t0速度速度：滞止声速：滞止声速a0无量纲速度无量纲速度：V=V/a0时间时间：t0=L/a0面积面积：喉道面积：喉道面积At无量纲面积无量纲面积：A=A/At121. 非守恒非守恒形式控制方程组的无量纲化形式控制方程组的无量纲化第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题p以连续方程（以连续方程（7-42）为例推导无量纲形式。）为例推导无量纲形式。p将连续方程按以下方式改写将连续方程按以下方式改写13第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题p换成无量纲量换成无量
12、纲量参考量都是常数参考量都是常数化简为化简为与有量纲形式完全相同与有量纲形式完全相同进一步化简为进一步化简为p同理可得动量方程和能量方程的无量纲形式。同理可得动量方程和能量方程的无量纲形式。14第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题p为简便计，去掉为简便计，去掉“ ”，仍用原符号表示，只要记住它们都是，仍用原符号表示，只要记住它们都是无量纲量即可。则适合时间推进求解的无量纲形式的控制方无量纲量即可。则适合时间推进求解的无量纲形式的控制方程组程组152. 守恒守恒形式控制方程组的无量纲化形式控制方程组的无量纲化第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷
13、管流动的三个物理问题p采用与非守恒控制方程组无量纲化相同的方法，可得无量纲采用与非守恒控制方程组无量纲化相同的方法，可得无量纲形式的守恒控制方程组，仍用原符号表示形式的守恒控制方程组，仍用原符号表示16第第7.1节节 拟一维喷管流动的三个物理问题拟一维喷管流动的三个物理问题p将守恒形式控制方程组写成矩阵形式将守恒形式控制方程组写成矩阵形式守恒变量守恒变量U与原始变量关系与原始变量关系17第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法解法（1）p沿沿x轴喷管的面积分布是二次函数轴喷管的面积分布是二次函数第第7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟一
14、、问题的提法一、问题的提法喉道位于喉道位于x=1.5处；处；截面积是以喉道面积为参考值的无量纲面积；截面积是以喉道面积为参考值的无量纲面积；轴向坐标轴向坐标x以喉道距进口距离以喉道距进口距离(而不是喷管全长而不是喷管全长)为参考值。为参考值。p本节使用本节使用MacCormack两步显式方法求解两步显式方法求解非守恒非守恒形式控制方形式控制方程组。程组。1. 喷管的形状喷管的形状182. 物理条件物理条件第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p气体为标准空气，气体为标准空气，=1.4；p进口给定总温和总压，无量纲值：进口给定总温和总压，无量纲
15、值：p0=1，T0=1，0=1p对出口达超声速的情况，环境反压必须很低。本计算不需要对出口达超声速的情况，环境反压必须很低。本计算不需要反压的具体值，它不会影响喷管内部的流动。反压的具体值，它不会影响喷管内部的流动。3. 精确解精确解p亚声速亚声速超声速喷管等熵流动为超声速喷管等熵流动为定常流动定常流动；p临界状态出现在喉道截面，喉道截面积就是临界面积临界状态出现在喉道截面，喉道截面积就是临界面积p已知已知p0、T0和和A*，对给定截面，对给定截面A，由，由(7-6)式解出式解出Ma，收敛段，收敛段取取Ma 1 ；用；用(7-7) (7-9a)式求出式求出p、T和和V流动参数分布完全取决于面积
16、分布；流动参数分布完全取决于面积分布；19第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）20第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）21第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p对一维问题，网格生成非常简单，只需代数方法即可；对一维问题，网格生成非常简单，只需代数方法即可； p有限差分方法使用均匀网格，取网格点数为有限差分方法使用均匀网格，取网格点数为N=31则网格间距为则网格间距为二、网格生成二、网格生成网格点坐标为网格点坐标为于是，每个网格点对应的
17、截面积可由于是，每个网格点对应的截面积可由(7-73)式计算出来。式计算出来。22第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p非守恒形式控制方程组由非守恒形式控制方程组由(7-46) 、(7-48) 和和(7-50)式组成式组成三、三、MacCormack显式两步格式算法显式两步格式算法1. 非守恒形式控制方程组非守恒形式控制方程组用非定常方程求解用非定常方程求解定常问题定常问题时间相关法时间相关法时间推进时间推进23第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p类似于（类似于（6-17）式，以一阶）
18、式，以一阶前差前差离散非守恒控制方程组，计离散非守恒控制方程组，计算时间偏导数算时间偏导数2. 预估步预估步24第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p由时间偏导数由时间偏导数计算预估值计算预估值3. 校正步校正步p以一阶以一阶后差后差离散非守恒控制方程组，用离散非守恒控制方程组，用预估值计算时间偏导预估值计算时间偏导数数25第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）26第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p计算时间偏导数的平均值计算时间偏导
19、数的平均值p计算校正值计算校正值27第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p稳定性条件规定了推进的时间步长，即稳定性条件规定了推进的时间步长，即四、稳定性条件四、稳定性条件均为无量纲量（有量纲和无量纲表均为无量纲量（有量纲和无量纲表示的稳定性条件形式相同）示的稳定性条件形式相同）无量纲声速无量纲声速a，计算公式与有量纲声速不同，计算公式与有量纲声速不同仍用仍用a表示无量纲声速表示无量纲声速按无量纲化所取参考值，有量纲声速按无量纲化所取参考值，有量纲声速(7-66)28第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解
20、法（解法（1）p因为每个网格点的流速和声速均不同，所以在每一个时间层，因为每个网格点的流速和声速均不同，所以在每一个时间层，计算出的时间步长随网格点不同而不同，例如计算出的时间步长随网格点不同而不同，例如柯朗数柯朗数稳定性条件稳定性条件C1(7-66)时间步长时间步长t29第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p因此，在一个具体的时间层，因此，在一个具体的时间层，t的选择有两种方法：的选择有两种方法：当地时间步长法当地时间步长法每个网格点均使用各自的每个网格点均使用各自的t时间步长，时间步长，按照各自的时间步长推进按照各自的时间步长推进各网格
21、点流动参数对应于不各网格点流动参数对应于不同时刻，流场呈现出人为扭曲，不能真实反映流动的瞬态同时刻，流场呈现出人为扭曲，不能真实反映流动的瞬态变化。显然，当地时间步长法不能用于非定常流动计算。变化。显然，当地时间步长法不能用于非定常流动计算。但对时间相关法的定常解时间精度是无关紧要的，并且往但对时间相关法的定常解时间精度是无关紧要的，并且往往能更快收敛；往能更快收敛；最小时间步长法最小时间步长法选取所有网格点对应时间步长的最小选取所有网格点对应时间步长的最小值值对所有网格点均使用这个时间步长，因而每个时间步所有对所有网格点均使用这个时间步长，因而每个时间步所有网格点流动参数都对应着同一时刻网格
22、点流动参数都对应着同一时刻最小时间步长法能最小时间步长法能够真实反映流动的瞬态变化，是本课程采用的方法。够真实反映流动的瞬态变化，是本课程采用的方法。30第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）五、初场布置五、初场布置p对应于非定常流动计算的初始条件，时间相关法计算也需要对应于非定常流动计算的初始条件，时间相关法计算也需要布置初场，即给出布置初场，即给出t=0时刻所有网格点上所有参数的值；时刻所有网格点上所有参数的值；p一般来说可以任意布置初场，对计算的影响主要是：一般来说可以任意布置初场，对计算的影响主要是：初场越接近定常解，收敛越快、越容易
23、；初场越接近定常解，收敛越快、越容易；若初场偏离最终定常解太远，在最初几个时间步中流动参若初场偏离最终定常解太远，在最初几个时间步中流动参数变化非常大，可能导致计算不稳定。数变化非常大，可能导致计算不稳定。p可以选择如下的简单线性函数计算初场参数，它们沿喷管轴可以选择如下的简单线性函数计算初场参数，它们沿喷管轴线的变化趋势与真实情况一致线的变化趋势与真实情况一致t=0时时31第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）六、边界条件分析六、边界条件分析p非守恒形式控制方程组非守恒形式控制方程组(7-46)式、式、(7-48)式和式和(7-50)式是双
24、曲型式是双曲型的，有三条特征线，在的，有三条特征线，在xt平面上分别为平面上分别为C0：斜率：斜率C：斜率：斜率C：斜率：斜率信息沿信息沿C0传播的速度为传播的速度为V信息沿信息沿C传播的速度为传播的速度为V +a信息沿信息沿C传播的速度为传播的速度为V-a传播传播32第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）pxt平面上的三条特征线平面上的三条特征线p流动信息是沿着特征线传播的，所以边界条件与特征线的方流动信息是沿着特征线传播的，所以边界条件与特征线的方向及沿特征线的传播有关；向及沿特征线的传播有关；33第第7.2节节 亚声速亚声速超声速喷管等
25、熵流动的超声速喷管等熵流动的CFD解法（解法（1）p边界点边界点P0处，处，C0和和C+斜率均为正，斜率均为正，C0和和C+指向流场内指向流场内C0和和C+将信息从流场外传递到流场内将信息从流场外传递到流场内在边界上要求在边界上要求C0和和C+所输所输运的量为已知值；运的量为已知值；pC斜率斜率V-a0 把流场内部信息传递到边界外把流场内部信息传递到边界外传递的信息传递的信息不能在边界上给定，只能从流场内部确定。不能在边界上给定，只能从流场内部确定。1. 亚声速流亚声速流进口进口34第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p边界点边界点P1处，
26、处，C0和和C+斜率均为正，斜率均为正，C0和和C+指向流场外指向流场外C0和和C+将信息从流场内传递到流场外将信息从流场内传递到流场外C0和和C+所输运的量在边界所输运的量在边界上不能指定，只能从流场内部确定；上不能指定，只能从流场内部确定；pC斜率斜率V-a0 把流场外部信息传递到把流场外部信息传递到流场流场内内在边界上要在边界上要求求C所输运的量为已知值。所输运的量为已知值。2. 亚声速流亚声速流出口出口35第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p边界点边界点P0处，处，C0、 C和和C+斜率均为正，三条特征线均指向斜率均为正，三条特征
27、线均指向流场内流场内将信息从流场外传递到流场内将信息从流场外传递到流场内三条特征线三条特征线所输运所输运的量在边界上的量在边界上均均为已知值。为已知值。3. 超声速流超声速流进口进口36第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p边界点边界点P1处，处，C0、 C和和C+斜率均为正，三条特征线均指向斜率均为正，三条特征线均指向流场外流场外将信息从流场内传递到流场外将信息从流场内传递到流场外三条特征线三条特征线所输运所输运的量在边界上的量在边界上不能指定，只能从流场内部确定不能指定，只能从流场内部确定。4. 超声速流超声速流出口出口37第第7.2节
28、节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p综上所述，一维流计算需指定的综上所述，一维流计算需指定的物理边界条件物理边界条件如下表所示。如下表所示。亚声速亚声速Ma1进口边界进口边界给定给定2个条件个条件给定给定3个条件个条件出口边界出口边界给定给定1个条件个条件无无p在实际计算中，仅有在实际计算中，仅有物理边界条件物理边界条件通常是不够的，多数情况通常是不够的，多数情况下，在边界上不能指定所有的变量，但数值求解前必须已知下，在边界上不能指定所有的变量，但数值求解前必须已知所有变量；所有变量；38第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵
29、流动的CFD解法（解法（1）p因此，当物理条件不够时，必须补充其它边界条件，称为因此，当物理条件不够时，必须补充其它边界条件，称为“数值边界条件数值边界条件”与那些由内部流动规定的变量相对应，与那些由内部流动规定的变量相对应，并依赖于尚未求出的内部流动；并依赖于尚未求出的内部流动； p施加数值边界条件的原则：施加数值边界条件的原则：应该与流动的物理特征相容；应该与流动的物理特征相容；不能影响物理边界条件。不能影响物理边界条件。p理论分析和数值实验均表明，对许多格式，数值边界条件的理论分析和数值实验均表明，对许多格式，数值边界条件的处理会对精度、稳定性和收敛速度产生重大影响。例如，很处理会对精度
30、、稳定性和收敛速度产生重大影响。例如，很多隐格式理论上都是无条件稳定的，但由于边界条件处理不多隐格式理论上都是无条件稳定的，但由于边界条件处理不当，会变成条件稳定。当，会变成条件稳定。39第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）七、边界条件施加七、边界条件施加p拟一维喷管流动的边界条件应遵循理论分析结果进行施加。拟一维喷管流动的边界条件应遵循理论分析结果进行施加。1. 喷管亚声速流进口边界喷管亚声速流进口边界p物理边界条件两个物理边界条件两个：给定总温：给定总温T0和滞止密度和滞止密度0（由热状态方（由热状态方程知，总压程知，总压p0确定）确定
31、）驻室内的参数，随时间固定不变；驻室内的参数，随时间固定不变；因为进口速度很小，滞止参数与静参数相差无几，所以进口边因为进口速度很小，滞止参数与静参数相差无几，所以进口边界网格点界网格点1处的无量纲温度和密度为处的无量纲温度和密度为40第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p数值边界条件一个数值边界条件一个：允许流速：允许流速V变化，从流场内部确定。变化，从流场内部确定。假设假设V沿沿x的梯度为常数，即的梯度为常数，即在边界附近离散为在边界附近离散为一阶外推法一阶外推法41第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的C
32、FD解法（解法（1）2. 喷管超声速流出口边界喷管超声速流出口边界p物理边界条件零个物理边界条件零个：不需要给定物理边界条件。：不需要给定物理边界条件。p数值边界条件三个数值边界条件三个：允许所有流动参数变化，从流场内部确：允许所有流动参数变化，从流场内部确定。定。采用一阶外推法设定边界条件采用一阶外推法设定边界条件42第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）八、计算结果与分析八、计算结果与分析1. 解的时间推进过程解的时间推进过程p喉道处各流动参数随时间推进步数的变化。喉道处各流动参数随时间推进步数的变化。43第第7.2节节 亚声速亚声速超声
33、速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p由图可见由图可见计算初期时间导数很大，且有振荡。这种振荡与初期的瞬态计算初期时间导数很大，且有振荡。这种振荡与初期的瞬态过程中各种非定常压缩波和稀疏波的传播有关；过程中各种非定常压缩波和稀疏波的传播有关；随着时间步数的增加，时间导数迅速减小；随着时间步数的增加，时间导数迅速减小；对于定常流动，理论上时间步数需要趋于无穷时间导数才会对于定常流动，理论上时间步数需要趋于无穷时间导数才会达到零值，这在计算中是难以实现的；达到零值，这在计算中是难以实现的；实际上，当时间导数小于预设的很小的值（精度要求）且不实际上，当时间导数小于预设的很小的值
34、（精度要求）且不再随时间步数变化时，就认为达到了定常解；再随时间步数变化时，就认为达到了定常解；从图上看，经过大约从图上看，经过大约500个时间步，解已经基本稳定了。个时间步，解已经基本稳定了。44第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）2. 收敛史收敛史p在定常计算中，某参数随时间推进步数的变化称为收敛史，在定常计算中，某参数随时间推进步数的变化称为收敛史，它反映了收敛过程中的变化和收敛的快慢；它反映了收敛过程中的变化和收敛的快慢；p喉道处密度和速度的收敛史喉道处密度和速度的收敛史喉道处密度和速度喉道处密度和速度时间导数的变化时间导数的变化p
35、图中时间导数实际是图中时间导数实际是(7-46)和和(7-48)式右端的数值。式右端的数值。当解趋于收敛时，这些数当解趋于收敛时，这些数值并非严格等于零，称为值并非严格等于零，称为残差，其大小和衰减速度残差，其大小和衰减速度是评价时间相关法的重要是评价时间相关法的重要指标。指标。45第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）3. 定常解与精确解的比较与误差分析定常解与精确解的比较与误差分析p密度和马赫数的定常解与精确解的对比情况密度和马赫数的定常解与精确解的对比情况p两者吻合令人满意；两者吻合令人满意；p比较具体数据比较具体数据(表表7-4)，数
36、值误差约在，数值误差约在0.33.29%之间。产生原之间。产生原因主要是因主要是进口边界条件误差进口边界条件误差、截断误差截断误差、柯朗数影响柯朗数影响等。等。46第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p进口边界条件误差进口边界条件误差进口边界上流动有速度，并非滞止条进口边界上流动有速度，并非滞止条件件(对应于无限大截面积对应于无限大截面积)。如密度解析解为。如密度解析解为/0=0.995，有，有0.5%的误差；的误差；p截断误差截断误差由有限大小网格间距引起。由有限大小网格间距引起。为观察网格间距对截断误差的影响，需要比较加密网格后的结为观
37、察网格间距对截断误差的影响，需要比较加密网格后的结果果网格无关性检验网格无关性检验，如表，如表7-5所示。所示。项目项目喉道处计算结果喉道处计算结果*/0T*/T0p*/p0Ma网格点数网格点数:31个个0.6390.8360.5340.999网格点数网格点数:61个个0.6380.8350.5331.000解析解解析解0.6340.8330.5281.000网格数加倍网格数加倍误误差减小，但幅度差减小，但幅度不大，继续加密不大，继续加密减小幅度更小。减小幅度更小。可以认为可以认为31个网个网格点已足够，计格点已足够，计算结果基本与网算结果基本与网格无关。格无关。47第第7.2节节 亚声速亚声
38、速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）p柯朗数影响柯朗数影响如前所述，为提高计算精度柯朗数应尽可能如前所述，为提高计算精度柯朗数应尽可能接近接近1，上述计算柯朗数为，上述计算柯朗数为C=0.5；表表7-6给出了改变柯朗数的计算结果给出了改变柯朗数的计算结果柯朗数柯朗数*/0T*/T0p*/p0Ma0.50.6390.8360.5340.9990.70.6390.8370.5350.9990.90.6390.8370.5350.9991.00.6400.8370.5350.9991.10.6400.8370.5350.9991.2计算变得不稳定，最终发散计算变得不稳定，
39、最终发散解析解解析解0.6340.8330.5281.00048第第7.2节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（1）可见，将柯朗数提高到可见，将柯朗数提高到0.9和和1.1所得结果不见得比所得结果不见得比0.5的好，的好，0.5的结果反而更接近精确解。这说明实的结果反而更接近精确解。这说明实际计算时如何选取柯朗数是非常经验性的；际计算时如何选取柯朗数是非常经验性的； C=1.1尽管已经破坏了尽管已经破坏了CFL准则，但计算仍是稳定准则，但计算仍是稳定的，原因是的，原因是CFL准则是从线性方程得来的，对非准则是从线性方程得来的，对非线性问题只有参考和指导意义
40、。线性问题只有参考和指导意义。 49第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法解法（2）p第第7.1节给出了无量纲守恒形式控制方程组节给出了无量纲守恒形式控制方程组(7-101)式式第第7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟一、守恒形式控制方程组一、守恒形式控制方程组p本节使用本节使用MacCormack两步显式方法求解两步显式方法求解守恒守恒形式形式控制方程组，并与非守恒形式控制方程组的结果进控制方程组，并与非守恒形式控制方程组的结果进行对比。行对比。 用非定常方程求解用非定常方程求解定常问题定常问题时间相关法时间相关法时间推进时间推进5
41、0第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）pU是守恒变量，它与原始变量的关系为是守恒变量，它与原始变量的关系为(7-102)式式(7-104)式，式，即即51第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）p观察守恒控制方程组的通量和源项表达式观察守恒控制方程组的通量和源项表达式可以看出，通量和源项仍然是由原始变量、而不是守恒变量可以看出，通量和源项仍然是由原始变量、而不是守恒变量表示的，这说明就自变量表示的，这说明就自变量因变量函数关系而言，因变量函数关系而言，守恒控制守恒控制方程组的表达还不彻底方程
42、组的表达还不彻底，还需要经过中间变量，还需要经过中间变量即原始变即原始变量的转换，表现在数值求解上就是需要进行变量的反复转换量的转换，表现在数值求解上就是需要进行变量的反复转换计算。经验表明，控制方程组表达的不彻底性可能会破坏时计算。经验表明，控制方程组表达的不彻底性可能会破坏时间推进求解的稳定性。间推进求解的稳定性。52第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）p解决的方法是将守恒控制方程组的通量和源项直接表示成守解决的方法是将守恒控制方程组的通量和源项直接表示成守恒变量的形式，即恒变量的形式，即方程组、通量和源项都是守恒变量的形式，方程组、通
43、量和源项都是守恒变量的形式，不涉及原始变量，求解稳定。不涉及原始变量，求解稳定。53第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）p喷管的形状仍然采用喷管的形状仍然采用(7-73)式给出的面积分布形式，网格数与式给出的面积分布形式，网格数与第第7.2节相同，即共划分节相同，即共划分31个网格点；个网格点； p数值格式仍是数值格式仍是MacCormack显式两步格式算法；显式两步格式算法；p亚声速进口边界条件：亚声速进口边界条件：数目与第数目与第7.2节完全相同，给法为节完全相同，给法为二、守恒控制方程组的数值求解二、守恒控制方程组的数值求解给定给定1
44、：U(2)一阶外推：一阶外推：由由(7-103)式计算进口式计算进口V1，是可变的：，是可变的：第三个守恒变量第三个守恒变量U(3)54第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）p超声速出口边界条件：超声速出口边界条件：同样不需要物理边界条件，数值边界同样不需要物理边界条件，数值边界条件的给法均为一阶外推条件的给法均为一阶外推首先按原始变量分段给出密度和温度初场首先按原始变量分段给出密度和温度初场p初场布置：初场布置：经验表明，守恒控制方程组的求解对初场更敏感，经验表明，守恒控制方程组的求解对初场更敏感，需要采用更接近于真实情况的初场。需要采用更
45、接近于真实情况的初场。0x0.555第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）0.5x1.51.5x3.0按喷管质量流率布置初始速度场按喷管质量流率布置初始速度场有了原始变量的初场，根据守恒变量与原始变量的关系计有了原始变量的初场，根据守恒变量与原始变量的关系计算出守恒变量的初场；算出守恒变量的初场；p稳定性条件（稳定性条件（CFL准则）和时间步长准则）和时间步长的确定和前面一样。的确定和前面一样。接近定常值接近定常值0.57956第第7.3节节 亚声速亚声速超声速喷管等熵流动的超声速喷管等熵流动的CFD解法（解法（2）p对比守恒型和非守恒型控制
46、方程组的计算结果，可以得出如对比守恒型和非守恒型控制方程组的计算结果，可以得出如下结论下结论三、计算结果与分析三、计算结果与分析守恒型方程组给出的质量流率更好些；守恒型方程组给出的质量流率更好些； 非守恒型方程组给出的残差较小；非守恒型方程组给出的残差较小； 两者在计算精度方面没有大的差别：非守恒型得到的原始两者在计算精度方面没有大的差别：非守恒型得到的原始变量更精确，而守恒型得到的守恒变量更精确；变量更精确，而守恒型得到的守恒变量更精确； 守恒型方程组的计算量比非守恒型的大。守恒型方程组的计算量比非守恒型的大。p总体而言，在喷管亚声速总体而言，在喷管亚声速-超声速等熵流动计算中，守恒型方超声
47、速等熵流动计算中，守恒型方程组并没有表现出优于非守恒型方程组的特别之处。程组并没有表现出优于非守恒型方程组的特别之处。57第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p对于给定的截面积分布，拉瓦尔喷管全亚声速等熵流动取决对于给定的截面积分布，拉瓦尔喷管全亚声速等熵流动取决于进口总压与出口环境压强（反压）的差值；于进口总压与出口环境压强（反压）的差值；第第7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟一、解析解一、解析解p对于全亚声速喷管等熵流动，使用对于全亚声速喷管等熵流动，使用MacCormack两两步显式方法求解守恒形式控制方程组步显式方法求解守恒形
48、式控制方程组(7-46)式、式、(7-48)式和式和(7-50)式。式。 58第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p假设初始反压假设初始反压pa=p0无压差无压差无流动无流动；ppa开始流动开始流动；p当降低到当降低到pa=(pa)a、(pa)b时，喷管内均为时，喷管内均为亚声速流亚声速流Ma1，如，如图图7-13所示；所示；p当降低到当降低到pa=(pa)c时，时，喉道达临界状态：喉喉道达临界状态：喉道道(Ma)t=1、压强、压强(p)t=0.528p0；59第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p在亚声速情况下，喉道
49、没有达到临界状态，因而临界流动面在亚声速情况下，喉道没有达到临界状态，因而临界流动面积积A*At，同时出口压强始终与反压相平衡，即，同时出口压强始终与反压相平衡，即不是喷管中实际存在的截面不是喷管中实际存在的截面A*Atp求出求出A*后，对于每一个待求截面，首先由面积比公式后，对于每一个待求截面，首先由面积比公式(7-6) 解解出出Ma，再用，再用(7-7)式式 (7-9)式求出其它参数。式求出其它参数。出口出口Mae1临界流动面积临界流动面积A*60第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p喷管的无量纲面积分布喷管的无量纲面积分布二、数值求解二、数值求解以喉
50、道距进口距离为参考值的无量纲以喉道距进口距离为参考值的无量纲轴向坐标，喉道位于轴向坐标，喉道位于x=1.5处。处。p网格划分仍为网格划分仍为31个网格点；个网格点；p初场布置初场布置61第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p边界条件边界条件进口亚声速边界条件进口亚声速边界条件与亚声速与亚声速-超声速流时一样；超声速流时一样；出口亚声速边界条件出口亚声速边界条件物理边界条件：物理边界条件：1个，给定反压值个，给定反压值pa，则，则热状态方程热状态方程数值边界条件数值边界条件：2个个速度一阶外推速度一阶外推温度一阶外推温度一阶外推热状态方程热状态方程或者密度用
51、一阶外推，由热状态方程确定温度。或者密度用一阶外推，由热状态方程确定温度。62第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p取压力比取压力比pa/p0=0.93，柯朗数为，柯朗数为0.5，不同时间步的计算结果如，不同时间步的计算结果如图图7-17所示；所示；p可见数值解与解析解吻合良好。实际上经过大约可见数值解与解析解吻合良好。实际上经过大约2500个时间个时间步计算就已收敛了。步计算就已收敛了。三、计算结果与分析三、计算结果与分析63第第7.4节节 全亚声速喷管等熵流动的全亚声速喷管等熵流动的CFD解法解法p可能出现的问题：可能出现的问题：对较大压对较大压差，如
52、差，如pa/p0=0.90，则在接近，则在接近收敛时计算可能变得不稳定、收敛时计算可能变得不稳定、最终发散。最终发散。1200t800t400t0t00.61.21.82.43.00.70.80.91.0x/Lp/p0pa/p0=0.9现象的原因：在较大压差现象的原因：在较大压差作用下，具有较大强度的作用下，具有较大强度的压缩波和稀疏波传播到出压缩波和稀疏波传播到出口时，被恒压边界阻挡而口时，被恒压边界阻挡而反射，波在喷管中的来回反射，波在喷管中的来回运动导致出口流动振荡，运动导致出口流动振荡，经过足够长时间后计算发经过足够长时间后计算发散散由较大压差和恒压由较大压差和恒压边界共同造成的数值现
53、象，边界共同造成的数值现象，而非真实物理过程；而非真实物理过程；解决方法解决方法采用更接近精确解的初场，采用更接近精确解的初场，缩短非定常流动发展过程；缩短非定常流动发展过程；用人工黏性抑制。用人工黏性抑制。64第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法第第7章章 拟一维喷管流动的数值模拟拟一维喷管流动的数值模拟p使用使用MacCormack两步显式方法求解含激波的喷管两步显式方法求解含激波的喷管流动，守恒形式控制方程组为流动，守恒形式控制方程组为(7-101)式、式、(7-108)式、式、(7-110)式、式、(7-112)式和式和(7-115)式；式；p激波自动捕捉。激
54、波自动捕捉。65一、问题的提法一、问题的提法第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p如图如图7-13所示，当反压降所示，当反压降低到低到(pa)c时，喉道达到临时，喉道达到临界状态界状态Ma=1；p若进一步若进一步降低降低反压：反压：pa(pa)c，则喉道，则喉道Ma不会变不会变化化质量流率也不会质量流率也不会变化，喷管变化，喷管壅塞壅塞了；了；p但是，当但是，当pa1临界面积临界面积A1*=At总压总压p02p01，熵，熵s2马赫数马赫数Ma21临界面积临界面积A2*69第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p解析解的解析解的关键关键是利用连续
55、方程是利用连续方程质量流率为常数：质量流率为常数：绝热流，绝热流，T0=Const。应用于激波前后应用于激波前后反压反压pa70第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p激波后仍为等熵流动，所以激波后仍为等熵流动，所以两式两式 相乘相乘71第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p激波前后的总压比激波前后的总压比=0.6882面积比公式面积比公式激波所在截面面积激波所在截面面积面积分布公式面积分布公式激波的轴向位置激波的轴向位置72三、数值求解三、数值求解第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p使用使用MacCormack两步
56、显式格式计算。两步显式格式计算。p密度和温度初场分三段布置：密度和温度初场分三段布置：1. 初场布置初场布置进口截面至喉道进口截面至喉道仍用仍用(7-120)式前两段布置初场式前两段布置初场0x0.50.5x1.573第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法喉道截面至激波，采用下式布置初场喉道截面至激波，采用下式布置初场1.5x2.1激波至出口截面，采用下式布置初场激波至出口截面，采用下式布置初场2.1x3.0p速度初场仍用速度初场仍用(7-121)式布置：式布置：74第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p亚声速亚声速进口进口边界：按使用守恒型控
57、制方程组求解亚声速边界：按使用守恒型控制方程组求解亚声速-超声超声速等熵流时的进口条件；速等熵流时的进口条件；p亚声速亚声速出口出口边界：与求解全亚声速流时的出口边界条件基本边界：与求解全亚声速流时的出口边界条件基本相同：相同：2. 边界条件边界条件物理边界条件物理边界条件1个：给定出口个：给定出口pN=0.6784，由，由U(3)定义定义热状态方程热状态方程75第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法数值边界条件：数值边界条件：2个，由一阶外推计算个，由一阶外推计算76第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p自动捕捉激波时，为得到光滑稳定解需要加
58、入人工黏性。自动捕捉激波时，为得到光滑稳定解需要加入人工黏性。p采用采用(6-58)式给出的人工黏性形式，即式给出的人工黏性形式，即3. 人工黏性人工黏性p预估步计算预估步计算77第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p校正步计算校正步计算以上各式中以上各式中78四、计算结果与分析四、计算结果与分析第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p为了提高对激波的分辨率，划分为了提高对激波的分辨率，划分61个网格点。个网格点。p为对比，首先不添加人工黏为对比，首先不添加人工黏性。计算到性。计算到1600时间步时如时间步时如图所示。图所示。p尽管激波位置准确
59、，但残差尽管激波位置准确，但残差仍在仍在10-1量级，尚未收敛，量级，尚未收敛，且存在振荡；且存在振荡；p继续求解到继续求解到2800时间步，虽时间步，虽未发散，但振荡明显增大，未发散，但振荡明显增大，残差增大到残差增大到101量级，不能接量级，不能接受。受。1. 无人工黏性计算无人工黏性计算79第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p添加人工黏性后，取添加人工黏性后，取Cx=0.2，柯朗数取，柯朗数取0.5，计算到，计算到1400时间时间步时收敛到定常状态。步时收敛到定常状态。2. 有人工黏性计算有人工黏性计算l解析解解析解l解析解解析解p由图可见由图可见人工黏性显著
60、抑制了振人工黏性显著抑制了振荡，完全消除了波前的荡，完全消除了波前的振荡；振荡；人工黏性通过抹平激波人工黏性通过抹平激波间断（激波跨过若干个间断（激波跨过若干个网格）得到了光滑的解；网格）得到了光滑的解；若增大人工黏性，则波若增大人工黏性，则波后振荡也可以抑制，但后振荡也可以抑制，但会使激波变得更宽。会使激波变得更宽。80第第7.5节节 含激波喷管流动的含激波喷管流动的CFD解法解法p质量流率的计算结果质量流率的计算结果无人工黏性的结果振荡显著，不能接受；无人工黏性的结果振荡显著，不能接受；有人工黏性时，波前质量流率与解析解吻合良好；跨过激有人工黏性时，波前质量流率与解析解吻合良好；跨过激波时
61、有大的跳跃，在出口稳定在比进口高约波时有大的跳跃，在出口稳定在比进口高约8.6%的数值；的数值；人工黏性直接加在人工黏性直接加在守恒变量上，而守守恒变量上，而守恒变量恒变量U(2)正是质正是质量流率量流率人工黏人工黏性相当于质量源，性相当于质量源，这是导致质量流率这是导致质量流率增大的原因。增大的原因。p一般来说，用人工黏性捕捉激波时，得到的原始变量是可以一般来说，用人工黏性捕捉激波时，得到的原始变量是可以接受的。接受的。无人工黏性无人工黏性人工黏性人工黏性Cx=0.2l解析解解析解81课程结束了！课程结束了！飞行器流动仿真飞行器流动仿真82p平时成绩平时成绩30%到课情况：到课情况：10%课
62、堂笔记：完整性课堂笔记：完整性10% +正确性正确性10%关于课程考试的说明关于课程考试的说明p考试成绩考试成绩70%：面试：面试知识点：知识点：20%编程计算：编程计算：50%83p编程要求与评分标准编程要求与评分标准考试内容：使用考试内容：使用MacCormack两步显式方法，编写拟两步显式方法，编写拟一维喷管流动数值模拟程序，可计算三个典型工况：一维喷管流动数值模拟程序，可计算三个典型工况：亚声速亚声速超声速等熵流；超声速等熵流；全亚声速等熵；全亚声速等熵；含含激波流动激波流动 。程序语言：可使用任何计算机语言程序语言：可使用任何计算机语言评分标准：计算条件修改与调试评分标准：计算条件修
63、改与调试10%+计算结果正确计算结果正确30%+正确处理数据正确处理数据10%。84程序编写举例程序编写举例p使用使用Fortran语言，编写语言，编写MacCormack两步显式方法的程序，两步显式方法的程序，数值求解拟一维喷管流动。数值求解拟一维喷管流动。 Program nozzleC 拟一维喷管计算，2013年7月25日修改C icase-计算工况：1-亚声速-超声速等熵流（非守恒型方程组）C 2-亚声速-超声速等熵流（守恒型方程组）C 3-全亚声速等熵流C 4-含激波流（守恒方程组）C nstop-最大时间步数；n-网格点数；deltx-轴向网格间距；cfl-CFL数；C dt-时间
64、步长；gamma=1.4-比热比；pe-反压；p0-进口总压C x()-轴向坐标；a()-面积；am()-马赫数；fm()-质量流率 C 原始变量：r()-密度；t()-温度；u()-速度；p()-压强C 守恒变量：U1()；U2()；U3()C 预测步计算的中间值：rm()-密度；um()-速度；tm()-温度；pm()-压强 C 一阶时间导数：drdtm-密度；dudtm-速度；dtdtm-温度C 校正步：一阶时间导数平均值drdtav-密度；dudtav-速度；dtdtav-温度 C 一阶时间导数drdt-密度；dudt-速度；dtdt-温度 C it-喉部网格点编号85程序编写举例程序
65、编写举例 parameter (n=31) dimension x(n),a(n),r(n),t(n),u(n),p(n),am(n),fm(n),rm(n),um(n),tm(n),pm(n) dimension U1(n),U2(n),U3(n),U1m(n),U2m(n),U3m(n)C常数计算 gamma=1.4 g0=gamma g1=gamma-1.0 g2=1.0/(gamma-1.0) g3=0.5*gamma g4=0.5*gamma*(gamma-1.0) g5=0.5*(3.0-gamma) g6=(gamma-1.0)/gamma g7=1.0/gammaC设置计算终止
66、条件 write(*,*) 输入最大时间步数： read (*,*) nstop86程序编写举例程序编写举例1000 write(*,*) 输入工况代号： write(*,*) 1=亚声速-超声速（非守恒方程）； write(*,*) 2=亚声速-超声速（守恒方程）； write(*,*) 3=全亚声速（非守恒方程）； write(*,*) 4=含激波流（守恒方程）。 read (*,*) icase if (icase.ge.1.and.icase.le.4) goto 1001 write(*,*)工况代号无效, 重新输入! goto 1000C布置网格1001 deltx=3.0/rea
67、l(n-1) it=1 x(1)=0.0 do 50 i=2,n x(i)=x(i-1)+deltx if (abs(x(i)-1.5).le.1.0e-5) it=I !寻找喉道网格编号50 continue write(*,*) 喉道网格编号it=,it87程序编写举例程序编写举例C布置初场 Call initial(icase,x,a,r,t,u,U1,U2,U3,gamma,pe) if (icase.eq.2.or.icase.eq.4) goto 2000 C非守恒型方程组计算 if (icase.eq.1) then ！打开数据文件，准备写入 open(12,file=非守恒方程
68、-亚超-喉道参数收敛史.dat) else open(12,file=非守恒方程-全亚-喉道参数收敛史.dat) end ifC喉部数据输出，排列顺序：时间步，密度，温度，压强，马赫数 write(12,22) 0,r(it),t(it),r(it)*t(it),u(it)/sqrt(t(it) do 500 j=1,nstop ！循环计算 CFL=0.5 dt=1.0 Do 110 i=1,n ！确定时间步长 as=sqrt(t(i) deltt=cfl*deltx/(as+u(i) if (dt.gt.deltt) dt=deltt ！找出最小的时间步长110 Continue88程序编写
69、举例程序编写举例 Do 300 i=2,n-1 ！MacCormack显式格式 drdt=(-r(i)*(u(i+1)-u(i)-u(i)*(r(i+1)-r(i)-r(i)*u(i)* + (log(a(i+1)-log(a(i)/deltx dudt=(-u(i)*(u(i+1)-u(i)-g7*(t(i+1)-t(i)+t(i)/r(i)*(r(i+1)-r(i)/deltx dtdt=(-u(i)*(t(i+1)-t(i)-g1*t(i)*(u(i+1)-u(i)+ + u(i)*(log(a(i+1)-log(a(i)/deltx rm(1)=r(1) ！计算预估步中间值 um(1)
70、=u(1) tm(1)=t(1) rm(i)=r(i)+drdt*dt um(i)=u(i)+dudt*dt tm(i)=t(i)+dtdt*dt drdtm=(-rm(i)*(um(i)-um(i-1)-um(i)*(rm(i)-rm(i-1) + -rm(i)*um(i)*(log(a(i)-log(a(i-1)/deltx dudtm=(-um(i)*(um(i)-um(i-1)-g7*(tm(i)-tm(i-1)+ + tm(i)/rm(i)*(rm(i)-rm(i-1)/deltx dtdtm=(-um(i)*(tm(i)-tm(i-1)-g1*tm(i)* + (um(i)-um(
71、i-1)+um(i)*(log(a(i)-log(a(i-1)/deltx 89程序编写举例程序编写举例 drdtav=(drdt+drdtm)/2.0 ！计算时间导数的算术平均值 dudtav=(dudt+dudtm)/2.0 dtdtav=(dtdt+dtdtm)/2.0 r(i)=r(i)+drdtav*dt ！计算校正步变量值 u(i)=u(i)+dudtav*dt t(i)=t(i)+dtdtav*dt p(i)=r(i)*t(i) am(i)=u(i)/sqrt(t(i) fm(i)=r(i)*a(i)*u(i) if (i.eq.it) write(12,22) j,r(i),t
72、(i),p(i),am(i)300 continueC进口边界条件i=1 u(1)=2.*u(2)-u(3) ！一阶外推 r(1)=1.0 ！物理边界条件 t(1)=1.0 p(1)=r(1)*t(1) am(1)=u(1)/sqrt(t(1) fm(1)=r(1)*a(1)*u(1)90程序编写举例程序编写举例C出口边界条件i=n if (icase.eq.1) then u(n)=2.0*u(n-1)-u(n-2) ！一阶外推 r(n)=2.0*r(n-1)-r(n-2) t(n)=2.0*t(n-1)-t(n-2) p(n)=r(n)*t(n) am(n)=u(n)/sqrt(t(n)
73、fm(n)=r(n)*a(n)*u(n) else if (icase.eq.3) then p(n)=pe ！物理边界条件 u(n)=2.0*u(n-1)-u(n-2) ！一阶外推 t(n)=2.0*t(n-1)-t(n-2) r(n)=p(n)/t(n) am(n)=u(n)/sqrt(t(n) fm(n)=r(n)*a(n)*u(n) end if 500 continue close(12) ！关闭数据文件91程序编写举例程序编写举例CC输出计算结果C if (icase.eq.1) then ！打开数据文件，准备写入 open(11,file=非守恒方程-亚超-参数分布.dat) e
74、lse open(11,file=非守恒方程-全亚-参数分布.dat) end ifC数据排列顺序：x，面积，密度，速度，温度，压强，马赫数，质量流率 do 901 i=1,n write(11,11) x(i),a(i),r(i),u(i),t(i),p(i),am(i),fm(i)901 continue close(11) ！关闭数据文件 stop92程序编写举例程序编写举例C守恒型方程组计算2000 if (icase.eq.2) then ！打开数据文件，准备写入 open(12,file=守恒方程-亚超-喉道参数收敛史.dat) else open(12,file=守恒方程-激波-
75、喉道参数收敛史.dat) end ifC喉部数据输出，排列顺序：时间步，密度，温度，压强，马赫数 write(12,22) 0,r(it),t(it),r(it)*t(it),u(it)/sqrt(t(it) CFL=0.5 Cx=0.2 ！人工黏性系数 do 600 j=1,nstop ！循环计算 dt=1.0 do 120 i=1,n ！确定时间步长 p(i)=r(i)*t(i) as=sqrt(t(i) deltt=cfl*deltx/(as+u(i) if (dt.gt.deltt) dt=deltt ！寻找最小时间步长120 continue93程序编写举例程序编写举例 do 400
76、 i=2,n-1 ！MacCormack显式格式 dU1dt=-(U2(i+1)-U2(i)/deltx cJ2=g7*p(i)*(a(i+1)-a(i)/deltx dU2dt=-(g5*(U2(i+1)*2/U1(i+1)-U2(i)*2/U1(i)+ + g6*(U3(i+1)-U3(i)/deltx+cJ2 dU3dt=-(g0*(U2(i+1)*U3(i+1)/U1(i+1)-U2(i)*U3(i)/U1(i)- + g4*(U2(i+1)*3/U1(i+1)*2-U2(i)*3/U1(i)*2)/deltxC if (icase.eq.2) then U1m(i)=U1(i)+dU
77、1dt*dt ！计算预估步中间值 U2m(i)=U2(i)+dU2dt*dt U3m(i)=U3(i)+dU3dt*dt else ！计算人工黏性 dU1=U1(i+1)-2.0*U1(i)+U1(i-1) dU2=U2(i+1)-2.0*U2(i)+U2(i-1) dU3=U3(i+1)-2.0*U3(i)+U3(i-1) dp=abs(p(i+1)-2.0*p(i)+p(i-1)/(p(i+1)+2.0*p(i)+p(i-1)94程序编写举例程序编写举例 s1=Cx*dp*dU1 s2=Cx*dp*dU2 s3=Cx*dp*dU3 U1m(i)=U1(i)+dU1dt*dt+s1 ！计算预
78、估步中间值（含人工黏性） U2m(i)=U2(i)+dU2dt*dt+s2 U3m(i)=U3(i)+dU3dt*dt+s3 end if rm(i)=U1m(i)/a(i) um(i)=U2m(i)/U1m(i) tm(i)=g1*U3m(i)/U1m(i)-g4*um(i)*2 pm(i)=rm(i)*tm(i)C U1m(1)=U1(1) U2m(1)=U2(1) U3m(1)=U3(1) pm(1)=p(1)95程序编写举例程序编写举例 U1m(n)=U1(n) U2m(n)=U2(n) U3m(n)=U3(n) pm(n)=p(n)C dU1dtm=-(U2m(i)-U2m(i-1)
79、/deltx cJ2=g7*pm(i)*(a(i)-a(i-1)/deltx dU2dtm=-(g5*(U2m(i)*2/U1m(i)-U2m(i-1)*2/U1m(i-1)+ + g6*(U3m(i)-U3m(i-1)/deltx+cJ2 dU3dtm=-(g0*(U2m(i)*U3m(i)/U1m(i)-U2m(i-1)*U3m(i-1)/U1m(i-1) + -g4*(U2m(i)*3/U1m(i)*2-U2m(i-1)*3/U1m(i-1)*2)/deltxC dU1dtav=(dU1dt+dU1dtm)/2.0 ！计算时间导数的算术平均值 dU2dtav=(dU2dt+dU2dtm)
80、/2.0 dU3dtav=(dU3dt+dU3dtm)/2.096程序编写举例程序编写举例 if (icase.eq.2) then U1(i)=U1(i)+dU1dtav*dt ！计算校正步变量值 U2(i)=U2(i)+dU2dtav*dt U3(i)=U3(i)+dU3dtav*dt else dU1m=U1m(i+1)-2.0*U1m(i)+U1m(i-1) dU2m=U2m(i+1)-2.0*U2m(i)+U2m(i-1) dU3m=U3m(i+1)-2.0*U3m(i)+U3m(i-1) dpm=abs(pm(i+1)-2.0*pm(i)+pm(i-1)/(pm(i+1)+2.0*
81、pm(i)+pm(i-1) s1m=Cx*dpm*dU1m s2m=Cx*dpm*dU2m s3m=Cx*dpm*dU3mC U1(i)=U1(i)+dU1dtav*dt+s1m ！计算校正步变量值（含人工黏性） U2(i)=U2(i)+dU2dtav*dt+s2m U3(i)=U3(i)+dU3dtav*dt+s3m end if97程序编写举例程序编写举例 r(i)=U1(i)/a(i) u(i)=U2(i)/U1(i) t(i)=g1*U3(i)/U1(i)-g4*u(i)*2 p(i)=r(i)*t(i) am(i)=u(i)/sqrt(t(i) fm(i)=r(i)*a(i)*u(i
82、) if (i.eq.it) write(12,22) j,r(i),t(i),p(i),am(i)400 continueC进口边界条件i=1 r(1)=1.0 ！物理边界条件 t(1)=1.0 U1(1)=r(1)*a(1) U2(1)=2.0*U2(2)-U2(3) ！一阶外推 u(1)=U2(1)/U1(1) U3(1)=U1(1)*(g2*t(1)+g3*u(1)*2) p(1)=r(1)*t(1) am(1)=u(1)/sqrt(t(1) fm(1)=r(1)*a(1)*u(1)98程序编写举例程序编写举例C出口边界条件i=n U1(n)=2.0*U1(n-1)-U1(n-2) ！
83、一阶外推 U2(n)=2.0*U2(n-1)-U2(n-2) r(n)=U1(n)/a(n) u(n)=U2(n)/U1(n) if (icase.eq.2) then U3(n)=2.0*U3(n-1)-U3(n-2) t(n)=g1*U3(n)/U1(n)-g4*u(n)*2 p(n)=r(n)*t(n) else if (icase.eq.4) then p(n)=pe ！物理边界条件 U3(n)=g2*p(n)*a(n)+g3*U2(n)*u(n) t(n)=g1*U3(n)/U1(n)-g4*u(n)*2 end if am(n)=u(n)/sqrt(t(n) fm(n)=r(n)*
84、a(n)*u(n)600 continue close(12)99程序编写举例程序编写举例C数据输出 if (icase.eq.2) then ！打开数据文件，准备写入 open(11,file=守恒方程-亚超-参数分布.dat) else open(11,file=守恒方程-激波-参数分布.dat) end ifC顺序：x,面积,密度,速度,温度,压强,马赫数,质量流率,三个守恒变量 do 902 i=1,n write(11,33) x(i),a(i),r(i),u(i),t(i),p(i),am(i),fm(i),U1(i),U2(i),U3(i)902 continue close(1
85、1) ！关闭数据文件C11 format(f5.2,7f9.5) ！数据输出格式22 format(i4,4f9.5)33 format(f5.2,10f9.5) stop end100程序编写举例程序编写举例C初始化子程序 Subroutine initial(icase,x,a,r,t,u,U1,U2,U3,ga,pe) parameter (n=31) dimension x(n),a(n),r(n),t(n),u(n),U1(n),U2(n),U3(n)C if (icase.eq.1) then do 100 i=1,n a(i)=1.0+2.2*(x(i)-1.5)*2 r(i)=
86、1.0-0.3146*x(i) t(i)=1.0-0.2314*x(i) u(i)=(0.1+1.09*x(i)*sqrt(t(i)100 continue else if (icase.eq.2) then do 200 i=1,n a(i)=1.0+2.2*(x(i)-1.5)*2 if (x(i).ge.0.0.and.x(i).le.0.5) then r(i)=1.0 t(i)=1.0101程序编写举例程序编写举例 else if (x(i).gt.0.5.and.x(i).le.1.5) then r(i)=1.0-0.366*(x(i)-0.5) t(i)=1.0-0.167*(
87、x(i)-0.5) else r(i)=0.634-0.3879*(x(i)-1.5) t(i)=0.833-0.3507*(x(i)-1.5) end if u(i)=0.59/(r(i)*a(i)C U1(i)=r(i)*a(i) U2(i)=U1(i)*u(i) U3(i)=r(i)*a(i)*(t(i)/(ga-1.0)+0.5*ga*u(i)*2)C200 continueelse if (icase.eq.3) then write(*,*) 输入反压值pa，数值范围为0.528至1.0： read(*,*) pe102程序编写举例程序编写举例 do 210 i=1,n if (x
88、(i).le.1.5) then a(i)=1.0+2.2*(x(i)-1.5)*2 elsea(i)=1.0+0.2223*(x(i)-1.5)*2 end if r(i)=1.0-0.023*x(i) t(i)=1.0-0.009333*x(i) u(i)=0.05+0.11*x(i)210 continueelse if (icase.eq.4) then pe=0.6784 write(*,*) 反压值已设定为：pa=,pe do 220 i=1,n a(i)=1.0+2.2*(x(i)-1.5)*2 if (x(i).ge.0.0.and.x(i).le.0.5) then r(i)
89、=1.0 t(i)=1.0103程序编写举例程序编写举例 else if (x(i).gt.0.5.and.x(i).le.1.5) then r(i)=1.0-0.366*(x(i)-0.5) t(i)=1.0-0.167*(x(i)-0.5) else if (x(i).gt.1.5.and.x(i).le.2.1) then r(i)=0.634-0.702*(x(i)-1.5) t(i)=0.833-0.4908*(x(i)-1.5) else r(i)=0.5892-0.10228*(x(i)-2.1) t(i)=0.93968-0.0622*(x(i)-2.1) end if u(i)=0.59/(r(i)*a(i) U1(i)=r(i)*a(i) U2(i)=U1(i)*u(i) U3(i)=U1(i)*(t(i)/(ga-1.0)+0.5*ga*u(i)*2)220 continue else end if return end104
