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

1、求解流體力學(xué)方程組的有限差分方法,第三講,1,控制方程,,2,貼體網(wǎng)格 H型C型O型,,,3,一、貼體坐標(biāo)中的基本方程,,,正變換,逆變換,假定已知上述變換的離散形式即已經(jīng)生成好貼體網(wǎng)格,4,導(dǎo)數(shù)的變換,一階導(dǎo)數(shù),,5,二階導(dǎo)數(shù),,,,,,6,二階導(dǎo)數(shù),,,,,,度量系數(shù)metrics,,,,,,7,度量系數(shù)及其計算方法,,困難,容易,,,,,,如果能夠找到 和

2、 之間的關(guān)系，就可以得到 等的計算方法,8,,,,,,,,,,,,,,,,Jacobi行列式（jacobian）,,,,9,,,,,,,,,,,,10,,,,11,,12,任意曲線坐標(biāo)系中流體力學(xué)方程組的守恒形式,,,13,,,14,,,,15,,,曲線坐標(biāo)系中Navier－Stokes方程的強(qiáng)守恒形式,16,二、守恒型Euler方程,在任意曲線坐

3、標(biāo)系中，二維守恒型Euler方程為：,,其中,,,,,17,二維守恒型Euler方程可以改寫為擬線性的形式：,其中,,令 或 ，則 可以寫為統(tǒng)一形式：,,,,18,其中， 為絕熱指數(shù)，對于完全氣體， ， 和 的定義為：,,,,,的特征值為：,,,其中 是音速。,,19,三、守恒型的Navier－Stokes方程,平面上Navier－Stokes方程的強(qiáng)守恒形

4、式為：,,,其中,,,,注意到，Euler方程可以看作Navier－Stokes方程的特例。一般而言，Euler方程中無粘通量 的離散方法，同樣可用于Navier－Stokes方程；而為了求解Navier－Stokes方程，我們還必須提供粘性通量 的離散方法。,,20,有限差分方法,,21,有限差分方法,對

5、 式直接進(jìn)行半離散差分近似,注意到 ，采用中心差分，則有：,,,,,,令，,則可以得到守恒格式：,22,,在有限差分方法中，在網(wǎng)格點上的度量系數(shù)(如 )也需要通過差分方法計算，因此要求計算網(wǎng)格是充分光滑的。,,,23,有限差分方法的主要特點是：1）有限差分方法

6、只需構(gòu)造偏導(dǎo)數(shù)的離散方法，這使得它比較容易推廣到高階精度，對于多維問題也是如此。2）有限差分方法對網(wǎng)格的光滑性有較高要求，不易在復(fù)雜形狀求解域上實施。3）在曲線坐標(biāo)系中，有限差分方法要對幾何量和物理量的確定組合進(jìn)行差分離散，這樣做的后果之一是有限差分方法可能產(chǎn)生所謂幾何誘導(dǎo)誤差（geometrically induced errors） 。例如，考慮一無界均勻流場，易知，該流動應(yīng)是定常的。即

7、 。在由一般曲線坐標(biāo)構(gòu)成的網(wǎng)格上采用有限差分方法計算此流動，可能有 （三維以及動網(wǎng)格中的有限差分方法更容易出現(xiàn)這種情況）。而采用有限體積方法，則恒有 。有限差分方法在某些情況下不能復(fù)現(xiàn)均勻流場，是幾何誘導(dǎo)誤差的表現(xiàn)之一。,,,,24,如何計算數(shù)值通量？,數(shù)值耗散數(shù)值色散,

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

9、場的計算問題。在數(shù)值求解Navier－Stokes方程時，雖然理論上激波附近物理量是連續(xù)的，但是在通常情況下，激波的厚度非常薄（為分子平均自由程的量級），遠(yuǎn)遠(yuǎn)小于計算網(wǎng)格的尺度；所以在計算中，我們?nèi)孕璋鸭げó?dāng)作間斷來處理。,在Euler方程的范疇內(nèi)如何表示含激波間斷的流場?,27,一、激波的形成,激波是一種非線性現(xiàn)象?？紤]線性對流方程的初值問題，如果其初始值是連續(xù)的，則解中不會自發(fā)產(chǎn)生間斷。但是，對于非線性問題，即使初值是連續(xù)

