計算流體課程-02

1、Copyright by Li Xinliang,1,計算流體力學(xué)講義 第九講 有限體積法（1）,知識點：,有限體積法的基本概念 重構(gòu)和反演 迎風(fēng)型有限體積法Riemann求解器；Roe格式的新理解：近似Riemann解 多維迎風(fēng)型有限體積法坐標(biāo)旋轉(zhuǎn),Copyright by Li Xinliang,2,知識回顧,1. 差分方法的基本概念：,差分格式、修正方程、相容性、收斂性、穩(wěn)定性、LAX等價定理,2. 精度分析、穩(wěn)定性分析與分辨率分析（修正波數(shù)）,Taylor分析,Fourier分析,修正波數(shù),激波捕捉格式 GVC, NND, Roe, Godnov, MUSCL, TVD, WENO,E

2、uler (N-S) 方程的通量分裂 逐點分裂、特征投影分裂 （建議使用Roe平均）,Copyright by Li Xinliang,3,5. 隱格式求解的LU-SGS方法,要點： a. 引入差量，方程線性化 b. 單邊差分，隱式代數(shù)方程顯式（推進）化,以一維為例，多維可直接推廣,方法1：直接隱式離散,直接求解,非線性方程組，計算量大,方法2,差量化,線性化,已知項,線化微分方程,Copyright by Li Xinliang,4,求解思路：如果直接離散，得到線性代數(shù)方程組，仍需求解，計算量大（多維情況）,如果能單側(cè)差分就好解了！,多對角方程組，不好解（多維情況）,中心（雙側(cè)）離散,如果單

3、側(cè)離散,單側(cè)離散，可推進求解，免受解方程組之苦。真簡單,Copyright by Li Xinliang,5,可是，A有正有負(fù)，無法單側(cè)差分化,還是個三對角的,奇思妙想：如果分成兩個子步，各自用單側(cè)值，就簡單多了,強行單側(cè)差分會不穩(wěn)定的,近似LU分解,Step 1:,近似LU分解,Step 2:,均為遞推求解 （兩次掃描），免受解方程組之苦,j -1 - j,j+1 j,以上描述適用于求解定常問題，求解非定常問題該過程可用于內(nèi)迭代。,迭代收斂后q趨于0，精度由右端項決定,Copyright by Li Xinliang,6, 9.1 有限體積法入門,有限體積法主要優(yōu)勢： 處理復(fù)雜網(wǎng)格,差分法處

4、理復(fù)雜外形 坐標(biāo)變換,坐標(biāo)變換函數(shù)必須足夠光滑 否則損失精度,實際問題： 外形復(fù)雜， 光滑的結(jié)構(gòu)網(wǎng)格生成困難,Copyright by Li Xinliang,7,9.1.1 有限體積法 的基本概念,實質(zhì)： 把幾何信息包含于離散過程中,雖然簡單，但有助于建立基本概念,j-1 j j+1,j-1/2 j+1/2,1. 全離散型過程,含義： f在j+1/2點的值 （注意與差分法的區(qū)別）,在控制體上積分原方程,定義：,空間平均,時間平均,精確推導(dǎo)，不含誤差,提示： 為區(qū)間內(nèi)的空間及時間平均值，如果把它們理解為某點的值，會產(chǎn)生誤差,Copyright by Li Xinliang,8,積分（精確）,重

5、構(gòu)（Reconstruction),有限差分法的離散：數(shù)值微分過程 有限體積法的離散：數(shù)值積分過程,積分方程,離散化,反演（evolution),(1) 重構(gòu)過程,A. 零階重構(gòu)，假設(shè)分片常數(shù),j-1,B. 線性重構(gòu)，假設(shè)分片線性函數(shù),零階重構(gòu)與一階重構(gòu)示意圖,j,j+1,or,or,或其他方法,C. 更高階的重構(gòu)例如: 分片二次函數(shù) （PPM）, WENO等,重構(gòu)是有限體積的空間離散化過程，有多種方法,Copyright by Li Xinliang,9,（2） 演化過程 （以線性方程為例）,需要得知時間演化信息，通常利用特征方程,若采用零階重構(gòu)：,則：,假設(shè)時間步長足夠小,則方程為：,等價

6、于一階迎風(fēng)差分,Riemann解,Copyright by Li Xinliang,10,若采用線性重構(gòu),若,Warming-Beam,Lax-Wendroff,0階重構(gòu) 1階精度 線性重構(gòu) 2階精度,一維均勻網(wǎng)格的有限體積法等價于有限差分法,Euler方程： 演化過程可通過Riemann解或近似Riemann解進行,Copyright by Li Xinliang,11,2. 半離散方法,全離散： 積分方程 代數(shù)方程 （守恒性好，但復(fù)雜） 半離散： 積分方程 常微分方程 （簡便，便于使用R-K等成熟方法）,僅空間積分,f 在j+1/2點的值，仍需要使用周圍點 進行插值,通常無法精確計算， 可

7、采用近似值 代替,等價于二階中心差分,半離散,j-1 j j+1,j-1/2 j+1/2,重構(gòu),Copyright by Li Xinliang,12,9.1.2 一維Euler方程的迎風(fēng)型有限體積法,j-1 j j+1,j-1/2 j+1/2,半離散,1. 重構(gòu),控制體積,j-1,j,j+1,左重構(gòu)值,右重構(gòu)值,選擇不同的模板會得到不同的重構(gòu)方案,向左偏的模板產(chǎn)生 向右偏的模板產(chǎn)生,差分法 同一點的導(dǎo)數(shù)可使用向前差分和向后差分，根據(jù)特征方向選擇之,例如： 0階重構(gòu) 1階單邊重構(gòu),根據(jù)特征方向，選擇左通量或右通量,途徑1： FVS,途徑2：FDS,Copyright by Li Xinlian

