高等計(jì)算流體力學(xué)-03

1、求解流體力學(xué)方程組的有限差分求解流體力學(xué)方程組的有限差分方法方法第三講1控制方程2 貼體貼體網(wǎng)格網(wǎng)格 H型 C型 O型( , )x y( , ) 3一、一、貼體坐標(biāo)中的基本方程貼體坐標(biāo)中的基本方程( , )( , )( , )( , )x yxxx yyy 正變換正變換逆變換逆變換假定已知上述變換的離散形式假定已知上述變換的離散形式即已經(jīng)生成好貼體網(wǎng)格即已經(jīng)生成好貼體網(wǎng)格4導(dǎo)數(shù)的變換導(dǎo)數(shù)的變換 一階導(dǎo)數(shù)xxyyxy5 二階導(dǎo)數(shù)222222222222222()()()()2()xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 222222222()2()yyyyyyyyy 2222

2、22()xyxyxyxyyxxyx y 6 二階導(dǎo)數(shù)22220 xy222222222()() 2()() ()()0 xyxxyyxyxxyyxxyy 度量系數(shù)度量系數(shù)metrics 7度量系數(shù)及其計(jì)算方法度量系數(shù)及其計(jì)算方法 1,1,1,1,1,1,1,1,)2)2)2)2ijiji ji ji ji jijiji ji ji ji jxxxxxxyyyyyy困難困難容易容易如果能夠找到如果能夠找到 和和 之間的關(guān)系，之間的關(guān)系，就可以得到就可以得到 等的計(jì)算方法等的計(jì)算方法 ,xyxy ,xxyy,xyxy 8xyxyddxdyddxdyxyxyddxddyxxdxdyydyddxx d

3、x ddyy dy d( , )( , )x yx y( , )( , )xxyy 1xyxyxxyy1xyxyyxyxJ1xyxyxxJx yx yyy1111,xxyyyyxxJJJJ Jacobi行列式行列式（jacobian） ,xyxy 910 xyxyxyxy11,xyyxJJ 22()()()()()()2()0()()()()()xyxyxyxxxyxyyyxyxxxyyyxyxyxyxyxxxyxyyyxxyxyxyxyxxyyxyxx yyxyxyxyxyxyxxyy()0yxxxyyyxyxyx xx yx yy yxy()()()()()()0()()()()()xyx

4、yxyxxxyxyyyxyxxxyyyxyxyxyxyxxxyxyyyxyxyxyxyxxyyxyx xx yx yy yxyxyxyxyxyxxyy 22()2()0 xyxxxyyyxyxyxx yyxy,xxxyyyxxxyyy102222()2()()()2()xyxxxyxyyyxyxyxx yyx xx yx yy yxyxx yyxy 12222()2()()()2()xyxxxyxyyyxyxyxx yyx xx yx yy yxyxx yyxy 11222222()2()1()()() ()2()xxxyyyxyxyxyJ yJy yJ yJx yJ x yx yJx yJx

5、 yx yJ xJx xJ xxyxyxy222222()2()1()()() ()2()xxxyyyxyxyxyJ yJy yJ yJx yJ x yx yJx yJx yx yJ xJx xJ xxyxyxy12任意曲線坐標(biāo)系中流體力學(xué)方程組的守恒形式任意曲線坐標(biāo)系中流體力學(xué)方程組的守恒形式()()0vvtxyFFGGU()()()()0vvvvxyxytFFGGFFGGU13()()()()0vvvxyxvyJJJJtJUFFGGFFGG()()() ()()()()()()()() ()()()()()(),vvxyvxvyvxvyvvxyvxvyvxvyJJttJJJJJJJJJJU

6、UFFGGFFGGFFGGFFGGFFGGFFGG14() ()() ()()()()()()()()0vxvyvxvyvxxvyyJJJtJJJJUFFGGFFGGFFGG()()()()0()()()()0 xxyyJJyyJJxx() ()() ()()0vxvyvxvyJtJJUFFGGFFGG15() ()() ()()0vxvyvxvyJtJJUFFGGFFGG( )()()0vvtUFFGGxyxyvvxvyvvxvyJJJJJUUFFGGFGFFGGFG曲線坐標(biāo)系中曲線坐標(biāo)系中NavierStokes方程的強(qiáng)守恒形式方程的強(qiáng)守恒形式 16二、守恒型Euler方程在任意曲線坐標(biāo)系

7、中，二維守恒型Euler方程為：( )( )()0tUFG其中xyxyJJJUUFFGGFGuvEU2()uupuvEp uF2()vvuvpEp vG17二維守恒型Euler方程可以改寫(xiě)為擬線性的形式：可以改寫(xiě)為擬線性的形式： ( )( )( )0ABtUUU其中xyxyxyxyAABBABFFGUUUGFGUUU令 或 ，則 可以寫(xiě)為統(tǒng)一形式：00(2)(1)(1)(1)(2)(1)(2)()(1)()(1)xyxxyxxyxyyyxyuuuvCvvuvEEuEv 18其中， 為絕熱指數(shù)，對(duì)于完全氣體， ， 和 的定義為：1.4221(),2xyuvuv,CACBCABxyCAxyCBFG

8、UU顯然時(shí)，時(shí)，。另外，由 的表達(dá)式，我們也容易導(dǎo)出直角坐標(biāo)系下Jacobi矩陣的表達(dá)式：當(dāng)取，而且時(shí)，。當(dāng)取且時(shí) 的特征值為：C1234xyxyxyuvauvuva其中 是音速。/ap19三、守恒型的NavierStokes方程平面上NavierStokes方程的強(qiáng)守恒形式為：( , ) ( )()()0vvtUFFGG其中xyxyvvxvyvvxvyJJJJJUUFFGGFGFFGGFG0 xxvxyxxxyTuvkxF0 xyvyyxyyyTuvkyG2(),2(),()xxxxyyyyxyxyyxyxuuvvuvuv注意到，Euler方程可以看作NavierStokes方程的特例。一般