10、的，在后續(xù)時間中，也可能發(fā)展出包含間斷的解。Bergers方程可以看作Euler方程的非線性模型方程；和Euler方程一樣，是雙曲型方程。下面，我們以Bergers方程為例，介紹間斷解或者激波的形成過程。,,,,,28,考慮Bergers方程,的初值問題，初始條件為 。根據(jù)特征理論，Bergers方程的解析解可以寫為下列形式：,,初始條件,,Bergers方程的解為

11、：,,29,圖1 Bergers方程初值問題的解,圖1顯示了解隨時間的發(fā)展過程。容易知道，當(dāng) 時，Bergers方程的解為,,,可見,對于非線性的Bergers方程,即使初始值是連續(xù)的,其解仍然可能出現(xiàn)間斷。,事實上，一般的非線性雙曲型守恒方程（組）或稱非線性雙曲型守恒律的解中都有可能存在間斷。通常，我們認(rèn)為偏微分方程的解是連續(xù)可微的。顯然，對于非線性雙曲型守恒律，這一點并不成立。因此，必須拓展雙曲型守恒律

12、解的概念。,30,二、廣義解,對于一般的一維雙曲型守恒律 的初值問題，我們引入廣義解（generalized solution）或稱弱解(weak solution)的概念，來表示包含間斷的解。,,定義：設(shè) 是分片連續(xù)可微的函數(shù)，在 的半平面，如果對于與 的間斷線只有有限個交點的任意分段光滑的閉曲線 ，

13、都有: 則稱 為方程 在初值 下的廣義解或弱解。,,,,,,,,,31,下面討論一下廣義解的意義。如果已知 是光滑的，設(shè) 圍成的 區(qū)域為 ，則由 式利用

14、Green公式知,,,,,由于閉曲線可以在光滑區(qū)內(nèi)任取，可得：,,在光滑區(qū)，弱解就是通常的連續(xù)可微解。,,,,,,,,,,,32,,,,,,如果 是由一條間斷線 分隔開的分片連續(xù)可微函數(shù)，取如圖所示的由 （ 是正數(shù)）圍成的閉曲線 ，在 上應(yīng)用

15、,有:,,,,,,,33,間斷關(guān)系的推導(dǎo),,令 ，并考慮到在 兩側(cè) 有間斷，則上式可簡化為：,,,,分別表示從左、右趨近間斷線時守恒變量的值。,,,令,,,在間斷兩側(cè)，顯然 。 考慮到 可以任意取值，可得：,,其中,Rankine-Hogoniot（R-H）關(guān)系,34,以一維Euler方程為例，R-H關(guān)系為：,

16、,注意到，當(dāng) 時，如果 ，則上面的關(guān)系成立，此時只有密度場存在間斷。這種類型的間斷稱為接觸間斷。不屬于接觸間斷的間斷稱為激波。,,,綜上所述，雙曲型守恒律的廣義解 是被有限個間斷線分開的分片光滑函數(shù)。在光滑區(qū)， 是滿足微分方程式 的連續(xù)可微解

17、；在間斷線的兩側(cè)， 滿足R-H關(guān)系。,,,,,,35,三、熵條件,廣義解推廣了偏微分方程初值問題連續(xù)可微解的概念，但其后果是導(dǎo)致了弱解不唯一。為了說明這一問題，我們舉一個例子：考慮Burgers方程在初值為 時的解,此時，初值在 處有一個間斷。,,把Bergers方程寫為守恒形式,,可知

18、 處的Rankine-Hogoniot條件為：,,,由上式知 。,,,在間斷處滿足Rankine-Hogoniot條件，在其他地方滿足微分方程，即 是Bergers方程的一個廣義解。,,36,另外，考慮,,分片連續(xù)可微函數(shù)，其導(dǎo)數(shù)在 處存在間斷(或者認(rèn)為函數(shù)值存在幅度為零的間斷)。,,根據(jù)廣義解的定義，容易驗證

