第四講 可壓縮流體力學(xué)方程組的離散方法 (下).

1、2013復(fù)旦大學(xué)暑期課程：復(fù)旦大學(xué)暑期課程：計算流體力學(xué)計算流體力學(xué) 李新亮李新亮Tel: 82543801力學(xué)所主樓力學(xué)所主樓219 1Copyright by Li Xinliang講義、課件上傳至講義、課件上傳至 （流體中文網(wǎng)）（流體中文網(wǎng)） - “流體論壇流體論壇” -“ CFD基礎(chǔ)理論基礎(chǔ)理論”下載地址：下載地址：https:/ 可壓縮流體力學(xué)方程組的離散方法可壓縮流體力學(xué)方程組的離散方法 （下）（下）1. 有限體積方法初步有限體積方法初步 1） 結(jié)構(gòu)網(wǎng)格有限體積法結(jié)構(gòu)網(wǎng)格有限體積法 2） 多塊結(jié)構(gòu)網(wǎng)格技術(shù)多塊結(jié)構(gòu)網(wǎng)格技術(shù) 3） 加速收斂技術(shù)簡介加速收斂技術(shù)簡介 2Copyright

2、 by Li Xinliang 1 結(jié)構(gòu)網(wǎng)格有限體積法結(jié)構(gòu)網(wǎng)格有限體積法有限體積法主要優(yōu)勢：有限體積法主要優(yōu)勢： 處理復(fù)雜網(wǎng)格處理復(fù)雜網(wǎng)格差分法處理復(fù)雜外形差分法處理復(fù)雜外形 坐標(biāo)變換坐標(biāo)變換),(),(),(zzyyxx321321VVVffftU)(32111fffJfzyx),(),(1zyxJ123123UfffVVVtxyzxyzxxxfffxf1111yyyfffyf2222zzzfffzf3333無粘項需要計算無粘項需要計算9次次導(dǎo)數(shù)，計算量大；導(dǎo)數(shù)，計算量大；守恒性被破壞守恒性被破壞111111111()()()0()()()0()()()0 xxxyyyzzzJJJJJJJJ

3、J（1）321321VVVffftU兩端同乘1J)(32111fffJfzyx12123()xyzfJfff13123()xyzfJfff無粘項僅需計算三次導(dǎo)數(shù)無粘項僅需計算三次導(dǎo)數(shù)Copyright by Li Xinliang4差分法差分法有限體積法有限體積法優(yōu)點優(yōu)點簡單、計算量小、易簡單、計算量小、易于提高精度于提高精度本身包含幾何信息，本身包含幾何信息，易處理復(fù)雜網(wǎng)格易處理復(fù)雜網(wǎng)格不足不足差分離散與幾何解耦，差分離散與幾何解耦，難以處理復(fù)雜網(wǎng)格難以處理復(fù)雜網(wǎng)格復(fù)雜、不易提高精度復(fù)雜、不易提高精度均勻流場數(shù)值解（誤差），未必為數(shù)值解（誤差），未必為0幾何誘導(dǎo)誤差幾何誘導(dǎo)誤差 ：對于均勻場

4、：對于均勻場0Ut312fff1111()()()0 xxxJJJf實際問題：實際問題： 外形復(fù)雜，外形復(fù)雜， 光滑的結(jié)構(gòu)網(wǎng)格生成困難光滑的結(jié)構(gòu)網(wǎng)格生成困難不易保障幾何守恒性不易保障幾何守恒性坐標(biāo)變換函數(shù)坐標(biāo)變換函數(shù)必須足夠光滑必須足夠光滑 否則損失精度否則損失精度Copyright by Li Xinliang51. 基本概念基本概念1） 控制體控制體節(jié)點（中心）型控制體與網(wǎng)格型控制體節(jié)點（中心）型控制體與網(wǎng)格型控制體Copyright by Li Xinliang62） 積分型控制方程積分型控制方程vvtxyxy1212F (U)F (U)F (U)F (U)U在控制體上積分110IJvI

5、JIJdsdsUF nF nt1IJIJUUdV物理含義：物理含義： 控制體內(nèi)總質(zhì)量控制體內(nèi)總質(zhì)量/動量動量/能量的增加能量的增加 = 穿過控制體邊界流入穿過控制體邊界流入的凈質(zhì)量的凈質(zhì)量/動量動量/能量能量無粘通量粘性通量Copyright by Li Xinliang7n2122()()()xyxxnnxnynuuvuupFnF nnnuvvpEp uEp vuuupnvupnEp uF nnF n物理含義：物理含義： 單位時間內(nèi)，無粘流動流過單位時間內(nèi)，無粘流動流過垂直于垂直于n方向的單位面積的質(zhì)量、動量和方向的單位面積的質(zhì)量、動量和能量能量nxyuunvn法向速度Copyright b

6、y Li Xinliang83）有限體積法中物理量的含義）有限體積法中物理量的含義1IJIJUUdV4） 殘差殘差11RIJIJvIJIJdsds UF nF nt殘差殘差=凈通量凈通量=右端項右端項含義：含義： 控制體內(nèi)的平均量控制體內(nèi)的平均量（平均質(zhì)量密度、平均動量密度、（平均質(zhì)量密度、平均動量密度、平均能量密度）平均能量密度）控制體控制體幾何中心處幾何中心處的當(dāng)?shù)孛芏取⒌漠?dāng)?shù)孛芏?、動量密度、能量密度動量密度、能量密度二階精度二階精度近似近似Copyright by Li Xinliang92. 無粘通量的計算無粘通量的計算dsF n常用方法常用方法 （流過（流過AB邊的通量）：邊的通量）

