基本流动方程.docx

基本流动方程本项目研究需要计算进气道内外流耦合的复杂三维粘性流动，采用的流场控制方程为三维的Navier-Stokes方程组，它是连续方程、动量方程、能量方程的联立方程组本章给出了积分形式和微分形式的N-S方程组，以及在有限体积法离散中需使用的坐标转换后N-S方程组，另外还给出了湍流计算使用的平均化后的湍流N-S方程1 积分形式N-S方程组在直角坐标系下，忽略重力做功和辐射传热的积分形式N-S方程组可写为如下的矢量形式：质量守恒方程 （2-1）动量守恒方程 （2-2）能量守恒方程 （2-3）这是一个关于时间的双曲型方程组，其中E代表总比能（内能和动能之和）；为能量流矢量（energy flux vector），这里假定能量流矢量仅仅表达分子能量的运输，它可以用Fourier定律来描述，即 （2-4）K为热传导系数，根据Prandtl数：为常数，及Sutherland公式：来确定，其中和 是与流体有关的常数，详见表2-1。

对于理想气体，总比能可表述为下述张量缩并的形式：并且总焓为：这里ｅ代表内能，Ｈ为总焓，γ为比热比方程（2-2）中是粘性应力张量（viscous stress tensor），对于牛顿流体，假定应力张量随着变形率张量连续变化（Stokes假设） （2-5）变形率张量表为 为粘性系数，层流的粘性系数与压力的关系可以忽略不计，于是可采用Sutherland公式 ，其中和 是与流体有关的常数，详见表2-1表2-1 气体的粘性系数和热传导系数气体种类 （N·S/㎡） （K） （W/m·K） （K）R （J/㎏·K）空气1.458700×10-6110.5552.5010830×10-3194.4287.07氮气N21.399312×10-6106.6672.3611054×10-3166.7296.80氧气O21.751720×10-6138.8892.6954795×10-3222.2259.902 微分形式N-S方程组在直角坐标系下，忽略重力做功和辐射传热的微分形式N-S方程组可以按守恒型变量写为如下的矢量形式： （2-6）其中，u是守恒变量，f、g、h 是对流通量，a、b、c 粘性通量，它们的表达式为： （2-7）粘性通量中，应力和热流量的表达式为： （2-8）上述公式中，t是时间，x， y， z是空间坐标，是密度，是三个坐标方向的速度分量，E 是单位质量气体的总能，p是压强，T是温度，为粘性系数，K为热传导系数。

对于量热完全气体，可给出下述形式的气体状态方程： （2-9）3 方程坐标变换N-S方程组可以从笛卡尔坐标系转换到任意曲线坐标系下，其中 （2-10）这里介绍的坐标变换是由Viviand[1]和Vinokur[2]发展的，变换的选取是如此进行，以便使曲线空间网格间隔是均匀的，且为单位长度这将导致计算域ξ、η 和ζ为矩形域且为均匀的规则网格所以在数值计算式中标准的非加权差分格式可以应用原笛卡尔空间将被称为物理域，在物理域的点与计算域中的点将形成一一对应的关系（奇异点例外）在奇异点区域，通常是一个物理点对应许多个计算点（这通常出现在计算边界上）按式（2-10）进行变换，对于某变量，有： 令上式的变换矩阵为J，则其行列式称为雅可比行列式，记作： 则有 变换后得守恒型N-S方程： （2-11）其中： 记 其中： 取： 上述方程中矩阵系数含有，，，，，，，，及J等参数。

在采用有限体积法离散时，变换后的各个等ξ、η、ζ面将计算区域划分为许多六面体（ 或四边形 ）的单元，即计算网格变换行列式J 就是坐标变换后网格单元的体积伸缩比，即：如设△ξ=1、△η=1、△ζ=1，则1/J 就是网格单元的体积令， 并将，，，，，，，，等换为网格分界面面积矢量（方向沿界面正法向，大小为面积值）各个分量来表示：， ， ， ， ， ， 则方程（2-11）可重新记为下述形式： （2-12）其中， 其中， 其中： 可取： 4 几种形式的Euler方程众所周知在N-S方程的求解中，困难最大的是对流项的差分解法，很多研究者努力寻求的高精度高分辨率格式也是为解决对流项模拟的问题对此所开展的研究首先都是针对只还有非定常项和对流项的Euler方程。

