CFDfudan可壓縮方程離散李新亮老師課件

2013復旦大學暑期課程：計算流體力學李新亮lixl@Tel:82543801力學所主樓219

1CopyrightbyLiXinliang講義、課件上傳至

（流體中文網(wǎng)）->“流體論壇”->“CFD基礎理論”下載地址：/?cid=1cc0dcbff560c149&id=1CC0DCBFF560C149%21128&sc=documents第四講可壓縮流體力學方程組的離散方法（上）

1.通量分裂技術:流通矢量分裂FVS2.Roe格式，通量差分分裂3.隱式時間推進：LU-SGS（下）

有限體積法初步2CopyrightbyLiXinliang§1通量分裂技術格式F+格式F-對流項：信息（波）從上游傳至下游——上游信息更重要——迎風差分擴散項：信息從中心向周圍擴散——不區(qū)分上、下游——中心差分迎風差分優(yōu)點：有效利用信息傳播的方向，增強穩(wěn)定性微分與差分方程的影響域N-S方程：單波方程：單波方程——一個波，容易判斷波傳播方向N-S對流項（Euler）——方程組：多波問題，復雜雙曲方程組的原則——特征分解，找到獨立傳播的波常系數(shù)矩陣A的情況——完全解耦，獨立求解變系數(shù)矩陣A的情況——局部討論3CopyrightbyLiXinliangCopyrightbyLiXinliang4Step1針對模型方程構造差分格式Step2將格式推廣到Euler方程方法1：流通矢量分裂(FVS)

特點：對流通矢量f(U)本身進行分裂

Steger-Warming,vanLeer,Lax-Friedrichs利用通量分裂技術，將模型方程推廣到Euler(或N-S)方程格式1格式2格式1格式2CopyrightbyLiXinliang5方法2：通量差分分裂（FDS）特點：對流通矢量f的導數(shù)進行分裂

Godnov,Roe,HLL,HLLC

方法3：AUSM類方法vanLeer分裂法+壓力項單獨處理

利用Riemann解利用差分表達式，計算求解Riemann問題，獲得通量1.Jacobian系數(shù)矩陣及其性質重要性質特點：A可以像常數(shù)一樣，和求導運算交換6CopyrightbyLiXinliang2流通矢量分裂(FVS)利用性質=+優(yōu)點：耗散小缺點：導數(shù)間斷方式A:特點：不必進行矩陣運算，計算量小Steger-Warming分裂A:Steger-Warming分裂7CopyrightbyLiXinliangSteger-Warming具體步驟（以一維為例）已知1）計算2）計算3）計算4）帶入（1）式得到5）利用不同的迎風格式，分別計算

（1）（后差，前差）6）計算7）時間推進8CopyrightbyLiXinliang二維問題的steger-Warming分裂令：則：具體使用步驟，以計算為例令

計算特征值

分裂特征值，計算

帶入左式，計算正、負流通矢量

計算計算設置，并注意對于曲線坐標系僅需令三維問題同樣處理

二維、三維具體公式見傅德薰等《計算空氣動力學》4.7節(jié)（158-162）書中公式有一定的排版錯誤，使用前務必重新仔細推導！9CopyrightbyLiXinliangB:Lax-Friedrichs(L-F)分裂特點：正特征值負特征值=+缺點：耗散偏大局部L-F分裂，每個點上計算全局L-F分裂，全局（一維）上計算足夠大數(shù)學性質（光滑性）最好，但耗散偏大常數(shù)與迎風格式結合，等價于人工粘性例如，可取10CopyrightbyLiXinliang方式很多=+S-W:L-F:=+VanLeer:=+11CopyrightbyLiXinliang

分裂后失去了A的性質（可以像常數(shù)一樣與求導交換）FVS分裂：

優(yōu)點：無需矩陣運算，計算量小

缺點：分裂后改變了特征方向，耗散大利用了性質一般情況下：變系數(shù)，

不能與導數(shù)交換實質：沒有做到解耦；只是把原變量重新組合，組合后波的傳播方向的保證f+向正向傳播，f-向負向傳播

缺點：由于未解耦，各變量的誤差會相互傳遞12CopyrightbyLiXinliang概念澄清：