8、g,13,2. 分裂方法 （1）: FVS方法 （流通矢量分裂 逐點分裂）,具體方法： Steger-Warming 分裂 Lax-Friedrichs分裂 Van Leer分裂： Liou-Steffen分裂： （壓力項與其他項分開， AUSM類格式的基礎(chǔ)）,根據(jù)當(dāng)?shù)豈ach數(shù)分裂,保證 的Jocabian陣特征值為正， 的為負(fù),正通量： 向左偏斜重構(gòu)； 負(fù)通量： 向右偏斜重構(gòu) 偏重向上游,與迎風(fēng)差分法類似： 網(wǎng)格基（或權(quán)重）偏重上游,差分、有限體積都可使用,Copyright by Li Xinliang,14,一個參數(shù)，反映全部特征,小知識： Liou-Steffen分裂,對流項,壓力項

9、,思路： 決定特征的關(guān)鍵參數(shù) 當(dāng)?shù)豈ach數(shù),超音速，x-方向,超音速，x+方向,因此，對Mach數(shù)進行分裂更為簡潔！,顯然：,參考文獻： Toro: Riemann Solvers and Numerical Methods for Fluid Dynamics, section 8.4.4 Liou: Ten Years in the making AUSM family, NASA TM-2001-210977,類似 Van Leer分裂，但壓力單獨處理,M,保證光滑過渡,M=1,Copyright by Li Xinliang,15,(3) FDS 方法 （通量差分分裂特征投影分裂）,

10、1. 利用精確Riemann解Godnov格式,目的：,j-1 j j+1,j-1/2 j+1/2,控制體積,j-1,j,j+1,左重構(gòu)值,右重構(gòu)值,1） 精確求解Riemann問題,2）,精度： 取決于重構(gòu)的精度 （原則上可任意階）,差分法：Godnov格式使用分片常數(shù)，精度1階 有限體積法：先重構(gòu)，再解Riemann問題，可高階,精確Riemann解（見本講座第2講）需迭代求解，計算量大 - 近似Riemann解,整體思路： 先重構(gòu)自變量（兩種方案得到 ）， 再求解Riemann問題（或用FVS）得到通量的方法通常稱為MUSCL方法。,Copyright by Li Xinliang,16

11、,差分法與有限體積法區(qū)別與聯(lián)系（二階迎風(fēng)+FVS為例）,差分、有限體積,差分（通常做法）： 直接插值通量fi+1/2,有限體積：先插值自變量U，然后計算通量f：,先插值自變量，再計算通量的方法，稱為MUSCL類方法。 是有限體積法的常用方法（差分法也可以用）,單側(cè)重構(gòu)，以避免跨過激波,還可使用FDS方法，重構(gòu)后求解Riemann問題,當(dāng)f=f(U) 連續(xù)時，對f插值與對U插值精度相同。,Copyright by Li Xinliang,17,（稱為數(shù)值流通量） 的含義,重要概念澄清： 重構(gòu)與插值,A. 有限差分法：,j+1/2,切線,j-1/2,j,j-1,注意： 與 f 在xj+1/2點的值

12、含義不同！,用周圍幾個點的值 計算 的過程稱為“重構(gòu)”，不能理解為用 來插值,記號 確實容易混淆，讓人容易聯(lián)想起 。記為 更好些,否則，最高只能達到2階精度了！,Copyright by Li Xinliang,18,是控制體內(nèi)的平均值,（稱為數(shù)值流通量） 的含義,重要概念澄清： 重構(gòu)與插值,B. 有限體積法：,j+1/2,j-1/2,確實為f在xj+1/2點的值 ！,通常做法： 1） 用 計算出 2）,u在xj+1/2點的值！,關(guān)鍵： 是用 計算 （稱為重構(gòu)） ，而不是用 計算 （是標(biāo)準(zhǔn)的插值）；否則最高也只能達到2階精度。,19,概念：MUSCL與非MUSC類方法,j+1/2,切線,j-1

13、/2,j-1,差分,有限體積,方法1 （非MUSCL類）： 直接利用周圍幾個點的函數(shù)值 或 ）直接計算 （或 ）,如何計算 或 ？,方法2 （MUSCL類）： 利用周圍幾個點的自變量值 （或 ）計算出 （或 ）； 然后再計算 （或 ）,當(dāng)f=f(u)是連函數(shù)時，二者精度相同,f的誤差與u的誤差同階,Copyright by Li Xinliang,20,2. 近似Riemann解 例： Roe格式,與差分法的Roe格式形式相同,理解： 近似Riemann解（Euler方程 常系數(shù)線性化解）,u,f(u),uL,uR,uRoe,利用Roe平均， 剛好是左右兩點間的平均增長率，實現(xiàn)了常系數(shù)線性化。

14、,常系數(shù)雙曲方程組，易解！,思路： 用平均增長率矩陣 取代瞬時增長率矩陣A,不但實現(xiàn)了線性化，而且實現(xiàn)了常系數(shù)化。 利用二次齊函數(shù)的性質(zhì)，可找到了Roe點（即Roe平均點），該點處的增長率剛好等于平均增長率。,Roe平均,常系數(shù)化,線性化,21,常系數(shù)線性單波方程的Riemann問題太簡單了,常系數(shù)方程組的Riemann問題,解耦了的單波方程，有精確解,初值,Copyright by Li Xinliang,22,解為,線性化條件 （并利用齊函數(shù)性質(zhì)）,與差分法的Roe格式相同,還有各種其他類型的近似Riemann解（今后介紹）,Copyright by Li Xinliang,23,9.1.

15、3 多維問題的有限體積法,二維問題,一維Riemann問題，坐標(biāo)選取不當(dāng)，變?yōu)?“二維”Riemann問題,x,y,差分法： 獨立計算 只考慮各自的特征方向,由于非線性， 實際（二維）特征方向并非x,y方向特征量的線性組合。 特征方向計算不嚴(yán)格，帶來誤差,差分方法： 多維情況，特征理論復(fù)雜，通常x,y方向獨立計算,轉(zhuǎn)化為 x方向與y方向的兩個一維問題,逐點分裂,特征投影分裂,完全按照一維情況獨立處理,局部坐標(biāo)旋轉(zhuǎn)？ 差分算法設(shè)計造成局部旋轉(zhuǎn)困難,差分法 的多維 處理方法,1. 小知識： 差分方法如何處理高維問題的 ？優(yōu)缺點？,優(yōu)點： 簡單 缺點：特征方向計算不準(zhǔn),Copyright by Li