7、： a. 利用周圍點的值，計算出利用周圍點的值，計算出(I+1/2,J) 點處的物理量；點處的物理量； b. 利用該處的物理量，計算出流過利用該處的物理量，計算出流過AB邊的流通量邊的流通量方法方法1： 中心型有限體積法中心型有限體積法1/2,1,12IJI JIJUUU人工粘性項人工粘性項1/2,1/2,()IJIJvisFF UCopyright by Li Xinliang10方法方法2： 迎風(fēng)型有限體積法迎風(fēng)型有限體積法Step 1: 利用（偏）左側(cè)點及（偏）右側(cè)點利用（偏）左側(cè)點及（偏）右側(cè)點的值，計算出的值，計算出I+1/2,J點的值點的值1/2,1,1,(,)LLIJIJI JI

8、JUgUUU1/2,1,2,(,)RRIJI JIJIJUgUUU計算方法：計算方法： 與差分法完全相同與差分法完全相同各種差分格式，均可直接使用各種差分格式，均可直接使用也稱為也稱為“差分格式差分格式”0uuatx1/21/21iiiuuuxx1/21/20:.0:.iiauau該過程稱為該過程稱為“重構(gòu)重構(gòu)”（很多文獻(xiàn)中稱為（很多文獻(xiàn)中稱為“插值插值”）有限體積與有限差有限體積與有限差分共通之處分共通之處, 可直可直接使用差分格式接使用差分格式Copyright by Li Xinliang11常見的差分格式：常見的差分格式：1/2111/211211/2minmod(,)1/2minmo

9、d(,)LIIIIIIRIIIIIIUUUU UUUUUU UU2階NND格式minmod(a,b) : a,b符號相反時取0， 符號相同時取絕對值小的1/2111/212(52)/6(25)/6LIIIIRIIIIUUUUUUUU 3階迎風(fēng)1/211111/4 (1/3)()(1/3)()LIIIIIIUUssUUsUU3階MUSCL格式111221 112()()()()IIIIIIIUUUUsUUUU1/21222121/ 4 (1/3)()(1/3)()RIIIIIIUUssUUsUU1212221212()()()()IIIIIIIIUUUUsUUUUTVD, WENO, GVC,