9、而言，Euler方程中無(wú)粘通量 的離散方法，同樣可用于NavierStokes方程；而為了求解NavierStokes方程，我們還必須提供粘性通量 的離散方法。,vvF G20有限差分方法21有限差分方法() () ()0 xyxyJJJtFGFGU對(duì) 式直接進(jìn)行半離散差分近似 注意到 ，采用中心差分，則有：1 ,1/2,1/2,1/2,1/2()1() () () () 0i jxyijxyiji jxyi jxyi jJJtJJJUFGFGFGFG,1,1,1,1()1() () 2() () 0i jxyijxyiji jxyi jxyi jJJtJJJUFGFGFGFG1/2,1,()

10、 () /2ijxyi jxyijJJHFGFG令，則可以得到守恒格式：221/2,1,() () /2ijxyi jxyijJJHFGFG在有限差分方法中，在網(wǎng)格點(diǎn)上的度量系數(shù)(如 )也需要通過(guò)差分方法計(jì)算，因此要求計(jì)算網(wǎng)格是充分光滑的。,()xi jJ,1/2,1/2,1/2,1/2,1()0i jijiji ji ji jtJUHHHH1/2,1,() () /2ijxyi jxyijJJHFGFG23有限差分方法的主要特點(diǎn)是：1）有限差分方法只需構(gòu)造偏導(dǎo)數(shù)的離散方法，這使得它比較容易推廣到高階精度，對(duì)于多維問(wèn)題也是如此。2）有限差分方法對(duì)網(wǎng)格的光滑性有較高要求，不易在復(fù)雜形狀求解域上實(shí)

11、施。3）在曲線坐標(biāo)系中，有限差分方法要對(duì)幾何量和物理量的確定組合進(jìn)行差分離散，這樣做的后果之一是有限差分方法可能產(chǎn)生所謂幾何誘導(dǎo)誤差（geometrically induced errors） 。例如，考慮一無(wú)界均勻流場(chǎng)，易知，該流動(dòng)應(yīng)是定常的。即 。在由一般曲線坐標(biāo)構(gòu)成的網(wǎng)格上采用有限差分方法計(jì)算此流動(dòng)，可能有 （三維以及動(dòng)網(wǎng)格中的有限差分方法更容易出現(xiàn)這種情況）。而采用有限體積方法，則恒有 。有限差分方法在某些情況下不能復(fù)現(xiàn)均勻流場(chǎng)，是幾何誘導(dǎo)誤差的表現(xiàn)之一。, ,0i ji jtU,0i jtU, ,0i ji jtU()()()()0()()()()0 xxyyJJyyJJxx24如何

12、計(jì)算數(shù)值通量？ 數(shù)值耗散 數(shù)值色散1/2,1,() () /2 ?ijxyi jxyijJJHFGFG線性: 分辨率, 精度非線性非線性: 激波捕捉激波捕捉分辨率, 相位畸變簡(jiǎn)單中心差分一般來(lái)講不是好的選擇簡(jiǎn)單中心差分一般來(lái)講不是好的選擇!25激波間斷和廣義解26可壓縮流動(dòng)中一個(gè)非常重要而又復(fù)雜的現(xiàn)象是激波現(xiàn)象。在NavierStokes方程的范疇內(nèi)，激波可以看作是一個(gè)連續(xù)但物理量的梯度非常大的有限厚度的結(jié)構(gòu)。在Euler方程的范疇內(nèi)，由于方程本身缺乏必要的耗散機(jī)制，激波厚度為零，激波兩側(cè)物理量存在間斷。在數(shù)值求解Euler方程時(shí)，必須解決包含間斷的流場(chǎng)的計(jì)算問(wèn)題。在數(shù)值求解NavierSto

13、kes方程時(shí)，雖然理論上激波附近物理量是連續(xù)的，但是在通常情況下，激波的厚度非常薄（為分子平均自由程的量級(jí)），遠(yuǎn)遠(yuǎn)小于計(jì)算網(wǎng)格的尺度；所以在計(jì)算中，我們?nèi)孕璋鸭げó?dāng)作間斷來(lái)處理。在Euler方程的范疇內(nèi)如何表示含激波間斷的流場(chǎng)?27一、激波的形成激波是一種非線性現(xiàn)象?？紤]線性對(duì)流方程的初值問(wèn)題，如果其初始值是連續(xù)的，則解中不會(huì)自發(fā)產(chǎn)生間斷。但是，對(duì)于非線性問(wèn)題，即使初值是連續(xù)的，在后續(xù)時(shí)間中，也可能發(fā)展出包含間斷的解。Bergers方程可以看作Euler方程的非線性模型方程；和Euler方程一樣，是雙曲型方程。下面，我們以Bergers方程為例，介紹間斷解或者激波的形成過(guò)程。0uuatx0uu

14、utx28考慮Bergers方程0uuutx的初值問(wèn)題，初始條件為 。根據(jù)特征理論，Bergers方程的解析解可以寫(xiě)為下列形式：0( ,0)( )u xux0( , )(,0)()u x tu xutuxut初始條件初始條件010( )10101xuxxxxBergers方程的解為：方程的解為：11( , )1101xtxu x ttxtx29ut=1/2t=1t=0 x圖圖1 Bergers方程初值問(wèn)題的解方程初值問(wèn)題的解圖1顯示了解隨時(shí)間的發(fā)展過(guò)程。容易知道，當(dāng) 時(shí)，Bergers方程的解為1t 111( , )01txu x tx可見(jiàn),對(duì)于非線性的Bergers方程,即使初始值是連續(xù)的,