16、 Xinliang,24,2. 二維有限體積方法的離散過程,在以某節(jié)點為中心的控制體上積分,i,j,k,非結(jié)構(gòu)網(wǎng)格的控制體,i+1,j,i-1,j,i,j+1,i,j-1,k3,k1,k2,k4,k5,結(jié)構(gòu)網(wǎng)格的控制體,x,y,n,體積平均,控制體邊界垂直于節(jié)點連線（也可選其他方式）,垂直平分線,n,1) 建立控制體,2) 在控制體上積分,離散方程,需要重構(gòu)： 由節(jié)點上平均值 給出函數(shù)分布，最終給出通量,表示第m個界面上的值,Copyright by Li Xinliang,25,1） 重構(gòu) 兩種不同的重構(gòu)方案，向左偏及向右偏。 給出兩種結(jié)果： 及,i,j,i+1,j,i-1,j,i,j+1,

17、i,j-1,n,左重構(gòu),右重構(gòu),2） 由左右重構(gòu)得到的自變量： 和 給出通量 方案A: FVS 方案B: 解Riemann問題 （常用）,3. 二維迎風(fēng)型有限體積法,例如： 0階重構(gòu)：,線性重構(gòu):,用i, i-1點的值 插i+1/2點的值 （網(wǎng)格劇烈變化時，應(yīng)當(dāng)用實際坐標(biāo)插值）,用i, i+1點的值 插i+1/2點的值,x,y,看似二維Riemann問題，其實是一維的，坐標(biāo)旋轉(zhuǎn)一下就行了,Copyright by Li Xinliang,26,x,y,x,y,（通常）進行坐標(biāo)旋轉(zhuǎn),旋轉(zhuǎn)q角后的坐標(biāo)系 (x,y),性質(zhì)： Euler方程的旋轉(zhuǎn)不變性,形式上完全不變 （僅需把u,v,x,y換成u,

18、v,x,y即可）,其中：旋轉(zhuǎn)矩陣,旋轉(zhuǎn) q 角,矩陣表示,Copyright by Li Xinliang,27,i,j,i+1,j,i-1,j,i,j+1,i,j-1,左重構(gòu),右重構(gòu),局部坐標(biāo)系,x,y,x,y,旋轉(zhuǎn)q角后的坐標(biāo)系 (x,y),習(xí)題： 設(shè) n 為平行x軸的向量，試證明：,證明：,坐標(biāo)旋轉(zhuǎn)，標(biāo)量不變 向量的模不變,Copyright by Li Xinliang,28,i,j,i+1,j,i-1,j,i,j+1,i,j-1,左重構(gòu),右重構(gòu),于是：,x,y,x,y,旋轉(zhuǎn)q角后的坐標(biāo)系 (x,y),其中下標(biāo)m表示控制體第m個面（線）， 表示該面的面積 （長度）,于是，問題轉(zhuǎn)化為求控