10、保單調(diào)格式保單調(diào)格式610Copyright by Li Xinliang12重構(gòu)方式：重構(gòu)方式： 原始變量、守恒變量及特征變量原始變量、守恒變量及特征變量以NND格式為例：1/2111/211211/2minmod(,)1/2minmod(,)LIIIIIIRIIIIIIUUUU UUUUUU UU守恒變量重構(gòu)uUvE1/2111/211211/2minmod(,)1/2minmod(,)LIIIIIIRIIIIIIqqqq qqqqqq qq原始變量重構(gòu)uqvp1/2111/211211/2minmod(,)1/2minmod(,)LIIIIIIRIIIIIIVVVV VVVVVV VV特

11、征變量重構(gòu)1/2(1,0,1,2)I kII kVSUk 先算出先算出UI+1/2 (可用可用UI和和UI+1的算術(shù)平均或的算術(shù)平均或Roe平均），再平均），再利用該值算出利用該值算出SI+1/211/21/21/211/21/21/2LLIIIRRIIIUSVUSV （稱為數(shù)值流通量）（稱為數(shù)值流通量） 的含義的含義Copyright by Li Xinliang13重要概念澄清：重構(gòu)與插值0)(xuftuA. 有限差分法：有限差分法：xffxfjjj2/12/1j+1/22/1jf 注意：注意： 與與 f 在在xj+1/2點的值含義不同！點的值含義不同！2/1jf用周圍幾個點的值用周圍幾個

12、點的值 計算計算 的過程稱為的過程稱為“重構(gòu)重構(gòu)”，不能理，不能理解為用解為用 來來插值插值2/1jf jf jf)(2/1jxf記號記號 確實容易混淆，讓人容易聯(lián)想起確實容易混淆，讓人容易聯(lián)想起 。記為。記為 更好些更好些2/1jf)(2/1jxf2/1jf否則，最高只能否則，最高只能達(dá)到達(dá)到2階精度了！階精度了！ 是控制體內(nèi)的平均值是控制體內(nèi)的平均值 （稱為數(shù)值流通量）（稱為數(shù)值流通量） 的含義的含義Copyright by Li Xinliang14重要概念澄清：重構(gòu)與插值0)(xuftuB. 有限體積法：有限體積法：02/12/1xfftujjjj+1/2j-1/2)(2/12/1jj

13、xff2/1jf2/1)(2/12/1jxxjjfxuff確實為確實為f在在xj+1/2點的值點的值 ！ 通常做法：通常做法： 1） 用用 計算出計算出 2） ju2/1ju)(2/12/1jjuff)(2/12/1jjxuuu在xj+1/2點的值！關(guān)鍵：關(guān)鍵： 是用是用 計算計算 （稱為（稱為重構(gòu)重構(gòu)） ，而不是用，而不是用 計算計算 （是標(biāo)準(zhǔn)的（是標(biāo)準(zhǔn)的插值插值）；否則最高也只能達(dá)到）；否則最高也只能達(dá)到2階精度。階精度。 ju2/1ju ju2/12/1)(1jjxxjdxxuxujuju1ju1ju2/1juCopyright by Li Xinliang15Step 2: 利用利用

14、， 計算出通量計算出通量 n1/2,LIJU1/2,LIJU 方法方法1：流通矢量分裂：流通矢量分裂 （FVS）1/21/21/2()()LRIIIUUfff方法方法2：通量差分分裂：通量差分分裂 (FDS)通過通過Riemann解，獲得通量解，獲得通量可利用近似可利用近似Riemann階（階（Roe, HLL, HLLC)11/2RLRLISS11ff(U )f(U )(UU )22Roe格式：格式：1111,pvu2222,pvuxyxy方法方法3： AUSM方法方法壓力項單獨處理，其余項類似壓力項單獨處理，其余項類似Van Leer 分裂分裂Copyright by Li Xinlian

