3-通量分裂

1§3.通量分裂3.1迎风原则先考虑线化的模型方程-对流方程0uuatx抖+=抖其中a为常数。1nnnjjjuuuxx--¶Ü禗空间导数用向后差分近似1nnnjjjuuuxx+-¶»禗空间导数用向前差分近似110nnnnjjjjuuuuatx+---+=DD向后差分格式110nnnnjjjjuuuuatx++--+=DD向前差分格式迎风原则的示意图23.2迎风格式向后差分格式和向前差分格式都可以看成是对中心差分格式的修改11111222nnnnnnnjjjjjjjuuuuuuuaatxx++-+----++=DDD（用于0a）11111222nnnnnnnjjjjjjjuuuuuuuaatxx++-+----++=-DDD（用于0a）不仅如此，经过这样的改写，我们发现它们还可以统一起来，写成11111222nnnnnnnjjjjjjjuuuuuuuaatxx++-+----++=DDD如果再经过进一步的改写11111022nnnnnnnnnnjjjjjjjjjjuuuuuuuuuuaatxx++-+---+---++-=DDD并重新组合格式中的各项，就得到111022nnnnnnjjjjjjaaaauuuuuutxx+-++----++=DDD记2aaa++=，2aaa--=，则上述格式可写成1110nnnnnnjjjjjjuuuuuuaatxx+-++----++=DDD注意到，无论0a还是0a，都有302aaa++=，02aaa--=而且22aaaaaaa+-+-+=+=于是我们可以对上面得到的格式做出新的解释：将a分解成了“正”和“负”两个部分aaa+-=+，0a+，0a-对流项也相应地分裂成了两项。然后根据迎风原则，对这两项分别用向后差分和向前差分来近似。这样得到的差分格式1110nnnnnnjjjjjjuuuuuuaatxx+-++----++=DDD称为迎风格式。3.3通量分裂考虑守恒律方程组()tx¶¶+=抖fww0此时，¶=¶fAw是一个矩阵，0A、0A、A等表达式没有意义，4因此迎风格式根本无法定义。要想解决上述困难，需要推广和拓展迎风格式的构造思路。前面曾对迎风格式的构造思路做过如下的解释：将a分解成了“正”和“负”两个部分，对流项也相应地分裂成了两项。然后根据迎风原则，对这两项分别用向后差分和向前差分来近似。对方程的守恒形式，则给出了守恒型格式的构造思路：对方程中通量的空间导数直接进行差分近似。现在，将这两个思路结合起来，就得出了守恒型迎风格式的构造思路：直接将通量分解成了“正”和“负”两个部分，对流项也相应地分裂成了两项。然后根据迎风原则，对这两项分别用向后差分和向前差分来近似。这种构造差分格式的方法就称为通量分裂。假设通量向量已经分解成“正”和“负”两个部分()()()+-=+fwfwfw则守恒律方程组就可分裂成()()txx+-抖¶++=抖?fwfww05于是就可以直接写出迎风格式的守恒形式()()()()111nnnnnnjjjjjjtxx++--+-+---++=DDDffffww0由此可见，现在的问题已经转化为如何实现通量向量的迎风分裂。为此，需要利用双曲型方程组的一些数学性质：雅克比矩阵A的特征值1、2、…、m都是实数雅克比矩阵A相似于一个对角矩阵。即：存在可逆矩阵R，使得121m-骣÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷==÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷÷ç桫RARDO通量具有齐次性，即：不仅有微分关系式dd=fAw，还有代数关系式=fAw。根据第一条性质，可以先将每一个特征值分裂成iii+-=+，02iii++=，02iii--=，（1,2,,im=L）就可以定义对角矩阵D的分裂6+-=+DDD，12m++++骣÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷=÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷÷ç桫DO，12m----骣÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷=÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷÷ç桫DO这里，+D的对角线元素都是大于等于零的，而-D的对角线元素都是小于等于零的。另外，还可以引进记号12m骣÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷=÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷÷ç桫DO则有2++=DDD，2--=DDD再根据第二条性质()1111-+--+---==+=+ARDRRDDRRDRRDR很自然地，可以定义雅克比矩阵A的分裂+-=+AAA，1++-=ARDR，1---=ARDR7以及1-=ARDR于是1111222--++--+++====DDRDRRDRAAARDRRR1111222---------====DDRDRRDRAAARDRRR最后根据第三条性质()+-+-==+=+fAwAAwAwAw最终可以定义通量f的分裂+-=+fff，++=fAw，--=fAw这样就解决了通量向量的迎风分裂问题。显然，如果雅可比矩阵A的特征值都是正的，则ii+=，0i-=（1,2,,im=L）+=DD，-=D0，+=AA，-=A0，+=ff，-=f0而如果雅可比矩阵A的特征值都是负的，则80i+=，ii-=（1,2,,im=L）+=D0，-=DD，+=A0，-=AA，+=f0，-=ff这实际上就是迎风原则。还要指出，上面给出的分解方案并不是惟一的。下面针对计算流体力学，以一维气体力学方程组为例，给出几种不同的分裂方式。3.4速度分裂通量雅可比矩阵的特征值为1uc=-，2u=，3uc=+注意到音速c总是大于零的，首先可以想到是，只需对速度u进行分解就行了，uuu+-=+，2uuu++=，2uuu--=于是就可以定义对角矩阵D的另一种分裂+-=+DDD，uuuc++++骣÷ç÷ç÷ç÷ç÷ç=÷ç÷ç÷÷ç÷ç÷+ç÷ç桫D，ucuu----骣-÷ç÷ç÷ç÷ç÷ç=÷ç÷ç÷÷ç÷ç÷ç÷ç桫D其余的就和前面一样了。93.5Steger-Warming分裂通量f可以改写成()2222222uupupuuEpEpu骣÷ç÷ç÷ç÷骣ç÷÷çç÷÷çç÷÷çç÷÷ç÷ç÷÷çç=+=+÷÷ç÷ç÷ç÷ç÷÷çç÷÷çç÷÷+çç÷÷ç桫+÷ç÷ç÷ç÷÷ç桫f利用通量雅可比矩阵的三个特征值1uc=-，2u=，3uc=+以及声速的表达式2pc=和关系式2222111112121112Epppupucu轾骣+÷ç犏÷=++=+ç÷犏ç÷ç--桫臌=+-有()()()()()21322122121uuuuucuc=-+=-+-++=-++()()()()()()()2222222132221222121puuucuuucucuucuc+=-++=-?-++=-+-++10()()()()()()()()()()()()()()()222332323332222211222121211326123111122312112Epuucuucuuucuucuucuccucucuuucuccucucucuc骣+÷ç÷=+=+ç÷ç÷ç--桫骣÷ç÷=-+-++ç÷ç÷ç-桫-轾轾=-+?+++-++犏犏臌犏-臌-轾=-?-++犏臌-轾+-?++?犏犏臌=()()()()()222221313311221ucucuc-轾-+++-++犏犏-臌这表明，通量向量f的每一个分量都可以写成三个特征值的线性函数112233112233112233aaabbbccc骣++÷ç÷ç÷ç÷ç÷ç=++÷ç÷ç÷÷ç÷ç÷++ç÷ç桫f（注）这些表达式实际上是拟线性关系式，其中的系数与u、c和有关，例如()2221122cuu-=?=。于是，就可以有如下的分裂方式，称为Steger-Warming分裂112233112233112233aaabbbccc++++++++++骣++÷ç÷ç÷ç÷ç÷ç=++÷ç÷ç÷÷ç÷ç÷++ç÷ç桫f，112233112233112233aaabbbccc----------骣++÷ç÷ç÷ç÷ç÷ç=++÷ç÷ç÷÷ç÷ç÷++ç÷ç桫f113.6vanLeer分裂利用马赫数的定义uMc=，可将通量向量f的第一个分量写成221122MMucMcmm+-轾骣骣+-犏鼢珑鼢==-=+珑犏鼢珑鼢珑桫桫犏臌式中212Mmc+骣+÷ç÷=ç÷ç÷ç桫，212Mmc-骣+÷ç÷=-ç÷ç÷ç桫于是，通量f的其余两个分量就可写成()()()()()()()2222222222222221111111111111112222211upuccMcMMcMMMMMMMcMcM+=+=+轾轾=-++=-?+犏犏臌臌禳轾轾镲骣骣骣骣镲+-+-犏犏鼢鼢镲珑珑鼢鼢=--++珑珑睚犏犏鼢鼢珑珑鼢鼢珑珑镲桫桫桫桫犏犏镲臌臌镲铪=-()222222211112222211111122222MMMMccMMMMcuc禳轾轾镲骣骣骣骣镲+-+-犏犏鼢鼢镲珑珑鼢鼢-++珑珑睚犏犏鼢鼢珑珑鼢鼢珑珑镲桫桫桫桫犏犏镲臌臌镲铪轾骣骣骣骣+-+-犏鼢鼢珑珑鼢?=--++珑珑犏鼢?珑珑鼢?珑珑桫桫桫桫犏臌()()()()2221111121222111212MMcuccucmucmuc+-禳轾镲镲犏镲÷睚犏÷÷镲犏镲臌镲铪骣骣+-鼢轾轾珑鼢=-+---珑鼢犏犏珑臌臌鼢珑桫桫轾轾=?++?-犏犏臌臌12()()()()()()()()()()()()22223223222232222322211122112121211121211122214211142121EpuEpuucuucMcMMucMMucMMMcMMM骣+÷ç÷+=?+ç÷ç÷ç-桫轾骣+犏÷ç÷=+=?ç犏÷ç÷ç--桫犏臌轾=?++犏臌-轾=?+-+-+犏犏臌-轾=?++-+犏犏臌-(){}()()()(){}()()()23222222322221114211211111422211121222cMMMMMMcMMMM轾=-++-+犏犏臌-ì轾ï骣骣ï+-轾犏鼢ï珑鼢=-+-犏珑í犏鼢珑鼢珑ï犏臌桫桫犏-ï臌ïîü轾ï骣骣ï+-犏鼢ï珑鼢+-?珑ý犏鼢珑鼢珑ï桫桫犏臌þ()()()()()()()223222222232111414221114142111222112McMMMMMMcMMïïìï骣+ï轾÷ïç÷=-+-+犏çí÷ç÷çï犏臌桫-ïïîüï骣-ï轾÷ïç÷----+犏çý÷ç÷çï犏臌桫ïïþìï骣+ï÷ï轾ç÷=-+çí÷犏ç臌÷çï桫-ïïîæ--()2212Müïöï÷ï轾ç÷--çý÷犏ç臌÷çï桫ïïþ13()()()()()()()()2222222222111222111122211112122121McMccMcMccmucmuc+-骣+÷轾ç÷=?+ç÷犏ç臌÷ç桫-骣-÷轾ç÷-?-ç÷犏ç臌÷ç桫-轾轾=?++?-犏犏臌臌--于是，又可定义如下的分裂方式()()()22111211221mucuc++骣÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷÷ç÷ç÷ç轾÷ç=-+÷ç犏÷臌ç÷ç÷÷ç÷ç÷ç÷ç÷ç轾÷-+ç÷犏ç臌÷-ç÷÷ç桫f，()()()22111211221mucuc--骣÷ç÷ç÷ç÷ç÷ç÷ç÷ç÷÷ç÷ç÷ç轾÷ç=--÷ç犏÷臌ç÷ç÷÷ç÷ç÷ç÷ç÷ç轾÷--ç÷犏ç臌÷-ç÷÷ç桫f称为vanLeer分裂。3.7Lax-Friedrichs分裂设雅克比矩阵A的特征值为1、2、…、m，则()()12max,,,mr==AL称为矩阵A的谱半径。根据线性代数的有关知识，矩阵()12r+=+AAI的特征值为14102r+³，202r+³，…，02mr+³而矩阵()12r