19、制面上的,這個量有兩個重構(gòu)方案,方法1： FVS: 方法2： 需要求解Riemann問題，旋轉(zhuǎn)后，轉(zhuǎn)化為“擴展的”1維Riemann問題,Copyright by Li Xinliang,29,解釋： “擴展的”一維Riemann問題,x,y,x,y,旋轉(zhuǎn)q角后的坐標(biāo)系 (x,y),問題本身是一維的 ： 所有變量都只沿著x方向分布，沿y方向均勻 允許有y方向的速度v (比純一維問題多一個變量）,v的存在對流動的一維性質(zhì)無任何影響,舉例： Sod激波管問題（一維）。 如果在沿y方向勻速運動的坐標(biāo)系中觀察，則方程為“擴展的一維問題”，但不影響其一維性質(zhì),坐標(biāo)系沿y方向勻速運動,x,y,可用精確Ri

20、emann解，也可用Roe等近似解,Copyright by Li Xinliang,30,二維迎風(fēng)型有限體積法求解步驟 1） 對n時刻的平均量 進行重構(gòu)，給出控制面上的左、右重構(gòu)值 ， 2）將以上值旋轉(zhuǎn)到（每個）控制面法向的局部坐標(biāo)系下： 3） 求解上述“擴展的”一維Riemann問題，給出后續(xù)時刻控制面上的值 4）利用積分型方程： 計算下一時刻的平均量,i,j,i+1,j,i-1,j,i,j+1,i,j-1,左重構(gòu),右重構(gòu),0階重構(gòu)：,1階重構(gòu) （線性重構(gòu)）：,更復(fù)雜的重構(gòu)（WENO等）,下標(biāo)m指的是第m個控制面上的值,Copyright by Li Xinliang,31,知識回顧： R

21、iemann問題精確解,Riemann問題,問題描述： 初始時刻，物理量分布存在單個間斷；間斷兩側(cè)物理量為常數(shù)。,求解思路： 采用積分方程 單個間斷，且間斷兩側(cè)物理量為常數(shù)情況下： 積分方程轉(zhuǎn)化為代數(shù)方程,代數(shù)方程： 質(zhì)量、動量、能量守恒,Copyright by Li Xinliang,32,計算出 ， 將 與這三個值進行比較，判斷會產(chǎn)生的情況。具體見下圖：,Riemann問題的具體計算步驟 （全流場）,1. 判斷可能會出現(xiàn)的情況（五種情形之一）,a. 定義函數(shù),b. 進行判斷,情況5,情況4,情況3,情況1,情況5,情況4,情況2,情況1,單調(diào)增函數(shù)，性質(zhì)很好,Copyright by L

22、i Xinliang,33,2. 求解中心區(qū)的壓力和速度,單未知數(shù)的代數(shù)方程，迭代求解（例如Newton法，F(xiàn)(p)性質(zhì)好，求解不困難）,3. 確定中心區(qū)接觸間斷兩側(cè)的密度 以及左、右波傳播的速度 a. 左波為激波的情況（情況1,3）,b. 左波為稀疏波的情況 （情況2，4，5）,中心區(qū) 接觸間斷左側(cè)的物理量,膨脹波的波頭及波尾速度,激波的傳播速度,對于情況（5），波尾速度為：,中心區(qū)為真空，音速 無定義，改由該式計算,Copyright by Li Xinliang,34,c. 右波為激波的情況（情況1，2）,中心區(qū) 接觸間斷右側(cè)的物理量,b. 右波為稀疏波的情況 （情況2，4，5）,4.

23、計算稀疏波區(qū)域的值（如果有稀疏波的話）,a. 左稀疏波 b. 右稀疏波,情況2，4,情況5：,注意： 教科書32頁c的公式有誤！,Copyright by Li Xinliang,35,有限體積法 “擴展的”Riemann問題的計算方法 （中心線x=0處）,迎風(fēng)型有限體積法，需要求解“擴展的”一維Riemann問題,x,y,物理問題分析： 所有物理量均沿x方向一維分布，沿y方向均勻分布。,僅需計算t時刻x=0處各物理量的值,v和w跟隨流體運動，相當(dāng)于“被動標(biāo)量”,穿過激波及稀疏波，切向速度不變,Copyright by Li Xinliang,36,求解t時刻x=0處物理量的具體步驟,Step

24、 1: 求解 得到中心區(qū)壓力 Step 2: 計算中心區(qū)的速度 Step 3: 根據(jù) 及 判斷會出現(xiàn)哪種情況（五種情況之一） Step 4: 根據(jù)具體情況（左、右波是激波還是膨脹波）計算出中心區(qū)接觸間斷兩側(cè)的密度 及 Step 5: 如果中心區(qū)x方向速度0, 則中心線（x=0）處的密度及切向速度為接觸間斷右側(cè)的值，否則為接觸間斷左側(cè)的值。,具體公式見本PPT33-34頁， 本頁中上標(biāo)“L”和“R”分別對應(yīng)原先的下標(biāo)“1”和“2”,x,t,v和w沒有給計算過程帶來任何麻煩,先無視v和w的存在，求解標(biāo)準(zhǔn)1維Riemann問題。再根據(jù)u的符號確定中心線的v和w,Copyright by Li Xin

25、liang,37,作業(yè)9.1 求解 “坐標(biāo)旋轉(zhuǎn)的Sod激波管問題”,x,y,物理問題描述：如右圖示，有一條與x軸夾角120的直線，左右兩側(cè)充滿同種介質(zhì)的無粘完全氣體。初始時刻，左右兩側(cè)的氣體狀態(tài)為： 左側(cè)： ，右側(cè) 試計算t=0.14時刻的流動分布。,計算要求： 1） 計算域 網(wǎng)格 2）初值設(shè)置如圖所示； 3）空間離散采用有限體積法計算。采用線性重構(gòu)（見上頁）。時間推進采用3階Runge-Kutta方法 4) 分別采用FVS方法及Roe-FDS （又稱Roe- Riemann近似解，見本PPT 22頁）兩種不同的方法計算通量。 5） 繪制出t=0.14時刻密度、速度u,v及壓力的二維分布。 6

26、） 繪制出t=0.14時刻x軸（垂直于初始間斷面，見右圖）上的密度、沿x方向的速度及壓力的一維分布，并同精確解進行比較。,x,y,計算流體力學(xué)講義 第十講 有限體積法（2）,知識點：,38,Copyright by Li Xinliang,近似Riemann求解器 HLL / HLLC 中心型有限體積法人工粘性 具體問題的計算 翼型繞流 (邊界條件、具體解法、粘性項處理),Copyright by Li Xinliang,39,知識回顧,在以某節(jié)點為中心的控制體上積分,i,j,k,非結(jié)構(gòu)網(wǎng)格的控制體,i+1,j,i-1,j,i,j+1,i,j-1,k3,k1,k2,k4,k5,結(jié)構(gòu)網(wǎng)格的控制體

27、,x,y,n,體積平均,控制體邊界垂直于節(jié)點連線（也可選其他方式）,垂直平分線,n,1) 建立控制體,2) 在控制體上積分,離散方程,重構(gòu)： 由節(jié)點上平均值 給出函數(shù)分布，最終給出通量,表示第m個界面上的值,1. 有限體積法的離散過程重構(gòu),1） 重構(gòu) 兩種不同的重構(gòu)方案，向左偏及向右偏。 給出兩種結(jié)果： 及,Copyright by Li Xinliang,40,i,j,i+1,j,i-1,j,i,j+1,i,j-1,n,左重構(gòu),右重構(gòu),2） 由左右重構(gòu)得到的自變量： 和 給出通量 方案A: FVS 方案B: 解Riemann問題 （常用）,例如： 0階重構(gòu)：,線性重構(gòu):,用i, i-1點的值

28、 插i+1/2點的值 （網(wǎng)格劇烈變化時，應(yīng)當(dāng)用實際坐標(biāo)插值）,用i, i+1點的值 插i+1/2點的值,x,y,看似二維Riemann問題，其實是一維的，坐標(biāo)旋轉(zhuǎn)一下就行了,2. 迎風(fēng)型有限體積法,（稱為數(shù)值流通量） 的含義,Copyright by Li Xinliang,41,重要概念澄清： 重構(gòu)與插值,A. 有限差分法：,j+1/2,切線,j-1/2,j,j-1,注意： 與 f 在xj+1/2點的值含義不同！,用周圍幾個點的值 計算 的過程稱為“重構(gòu)”，不能理解為用 來插值,記號 確實容易混淆，讓人容易聯(lián)想起 。記為 更好些,否則，最高只能達到2階精度了！,是控制體內(nèi)的平均值,（稱為數(shù)值