15、gn小技巧：小技巧： 利用局部坐標(biāo)系，計算通量時簡化利用局部坐標(biāo)系，計算通量時簡化為為1維問題維問題ABx yxy1111,pvu2222,pvuyxyuvpuvp x y1111,u v p2222,u vp“擴展的”一維問題xx,y坐標(biāo)系下的Riemann問題（二維）x,y 坐標(biāo)系下的Riemann問題（一維）x, y 坐標(biāo)系下，切向速度坐標(biāo)系下，切向速度v表現(xiàn)相當(dāng)于表現(xiàn)相當(dāng)于被動標(biāo)量被動標(biāo)量Copyright by Li Xinliang173. 粘性通量的計算粘性通量的計算110IJvIJIJdsdsUF nF nt1111211120( )PrRevpUCTuvx F12222122

16、20( )PrRevpUCTuvyFjidivVxujixuxuiiijjiij),322(),(關(guān)鍵問題：關(guān)鍵問題： 計算計算yx , ,u v T 對對 的導(dǎo)數(shù)同樣計算的導(dǎo)數(shù)同樣計算Copyright by Li Xinliang18方法方法1. 利用利用Jacobian變換變換yx ,1/2,1/2,xxIJIJx1,1/2,IJI JIJ1/2,11/2,11/2,11,1,11,11214IJIJIJI JIJI JIJxyxyJyJxJyJxJx yy x , x y, 1/2,1/2,yyIJIJy方法方法2. 利用利用Green積分公式計算積分公式計算1/2,1/2,IJIJd

17、sdn1/2,1/2,1IJIJdsn最終公式與方法最終公式與方法1相同相同Copyright by Li Xinliang19 2 邊界條件及網(wǎng)格分塊邊界條件及網(wǎng)格分塊1. 邊界條件的處理方法邊界條件的處理方法常用方法：常用方法： 虛網(wǎng)格虛網(wǎng)格設(shè)置設(shè)置1或或2層虛網(wǎng)格（如果采用高層虛網(wǎng)格（如果采用高精度格式，需設(shè)置更多重虛網(wǎng)格）精度格式，需設(shè)置更多重虛網(wǎng)格）虛網(wǎng)格點上的值根據(jù)邊界條件給定虛網(wǎng)格點上的值根據(jù)邊界條件給定 絕熱固壁邊界條件絕熱固壁邊界條件11111111ppuuvv 超聲速入口邊界條件超聲速入口邊界條件1111uuvvpp超聲速出口邊界條件超聲速出口邊界條件11111111ppu

18、uvv 0111()2qqq( , , , )Tqu v pCopyright by Li Xinliang20亞聲速入口邊界條件亞聲速入口邊界條件外部計算域內(nèi)部(,)Tquvp1q0q1q0111()2qqq( , , , )Tqu v p邊界特征分析： 內(nèi)點提供一個邊界條件，外部提供其余邊界條件n0111()2pppcVVn調(diào)節(jié)作用， 內(nèi)點速度降低時，增加入口壓力； 內(nèi)點速度增加時，降低入口壓力；2000000()/()/()/xyppcuunppcvvnppcCopyright by Li Xinliang21亞聲速出口邊界條件亞聲速出口邊界條件外部計算域內(nèi)部1q0q1q0111()2q

19、qq( , , , )Tqu v p邊界特征分析：特征分析： 外部提供一個邊界條件；外部提供一個邊界條件； 通常情況下，指定背壓通常情況下，指定背壓n0pp如果邊界處壓力比內(nèi)點高，則邊界如果邊界處壓力比內(nèi)點高，則邊界處密度增加，速度降低；處密度增加，速度降低； 如果邊界壓力比內(nèi)點低，則邊界處如果邊界壓力比內(nèi)點低，則邊界處密度降低，速度升高；密度降低，速度升高；2010101100110()/()/()/xyppcuunppcvvnppcpCopyright by Li Xinliang222. 多塊網(wǎng)格多塊網(wǎng)格復(fù)雜外形情況下，單塊結(jié)構(gòu)網(wǎng)格很難實現(xiàn)，需要由多塊網(wǎng)格實現(xiàn)復(fù)雜外形情況下，單塊結(jié)構(gòu)網(wǎng)格

20、很難實現(xiàn)，需要由多塊網(wǎng)格實現(xiàn)Copyright by Li Xinliang231） 對接網(wǎng)格對接網(wǎng)格對接網(wǎng)格示意圖對接網(wǎng)格示意圖相鄰的網(wǎng)格塊共享交界網(wǎng)格點相鄰的網(wǎng)格塊共享交界網(wǎng)格點處理方法：處理方法： 利用虛網(wǎng)格傳遞信息；利用虛網(wǎng)格傳遞信息； 利用連接信息，將內(nèi)點的物理量利用連接信息，將內(nèi)點的物理量復(fù)制給對應(yīng)的虛網(wǎng)格點；復(fù)制給對應(yīng)的虛網(wǎng)格點；Block2 的值的值 Block 1虛網(wǎng)格的值虛網(wǎng)格的值Copyright by Li Xinliang242） 重疊網(wǎng)格重疊網(wǎng)格重疊網(wǎng)格示意圖重疊網(wǎng)格示意圖處理方法：處理方法： 使用虛網(wǎng)格交換信息；使用虛網(wǎng)格交換信息； 利用插值，用網(wǎng)格內(nèi)點利用插值，

21、用網(wǎng)格內(nèi)點的信息計算出另一套網(wǎng)格的信息計算出另一套網(wǎng)格虛網(wǎng)格上的信息。虛網(wǎng)格上的信息。Copyright by Li Xinliang25 3 加速收斂技術(shù)加速收斂技術(shù)1. 當(dāng)?shù)貢r間步長法當(dāng)?shù)貢r間步長法時間步長時間步長：*/minxCFLt受制于最小空間步長受制于最小空間步長邊界層近壁空間網(wǎng)格邊界層近壁空間網(wǎng)格1wy641010 x時間步長小，計算速度慢時間步長小，計算速度慢當(dāng)?shù)貢r間步長：當(dāng)?shù)貢r間步長： 每個點采用不同的時間步長推進每個點采用不同的時間步長推進IJIJRtUIJnIJnIJRtUU1IJIJnIJnIJRtUU1不同的點采用不同的點采用不同時間步長不同時間步長對于定常問題，收斂