19、此函數(shù)也是Bergers方程的一個廣義解。所以廣義解一般不唯一，但是對于有明確物理意義的守恒律，其中只有一個解是有物理意義的，我們稱之為物理解。為了得到我們關(guān)心的物理解，廣義解除了必須滿足式外，還必須滿足附加的條件，這個條件因為與熱力學(xué)第二定律所起的作用相同，被稱為熵條件（entropy condition）。,,,,,,37,如何在不唯一的廣義解中挑選所謂物理解呢？最直觀的方法是認(rèn)為物理解是一個光滑的粘性問題的解在粘性趨于零

20、時的極限。也就是說：,方程 的解如果當(dāng) 時，幾乎處處有界的收斂到分片連續(xù)可微函數(shù) ，則 是 的物理解。,這種方法的物理基礎(chǔ)是：在真實的物理過程中，一般會包含某種耗散機(jī)制。而雙曲型守恒律往往是在建立數(shù)學(xué)模型

21、時忽略這種耗散機(jī)制的結(jié)果。例如，我們前面說過，對Navier-Stokes方程而言，激波是有一定厚度的連續(xù)結(jié)構(gòu)，但忽略粘性效應(yīng)以后的Euler方程，解中則可能出現(xiàn)理想的間斷。很明顯，Euler方程的物理解應(yīng)該是Navier-Stokes方程在粘性趨于零時得解。,38,這種選擇物理解的方法，雖然非常直觀，但并不實用，因為粘性問題的解一般也是未知的。對于標(biāo)量非線性雙曲型守恒律，Oleinik 提出了下面的熵條件：,熵條件：設(shè) 是定

22、義在t≥0上半平面上，存在有限條光滑間斷線的分片連續(xù)可微函數(shù);且是標(biāo)量守恒律： 的弱解,,若 在間斷線附近滿足,,,其中 。 則弱解是唯一的，并且就是物理解。,,39,40,,,,41,含激波的流場的計算方法,,42,,,包含激波的流場的數(shù)值計算方法可以

23、大致分為兩類，即激波裝配（shock fitting）方法和激波捕捉（shock capturing）方法,43,,激波裝配方法中，激波把求解域分為兩個或者多個子域，激波作為各子域邊界的一部分，激波的兩側(cè)通過R-H關(guān)系相聯(lián)系。在每個子域內(nèi)，流場是光滑的，可以用任何計算可壓縮流動的數(shù)值方法求解；在激波處，用R-H關(guān)系和必要的補(bǔ)充條件計算出激波運動的速度，從而可以在計算中不斷更新激波的位置。激波裝配方法的優(yōu)點是：可以比較準(zhǔn)確地計算出

24、激波位置和光滑區(qū)的流場。但缺點是計算過程比較復(fù)雜，不容易設(shè)計通用的計算軟件。顯然，如果在各個子域分別生成計算網(wǎng)格，則網(wǎng)格一般是隨激波運動的；如果各個子域共用一套靜止的網(wǎng)格（稱為浮動激波裝配方法），則激波附近的計算格式要進(jìn)行特殊處理。另外，當(dāng)流場中激波的結(jié)構(gòu)比較復(fù)雜，或者存在激波的產(chǎn)生和消失等過程時，激波裝配方法將更難應(yīng)用。,,44,激波捕捉方法的基本思想是：在計算包含激波的流場時，采用統(tǒng)一的計算格式，不對激波進(jìn)行任何特殊處理。當(dāng)

25、計算格式滿足一定要求時，可以自動計算出流場中的間斷。這種方法的優(yōu)點是：編程計算比較簡單，且可以適用于具有任意激波結(jié)構(gòu)（包括激波的產(chǎn)生、消失、運動等各種情況）的流場計算。其缺點是：此時激波不再是理想的間斷，而是厚度為一至數(shù)個網(wǎng)格寬度的連續(xù)變化的結(jié)構(gòu)；另外，由于激波內(nèi)部流場的梯度較大，常規(guī)的計算格式或者會把激波“抹平”，或者在激波附近產(chǎn)生虛假數(shù)值振蕩。因此，為了高分辨率、無振蕩的計算激波，對數(shù)值計算方法提出了很高的要求。事實上，