15、其解仍然可能出現(xiàn)間斷。 事實(shí)上，一般的非線性雙曲型守恒方程（組）或稱非線性雙曲型守恒律的解中都有可能存在間斷。通常，我們認(rèn)為偏微分方程的解是連續(xù)可微的。顯然，對(duì)于非線性雙曲型守恒律，這一點(diǎn)并不成立。因此，必須拓展雙曲型守恒律解的概念。30二、廣義解對(duì)于一般的一維雙曲型守恒律 的初值問(wèn)題，我們引入廣義解（generalized solution）或稱弱解(weak solution)的概念，來(lái)表示包含間斷的解。0txUF定義：設(shè) 是分片連續(xù)可微的函數(shù)，在 的半平面，如果對(duì)于與 的間斷線只有有限個(gè)交點(diǎn)的任意分段光滑的閉曲線 ，都有: 則稱 為方程 在初值 下的廣義解或弱解。( , )x tU0t

16、( , )x tU( )0dtdxF UU( , )x tU0txUF0( ,0)( )xxUUx 31 下面討論一下廣義解的意義。如果已知 是光滑的，設(shè) 圍成的 區(qū)域?yàn)?，則由 式利用Green公式知( , )x tU( )0dtdxF UU()0dxdtdtdxtxUFFU由于閉曲線可以在光滑區(qū)內(nèi)任取，可得：0txUF在光滑區(qū)，弱解就是通常的連續(xù)可微解。在光滑區(qū)，弱解就是通常的連續(xù)可微解。32如果 是由一條間斷線 分隔開(kāi)的分片連續(xù)可微函數(shù)，取如圖所示的由 （ 是正數(shù)）圍成的閉曲線 ，在 上應(yīng)用 ,有:( , )x tU xx t 12,xx ttt tt( )0dtdxF UU 2212(

17、)2()( ),( , )tx ttx tx x tdxx ttdtx ttdtx t dxdtF UUU 1121( )1( )( ),( , )0tx ttx tx x tdxx ttdtx ttdtx t dxdtF UUUx=x(t)xtt=t2t=t1P xx t xx t33間斷關(guān)系的推導(dǎo)間斷關(guān)系的推導(dǎo) 21( )( )0,0,00,0,tx x ttx x tdxx ttx ttdtdtdxx ttx ttdtF UUF UU令 ，并考慮到在 兩側(cè) 有間斷，則上式可簡(jiǎn)化為：0 xx t( , )x tU 分別表示從左、右趨近間斷線時(shí)守恒變量的值。 0,x ttU 0,x ttU令

18、令 0,x ttUU 0,x ttUU( )x x tdxDdt在間斷兩側(cè)，顯然 。 考慮到 可以任意取值，可得：UU12,t tDFU其中其中 FF UF UUUURankine-Hogoniot（R-H）關(guān)系）關(guān)系x=x(t)xtt=t2t=t1P xx t xx t34以一維Euler方程為例，R-H關(guān)系為：2()uupDuEp uE注意到，當(dāng) 時(shí)，如果 ，則上面的關(guān)系成立，此時(shí)只有密度場(chǎng)存在間斷。這種類型的間斷稱為接觸間斷。不屬于接觸間斷的間斷稱為激波。Duupp綜上所述，雙曲型守恒律的綜上所述，雙曲型守恒律的廣義解廣義解 是被有限個(gè)間斷線分開(kāi)的分片光滑是被有限個(gè)間斷線分開(kāi)的分片光滑函

19、數(shù)。在光滑區(qū)函數(shù)。在光滑區(qū)， 是滿足微分方程是滿足微分方程式式 的連續(xù)可微解；在間的連續(xù)可微解；在間斷線的兩側(cè)斷線的兩側(cè)， 滿足滿足R-H關(guān)系。關(guān)系。, x tU, x tU0txUF, x tU35三、熵條件廣義解推廣了偏微分方程初值問(wèn)題連續(xù)可微解的概念，但其后果是廣義解推廣了偏微分方程初值問(wèn)題連續(xù)可微解的概念，但其后果是導(dǎo)致了弱解不唯一導(dǎo)致了弱解不唯一。為了說(shuō)明這一問(wèn)題，我們舉一個(gè)例子：考慮Burgers方程在初值為 時(shí)的解010( )10 xuxx此時(shí)，初值在 處有一個(gè)間斷。0 x 把Bergers方程寫(xiě)為守恒形式2(/2)0txuu可知 處的Rankine-Hogoniot條件為：0

20、x 220000(/2)|(/2)|( |)xxxxuuD uu 由上式知 。0D 010( , )( )10 xu x tuxx在間斷處滿足Rankine-Hogoniot條件，在其他地方滿足微分方程，即 是Bergers方程的一個(gè)廣義解。0( , )( )u x tux36另外，考慮1( , )/1xtu x tx ttxtxt 分片連續(xù)可微函數(shù)，其導(dǎo)數(shù)在 處存在間斷(或者認(rèn)為函數(shù)值存在幅度為零的間斷)。xt 根據(jù)廣義解的定義，容易驗(yàn)證此函數(shù)也是Bergers方程的一個(gè)廣義解。所以廣義解一般不唯一，但是對(duì)于有明確物理意義的守恒律，其中只有一個(gè)解是有物理意義的，我們稱之為物理解。為了得到我們

