北航CFD讲义第18课

迁移方程0t x u au += （10-1）若a 为大于零的正数，则Euler 后差格式110n n n n i i i i u u u u a t x+---+=∆∆ （10-2） 稳定。

若a 为小于零的负数，则Euler 前差格式110n n n ni i i i u u u u a t x++--+=∆∆ （10-3） 稳定。

若a 是可正可负的变量，则需将a 进行分裂，令a a a +-=+式中，()12a a a +=+ ， ()12a a a -=- 于是有，0a +≥ ， 0a -≤迁移方程可写为，0t x x u a u a u +-++= （10-4）迎风格式，1110n n n n n ni i i i i iu u u u u u a a t x x++--+---++=∆∆∆ （10-5） 总是稳定的。

当0a ≥时，式（10-5）与（10-2）相当，当0a ≤时，式（10-5）与（10-3）相当。

对于一维Euler 方程0t x U F += （10-6）式中，123U U u U e U ρρ⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ ()1223u F F u p F e p u F ρρ⎡⎤⎡⎤⎢⎥⎢⎥=+=⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎣⎦补充状态方程，2112p e u ργ=+- （10-7） 对于方程（10-6），若设F A U∂=∂ 则方程（10-6）成为0t x U AU += （10-8）如何分裂A ？两个步骤： 1，求出A2，求出A 的特征值。

可将矢通量F 写成矢恒量U 的显函数形式，即123U U U U ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦ ()2223122231131212U U F U U U U U U U γγγγ⎡⎤⎢⎥⎢⎥⎢⎥-=+-⎢⎥⎢⎥⎢⎥⎛⎫-⎢⎥- ⎪⎢⎥⎝⎭⎣⎦第1个步骤：()()111123123222123123333123,,,,F F F U U U F F F F F F F A U U U U U U U F F F U U U ⎡⎤∂∂∂⎢⎥∂∂∂⎢⎥∂⎢⎥∂∂∂∂===⎢⎥∂∂∂∂∂⎢⎥⎢⎥∂∂∂⎢⎥∂∂∂⎣⎦（10-9） ()()()23201033123112F A u u U eu e u u u γγγγγγγγρρ⎡⎤⎢⎥⎢⎥∂-⎢⎥==---⎢⎥∂⎢⎥⎢⎥-+---⎢⎥⎣⎦（10-10） 第2个步骤，求出A 的特征值。

对A 作相似变换1A P AP -=Λ （10-11）式中，A u u au a ⎡⎤⎢⎥Λ=+⎢⎥-⎢⎥⎣⎦对A Λ可以进行分裂，设A A A +-Λ=Λ+Λ （10-12）其中，222A u uu a u au a u a ++⎡⎤⎢⎥⎢⎥+++⎢⎥Λ=⎢⎥⎢⎥-+-⎢⎥⎢⎥⎣⎦由式（10-11）可得，1A A P P -=Λ （10-13）将（10-12）代入（10-13）得()111A A A A A P P P P P P +--+---=Λ+Λ=Λ+Λ （10-14）设 1A A P P ++-=Λ 1A A P P ---=Λ 于是式（10-14）成为A A A +-=+从而完成了第2个步骤，分裂了A ，于是Euler 方程（10-8）可写成0t x x U A U A U +-++=迎风格式可写成1110n n n n n ni i i i i iU U U U U U A A t x x++--+---++=∆∆∆ （10-16） 这是非守恒型的迎风格式，并不常用。

222Au u u a u au a u a --⎡⎤⎢⎥⎢⎥+-+⎢⎥Λ=⎢⎥⎢⎥---⎢⎥⎢⎥⎣⎦（10-15）一．对守恒型Euler 方程采用矢通量分裂措施应当直接分裂Euler 方程中的F ，即对Euler 方程（10-6）0t x U F +=实施迎风格式，这就是所谓矢通量分裂格式。

定理：若函数()F U 恒等地满足下列关系()()k F lU l F U =则称()F U 是一个k 次的齐次函数。

对于这种函数，只要它可微，就有F kF U U∂=∂上式称为齐次函数的Euler 公式。

在Euler 方程中，F 是关于U 的一次齐次函数，即1k =。

根据微分学中齐次函数的Euler 公式，有F AU = （10-17）由于对A 已进行了分裂（A A A +-=+），于是有()F A A U +-=+令F A UF A U++--⎧=⎪⎨=⎪⎩ （10-18） 于是，完成了对F 的分裂F F F +-=+Euler 方程（10-6）成为0t x x U F F +-++=矢通量分裂格式可写为：()()()()1110n n n n n n iiii i iF F F FUU txx++--+-+---++=∆∆∆归纳一下，分裂F 的步骤：1．求出A （FA U∂=∂）2．找到矩阵P ，使1A P AP -=Λ 3．求出A +Λ，A -Λ（()12A A A +Λ=Λ+Λ，()12AA A -Λ=Λ-Λ） 4．求出A +，A -（1A A P P ++-=Λ，1A A P P ---=Λ） 5．求得F +，F -（F A U ++=，F A U --=） 对于二维和三维情况，可以依次类推。