26、能有效計算激波的高分辨率格式，是最近30年以來，計算流體力學(xué)研究的熱點問題之一。隨著高分辨率格式研究的進(jìn)展，目前，激波捕捉方法已經(jīng)成為計算包含激波的可壓縮流場的主流方法。,45,一、守恒格式和Lax-Wendroff定理,在求解守恒律,,我們可以定義一類特殊的格式，稱為守恒格式。,的有限差分和有限體積格式中,,定義：一維守恒律 式的差分或者有限體積格式,,稱為守恒型差分格

27、式，如果,,且,,等稱為數(shù)值通量。,,46,守恒格式的特點是： 對于任意的 ，數(shù)值通量用相同的函數(shù)關(guān)系計算，且函數(shù)涉及的自變量（離散點上的數(shù)值解）與 的相對位置關(guān)系不變。這一特點使得相鄰網(wǎng)格點，如 和 點上的計算格式在公共界面 處的數(shù)值通量是相同的。稍后我們將說明，正是這個特點使守恒格式具有自動捕捉間斷的能

28、力。,,,,,,,,,,47,守恒格式的相容性條件,48,,,,49,因此，,,50,整理上式，可得：,,當(dāng)數(shù)值通量滿足條件,,數(shù)值通量的相容條件,,,,51,相容的守恒格式具有自動捕捉激波和接觸間斷的能力。,52,,,,,,,,由于 的任意性， 可以代表 平面的任意閉曲線；這樣，當(dāng) 時，相容的守

29、恒格式的數(shù)值解滿足弱解的定義。即當(dāng) 時，相容的守恒格式的解就是守恒律的廣義解或弱解。這一結(jié)論可以概括為下面的定理：,由于守恒格式的數(shù)值通量的特點，內(nèi)部的數(shù)值通量可以相互抵消，從而有：,53,,,,,,,Lax-Wendroff 定理:如果 是守恒律 初值問題的相容守恒格式的離散的數(shù)值解，且當(dāng)

30、 時， 在某種范數(shù)的意義下趨于 ，則 是守恒律的一個弱解。,提示（1）Lax-Wendroff 定理并不保證弱解的存在性和相容的守恒格式的收斂性。同時，該定理也并不保證得到的弱解滿足熵條件。但是，通過該定理我們可以確信：用相容的守恒格式可以自動捕捉激波和接觸間斷。（2）Lax-Wendrof

31、f 定理給出了用數(shù)值方法計算弱解的一個充分條件，也就是說，該定理并未排除用非守恒型格式計算出弱解的可能性。事實上，在某些情況下，非守恒格式也可以得到弱解的良好近似。但是，非守恒型格式是否可以計算弱解，對具體的格式要進(jìn)行具體分析，而這種分析通常是復(fù)雜的。所以，我們一般用守恒格式計算包含間斷的流場。,54,二、人工粘性和格式粘性,在計算包含激波的流場時，考察計算格式的“數(shù)值粘性”特征是非常重要的。這主要出于兩個原因：（1）物理解是所謂“粘

32、性消失解”，所以，在計算包含激波間斷的流場時，計算格式包含足夠的數(shù)值粘性是必要的；（2）某些計算格式在間斷附近會產(chǎn)生非物理振蕩，合理地調(diào)整計算方法的數(shù)值粘性，對于消除計算結(jié)果在間斷附近的數(shù)值振蕩是有益的。,55,為了保證格式是相容的，應(yīng)該有：,在能夠保證計算穩(wěn)定和抑制間斷附近振蕩的前提下, 人工粘性和格式粘性越小越好!,56,,,,,中心型格式需要添加人工粘性;迎風(fēng)型格式自帶格式粘性!,先考慮中心格式……,57,考慮Euler顯