在本文的算法研究中需要用到当地线性化的守恒型、原始变量型、特征变量型的Euler方程，因此这一节给出了这几种形式的Euler方程4.1守恒形式的Euler方程将N-S方程（2-12）中的粘性项A、B、C去掉，就得到了无粘流动的Euler方程通量F、G、H是守恒变量u的非线性函数，假设F、G、H在空间某点（ξ0，η0，ζ0）的某个邻域近似线性分布，则方程（2-12）可当地线性化为： （2-13）式中，矩阵Ac、Bc、Cc是方程的雅可比系数矩阵，， ， 以k表示某个坐标方向（ξ、η、ζ中的任意一个），令：，， 则可以用一个形式统一的系数矩阵Kc来表示上述三个矩阵：（2-14）其中，4.2原始变量型Euler方程在算法研究中，将会用到以求解原始变量Up的方程，称为原始变量型方程 （2-15）设有：，则方程（ 2-13 ）可写为 （2-16）将式（ 2-16 ）左乘以M即可得： （2-17）若令：， ， 则得到了当地线性化的任意曲线坐标系下原始变量型Euler方程： （2-18）与Ac、Bc、Cc类似，系数Ap、Bp、Cp矩阵也可以用统一的形式Kp来表示。

M、M –1及Kp的表达式为： （2-19）式中a是音速，可以看出矩阵Kp与矩阵Kc相似，所以它们有相同的特征值而矩阵Kp具有较多的零元素，求解它的特征值要容易得多由式（2-19）很容易求得矩阵Kp的特征值： （2-20）这些特征值的量纲都是速度乘以，但并不影响依据特征波进行的分析，所以有时又把这些特征值都叫做特征波速4.3特征变量型Euler方程在求得特征值之后，可以将矩阵Kp进行对角化： （2-21）导出矩阵Pk并不困难，本文在推导中注意到：Kp有三个特征值是相同的，所以矩阵Pk的形式不唯一如果采用了不恰当的形式，则在网格界面在某些特殊位置时，矩阵Pk是不可逆的，式（2-21）将不再成立，因此推导建立算法时需要特别小心，本研究给出了在任意三维情况下均适用的特征矢量矩阵Pk及其逆矩阵Pk-1（详见第三章，节）由特征矢量矩阵可以导出特征变量形式的Euler方程。

在方程（2-18）中应用式（2-21），并且对方程左乘以Pk-1，得：（2-22）在这里，k指的是ξ、η、ζ三个坐标方向中任意一个，m、n是其余两个坐标方向合并上式中的第三、四两项，令：则，方程（2-22）可改写为：定义矢量Uch，使得：可以在网格k向，得到特征变量Uch的特征变量型Euler方程： （2-23）对于一维流动，方程中第三项不存在，此时方程成为对角化形式，组成上式的各分量方程互相解耦，这对于各种数值处理来说都是非常方便的当流动是多维时，方程第三项一般不为零，只有在满足当地一维流动条件：方程（2-23）才能成为对角化形式守恒型方程（2-13）也能转化为特征变量型方程因：代入式（2-21）得： （2-24）令：， 可得： （2-25）与原始变量方程类似，在方程（2-13）上左乘以矩阵Tk-1，并结合式（2-25）进行化简，同样能得到特征变量方程（2-23）。

5 湍流平均N-S方程组在湍流的数值模拟中目前出现的方法有三类，分别是：1）使用湍流模型（Turbulence Model）的平均N-S方程组；2）大尺度涡模拟（Large Eddy Simulation，LES）加小涡模型；3）湍流直接数值模拟（Direct Numerical Simulation，DNS）使用DNS计算在理论上说是最直接的算法，它求解湍流中所有时刻和空间尺度参数的变化分布，这样一个问题所需的努力将与雷诺数的三次方成正比增加既使借助于当代的或未来的超级计算。