22、后不影對于定常問題，收斂后不影響計算精度；響計算精度；可大幅加速收斂；可大幅加速收斂；See J. Blazek: 6.1.40yxt(U)f(U)fU21011yxtnn)(Uf)(UfUU21n1nn1nnUUQ)()(21212111yxtyxtnnnnnnffffffQnRHSBQAQQntyxtnn1,11,1ijijijijiji jijijiji jijabcdeQQQQQRRHSQQAQQAQtBytxtBytxtytxtjijijijijijijijiBAij1 1,1, 1, 11,1, 1, 1*2. 隱格式隱格式 i,ji-1,ji+1,ji,j+1i,j+1Copyr

23、ight by Li Xinliang27LU-ADI1,11,1ijijijijiji jijijiji jijabcdeQQQQQRijjiijjiijijijRQdQbQa, 1, 1i,ji-1,ji+1,ji,j+1i,j+1Step 1,1,1,1ijijiji jiji jiji ja Qc Qe Qa Qi,ji-1,ji+1,ji,j+1Step 2三對角方程，追趕法Copyright by Li Xinliang281,11,1ijijijijiji jijijiji jijabcdeQQQQQRStep 1: 求解1,1ijijijijiji jijadeQQQR1,1(

24、)/ijijijijiji jijdeaQRQQ111111/ aQRStep 2: 求解1,1ijijijijiji jijijabcaQQQQ1,1()/ijijijijijiji jijabcaQQQQ/mnmnmnaQRLU-SGSi,ji-1,ji+1,ji,j+1i,ji-1,ji+1,ji,j+1i,j+1i,j+1ULDLUD1)()(1UDDLDCopyright by Li Xinliang293. 多重網(wǎng)格多重網(wǎng)格細(xì)網(wǎng)格：細(xì)網(wǎng)格： 精度高精度高粗網(wǎng)格：粗網(wǎng)格： 收斂快收斂快利用粗網(wǎng)格加速收斂利用粗網(wǎng)格加速收斂需克服粗網(wǎng)格降低精度的缺點，需克服粗網(wǎng)格降低精度的缺點，不能讓粗

25、網(wǎng)格的結(jié)果不能讓粗網(wǎng)格的結(jié)果“污染污染”細(xì)網(wǎng)格細(xì)網(wǎng)格1） 不能讓粗網(wǎng)格獨立計算，細(xì)網(wǎng)格的結(jié)果不能讓粗網(wǎng)格獨立計算，細(xì)網(wǎng)格的結(jié)果必須必須控制控制粗網(wǎng)格；粗網(wǎng)格；2） 不能將粗網(wǎng)格的計算直接插值到細(xì)網(wǎng)不能將粗網(wǎng)格的計算直接插值到細(xì)網(wǎng)格，而是把格，而是把粗網(wǎng)格的差量粗網(wǎng)格的差量插值到細(xì)網(wǎng)格插值到細(xì)網(wǎng)格算法步驟算法步驟1hhhddt WR守恒變量守恒變量殘差殘差下標(biāo)下標(biāo)“h”細(xì)網(wǎng)格；細(xì)網(wǎng)格； “2h”粗網(wǎng)格；粗網(wǎng)格； 須注意：須注意：Copyright by Li Xinliang30Step 1: nhW1nhW1nhR細(xì)網(wǎng)格上細(xì)網(wǎng)格上推進推進1步步， 并計算出殘差并計算出殘差Step 2: 將細(xì)網(wǎng)

26、格上的將細(xì)網(wǎng)格上的物理量插值到粗網(wǎng)格物理量插值到粗網(wǎng)格 （作為粗網(wǎng)格的初值）（作為粗網(wǎng)格的初值）12)0(2nhhhhIWWStep3: 計算計算強迫函數(shù)強迫函數(shù)21(0)22()hnFhhhhQIRR細(xì)網(wǎng)格的殘差細(xì)網(wǎng)格的殘差（插值到粗網(wǎng)格）（插值到粗網(wǎng)格）粗網(wǎng)格的殘差粗網(wǎng)格的殘差細(xì)網(wǎng)格對粗網(wǎng)細(xì)網(wǎng)格對粗網(wǎng)格的修正格的修正1步Euler3步步RKLU-SGS2221() )hhFhhdQdt WRStep4: 粗網(wǎng)格上進行時間推進粗網(wǎng)格上進行時間推進使用多步法，推進使用多步法，推進1步即可；步即可；使用單步法，推進使用單步法，推進2-3步，步，強迫函數(shù)保持不變強迫函數(shù)保持不變Step 5: 計算粗網(wǎng)格上的計算粗網(wǎng)格上的差量差量1(0)222nhhhWWWStep 6: 將差量將差量插值到細(xì)網(wǎng)格插值到細(xì)網(wǎng)格，對細(xì)網(wǎng)格進行，對細(xì)網(wǎng)格進行修正修正，得到細(xì)網(wǎng)格上的物理量，得到細(xì)網(wǎng)格上的物理量hhhnhnhWIWW2211兩重網(wǎng)格的計算步驟兩重網(wǎng)格的計算步驟細(xì)細(xì)網(wǎng)網(wǎng)格格粗粗網(wǎng)網(wǎng)格格細(xì)細(xì)網(wǎng)網(wǎng)格格Copyright by Li Xinliang31粗、細(xì)