33、式格式,,其修正方程為：,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,這是一個一階精度的格式，二階耗散項對數(shù)值解的幅值有主要影響。該格式二階耗散項是負(fù)的，即格式粘性為負(fù)，因而是無條件不穩(wěn)定的。為了使計算穩(wěn)定，我們可以在格式的右端添加一個二階導(dǎo)數(shù)項，如,,,,,,,,,,,,,此格式就是前面介紹過的Lax格式，其修正方程為：,,則當(dāng) 時，格式是穩(wěn)定的。,,,58,注意到：,而且上述兩個格式修正方

34、程之間的,，即,差的主項正好為,顯然在Euler顯式格式的右端增加,,相當(dāng)于增加了Euler顯式,格式的耗散，所以這一項通常稱為,Euler顯式格式的“人工粘性”。,由于添加了這種人工粘性，使得格式是穩(wěn)定的，而且間斷會被抹平。,推廣到Euler方程，這種添加了人工粘性的Euler顯式格式（實際上就是Lax格式），可以寫為,,上述格式可以改寫為：,其中,,所以，Lax格式是守恒格式，可以自動捕捉流場中的間斷。,59,線性波動方程的L

35、ax-Wendroff格式為,,其修正方程為：,可以看到，這是一個二階格式，其格式粘性是四階的。四階格式粘性通常不足以抑制間斷附近的振蕩。為了解決這一問題，我們也可以采取添加人工粘性的方法。,例如，可以把差分格式改寫為：,,其中 為二階粘性系數(shù)。如果 ，雖然可以有效的消除間斷附近的振蕩，但是，間斷被嚴(yán)重抹平，而且格式變?yōu)橐浑A精度，在光滑區(qū)的精度也大為降低。因此，這種方法顯然不是一個好的選

36、擇。,,,60,比較好的取法是取,,其中 為量級為 的正數(shù)， 是一個小的正數(shù)，如 ，以防止出現(xiàn)分母為零的情況。對于這種取法，在光滑區(qū)，有 ，則人工粘性項為 的量級，從而不會影響光滑區(qū)格式的精度。在間斷附近， ，人工粘性為

37、 的量級， 激波附近的振蕩可以得到有效控制。注意到，在人工粘性項中一個自由參數(shù) ，一般可取1/2，也可以根據(jù)試算調(diào)整 ，以達(dá)到更好的效果。但是， 對于不同的問題，最佳的 值是不同的。,,,,,,在添加上述人工粘性項以后，我們并不能從理論上保證間斷附近解沒有振蕩。而且由于可調(diào)參數(shù) 的存在，使得這種方法的實施在一定程度上依賴于經(jīng)驗。,,61,,,,推廣到Euler方程，有,其中，,,是Jacobi

38、矩陣， 是壓力， 為Jacobi矩陣 中所有特征值中絕對值最大者。注意到，為了保證格式是守恒格式，Euler方程人工粘性的形式與線性波動方程略有不同。,,62,,,,,中心型格式需要添加人工粘性;迎風(fēng)型格式自帶格式粘性!,再考慮迎風(fēng)格式格式……,63,迎風(fēng)格式如何應(yīng)用于方程組?,64,四、一維線性波動方程組的迎風(fēng)格式,考慮線性雙曲型方程組,由于Jacobi矩陣是常數(shù)矩陣，上式也可以寫為守恒形式,,由于雙曲性假設(shè)，矩

39、陣 可以對角化：,,,,其中 是 的特征值， 是由特征值組成的對角陣， 和 分別是 的左特征向量和右特征向量組成的矩陣，且二者互為逆矩陣：,,,,,,,,65,令,其中W 稱為特征變量，則有,寫成分量形式,,其中,為W的第i個分量。,,在形式上與線性波動方程相同,因此，線性波動方程的迎風(fēng)格式可以直接應(yīng)用于,以一階迎風(fēng)格式為例，有：,,66,寫成向量形式，有,,其中,,是