29、流通量） 的含義,Copyright by Li Xinliang,42,重要概念澄清： 重構(gòu)與插值,B. 有限體積法：,j+1/2,j-1/2,確實為f在xj+1/2點的值 ！,通常做法： 1） 用 計算出 2）,u在xj+1/2點的值！,關(guān)鍵： 是用 計算 （稱為重構(gòu)） ，而不是用 計算 （是標(biāo)準(zhǔn)的插值）；否則最高也只能達到2階精度。,43,概念：MUSCL與非MUSC類方法,j+1/2,切線,j-1/2,j-1,差分,有限體積,方法1 （非MUSCL類）： 直接利用周圍幾個點的函數(shù)值 或 ）直接計算 （或 ）,如何計算 或 ？,方法2 （MUSCL類）： 利用周圍幾個點的自變量值 （或

30、）計算出 （或 ）； 然后再計算 （或 ）,當(dāng)f=f(u)是連函數(shù)時，二者精度相同,f的誤差與u的誤差同階,Copyright by Li Xinliang,44,10.1 近似Riemann解,迎風(fēng)型有限體積法通常需要求解Riemann問題； 精確Riemann解計算量大,10.1.1 HLL Riemann近似求解器 （Harten, Lax 其他壁面u=v=0; 流函數(shù)的邊界條件：,邊界是一條流線， 流線是流函數(shù)的等值線,渦量的邊界條件： 由速度給出,可用更高階的格式,Copyright by Li Xinliang,122,12.4 SIMPLE方法,基本思想： 與（離散型）投影法類似

31、， 但速度推進是隱式的；,1） 已知預(yù)估壓力 計算速度,采用隱式離散,2） 壓力及速度修正,重要簡化,類似“人工壓縮性方法”,修正方程對角化 （顯式化）,已知,聯(lián)立求解,Copyright by Li Xinliang,123,帶入離散的連續(xù)性方程：,得到離散的壓力Poisson方程：,求解后，得到壓力修正值：,（4）,帶入（4）時得到n+1時刻的速度,具體步驟： 1） 已知n時刻的速度壓力 2） 預(yù)估壓力 （可取為n時刻的壓力） 3） 帶入（1）（2）式，解出 （隱格式，需迭代求解） 4） 求解壓力的修正方程 （5）得到修正壓力 5） 帶入（4）式，得到n+1時刻的速度及壓力 6） 推進求解

32、直到給定時刻（或收斂）,如該步改用顯格式，則為（離散型）投影法,提示： 對于定常問題，內(nèi)迭代無需收斂，最終時刻收斂即可,習(xí)題 12.1 求解方腔問題,問題描述： 如圖示邊長為L的方腔，上表面流體以常速度U運動，求解里面的流場（假設(shè)流動定常）。 考慮 三種情況,要求： 數(shù)值方法不限 （人工壓縮性方法、投影法、渦量-流函數(shù)方法及SIMPLE方法均可）； 空間離散采用差分法，建議采用較高階精度的方法。 繪制出定常解的流線圖。 請詳細(xì)寫明方程及公式的推導(dǎo)過程及計算流程，切勿只上交計算結(jié)果。,計算流體力學(xué)講義 第十三講 湍流與轉(zhuǎn)捩 （1）,知識點：,125,Copyright by Li Xinlian