例如，对于二维Euler 方程，0t x y U F G ++=同样可进行矢通量分裂，即F F F +-=+G G G +-=+于是，Euler 方程就成为0t x x y y U F F G G +-+-++++=那么矢通量分裂格式可写成：()()()()()()()()1,,,1,1,,,,1,1, 0n n n n n n i ji ji ji ji ji jn n n n i ji j i j i jF F F F UUtx xG G G Gyy++--+-+++---+---++∆∆∆--++=∆∆二．矢通量分裂格式（采用矢通量分裂措施的差分格式）1．一阶精度显式格式()()()()1110n n n n n n iiii i iF F F FUU txx++--+-+---++=∆∆∆2．空间精度为二阶的格式()()()()()()1122113421 4302n n n n n i i i i i n n n i i i U U F F F t x F F F x ++++-----++-⎡⎤+-+⎢⎥⎣⎦∆∆⎡⎤+-+-=⎢⎥⎣⎦∆3．MacCormack 格式()()()()()()()()1111111111112n n n n n n i i i i i i n n n n n n n i i i i i i i t t U U F F F F x x t t U U U F F F F x x +++---+++++++++---+∆∆⎧⎡⎤⎡⎤=----⎪⎢⎥⎢⎥⎣⎦⎣⎦∆∆⎪⎨∆∆⎧⎫⎡⎤⎡⎤⎪=+----⎨⎬⎢⎥⎢⎥∆⎣⎦∆⎣⎦⎪⎩⎭⎩由于采用了一侧差分，破坏了格式的对称性，只有一阶精度。

4．AF 格式（C-N 格式）22n n n n nb f b f t A t A F F I U t x x x x +-+-⎡⎤⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫∆∂∆∂∂∂++∆=-∆+⎢⎥⎢⎥ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎢⎥⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎣⎦式中，下标b 表示后差，f 表示前差。

对上式作近似分解：22n n n n nb f b f t A t A F F I I U t x x x x +-+-⎡⎤⎡⎤⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫∆∂∆∂∂∂++∆=-∆+⎢⎥⎢⎥⎢⎥ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎢⎥⎢⎥⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎣⎦⎣⎦以上四个格式都是针对一维Euler 方程的，可以直接推广到二维和三维。

三．人工粘性MC 格式和AF 格式由于采用中心差分使其等价微分方程中的二阶导数项都消失了，引起非物理的数值振荡，所以必须人为地把这些消失了的二阶项再加到格式中去，这些二阶项称为人工粘性。

对于MC 格式（以一维Euler 方程为例），有()()11111111112n n n n n n i i i i xx xxxx n n n n n n n i i i i i xx xxxx t U U F F tU tU x t U U U F F tU tU x αβαβ+-+++++++∆⎧=--+∆-∆⎪∆⎪⎨∆⎡⎤⎪=+--+∆-∆⎢⎥∆⎪⎣⎦⎩对于AF 格式（以二维Euler 方程为例）()()()22422 n nn x y n n t A t B I D I D U x y F G t D x y ⎛⎫⎛⎫∆∂∆∂++++∆= ⎪⎪∂∂⎝⎭⎝⎭⎡⎤⎛⎫∂∂⎛⎫-∆++⎢⎥ ⎪ ⎪∂∂⎝⎭⎝⎭⎢⎥⎣⎦式中，()222xx D x ε∂=-∂ ， ()222y y D y ε∂=-∂ ， ()44444n D U x y ε⎛⎫∂∂=+ ⎪∂∂⎝⎭在以上格式中，α、β、x ε、y ε和ε均为待定系数。

如何确定这些系数？四．人工粘性与矢通量分裂格式的关系对于迁移方程0t x u au +=迎风格式1110n n n n n ni i i i i i u u u u u u a a t x x++--+---++=∆∆∆ （10-5） 式中，()12a a a +=+，()12a a a -=-，代入上式，整理得， 111112222n n n n n n n i i i i i i i u u u u x u u u a a t x x++-+---∆-++=∆∆∆ 上式表明，迎风格式相当于Euler 中心差分格式加上二阶粘性。

对于一维Euler 方程为例0t x x U F F +-++= （10-19）简记 1i i x f f f x--∇=∆ —— 后差 1i ix f f f x+-∆=∆ —— 前差112i i x f f f xδ+--=∆ —— 中心差 则一阶精度矢通量分裂格式可写为：10n ni i x x U U F F t++--+∇+∆=∆ （10-20） 我们已知 F F F +-=+而F A U ++= ， F A U --=（10-21）而1A A P P ++-=Λ ， 1A A P P ---=Λ （10-22）而 ()12AA A +Λ=Λ+Λ ，()12A A A -Λ=Λ-Λ （10-23） 于是，将式（10-23）代入式（10-22），有()()1111122A A A AA PP P P P P +---=Λ+Λ=Λ+Λ 显然 1A A P P -=Λ ，令 1A A P P -=Λ于是 ()12A A A +=+， 同理 ()12A A A -=-， 代式（10-24）入式（10-21），有()12F A A U +=+ ， ()12F A A U -=- 于是()()2212x x x x x x x x A A A A F F U U AU A U +-+-⎛⎫⎛⎫∇+∆=∇+∆ ⎪ ⎪⎝⎭⎝⎭=∇+∆+∇-∆⎡⎤⎣⎦ （10-25） 根据定义：2x x x δ∇+∆= x x x x x ∇-∆=-∆∇∆代入式（10-25），得()2x x x x x xF F F AU δ+-∆∇+∆=-∇∆ （10-26）代入一阶精度矢通量分裂格式（10-20），得（10-24）()12n n i i x x x U U x F A U t δ+-∆+=∇∆∆ （10-27） 可见，一阶精度的矢通量分裂格式相当于中心差分格式加上二阶粘性。