21、關(guān)心的物理解，廣義解除了必須滿足式外，還必須滿足附加的條件，這個(gè)條件因?yàn)榕c熱力學(xué)第二定律所起的作用相同，被稱為熵條件（entropy condition）。37如何在不唯一的廣義解中挑選所謂物理解呢？最直觀的方法是認(rèn)為物理解是一個(gè)光滑的粘性問(wèn)題的解在粘性趨于零時(shí)的極限。也就是說(shuō)：方程方程 的的解如果解如果當(dāng)當(dāng) 時(shí)，幾乎處處有界的收斂到時(shí)，幾乎處處有界的收斂到分片連續(xù)可微函數(shù)分片連續(xù)可微函數(shù) ，則則 是是 的物理解。的物理解。22(0)txxUFU0, x tU, x tU0txUF這種方法的物理基礎(chǔ)是：在真實(shí)的物理過(guò)程中，一般會(huì)包含某種耗散機(jī)制。而雙曲型守恒律往往是在建立數(shù)學(xué)模型時(shí)忽略這種耗散

22、機(jī)制的結(jié)果。例如，我們前面說(shuō)過(guò)，對(duì)Navier-Stokes方程而言，激波是有一定厚度的連續(xù)結(jié)構(gòu)，但忽略粘性效應(yīng)以后的Euler方程，解中則可能出現(xiàn)理想的間斷。很明顯，Euler方程的物理解應(yīng)該是Navier-Stokes方程在粘性趨于零時(shí)得解。38這種選擇物理解的方法，雖然非常直觀，但并不實(shí)用，因?yàn)檎承詥?wèn)題的解一般也是未知的。對(duì)于標(biāo)量非線性雙曲型守恒律，Oleinik 提出了下面的熵條件：熵條件：設(shè) 是定義在t0上半平面上，存在有限條光滑間斷線的分片連續(xù)可微函數(shù);且是標(biāo)量守恒律： 的弱解0( .0)( )uftxu xx若 在間斷線附近滿足( , )u x t()( )()()()( )f

23、uf wf uf uf uf wuwuuuw其中 。 則弱解是唯一的，并且就是物理解。(min,max,)wIuuuu39()()()( )lim()()()()lim()()()()wuwuf uf udxDuudtf uf wwuuwf uf uf wuwuuuwf uuuuDu在熵條件中，是間斷面的傳播速度。在間斷左側(cè)，若令，則表示間斷左側(cè)的特征斜率。在右側(cè)，令，則表示間斷右側(cè)的特征線。所以對(duì)于標(biāo)量守恒律，Oleinik熵條件蘊(yùn)涵了40用這個(gè)條件考察Bergers方程的上式兩個(gè)弱解010( )10 xuxx1( , )/1xtu x tx ttxtxt ()101()uDu 在導(dǎo)數(shù)的間斷

24、線兩側(cè)滿足熵條件，因而是物理解！因此不是物理解熵條件有多種表達(dá)形式。對(duì)于標(biāo)量雙曲型守恒律，熵條件的研究已經(jīng)比較成熟。對(duì)于雙曲型守恒方程組，熵條件的研究還不夠充分。41含激波的流場(chǎng)的計(jì)算方法42包含激波的流場(chǎng)的數(shù)值計(jì)算方法可以大致分為兩類，即激波裝配（shock fitting）方法和激波捕捉（shock capturing）方法43激波裝配方法中，激波把求解域分為兩個(gè)或者多個(gè)子域，激波作為各子域激波裝配方法中，激波把求解域分為兩個(gè)或者多個(gè)子域，激波作為各子域邊界的一部分，激波的兩側(cè)通過(guò)邊界的一部分，激波的兩側(cè)通過(guò)R-H關(guān)系相聯(lián)系。關(guān)系相聯(lián)系。在每個(gè)子域內(nèi)，流場(chǎng)是光滑的，可以用任何計(jì)算可壓縮流動(dòng)

25、的數(shù)值方法求解；在激波處，用R-H關(guān)系和必要的補(bǔ)充條件計(jì)算出激波運(yùn)動(dòng)的速度，從而可以在計(jì)算中不斷更新激波的位置。激波裝配方法的優(yōu)點(diǎn)是：可以比較準(zhǔn)確地計(jì)算出激波位置和光滑區(qū)的流場(chǎng)。但缺激波裝配方法的優(yōu)點(diǎn)是：可以比較準(zhǔn)確地計(jì)算出激波位置和光滑區(qū)的流場(chǎng)。但缺點(diǎn)是計(jì)算過(guò)程比較復(fù)雜，不容易設(shè)計(jì)通用的計(jì)算軟件。點(diǎn)是計(jì)算過(guò)程比較復(fù)雜，不容易設(shè)計(jì)通用的計(jì)算軟件。顯然，如果在各個(gè)子域分別生成計(jì)算網(wǎng)格，則網(wǎng)格一般是隨激波運(yùn)動(dòng)的；如果各個(gè)子域共用一套靜止的網(wǎng)格（稱為浮動(dòng)激波裝配方法），則激波附近的計(jì)算格式要進(jìn)行特殊處理。另外，當(dāng)流場(chǎng)中激波的結(jié)構(gòu)比較復(fù)雜，或者存在激波的產(chǎn)生和消失等過(guò)程時(shí)，激波裝配方法將更難應(yīng)用。44

26、激波捕捉方法的基本思想是：在計(jì)算包含激波的流場(chǎng)時(shí)，采用統(tǒng)一的計(jì)算激波捕捉方法的基本思想是：在計(jì)算包含激波的流場(chǎng)時(shí)，采用統(tǒng)一的計(jì)算格式，不對(duì)激波進(jìn)行任何特殊處理。格式，不對(duì)激波進(jìn)行任何特殊處理。當(dāng)計(jì)算格式滿足一定要求時(shí)，可以自動(dòng)計(jì)算出流場(chǎng)中的間斷。當(dāng)計(jì)算格式滿足一定要求時(shí)，可以自動(dòng)計(jì)算出流場(chǎng)中的間斷。這種方法的優(yōu)點(diǎn)是：這種方法的優(yōu)點(diǎn)是：編程計(jì)算比較簡(jiǎn)單，且可以適用于具有任意激波結(jié)構(gòu)（包括激波的產(chǎn)生、消失、運(yùn)動(dòng)等各種情況）的流場(chǎng)計(jì)算。其缺點(diǎn)是：其缺點(diǎn)是：此時(shí)激波不再是理想的間斷，而是厚度為一至數(shù)個(gè)網(wǎng)格寬度的連續(xù)變化的結(jié)構(gòu)；另外，由于激波內(nèi)部流場(chǎng)的梯度較大，常規(guī)的計(jì)算格式或者會(huì)把激波“抹平”，或者在

