可压缩边界层的抛物化稳定性方程 及其在某些问题中的应用张永明 2008年10月摘要推导出可压缩边界层的抛物化稳定性方程（PSE），验证了PSE方法的可靠性，再将其应用于非平行效应、二次失稳和转捩预测三个具体问题中。发现了边界层非平行性对中性曲线的影响主要集中在临界雷诺数处；验证了二次失稳机制的存在，得到了二次失稳问题中三维亚谐波的放大率随其展向波数和二维基本波幅值的变化关系；得到了准确的转捩起始位置，并再现了转捩中breakdown过程的机理。PSE方法是研究有关边界层内扰动演化的多种问题的有效工具。1 引言PSE是研究基本流中扰动演化的比较新的一种方法，此外常用的方法还有线性稳定定性理论（LST）和直接数值模拟（DNS）。与另外二者相比，PSE的优点在于同时具备以下特点：可以考虑基本流的非平行性，可以用于非线性问题，计算量不大。因此，从提出PSE至今不长的时间里，PSE方法得到了充分的发展，并应用到多个问题的研究中。Herbert和Bertolotti（1987）首先提出了抛物化稳定性方程。他们研究的是不可压边界层的线性PSE，这是一种比较简单的情况。Bertolotti等（1992）用线性PSE在大雷诺数处得到的小扰动演化与线性稳定性理论和直接数值模拟的结果一致，验证了线性PSE计算的可靠性，如图1-1所示。Bertolotti等（1992）用PSE研究不可压边界层的非平行性对其稳定性的影响，得到的中性曲线比用LST得到的更符合Schubauer和Skramstad（1943，1947）、Ross等（1970）、Strazisar等（1977）和Kachanov等（1977）的实验结果，也与Gaster（1974）、于秀阳和周恒（1986）用另外一种方法得到的计算结果一致，如图1-2所示。Bertolotti等（1992）提出了不可压边界层的非线性PSE，所得扰动演化结果与直接数值模拟所得一致，如图1-3所示，用计算的方法验证了非线性PSE的可靠性。而前者计算花费的时间比后者少很多，这显示出PSE在研究某些非线性问题时的优势：准确且快速。Joslin等（1992）和Esfahanian等（2001）也做了同样的验证。Malik等（1999）将非线性PSE计算所得与实验测量结果对比，发现二者一致，又从实验上验证了其可靠性。在此基础之上，Malik等（1999）开始用PSE作为可靠的工具来研究不可压边界层的稳定性和转捩问题。Herbert（1997）对不可压缩流PSE的研究进行了总结。图1-1 不同方法得到的大雷诺数处小扰动增长率。实线为Bertolotti等（1992）的线性PSE结果，虚线为LST结果，圆形标记为DNS结果。图1-2 不同方法得到的中性曲线。实线为Bertolotti等（1992）的线性PSE结果，虚线为LST结果，圆形标记为Gaster（1974）的结果。图1-3 Bertolotti等（1992）的非线性PSE计算与空间模式DNS的对比，包括平均流修正，一阶扰动，二阶扰动和三阶扰动。实线为NPSE结果，方形标记为DNS结果。相对不可压边界层，可压缩边界层的PSE更为复杂，其结果也较少。但由于航空航天技术发展的需求，这方面的研究逐渐多了起来。Bertolotti和Herbert（1991），Bertolotti（1991）提出了可压缩平板边界层的线性PSE，并发现用线性PSE在大雷诺数处得到的小扰动演化与线性稳定性理论的结果一致，验证了线性PSE计算的可靠性。可压缩边界层的非线性PSE研究工作比较少，目前见到的工作主要来自美国国家航空航天局（NASA）Langley研究中心的Chang等人。Chang等（1991）提出了可压缩边界层的非线性PSE。为了验证其可靠性，Chang等（1991，1993）、Pruett等（1995）、Pruett和Chang（1995）、Jiang等（2003，2004）分别将非线性PSE结果与DNS的对比，发现二者一致，但PSE计算花费时间少得多。遗憾的是，上述关于可压缩边界层的线性和非线性PSE的验证工作都不是足够充分、完整和严格的，其原因将在本文第三部分做具体阐述。本文以可压缩平板边界层为研究对象。先从完全的NavierStokes方程出发，推导出线性和非线性的抛物化稳定性方程（见第二部分），再介绍具体的数值计算方法（见第二部分）。然后对线性PSE和非线性PSE计算方法做充分、完整和严格的检验，验证其可靠性（见第三部分）。此后则以PSE作为工具研究具体问题。应用线性PSE研究基本流的非平行性对可压缩边界层中性曲线的影响（见第四部分）。应用非线性PSE研究超音速边界层的二次失稳问题（见第五部分）。应用非线性PSE研究转捩可压缩边界层的转捩预测问题（见第六部分）。2 可压缩边界层PSE的控制方程和数值求解方法这一部分先介绍可压缩边界层的抛物化稳定性方程的推导过程，以及求解方程所用的数值方法。这些是前人做过的工作，但都是PSE方法的基础，所以有必要做详细的介绍。首先推导出线性和非线性PSE的控制方程（见2.1）；然后说明数值计算中所用的基本流和入口条件（见2.2）；最后阐述计算网格和差分格式（见2.3）。2.1 控制方程推导抛物化稳定性方程需经四个步骤：（1）先由有量纲的完全的N-S方程出发，利用特征量得到无量纲的方程，（2）再推导出完全的扰动方程，（3）然后针对扰动的特点推出稳定性方程，（4）最后利用边界层厚度增长缓慢的性质将其抛物化而得到抛物化稳定性方程。下面作具体的介绍。用上标表示有量纲的量，在三维笛卡尔坐标系下用x，y，z表示流向、法向和展向，则有量纲的完全的N-S方程为 （2.1-1）其中是速度矢量，是时间，是密度，是压力，是温度，是定压比热，是热传导系数，是第一动力粘性系数，是第二动力粘性系数，为普适气体常数。为粘性耗散函数，表示为 ， （2.1-2）其中为应变率张量，表达式为。 （2.1-3）方程（2.1-1）中第一到第四个方程分别为连续性方程、动量方程、内能方程和状态方程。本文假设可压缩流体为完全气体，使用完全气体的状态方程。下面要将N-S方程无量纲化，其中用到的特征量如下。取特征长度为， （2.1-4）其中为平板前缘到计算域入口的距离，为运动粘性系数，运动粘性系数与第一动力粘性系数关系为，为流向速度，下标e表示边界层外的自由来流中的量。其它的原始特征量为：自由流速度、自由流温度、自由流密度，自由流第一动力粘性系数，自由流第二动力粘性系数，自由流热传导系数。导出的特征量为：特征压力为，特征时间为。用有量纲量除以特征量即为无量纲量，如。（2.1-5）将上式代入有量纲NS方程（2.1-1），即得无量纲NS方程。无量纲N-S方程会用到几个无量纲参数，即入口出雷诺数、普朗特数和马赫数，下面分别介绍。平板前缘到入口处的无量纲距离为，（1）而入口处的雷诺数的表达式则为。（2由一直接化简即可）平板前缘到入口下游某处的距离为，显见。由此得到各变量对x的一阶、二阶导数的量级分别为。 （2.1-6）这将是后来对方程进行抛物化的依据。普朗特数表达式为。马赫数表达式为，其中为自由来流中的声速，其表达式为，其中为比热比，其表达式为，其中为定容比热。在整理无量纲N-S方程时，还要注意粘性系数和热传导系数的处理。本文使用了Stokes假设（见流体力学上册，68页）（周光炯2000），第二动力粘性系数可用第一动力粘性系数表示。 （2.1-7）而无量纲的第一动力粘性系数则使用Sutherland公式，假设其随温度的变化满足如下关系式， （2.1-8）其中C1110.4K为常数。本文中无量纲的热传导系数也使用Sutherland公式，形式为。 （2.1-9）一般来讲可以取C2C1，则无量纲的热传导系数与第一动力粘性系数相同， （2.1-10）以后的公式中可以都用粘性系数表示。整理之后的无量纲N-S方程形式如下(此种写法挺好) （a）（b）（c）（d）（e） （f）（2.1-11）其中（a）为连续性方程，（b）、（c）和（d）分别为流向、法向和展向的动量方程，（e）为内能方程，（f）为状态方程。将（2.1-11）中的状态方程代入连续性方程、动量方程和内能方程，消去压力p，则无量纲方程化简为5个方程，其中的变量也是5个，即： （a）（b）（c）（d）（e）（2.1-12）方程中有一点需要注意，无量纲长度x、y、z、时间t、入口雷诺数Re0等用到特征长度的量，都与常见的直接数值模拟中的对应量不同。边界层的直接数值模拟常用入口处位移排移厚度作为特征长度，（常用有两种，这是其中一种）其表达式为， （2.1-13）其中为马赫数M的函数，其表达式较复杂，详见文献（黄章峰2003第17页，黄章峰2006第26页）。对比式（2.1-4）可见，与本文所用特征长度相差倍。那么无量纲量x、y、z、t、Re0等也相差倍。本文后面还会提到一些用特征长度无量纲化的量，也会与DNS中的对应量相差倍，这一点以后要注意。将方程（2.1-12）中的瞬时量设为定常基本流与扰动量之和 （2.1-14）其中上标表示基本流，上角标表示扰动量。将（2.1-14）式代入NS方程（2.1-12），并减去定常基本流满足的方程，得到扰动方程（方便起见，较复杂的方程用矩阵形式表达） （2.1-15）其中为扰动矢量。系数矩阵Vxx，Vyy，Vzz，Vxy，Vxz，Vyz，A，B，C，D都可分解为线性部分与非线性部分，其中线性部分只包含定常基本流的量，非线性部分则包含扰动量。二者分别用上标l和n表示，则这些矩阵可记为。通过量级分析可以发现，对流项和粘性项的量级分别为1和1/Re0。系数矩阵，A，B，C，D同时包含对流项和粘性项，其总的量级为1，系数矩阵Vxx，Vyy，Vzz，Vxy，Vxz，Vyz只包含粘性项，其量级为1/Re0。将线性和非线性部分分别放在方程两边，则扰动方程（2.1-15）整理为如下形式 （2.1-16）其中右边的矢量函数表示所有的非线性项，其表达式如下（2.1-17）非线性扰动方程（2.1-16）的系数矩阵和非线性项矢量的元素详见附录1。其中需要注意的是的处理。本文利用Sutherland公式（2.1-8）将表示为和的函数，下面介绍详细的推导过程。由Sutherland公式可知是的函数， （2.1-18）相对来说是很小的量，可用的微分近似， （2.1-19）考虑到，代入上式，则得到的近似表达式， （2.1-20）或简写为， （2.1-21）其中，是基本流温度的函数，即。 （2.1-22）（第一部分工作是方程无量纲化，变为关于五个未知量的五个方程，注意特征长度的选择和以前略有不同，再就是将瞬时量设为定常基