高等计算流体力学讲义(3)

1高等计算流体力学讲义（3）§2Riemann问题1．预备知识：Euler方程解的结构我们讨论Euler方程解的结构。在上一节，我们已经得到，在均熵流动条件下，有constR，沿audtdx（1）其中auR12。且全场Sconst。（2）在这种情况下，Euler方程的光滑解有如下几种可能。1）在求解域中，Riemann不变量auR12均不为常数。这是最一般的情况，Euler方程的解比较复杂，通常无解析解。2）均匀流：Riemann不变量auR12均为常数。此时，令0RR，有：0000()/21()4uRRaRR，可见，此时流动是均匀的。3）简单波：有一个Riemann不变量在某区域内为常数（00RRorRR）。以0RR的情况为例。此时021RuaR。（3）且沿dxuadt，有21uaconst。这个常数具体的数值与特征线的起点有关。由此我们知道，沿dxuadt，有200()/21()4uRconstaRconst。这说明，沿dxuadt，u和a均为常数，即特征线是直线。由均熵条件，密度和压力p沿特征线dxuadt也为常数。参见上图，由于uau，所以流线dxudt（或流体质点）从左侧穿过特征线dxuadt，这种简单波称为左简单波或向后简单波。简单波可以分为压缩波和稀疏波（膨胀波）两类。设流线与dxuadt交点处，流线的切线方向为。把（3）式沿求方向导数，得：201ua当0u，有()0,0,0,0apuc。此时，压力密度沿流线减小，且特征线dxuadt是发散的。这种简单波称为稀疏波。当0u，有()0,0,0,0apuc。txdxudtdxuadt3此时，压力密度沿流线增加，且特征线dxuadt是收敛的的。这种简单波称为压缩波。对于0RR的情况可类似讨论。4）中心稀疏波中心稀疏波是一种特殊的简单波。以向左中心稀疏波为例R对应的特征线dxuadt是通过某一点的中心直线族。设这一点的坐标为(0,0)，则特征线的方程为：xuat（4）设左中心稀疏波的左边界（波头）的特征线为11xuat，右边界（波尾）特征线为22xuat。另外，由简单波定义，还有021RuaR。（5）显然，011222211Ruaua，（6）这就是波头波尾状态之间的关系。此外波头波尾之间还满足等熵关系12SS。（7）由（4）、（5）式可得左中心稀疏波的状态分布00112212[]111[],1xuRtxxaRuauatt压力和密度可由等熵条件得出。右中心稀疏波的讨论类似。txdxudtdxuadtdxuadtdxudt稀疏压缩45）激波和接触间断当流场中存在间断时，Rankine-Hogoniot关系为：2(,,)(,,)TTFDUUuEFuupuH。（8）假定间断两侧速度是连续的，则必有Duu。这种间断称为接触间断。由Rankine-Hogoniot关系易知，接触间断两侧压力也是连续的，发生间断的只有密度。不是接触间断的间断称为激波。对于激波而言，必有()()0uDuD。通过上面的讨论，我们可以得到一个重要的结论：与均匀流区相邻的区域，一定是简单波区或者间断线。间断显然可以和均匀流区相邻。如不存在间断，必有一族特征线从均匀流区进入这一区域。这样，相应的Riemann不变量不仅沿特征线是常数，而且在这个区域内保持常数，因此这一区域必然是简单波区。2．Riemann问题所谓Riemann问题就是求解Euler方程0xFtU(9)在初值00)0,(xUxUxURL(10)下的解。Riemann问题存在解析解。那么它的解是什么结构呢？在一般情况下，,LRUU之间不满足Rankine-Hogoniot关系，所以，随着时间的变化，初始间断会分解为一些特定的流动结构，而且这些结构均发源于初始间断处；而在远离这些结构的地方，流动仍保持为初始值,LRUU。因此，我们可以预期，Riemann问题解的结构是两侧的均匀流区和中部的间断影响区。见下图。5与均匀流区相邻的只能是简单波区或者间断。如果是简单波区，容易知道，其必为中心稀疏波。所以，我们知道，中心的初始间断影响区和均匀流区通过左波（激波或中心稀疏波）和右波（激波或中心稀疏波）分开。但是，是不是只有这两个波，或者说，Riemann问题的解是所谓“双波结构”呢？容易证明，双波结构是超定的，所以，初始间断影响区必然包含其他流动结构。这个结构就是上文所说的接触间断。它处于左波和右波之间。此时，未知量的个数和所能提供的独立方程数相同，我们可以唯一确定Riemann问题的解。4波或4波以上的解都不可能是稳定的物理解。因此，Riemann问题的解是一种“三波结构”：左波和右波可以是激波或膨胀波，中间的波为接触间断（接触间断两侧速度、压力是连续的，但密度有间断）。根据上述分析，Riemann问题解的结构如下图所示：在求解Riemann问题过程中，比较方便的是采用原始变量TpuW),,(所以Riemann问题的初始条件也可写为xt均匀均匀?txtxxttx600)0,(xWxWxWRL(11)各个区域内Riemann问题的解可以记为如下图的形式：tW*LW*RWLWRx注意LW*和LW*对应的压力和速度相同，密度不一定相同。即：RLRLRLpppuuu********,,3.Riemann问题的求解具体的求解过程可参见Toro的书。压力*P的计算压力*P是下述方程的解：***(,,)(,)(,)0LRLLRRfPWWfPWfPWu(12)其中uRu-Lu，12***1*2**,211LLLLLLLLLAPPifPPPBfPWaPifPPP712***1*2**,211RRRRRRRRRAPPifPPPBfPWaPifPPP且LA=L12LB=LP11RA=R12RB=RP11速度*u的计算：*u=21RLuu+21**PfPfLR(13)(12)式不能得到显式的解析解，一般要通过迭代法求解。在计算出*P,*u后，就可以得到Riemann问题的完整解。以左波附近的解为例，简述如下：如果已知*P，则：①当*PLP时，*P和LP之间有一道激波。由R-H关系，LLLLLUUDFF**（14）由于F是U的函数，所以上式中有四个未知量，三个方程。在得到*P后，(14)式封闭，我们可以通过(14)式求出*u(已算出),L*,LD。则其解如下图所示即x/*utx/LDttLW*LW8**,LxutLWWxtW*LLxDtxDut(15)②当*LPP时，*P和LP之间是中心膨胀波。此时左波附近解的结构如下图所示：其中,HLTLDD分别是膨胀波波头和波尾的速度：**,HLLLTLLDuaDua(16)*La的计算方法为，利用等熵条件1/**()LLLpp可得(1)/(2)**()LLLpaap(17)此时，左波附近的解为**/,//LHLxLfanHLTLutLTLWifxtDWxtWifDxtDWifxtD(18)其中，LfanW为中心稀疏波内的解，根据特征关系，可得：2/(1)2/(1)2(1)[(/)](1)(1)2(1)[/](1)22(1)[(/)](1)(1)LLLLfanLLLLLuxtaWuauxtppuxta(19)/HLxtDx/*utx/TLtDtLW*LW9右波附近的解可以用类似方法确定。根据上面的分析，我们可以看出，Riemann问题的解有相似性，即(,)(/)WxtWxt§3Godunov格式考虑Euler方程0UFtx(1)我们采用有限体积方法离散上述方程。在任意控制体1/21/21[,][,]iinnxxtt上积分上式,有:1/211/211/21/2[(,)(,)][((,))((,))]0ininxtnniixtUxtUxtdxFUxtFUxtdt(2)令ixxnnixdxtxUUii2121,（3）其中2121iiixxx，则（2）式可以写为：0,,11121211nnnnttttiininidttxUFdttxUFtxtUU（4）到目前为止整个过程是精确的。假定有某种方法计算出控制体边界12ix处的通量积分1,21nnttidttxUF，则通过（4）式可以求解1niU。可见有限体积法得到的是控制体内物理量在空间的平均值。即我们直接求得的只是IiUUUUniniii,,,,,,121（I是求解域的空间指标的集合）。但是在计算通量1,21nnttidttxUF时，我们需要知道U在21ix处的值，即U沿x的分布。因此，在有限体积格式中，要解决的一个重要问题就是如何由xU在单元内的平均值,,,,11iiiUUU求出一个函数xU~，使得xU~近似xU到一定精度。这个过程称为重构（Reconstruction）。最简单的重构是零阶重构：iixxxUxUxUii~~2121,（5）即假定xU的近似值xU~在每个控制体内均为常数，显然此时UxUxOx。可以证明，这种重构方法得到的差分格式最多有空间一阶精度。一阶重构可以写为：10iiiixxxxxDUxUxUii~~2121,其系数iD的确定方法以及更高阶的数据重构方法见后续内容。这一节和后面的几节，我们均采用零阶重构，因此所得到的格式只有一阶精度。以后我们还会详述高阶精度格式的构造方法。在采用零阶重构时，从n时刻物理量的平均值niU得到的U的分布为：1/21/2[,](,)(,)iinnxxxiniUxtUxtU（6）（此时),(ntxU不再是精确的，但为了简单，我们略去了上标“～”）。同理niniUtxU11),(。所以在21ix附近U(x,tn)的分布为),(),(),(1ninintxUtxUtxU2/12/1iixxxx（7）令niniiLininiiRiUtxUUUtxUU),(),(2121121121（8）则（7）式也可以写为LiRinUUtxU2121),(2121iixxxx（9）因此在x=xi+1/2处，U(x,tn)的分布有一间断。在解决了数据重构问题后，我们必须解决通量1)),((2/1nnttidttxUF的计算问题。具体说是要解决tn+1ttn时，U(xi+1/2,t)的计算问题。Godunov通过在21ix附近求解以（9）为初值的Riemann问题来计算U(xi+1/2,t)。通过求解Riemann问题可以得到其精确解为)(),(2/1nittxxWtxW（10）由此容易得到对应的守恒变量解为：)(),(2/1nittxxUtxU（11）ttn时，)0(),(2/1UtxUi（12）因此，在xi+1/2处，Riemann问题的解为在tnttn+