27、激波附近產(chǎn)生虛假數(shù)值振蕩。因此，為了高分辨率、無(wú)振蕩的計(jì)算激波，對(duì)數(shù)值計(jì)算方法提出了很高的要求。事實(shí)上，能有效計(jì)算激波的高分辨率格式，是最近30年以來(lái)，計(jì)算流體力學(xué)研究的熱點(diǎn)問(wèn)題之一。隨著高分辨率格式研究的進(jìn)展，目前，激波捕捉方法已經(jīng)成為計(jì)算包含激波的可壓縮流場(chǎng)的主流方法。45一、守恒格式和Lax-Wendroff定理在求解守恒律0txUF 我們可以定義一類特殊的格式，稱為守恒格式。的有限差分和有限體積格式中, 定義：一維守恒律 式的差分或者有限體積格式0txUF11(,)nnnnjj pj pj q UG UUU稱為守恒型差分格式，如果11/21/2(,)()nnnnnnj pj pj qj

28、jjtx G UUUUFF且且1/211/21(,)(,)nnnjj pj qnnnjj pj q FF UUFF UU 等稱為數(shù)值通量。1/2njF46守恒格式的特點(diǎn)是： 對(duì)于任意的 ，數(shù)值通量用相同的函數(shù)關(guān)系計(jì)算，且函數(shù)涉及的自變量（離散點(diǎn)上的數(shù)值解）與 的相對(duì)位置關(guān)系不變。這一特點(diǎn)使得相鄰網(wǎng)格點(diǎn)，如 和 點(diǎn)上的計(jì)算格式在公共界面 處的數(shù)值通量是相同的。稍后我們將說(shuō)明，正是這個(gè)特點(diǎn)使守恒格式具有自動(dòng)捕捉間斷的能力。1/2j1/2jj1j1/2j47守恒格式的相容性條件守恒格式的相容性條件對(duì)于由 三個(gè)模板點(diǎn)構(gòu)成的三點(diǎn)格式，數(shù)值通量為：(1, ,1)jj j1/211/21(,)(,)nnnj

29、jjnnnjjjFF U UFF UU我們要求差分或者有限體積格式 與微分方程是相容的。下面我們以上面的三點(diǎn)格式為例，討論守恒格式相容的條件。顯然，三點(diǎn)格式可以寫(xiě)為：11(,)nnnnjj pj pj q UG UUU11/21/2(,)()nnnnnnj pj pj qjjjtx G UUUUFF111 (,)(,)0nnnnnnjjjjjjtxUUF U UF UU48則則111211 (,)(,)()()(,) (,)()(,) (,)()nnnnnnjjjjjjnnnjjjnnjjnnnnjjjjnjRnnjjnnnnjjjjnjLtxtOtttxtx UUF U UF UUUUUF

30、U UF U UUUUF U UF U UUUU其中 分別表示對(duì) 表達(dá)式中第一個(gè)變量和第二個(gè)變量的導(dǎo)數(shù)。(,)(,),nnnnjjjjnnjLjRF UUF UUUU49因此，因此，111222 (,)(,)()()(,)()()(,)()()(,)()(nnnnnnjjjjjjnjnnjjnnnjxjjnjRnnjjnnnjxjjnjLnnjjnjxnjRtxtOtttxOxxtxOxxtttx UUF U UF UUUF U UUUUUF U UUUUUF U UUUU222(,)() ()()( ,)()()()( ,)()()nnjjnnjxjnjLnnjjnjxOtOx tttxtx

31、xOtOx ttOtOx ttx F U UUUUF U UUF U U50整理上式，可得：整理上式，可得：1111( ,) (,)(,)()()nnnjjnnnnjjjjjOtOxtxtxUUUF U UF U UF UU當(dāng)數(shù)值通量滿足條件( ,)()F U UF U數(shù)值通量的相容條件數(shù)值通量的相容條件 1110,01lim (,)(,)0nnnjjnnnnjjjjtxjtxtx UUUFF U UF UU對(duì)上式兩側(cè)取極限，有對(duì)上式兩側(cè)取極限，有 即三點(diǎn)格式與守恒方程式即三點(diǎn)格式與守恒方程式 相容相容 0txUF在構(gòu)造數(shù)值通量的計(jì)算方法時(shí)，一般應(yīng)滿足這一條件。對(duì)于一般的守恒在構(gòu)造數(shù)值通量的計(jì)

32、算方法時(shí)，一般應(yīng)滿足這一條件。對(duì)于一般的守恒格式，數(shù)值通量的相容條件為：格式，數(shù)值通量的相容條件為： ( ,)( )F UUF U51任取一與間斷只有有限個(gè)交點(diǎn)并包含在, , LRbtxxt t之內(nèi)的封閉凸區(qū)域12( , ),( )( ),btx t x txx tttt 這個(gè)區(qū)域也可以等價(jià)地表示為12( , ), ( )( ),LRx tt xttx xxx 區(qū)域的邊界記為我們把區(qū)域, , LRbtxxt t剖分為矩形控制體， 在內(nèi)部的控制體的集合記為,12( , ),nnj nbtj njjjnnn或者等價(jià)地,12( , ),jjj nLRj n nnnjjj我們把11/21/2()()n

33、nnnjjjjxt UUFF在 上求和。, j n相容的守恒格式具有自動(dòng)捕捉激波和接觸間斷的能力。相容的守恒格式具有自動(dòng)捕捉激波和接觸間斷的能力。52221111/21/2()()jtRnnjLbnnjnnnnjjjjjjn nn nxt UUFF化簡(jiǎn)為為 212111/21/2()()tRjjnnLbnjnnnnjjjjjjn nxt UUFF0,0 xt 當(dāng)()RLjjx()tbnnt（但（但和和不變）時(shí)，不變）時(shí)， , j n區(qū)域區(qū)域的邊界趨近于的邊界趨近于。考慮到相容條件?？紤]到相容條件 ,可知上式趨于：( )0dtdxF UU由于 的任意性， 可以代表 平面的任意閉曲線；這樣，當(dāng) 時(shí)