流通矢量分裂本身不帶來耗散，但其會影響到差分的耗散;舉例：分裂過程耗散如果差分格式無耗散（例如都用中心差分），則通量分裂不帶來耗散。=+向上平移向下平移分裂差分格式耗散分裂后的流場越偏離原先流場，則總體耗散越大精確滿足，不引入誤差！如使用低精差分度格式，則對分裂形式敏感（推薦使用特征分裂）如使用高精度格式（低耗散），則對分裂形式不敏感（可使用逐點分裂）13CopyrightbyLiXinliang3.特征重構方法常系數(shù)方程組：完全解耦變系數(shù)情況——局部凍結系數(shù)…j-2j-1jj+1…在基架點上系數(shù)不變計算：在差分基架點上Aj

不變，可按常矩陣處理局部凍結系數(shù)分別采用后差和前差優(yōu)點：嚴格保證（局部）特征方向，數(shù)值解質量好；缺點：大量矩陣運算，計算量大。14CopyrightbyLiXinliang通常寫成守恒型差分，計算…j-2j-1jj+1…在基架點上系數(shù)不變具體步驟：

假設已知U，且針對模型方程（線性單波方程）

已構造出差分格式（1）1）計算出教材130頁的公式(6.1.11-6.1.13),式中用到各變量在j+1/2的值（例如）

可使用j,j+1點值的算術平均（如）或Roe平均（教材6.4節(jié)）；由計算；方法很多，例如前面介紹的或15CopyrightbyLiXinliang

均可2）

在網(wǎng)格基上計算…j-2j-1jj+1…計算fj+1/2用到的點注意，在該網(wǎng)格基上（例如k=j-1,j,j+1）保持不變例如：3)利用已構造好的差分格式，計算通量4）得到總通量

5）計算差分（j點處）步驟的算法描述（注意：實際上是兩重循環(huán)）doj=1,Ndok=j-1,j+1(網(wǎng)格基，可以是更多或更少點)

enddoenddodoj=1,N

enddo需要多次矩陣運算，計算量大

守恒性好，耗散小，數(shù)值解質量好16CopyrightbyLiXinliangCopyrightbyLiXinliang171.單方程的Roe格式線性化，用平均變化率代替(j,j+1)之間的變化率a(u)“平均斜率”，不等于“斜率的平均值”，也不等于中點處的斜率§2Roe格式——FVS分裂技術非線性情況根據(jù)Langrage中值定理，[uL,uR]之間必有一點uRoe,該點處的斜率為平均斜率;二次函數(shù)f(u)=u2中點處的斜率=平均斜率CopyrightbyLiXinliang182.方程組的情況平均斜率線性化，以平均增長率代替瞬時增長率[j,j+1]區(qū)間內連續(xù)，且可通過相似變換對角化

應當具有的性質常系數(shù)方程的Riemann解CopyrightbyLiXinliang19平均斜率xj+1/2常系數(shù)單波方程的Riemann解Roe格式：微分型近似Riemann解CopyrightbyLiXinliang203.矩陣的構造關鍵：“向量除以向量”？直接求平均增長率：uf(u)uLuRuRoeRoe點的斜率為平均斜率（根據(jù)拉格朗日中值定理，[UL,UR]區(qū)間內肯定存在Roe點）思路1：在UL與UR之間尋找一個點URoe,該點處的增長率為平均增長率f(u)=u2u二次函數(shù)——Roe點與中點重合標量函數(shù)的啟示：Roe點肯定存在（Langrage中值定理）

二次函數(shù)的中點即為Roe點思路2：進行坐標變換，得到一個二次（齊）函數(shù)引入如果

是二次（齊）函數(shù)，則其中點即為Roe點重要啟示更準確地講，應當是要求為W的線性函數(shù)，

即增長率為線性函數(shù)（中點處的增長率剛好為平均增長率）CopyrightbyLiXinliang21針對Euler方程的具體構造引入新變量：則：目的：使得F(w)是W二次齊函數(shù)（增長率為線性函數(shù)）f(U)不是U的二次齊函數(shù)二次齊函數(shù)！中點處的斜率即為平均斜率。Roe點Roe點為：增長率為線性函數(shù)！CopyrightbyLiXinliang22最終：其中如下計算：平均增長率（矩陣）含義：左、右兩個狀態(tài)點的某種平均（稱為Roe平均，為密度加權平均）