33、g,線性穩(wěn)定性理論 轉(zhuǎn)捩的預(yù)測方法 壁湍流轉(zhuǎn)捩的渦動力學(xué)機制,Copyright by Li Xinliang,126, 13.1 線性穩(wěn)定性理論,一、 穩(wěn)定性基本概念,常識：流體中的不穩(wěn)定性,K-H不穩(wěn)定性,A. K-H (Kelvin-Helmholtz）不穩(wěn)定性 自由剪切流的無粘不穩(wěn)定性,混合層 K-H不穩(wěn)定性,K-H不穩(wěn)定性的關(guān)鍵： 速度剖面有拐點,Lee-Lin: 速度剖面的拐點是無粘不穩(wěn)定性的必要條件,流體不禁搓，一搓搓出渦,已知某運動狀態(tài)； 在此基礎(chǔ)上施加微小擾動； 如擾動隨時間（或空間）衰減，則稱系統(tǒng)穩(wěn)定，否則為不穩(wěn)定,Copyright by Li Xinliang,127,

34、自然界中 K-H不穩(wěn)定性圖片,智利塞爾扣克島的卡門渦街,澳大利亞Duval山上空的云,KelvinHelmholtz instability clouds in San Francisco,佛蘭格爾島周圍的卡門渦街,高速流,低速流,自由剪切層受到擾動界面變形后的情況 K-H不穩(wěn)定性的產(chǎn)生機理,受阻減速，壓力升高，產(chǎn)生高壓區(qū),高壓導(dǎo)致變形加劇,Copyright by Li Xinliang,128,B. T-S (Tollmien-Schlichting) 不穩(wěn)定性不可壓 壁面剪切流的粘性不穩(wěn)定性,Mack 不穩(wěn)定性 超聲速壁面剪切流的不穩(wěn)定性,不可壓邊界層速度剖面 （Blasius解） 無拐

35、點,可壓縮情況 Mach數(shù)足夠高時會出現(xiàn)廣義拐點 出現(xiàn)無粘不穩(wěn)定性,不可壓縮無粘不穩(wěn)定性 需存在拐點 可壓縮無粘不穩(wěn)定性需存在廣義拐點,Mach 6 鈍錐（1攻角） 不同子午面 的分布,超音速平板邊界層的不穩(wěn)定波,第1模態(tài)（T-S波）,第2模態(tài) （Mack模態(tài)）,Copyright by Li Xinliang,129,激波,密度界面,R-M (Richtmyer-Meshkov)不穩(wěn)定性 激波與密度界面作用的斜壓效應(yīng),慣性約束聚變（ICF）示意圖,小知識 渦的產(chǎn)生機制： 粘性、 斜壓、有旋的外力,激波,密度界面,斜壓項,Copyright by Li Xinliang,130,D. R-T

36、(Reyleigh-Taylor)不穩(wěn)定性 重力帶來的不穩(wěn)定性,R-T (Reyleigh-Taylor)不穩(wěn)定性,重介質(zhì),輕介質(zhì),Copyright by Li Xinliang,131,E Barnard熱對流不穩(wěn)定性,其他學(xué)科的不穩(wěn)定性：,Euler壓桿的不穩(wěn)定性,Barnard 熱對流的胞格結(jié)構(gòu),板殼的不穩(wěn)定性,Copyright by Li Xinliang,132,二、 穩(wěn)定性問題的常用數(shù)學(xué)方法 線性穩(wěn)定性分析,Step 1: 得到線性化的擾動方程,控制方程為：,已知其具有解,最好是精確解，也可用高精度的數(shù)值解,令：,舍棄高階小量，得到線性化的擾動方程,（1）,例如： 平板的Bla

37、sius解，槽道的Poiseuille 解,線性方程,Copyright by Li Xinliang,133,Step 2: 求解 的特征值問題,什么條件下具有非零解，非零解如何？,通常假設(shè)在某些方向具有周期性，轉(zhuǎn)化為一維問題,數(shù)值方法： 將 （1） 離散代數(shù)方程 何時有非零解， 非零解如何？ 特征值問題,什么條件下有非零解？,特征值問題計算量巨大，目前通常只能處理一維問題,Copyright by Li Xinliang,134,三、 穩(wěn)定性問題示例 不可壓縮槽道流動的線性穩(wěn)定性（LST）理論 (以二維為例),Step 1: 獲得線性化擾動方程,令：,Poiseuille解：,(2),代入

38、方程（2），并舍去高階小量得到線性化的擾動方程,（3）,1） 控制方程及邊界條件,Copyright by Li Xinliang,135,研究擾動發(fā)展的空間模式和時間模式,擾動源,空間模式： 任一點的擾動具有時間周期性 符合物理條件,假設(shè)擾動具有如下形式：,沿流向及時間方向具有波動特性 稱為Tolmien-Schlichting（T-S）波,任意擾動可分解為正弦波的疊加 線性系統(tǒng)各成分無相互作用 可獨立研究,為實數(shù),為復(fù)數(shù),擾動波的振幅沿流向指數(shù)變化,空間增長率,時間模式： 擾動具有流向的周期性 假設(shè)一窗口沿流向運動，研究窗口內(nèi)擾動的演化,為實數(shù),為復(fù)數(shù),擾動波的振幅雖時間變化,時間增長率,

39、Copyright by Li Xinliang,136,以時間模式為例：,（4）,（5）,（6）,線性偏微方程（3）轉(zhuǎn)化成為含參數(shù)的線性常微方程組（4）-（6）,譜方法的常規(guī)做法,通過消元法，轉(zhuǎn)化為更高階的常微方程 （不是必須的）,常用做法，通常還可以反向為之： 高階方程轉(zhuǎn)化為低階方程組,消去,Orr-Sommerfeld（O-S) 方程,其中：,最終，控制方程為O-S方程：,Copyright by Li Xinliang,137,邊界條件：,y=1 （固壁）:,y=0 （中心線，對稱）：,可以取計算域-1,1,使用固壁邊界條件； 也可以取計算域-1,0,使用固壁及對稱邊界條件,流函數(shù)形式

40、的O-S方程,引入流函數(shù)，使得：,計算出 后，利用公式,計算其他兩個量,則：,令：,常數(shù)倍,滿足的方程及邊界條件與 完全相同。,如果 恒大于（或恒小于0），則必有,Copyright by Li Xinliang,138,小知識： 關(guān)于O-S方程,1） O-S方程適用于不可壓平行流的穩(wěn)定性問題 （不僅槽道流） 2） 準(zhǔn)平行流 （流線沿x方向接近平行）也可使用（例如邊界層流動） 3）如果舍去粘性（左端）項，則方程稱為Rayleigh方程,Rayleigh拐點定理： Rayleigh方程存在不穩(wěn)定解的必要條件是速度型存在拐點。即存在某點 使得,若存在無粘不穩(wěn)定性，該項必有0點。,分部積分，并取虛部

41、，得：,不存在非穩(wěn)定解,Copyright by Li Xinliang,139,2） O-S方程的解法,數(shù)學(xué)表述 奇性（特征值）問題： 參數(shù) 為何值時，方程有非零解？ 非零解如何？,時間發(fā)展槽道湍流： (通常) 給定Re及a，問 w取何值時，O-S方程有非零解？,增長率,求解步驟： 1） 將O-S方程離散，得到線性代數(shù)方程組 離散方法： 差分法、有限元法、譜方法、打靶法 2） 求w，使得該方程有非零解（奇性或特征值問題）.,求出w,局部法：只求出一個w 全局法：計算出全部的w,試計算 的時間發(fā)展槽流中 （即波長2p）T-S波的頻率及增長率,Copyright by Li Xinliang,1

42、40,四: 例題,Step 1: 離散 （差分法）,一維問題,網(wǎng)格： 均勻網(wǎng)格 （簡單，但需要較多網(wǎng)格點 N 300） 非均勻網(wǎng)格,差分離散：,2階格式,4階格式：自行推導(dǎo) （利用小程序）,非均勻網(wǎng)格：,問題： 會產(chǎn)生大量非物理解 （例： 1000個網(wǎng)格點算出1000個特征值） ； 可通過不同網(wǎng)格的對比，進行篩選,Copyright by Li Xinliang,141,離散化后得：,Step 2: 求廣義解特征值問題（1） 即w為何值時（1）有非零解,(1),方法1： 全局法 一次計算出全部特征值,常用Q-Z分解法; 其他： 冪法、反冪法、Jacobi法，Householder ),狹義特征