34、，相容的守恒格式的數(shù)值解滿足弱解的定義。即當(dāng) 時(shí)，相容的守恒格式的解就是守恒律的廣義解或弱解。這一結(jié)論可以概括為下面的定理：, ,LRbtxxt t( , )x t0,0 xt 0,0 xt 由于守恒格式的數(shù)值通量的特點(diǎn)，內(nèi)部的數(shù)值通量可以相互抵消，從而有： 53Lax-Wendroff 定理定理:如果如果 是是守恒律守恒律 初值問(wèn)題的相容守恒格初值問(wèn)題的相容守恒格式的離散的數(shù)值解，且式的離散的數(shù)值解，且當(dāng)當(dāng) 時(shí)，時(shí)， 在在某種范數(shù)的意義下趨于某種范數(shù)的意義下趨于 ，則則 是守恒律的一個(gè)弱解。是守恒律的一個(gè)弱解。njU0txUF0,0 xt njU( , )x tV( , )x tVV提示（1

35、）Lax-Wendroff 定理并不保證弱解的存在性和相容的守恒格式的收斂性。同時(shí)，該定理也并不保證得到的弱解滿足熵條件。但是，通過(guò)該定理我們可以確信：用相容的守恒格式可以自動(dòng)捕捉激波和接觸間斷。（2）Lax-Wendroff 定理給出了用數(shù)值方法計(jì)算弱解的一個(gè)充分條件，也就是說(shuō)，該定理并未排除用非守恒型格式計(jì)算出弱解的可能性。事實(shí)上，在某些情況下，非守恒格式也可以得到弱解的良好近似。但是，非守恒型格式是否可以計(jì)算弱解，對(duì)具體的格式要進(jìn)行具體分析，而這種分析通常是復(fù)雜的。所以，我們一般用守恒格式計(jì)算包含間斷的流場(chǎng)。54二、人工粘性和格式粘性在計(jì)算包含激波的流場(chǎng)時(shí)，考察計(jì)算格式的“數(shù)值粘性”特征

36、是非常重要的。這主要出于兩個(gè)原因：（1）物理解是所謂“粘性消失解”，所以，在計(jì)算包含激波間斷的流場(chǎng)時(shí)，計(jì)算格式包含足夠的數(shù)值粘性是必要的；（2）某些計(jì)算格式在間斷附近會(huì)產(chǎn)生非物理振蕩，合理地調(diào)整計(jì)算方法的數(shù)值粘性，對(duì)于消除計(jì)算結(jié)果在間斷附近的數(shù)值振蕩是有益的。55 在我們已經(jīng)介紹過(guò)線性波動(dòng)方程 的修正方程。修正方程的偶數(shù)階導(dǎo)數(shù)項(xiàng)影響數(shù)值解的振幅，當(dāng) 等關(guān)系成立時(shí)，偶數(shù)階導(dǎo)數(shù)項(xiàng)起耗散作用，相當(dāng)于計(jì)算格式中的數(shù)值粘性。某些差分格式本身就帶有數(shù)值粘性這部分?jǐn)?shù)值粘性稱為“格式粘性格式粘性”。某些格式缺乏必要的格式粘性（如中心差分格式）或者格式粘性不足以抑制間斷附近的非物理振蕩時(shí)，常常需要在計(jì)算格式中添

37、加顯式的數(shù)值粘性項(xiàng) ，稱為“人工粘性人工粘性”。根據(jù)修正方程的性質(zhì)，人工粘性項(xiàng)一般可以寫(xiě)成 或者 的離散形式。 等稱為人工粘性系數(shù)。uuaotx240,0222ux444ux240,0為了保證格式是相容的，應(yīng)該有：240,00,0lim0,lim0 xtxt 在能夠保證計(jì)算穩(wěn)定和抑制間斷附近振蕩的前提下在能夠保證計(jì)算穩(wěn)定和抑制間斷附近振蕩的前提下, 人工粘性人工粘性和格式粘性越小越好和格式粘性越小越好!56中心型格式需要添加人工粘性中心型格式需要添加人工粘性;迎風(fēng)型格式自帶格式粘性迎風(fēng)型格式自帶格式粘性!先考慮中心格式先考慮中心格式57考慮Euler顯式格式11102nnnnjjjjuuuua

38、tx其修正方程為：這是一個(gè)一階精度的格式，二階耗散項(xiàng)對(duì)數(shù)值解的幅值有主要影響。該格式二階耗散項(xiàng)是負(fù)的，即格式粘性為負(fù)，因而是無(wú)條件不穩(wěn)定的。為了使計(jì)算穩(wěn)定，我們可以在格式的右端添加一個(gè)二階導(dǎo)數(shù)項(xiàng)，如22( 231)26txxxxxxa xca xuauuccu 11111222nnnnnnnjjjjjjjuuuuuuuatxt此格式就是前面介紹過(guò)的Lax格式，其修正方程為：221()( 231)26txxxxxxa xa xuauc uccuc 則當(dāng) 時(shí)，格式是穩(wěn)定的。1c 5811222nnnjjjxxuuua xutc注意到： 而且上述兩個(gè)格式修正方程之間的2xxa xuc，即 差的主項(xiàng)正