該狀態(tài)點對應的增長率（矩陣）為平均增長率（矩陣）

實際上是一種“等效平均”。效果優(yōu)于簡單的算數(shù)（或幾何）平均。

三維情況下，還有其他量（如壓力、溫度、音速等）用這三個量計算（5）簡單易記：CopyrightbyLiXinliang23Roe格式的計算步驟（半離散）已知n時刻所有網(wǎng)格點上的物理量，對于j點：1）利用差分格式計算UR,UL

2）采用Roe平均公式（5）計算Roe平均值3）將Jacobian矩陣進行特征分解:

計算4)計算5）計算6）計算空間導數(shù)7）時間推進，計算下一時間步的值。j-1jj+1與前文（第3，4講）的形式相同，僅需把式中的密度、壓力、速度等換成經(jīng)過Roe平均的密度、壓力、速度即可其中：CopyrightbyLiXinliang24可能出現(xiàn)導數(shù)不連續(xù)，可能引起數(shù)值振蕩實際使用時可用如下函數(shù)代替——所謂“熵修正”實際上是在特征值0點周圍增加了耗散Roe格式的優(yōu)點：1）保持守恒性的同時，嚴格保證了特征方向2）便于推廣到高精度格式——特征投影分裂中使用Roe平均即可（見本PPT第5頁）。推廣到高階后，雖不再保證嚴格的特征方向，但仍優(yōu)于采用算數(shù)平均方法。Roe格式的不足：

本身精度只有一階；

推廣到高階后，特征方向無法嚴格保證；

推廣到二維或三維后，特征方向無法嚴格保證，出現(xiàn)振蕩。CopyrightbyLiXinliang25關于f(U)與f(W)深入討論新變量雖然是“一次齊函數(shù)”但有變量在分母上干凈的二次齊函數(shù)自變量W的線性函數(shù)！實質區(qū)別：是自變量的線性函數(shù)，而是自變量的非線性函數(shù)！自變量U的非線性函數(shù)使用U做自變量的優(yōu)點：物理意義鮮明（質量密度、動量密度和能量密度），守恒性好使用W做自變量的優(yōu)點：Jocabian矩陣為線性矩陣思考：如果在CFD計算中，使用W替換U做自變量會怎樣？4階3階（TVD型）2階CopyrightbyLiXinliang26

§

3時間推進方法

1.顯格式推薦方法：Runge-Kutta法更高階……CopyrightbyLiXinliang272.隱格式算法簡介（以二維Euler方程為例）

1）原理介紹（1）時間離散后方案A：直接將（1）進行空間離散，得到Un+1的代數(shù)方程組困難：大型非線性方程組，求解困難方案B:設計一種迭代算法令兩端同時添加顯式項右端項，已知是已知項，可采用某種差分方法顯式計算得到（對算法無限制）CopyrightbyLiXinliang28(2)已知，（2）為線性方程。（1）(2)是一個線性化過程

含義：先用顯格式計算，再用隱格式計算修正量線性方程，離散求解離散后為大型帶狀方程組，求解計算量大LU-SGS原理：將矩陣分解為上、下三角陣，避免矩陣求逆運算隱式問題顯式化未知量：顯式部分隱式修正若隱式修正為0，則為顯格式CopyrightbyLiXinliang29Step1:求解Step2:求解LU-SGSi,ji-1,ji+1,ji,j+1i,ji-1,ji+1,ji,j+1i,j+1i,j+12）LU-SGS方法原理簡介：忽略右、上未知量忽略左、下未知量i,ji-1,ji+1,ji,j+1i,j+1原則上：（1）（2）兩步反復迭代，直至收斂CopyrightbyLiXinliang30LU-SGS方法求解Euler方程a)將矩陣A,B分裂為了簡化計算，通常采用L-F分裂矩陣分裂1階迎風格式離散整理迭代收斂后，