43、值問題,廣義特征值問題,求解特征值是計算數(shù)學(xué)的主要研究方向，有大量成熟的方法,可借助軟件包或Lapack庫等 （自行到網(wǎng)上搜索）,通常關(guān)心最不穩(wěn)定的擾動波,最大的那個,Copyright by Li Xinliang,142,方法2： 局部法,令：,方法： Newton法， 弦位法，拋物線（Muler）法,Newton 法：,弦位法：,差分化,拋物線法：,已知：,可連成一條拋物線,令：,求出新的值,消元法計算行列式（利用5對角特征）,計算出 后，求解方程（1）就可得到特征向量。,（1）,方程有奇性，可補充一個條件，例如給定某個點的值,Copyright by Li Xinliang,143,效

44、果更好的方法 Malik提出的緊致差分格式求解 參考文獻: 周恒等： 流動穩(wěn)定性 （p. 10-13），國防工業(yè)出版社,非均勻網(wǎng)格的超緊致格式（4階精度）,優(yōu)點： 緊致網(wǎng)格基僅2點 精度高 4階 直接適用于非均勻網(wǎng)格無需坐標(biāo)變換,原方程組：,形式變換 變?yōu)橐浑A方程組,通常的做法,推導(dǎo)倉促 可能有誤 請仔細(xì)推導(dǎo),與y無關(guān),Copyright by Li Xinliang,144,令：,帶入差分格式：,得：,Copyright by Li Xinliang,145,令：,邊界處表達式 (C1,D1,CN,DN) 根據(jù)邊界條件而定,內(nèi)點的表達式,特點： 離散形成的代數(shù)方程組呈 塊兩對角特征,求行列式

45、，特征值都非常便利！,其余步驟與前文相同 可用矩陣廣義特征值理論計算全部模態(tài) 也可利用 計算單個模態(tài),計算域可以取為-1,1。 也可取為-1,0， 在中心線給對稱（或反對稱）邊界條件,Copyright by Li Xinliang,146,邊界條件,1） 如果取完整計算域 -1,1,2) 如果取一半計算域-1,0,處可設(shè)定對稱或反對稱邊界條件,如設(shè)定對稱條件，只能計算出對稱擾動模態(tài) 如設(shè)定反對稱條件，只能計算出反對稱模態(tài),對于槽道流，最不穩(wěn)定模態(tài)是對稱模態(tài),但有些情況下，穩(wěn)定模態(tài)在轉(zhuǎn)捩過程中也發(fā)揮作用（例如感受性過程，見Zhong et al. JFM 556,55-103,2006）,按照

46、內(nèi)點方法計算,先根據(jù)處理內(nèi)點的方法，求出所有點上的C,D值，再對邊界點進行特殊處理（D的前兩行設(shè)為0， C的前兩行設(shè)為(1,0,0,0)及（0，1，0，0）,Copyright by Li Xinliang,147,具體解法（局部法）,求解,顯然,因此只需計算每個4*4矩陣的行列式即可 （可直接寫出表達式，也可用消元法計算）,編制好計算 的子程序后，可利用Newton法，弦位法，拋物線法等求解 ，得到復(fù)增長率,Copyright by Li Xinliang,148,計算出 后，利用消元法求解方程 ，即可得到特征向量,消元過程中，充分利用兩對角塊矩陣的性質(zhì)，可減少計算量,獨立消元,注： 由于方

47、程有奇性，消元后，最后一個方程為0=0, 舍棄最后一個方程，并令最后一個未知數(shù)為1 （ ），即可解出 顯然 （a為任意常數(shù)）也是原方程的解，因此 可用某一值（例如 ）歸一化。,最終，得到 條件下，波數(shù)為 的最不穩(wěn)定擾動波：,擾動波型函數(shù)的分布,Copyright by Li Xinliang,149, 13.2 邊界層轉(zhuǎn)捩的預(yù)測方法,1. 經(jīng)驗公式法,轉(zhuǎn)捩位置,Mach 6 鈍錐邊界層表面的摩擦系數(shù)分布 (Li et al. Phys. Fluid. 22， 025105， 2010. Li et al. AIAA J. 46（11），2899-2913，2008),x,摩阻或熱流,轉(zhuǎn)捩起始點

48、（transition onset),轉(zhuǎn)捩峰（transition peak),充分發(fā)展湍流區(qū),球錐的轉(zhuǎn)捩Reynolds數(shù),邊界層外緣的Mach數(shù),動量厚度定義的轉(zhuǎn)捩Reynolds數(shù),Copyright by Li Xinliang,150,2. eN 方法,LST理論：,積分起始點,擾動波進入中性曲線后，開始增長，局部增長率為,eN理論：擾動波增長到eN倍，即發(fā)生轉(zhuǎn)捩,N值需要由實驗（或經(jīng)驗）確定，通常為811, 即擾動增長4個量級（10000倍）左右。 在不可壓縮及航空領(lǐng)域（亞、跨、超）較為成熟。 在航天領(lǐng)域（高超聲速），還有待檢驗。,不足之處： 未考慮擾動波進入中性曲線前的衰減過程，

49、沒考慮感受性過程。,他人的改進： 蘇彩虹，周恒等 考慮衰減過程 C. H. Su, and H. Zhou, Science in China G, 52 (1):115-123 (2009).,Copyright by Li Xinliang,151,3. PSE （拋物化擾動方程）法,優(yōu)點： 1） 無需平行流假設(shè)； 2） 可處理非線性(N-PSE),Step 1: 得到擾動的控制方程,已知解,線性化,L-PSE,N-PSE,Step 2: 假設(shè)擾動具有波動形式,振幅，沿x方向是個緩變量 （相對y方向而言）,Step 3: 帶入擾動方程，得到振幅 的控制方程,“緩變量”是個很有用的概念，可用