39、好為22()()2LaxxxEulerxxxxa xuuuc顯然在Euler顯式格式的右端增加 1122nnnjjjuuut相當(dāng)于增加了Euler顯式格式的耗散，所以這一項(xiàng)通常稱為 Euler顯式格式的“人工粘性”。 由于添加了這種人工粘性，使得格式是穩(wěn)定的，而且間斷會(huì)被抹平。推廣到Euler方程，這種添加了人工粘性的Euler顯式格式（實(shí)際上就是Lax格式），可以寫(xiě)為11111222nnnnnnnjjjjjjjtxtUUFFUUU上述格式可以改寫(xiě)為：11/21/2()nnnnjjjjtxUUFF其中1/2111/2111()()221()()22nnnnnjjjjjnnnnnjjjjjxtx

40、tFFFUUFFFUU所以，Lax格式是守恒格式，可以自動(dòng)捕捉流場(chǎng)中的間斷。59線性波動(dòng)方程的Lax-Wendroff格式為1211112222nnnnjjjjnnnjjjuuuuatauuutxx。其修正方程為：2322(1)(1)68txxxxxxxxa xa xuauc ucc u 可以看到，這是一個(gè)二階格式，其格式粘性是四階的。四階格式粘性通常不足以抑制間斷附近的振蕩。為了解決這一問(wèn)題，我們也可以采取添加人工粘性的方法。例如，可以把差分格式改寫(xiě)為：12112112()222nnnnjjjjnnnjjjuuuuatauuutxx。其中 為二階粘性系數(shù)。如果 ，雖然可以有效的消除間斷附近的

41、振蕩，但是，間斷被嚴(yán)重抹平，而且格式變?yōu)橐浑A精度，在光滑區(qū)的精度也大為降低。因此，這種方法顯然不是一個(gè)好的選擇。2212 t60比較好的取法是取112112nnnjjjnnnnjjjjuuuax uuuu其中 為量級(jí)為 的正數(shù)， 是一個(gè)小的正數(shù)，如 ，以防止出現(xiàn)分母為零的情況。對(duì)于這種取法，在光滑區(qū)，有 ，則人工粘性項(xiàng)為 的量級(jí)，從而不會(huì)影響光滑區(qū)格式的精度。在間斷附近， ，人工粘性為 的量級(jí)， 激波附近的振蕩可以得到有效控制。注意到，在人工粘性項(xiàng)中一個(gè)自由參數(shù) ，一般可取1/2，也可以根據(jù)試算調(diào)整 ，以達(dá)到更好的效果。但是， 對(duì)于不同的問(wèn)題，最佳的 值是不同的。(1)O6102(1)O222

42、()uOxx2(1/)Ox22()uOxx 在添加上述人工粘性項(xiàng)以后，我們并不能從理論上保證間斷附近解沒(méi)有振蕩。而且由于可調(diào)參數(shù) 的存在，使得這種方法的實(shí)施在一定程度上依賴于經(jīng)驗(yàn)。 61 推廣到Euler方程，有111221/21/211/21/21222()()()()22nnnnjjjjnnnnnnjjjjjjjjtxttAAxxUUFFUUUU11/21/2()nnnnjjjjtxUUFF其中，21/211/21/2121/211/21/211()()()221()()()22nnnnnnjjjjjjjnnnnnnjjjjjjjtAxxtAxxFFFUUFFFUU11111/212max

43、(,)nnnjjjjjnnnnjjjjjjjpppx pppp 是Jacobi矩陣， 是壓力， 為Jacobi矩陣 中所有特征值中絕對(duì)值最大者。注意到，為了保證格式是守恒格式，Euler方程人工粘性的形式與線性波動(dòng)方程略有不同。ApFU62中心型格式需要添加人工粘性中心型格式需要添加人工粘性;迎風(fēng)型格式自帶格式粘性迎風(fēng)型格式自帶格式粘性!再考慮迎風(fēng)格式格式再考慮迎風(fēng)格式格式630uuatx(0)a 110nnnnjjjjuuuuatx110nnnnjjjjuuuuatx0a 一階迎風(fēng)格式一階迎風(fēng)格式修正方程: 22(1)( 231)26txxxxxxa xa xuauc uccu(0)a 22

44、(1)( 231)26txxxxxxa xa xuauc uccu0a 11111222nnnnnnnjjjjjjjuuuuuuuatxt LAX格式格式221()( 231)26txxxxxxa xa xuauc uccuc 修正方程: 兩種格式的數(shù)值粘性都足以抑制間斷附近的振蕩兩種格式的數(shù)值粘性都足以抑制間斷附近的振蕩, 但一階迎風(fēng)格式的數(shù)值粘但一階迎風(fēng)格式的數(shù)值粘性較小性較小!迎風(fēng)格式如何應(yīng)用于方程組迎風(fēng)格式如何應(yīng)用于方程組?64四、一維線性波動(dòng)方程組的迎風(fēng)格式考慮線性雙曲型方程組0AAconsttxUU由于Jacobi矩陣是常數(shù)矩陣，上式也可以寫(xiě)為守恒形式0,AtxUFFU由于雙曲性假