40、特征值的絕對值組成的對角陣。,,,左右同乘以,,寫成守恒變量的形式，有：,,67,,,,,,,68,,,,,,,,69,注意到，對于常系數(shù)的雙曲型方程組，有,70,,,SCM：Splitting Coefficient Matrix,71,五、 Euler方程的迎風(fēng)型有限差分格式,考慮一維Euler方程,Euler方程是非線性的方程，如何把上面所述迎風(fēng)格式推廣到求解非線性雙曲型方程組特別是Euler方程呢？我們下面介紹幾種方法。,72,

41、1. SCM方法,和線性情況不同，上式中的矩陣A不是常數(shù)，而是U的函數(shù)。但是，我們?nèi)匀豢梢越梃b線性雙曲型方程組迎風(fēng)格式的思路，把矩陣A進(jìn)行分裂。其做法與線性的情況完全類似，唯一的區(qū)別是由于A在各個網(wǎng)格點是不同的，所以，在每個網(wǎng)格點均需進(jìn)行矩陣分裂。這種方法稱為分裂系數(shù)矩陣方法（Split Coefficient Matrix Scheme），簡稱SCM方法。,把Euler方程寫成擬線性形式,應(yīng)用SCM方法的一階迎風(fēng)格式為：,，,,73,

42、2. FVS方法FVS是矢通量分裂（Flux Vector Splitting） 的簡稱，最早由Steger-Warming提出。我們注意到，雖然SCM方法可以看作線性雙曲型方程組迎風(fēng)格式的直接推廣，但是卻不是守恒格式。為了構(gòu)造守恒的迎風(fēng)格式，Steger 和Warming注意到Euler方程的通量F是守恒變量U的一次齊次函數(shù)，即滿足,在這種條件下，有,這一關(guān)系可證明如下：,則,74,利用這一性質(zhì)，我們可以把通量分解為：,其中,所以，

43、Euler方程可以改寫為：,同線性雙曲型方程組迎風(fēng)格式一樣，,用后差離散，而,用前差離散。,,一階迎風(fēng)格式為,75,76,3. Roe方法,Roe方法是由Roe在1981年提出的。 的守恒型差分格式為：,,下面討論通量 的算法。首先，我們通過 和 給出 的估計值

44、 。 代表某種平均算子，其具體形式下面討論。則我們可以算出與 對應(yīng)的特征值 以及左右特征向量組成的矩陣：,,,,,,,,,77,其次，我們記 對應(yīng)的特征通量為 ，其分量形式為,,,,顯然 ，

45、 與特征值 對應(yīng)。當(dāng) 大于零時，代表波動從左向右傳播，對于一階精度的格式，可以近似認(rèn)為： ；反之當(dāng) 小于零時，有 。寫成統(tǒng)一形式，有：,,,,,,,或者,,78,最后，我們由特征空間的通量反推出物理空間的通量。,或者,,關(guān)于

46、 的確定，定義,,Roe要求 滿足下面的條件：,,（A） 雙曲性 即 的特征值均為實數(shù)，且存在完備的左右特征向量；,,（B） 相容性,,（C） 守恒性,,,,79,利用上述三個條件，可以推導(dǎo)出

47、的確定方法為：,,,上述表達(dá)式是對左右狀態(tài)的一種特殊平均，一般稱為Roe平均。,80,4. 二階精度的迎風(fēng)格式,上面的討論主要針對一階迎風(fēng)格式。事實上，依照上面的基本思路，我們也可以構(gòu)造二階和更高階精度的格式。舉例說明如下：,1）基于FVS方法的二階格式,,,,81,時間積分: Runge-Kutta方法,2）二階精度的Roe格式,考慮一維Euler的半離散型守恒格式,對應(yīng)的特征變量

48、 的分量形式為：,,,為了達(dá)到空間二階精度，上式中的 采用下列近似：,,,82,即,,因此，數(shù)值通量為,,時間方向可以用Runge-Kutta等方法求解。,83,六、曲線坐標(biāo)中多維問題的差分格式,84,85,矢通量分裂,86,Roe,Development of Upwind Schemes for the Euler Equations Sukumar R., Chakravarthy, NASA CR-