50、來簡化方程,Plantdl邊界層理論就是利用“緩變量”的概念進行簡化的。,LST的 方程是一維的 PSE的 方程是二維的,Step 4: 利用緩變量性質(zhì)，舍棄方程中的橢圓項（為高階小量），得到拋物化的擾動方程,沿x方向推進求解 （類似時間方向的處理），計算量相當(dāng)于一維問題。 （“拋物化”的優(yōu)勢）,非線性項的處理方法與譜方法相似,Copyright by Li Xinliang,152,4. 轉(zhuǎn)捩模型法 （間歇因子模型）,實際粘性系數(shù),層流粘性系數(shù),湍流粘性系數(shù) （由湍流模型給定）,湍流間歇因子 （0表示純層流，1表示純湍流）,方法 1） 根據(jù)經(jīng)驗公式，給定 沿流向的分布 方法2） 給出 的發(fā)展

51、方程，進行求解,Copyright by Li Xinliang,153,常識： 湍流的間歇性,外間歇性： 層流及湍流交替出現(xiàn)的現(xiàn)象,Mach 6 鈍錐邊界層的密度分布 （Li et al. PoF 22, 2010),邊界層有清晰銳利的界面（層流區(qū)、湍流區(qū)“涇渭分明”）,湍流信號,層流信號,內(nèi)間歇性： 湍流脈動的概率密度分布偏離Gauss分布（隨機分布）,概率論（中心極限定理）獨立隨機事件滿足Gauss分布 偏離Gauss分布 湍流不是完全隨機的。,既非確定，又非隨機 湍流的復(fù)雜性,邊界層湍流DNS, Ma=0.7, 2.25,6 擾動: 二維不穩(wěn)定波 + 一對三維斜波 (自然轉(zhuǎn)捩） 壁面吹

52、吸擾動 (Bypass),計算模型：平板邊界層, 13.3 平板邊界層轉(zhuǎn)捩過程的渦動力學(xué)機制,動畫演示： 平板邊界層擬序結(jié)構(gòu)的形成及演化,Q=10的等值面 （流場中的渦結(jié)構(gòu)）,發(fā)卡渦形成及發(fā)展的渦動力學(xué)機制,Copyright by Li Xinliang,157,作業(yè)題 13.1,以不可壓縮槽道湍流為例，推導(dǎo)線性穩(wěn)定性理論的控制方程及邊界條件。要求： 1） 給出擾動量 滿足的線性化控制方程及邊界條件（必須有推導(dǎo)步驟） 2） 推導(dǎo)擾動量振幅滿足的O-S方程及邊界條件 （必須有詳細(xì)推導(dǎo)步驟）,作業(yè)題 13. 2,利用差分法 （最好是MaliK的方法，見本PPT19-24）計算不可壓縮槽道流 (R

53、e=7500) 波長為2p的最不穩(wěn)定T-S波的頻率及時間增長率 （時間模式）。要求給出復(fù)頻率 及擾動波型函數(shù) 的分布曲線。,要求： 必須寫出詳細(xì)的計算過程 （例如類似本PPT 19-21頁，推導(dǎo)出矩陣A,B,C,D的表達式并寫出來） 。本PPT中的公式推導(dǎo)倉促，難免出現(xiàn)問題，切勿照搬使用，請務(wù)必推導(dǎo)！,可使用局部法或全局法，如使用全局法建議使用Q-Z分解法，可使用線性代數(shù)庫或自行編制程序。 使用局部法可使用初值：,建議使用非均勻網(wǎng)格 ， 例如 201網(wǎng)格點效果即可。,計算流體力學(xué)講義 第十四講 湍流與轉(zhuǎn)捩 （2） 李新亮 ；力學(xué)所主樓219； 82543801,知識點：,158,講義、課件上傳

54、至 (流體中文網(wǎng)） - “流體論壇” -“ CFD基礎(chǔ)理論 ” 講課錄像及講義上傳至網(wǎng)盤,Copyright by Li Xinliang,湍流的模式理論（ RANS） 湍流的大渦模擬(LES),Copyright by Li Xinliang,159,知識回顧,1. 流體力學(xué)中的不穩(wěn)定性,Kelvin-Helmholtz; Tollmien-Schlichting; Mack ; Richtmyer-Meshkov ; Reyleigh-Taylor ; Barnard ,擾動的線性化控制方程 不可壓平行流 Orr-Sommerfeld 方程,O-S方程的求解 差分法 （Malik的緊致差分

55、法）,全局法： 解出全部特征值,局部法,Copyright by Li Xinliang,160, 14.1 湍流的工程模式理論RANS,1. 為什么用湍流模型,N-S方程適用于湍流，但其解過于復(fù)雜 如果網(wǎng)格分辨率不夠，數(shù)值解誤差較大,常用方法 進行平均，求解平均量滿足的方程,以不可壓縮為例研究 - 推廣的可壓縮情況,壓縮折角流動,例1： 壓縮折角流動： 如果網(wǎng)格分辨率不足，且不用湍流模型，則分離區(qū)過大 例2： 有攻角機翼流動，如果分辨率不足，且不用湍流模型，則造成“非物理分離”,翼型繞流,平進行均,Copyright by Li Xinliang,161,2. Reynolds平均的N-S方

56、程,以不可壓縮N-S方程為例,時間平均； 空間平均； 系綜平均,脈動,引入平均：,RANS,比N-S方程多了該項,稱為Reynolds應(yīng)力,Copyright by Li Xinliang,162,3. Reynolds平均N-S方程的求解,未知量，必須用已知量表示才能求解,湍流模型,方法 1） Boussinesq 渦粘假設(shè) （常用）,與原先方程的唯一區(qū)別： 改變了粘性系數(shù) 程序?qū)崿F(xiàn)方便,的計算模型： 0方程（代數(shù)模型）： B-L 1方程模型 : S-A 2 方程模型 ：,SST,方法2） Reynolds應(yīng)力模型,給出,的控制方程,可并入壓力項中,Copyright by Li Xinliang,163, 14.2 湍流邊界層的結(jié)構(gòu)及平均速度剖面,內(nèi)層