45、設(shè)，矩陣 可以對(duì)角化：AAR Li jiij其中 是 的特征值， 是由特征值組成的對(duì)角陣， 和 分別是 的左特征向量和右特征向量組成的矩陣，且二者互為逆矩陣：iALRALRI65令LWU其中W 稱為特征變量，則有0txWW寫(xiě)成分量形式0iiiwwtx其中iw為W的第i個(gè)分量。 0iiiwwtx在形式上與線性波動(dòng)方程相同 因此，線性波動(dòng)方程的迎風(fēng)格式可以直接應(yīng)用于 0iiiwwtx以一階迎風(fēng)格式為例，有： 111()()()()()()022nnnnnnijijijijijijiiiiwwwwwwtxx66寫(xiě)成向量形式，有111()()nnnnnnjjjjjjtxWWWWWW其中22 是特征值的

46、絕對(duì)值組成的對(duì)角陣。 ,ii j 111()()nnnnnnjjjjjjtxWWWWWW左右同乘以R 寫(xiě)成守恒變量的形式，有： 111()()nnnnnnjjjjjjtAAxUUUUUU67其中22AAARLAAARLARL111()()nnnnnnjjjjjjtAAxUUUUUU即為線性雙曲型方程組式 的一階迎風(fēng)格式。 0AAconsttxUU 注意到 的特征值非負(fù)， 的特征值非正，且 。AAAAA68所以，線性雙曲型方程組的迎風(fēng)格式的構(gòu)造方法可以歸納為： AAAAAA把系數(shù)矩陣分裂為，其中的特征值非負(fù)，的特征值非正。 把方程 改寫(xiě)為0AAconsttxUU0AAtxxUUU與 對(duì)應(yīng)的 后差

47、離散， 對(duì)應(yīng)的 用前差離散。當(dāng) 依照上述原則采用高階差分近似時(shí)，可以得到高階的迎風(fēng)格式。AxUAxUxU690txxUFFxFxF其中用后差離散，而用前差離散，則可得到雙曲型守恒方程組的迎風(fēng)格式。注意到，對(duì)于常系數(shù)的雙曲型方程組，有AFUAAAFUUUFF70111/21/21/2111111100022()()221()()22jjjjjjjjjjjjjjjjjjjjjAAtxxAAtxxtxAAFVSAAAAAAARoeUUUUUUUUFFUFUUFFUUUUUUFFUUAAAFUUUFF22AAARLAAARLARLSCM：Splitting Coefficient Matrix71五、

48、五、 Euler方程的迎風(fēng)型有限差分格式方程的迎風(fēng)型有限差分格式考慮一維Euler方程0txUFEuler方程是非線性的方程，如何把上面所述迎風(fēng)格式推廣到求解非線性雙曲型方程組特別是Euler方程呢？我們下面介紹幾種方法。721. SCM方法0AtxUU 和線性情況不同，上式中的矩陣A不是常數(shù)，而是U的函數(shù)。但是，我們?nèi)匀豢梢越梃b線性雙曲型方程組迎風(fēng)格式的思路，把矩陣A進(jìn)行分裂。其做法與線性的情況完全類似，唯一的區(qū)別是由于A在各個(gè)網(wǎng)格點(diǎn)是不同的，所以，在每個(gè)網(wǎng)格點(diǎn)均需進(jìn)行矩陣分裂。這種方法稱為分裂系數(shù)矩陣方法（Split Coefficient Matrix Scheme），簡(jiǎn)稱SCM方法。

49、把Euler方程寫(xiě)成擬線性形式應(yīng)用SCM方法的一階迎風(fēng)格式為：111()()nnnnnnjjjjjjjjtAAxUUUUUU， ,22nnnnjjjjnnnjjjjjAAAAAAARL73()( )FUF UAFU2. FVS方法FVS是矢通量分裂（Flux Vector Splitting） 的簡(jiǎn)稱，最早由Steger-Warming提出。我們注意到，雖然SCM方法可以看作線性雙曲型方程組迎風(fēng)格式的直接推廣，但是卻不是守恒格式。為了構(gòu)造守恒的迎風(fēng)格式，Steger 和Warming注意到Euler方程的通量F是守恒變量U的一次齊次函數(shù)，即滿足在這種條件下，有()( )()()()AFUF U

50、FFUFUUFUUUUUAFU這一關(guān)系可證明如下：則74AAAFUUUFF,22AAAAAAAAARLFU FU0txxUFF利用這一性質(zhì)，我們可以把通量分解為：其中所以，Euler方程可以改寫(xiě)為：xFxF111()()nnnnnnjjjjjjtxUUFFFF同線性雙曲型方程組迎風(fēng)格式一樣，用后差離散，而用前差離散。 一階迎風(fēng)格式為7511/21/20nnjjjjtxUUFF1/211/21,nnnnnnjjjjjjFFFFFF1/2j jF1jF1/2j 寫(xiě)成守恒格式，其數(shù)值通量為：也就是說(shuō)，考慮到波的傳播方向，左側(cè)的和右側(cè)的對(duì)處的數(shù)值通量有貢獻(xiàn)。 maxmax,22AAAIAAIAFU FUmaxI注意到，如果通量F不是守恒變量U的一次齊次函數(shù)，則不能構(gòu)造相應(yīng)的FVS方法。另外，通量分裂的方法不是唯一的。除了式所示的分裂方法外，另一種簡(jiǎn)單的分裂方法為其中：為特征值絕對(duì)值的最大值，是單位矩陣。763. Roe方法Roe方法是由Roe在1981年提出的。 的守恒型差分格式為：0txUF11/21/20nnjjjjtxUUFF下面討論通量 的算法。首先，我們通過(guò) 和 給出 的估計(jì)值 。 代表某種平均算子，其具體形式下面討論。則我們可以算出與 對(duì)應(yīng)的特征值 以及左右特征向量組成的矩陣：1/2njF1njUnjU1/2njU1/21(,)nnnjjjBUUUB1