計算流體力學(xué)課件：有限體積法2

有限體積法（2）知識點：

1通量分裂技術(shù)簡述加速收斂技術(shù)DG方法2知識回顧：有限體積法基本流程無粘項常用方法（流過AB邊的通量）：a.利用周圍點的值，計算出(I+1/2,J)點處的物理量；

直接利用“差分格式”b.利用該處的物理量，計算出流過AB邊的流通量

迎風(fēng)型方法需利用“通量分裂技術(shù)”FVS類：FDS類：利用Riemann解Reimann解：Godunov,Roe,HLL,HLLC利用坐標(biāo)變換，轉(zhuǎn)化為一維Riemann問題§8.1通量分裂技術(shù)簡述多維問題可通過坐標(biāo)旋轉(zhuǎn)，變?yōu)榫植繑U展一維問題該情況下，v,w表現(xiàn)為被動標(biāo)量新坐標(biāo)系下的控制方程：含（一個或多個）被動標(biāo)量的一維Euler方程被動標(biāo)量與一維問題相比，增加了（一個或多個）被動標(biāo)量，對方程性質(zhì)沒有影響。被動標(biāo)量方程分裂形成簡單41.流通矢量分裂(FVS)1)Steger-Warming分裂1維：3維：52)Lax-Friedrichs(L-F)分裂=+優(yōu)點：簡單，計算量小缺點：耗散偏大足夠大L-F分裂=+優(yōu)點：耗散小缺點：導(dǎo)數(shù)間斷S-W分裂63）VanLeer分裂及Liou-Steffen分裂（AUSM類方法）根據(jù)當(dāng)?shù)豈ach數(shù)分裂亞聲速情況下，均勻過渡決定特征傳播方向的關(guān)鍵參數(shù)：當(dāng)?shù)豈ach數(shù)方法1：Liou-Stenffen分裂壓力項單獨處理-117特點：連續(xù)、光滑、無可調(diào)參數(shù)參考文獻：Toro:RiemannSolversandNumericalMethodsforFluidDynamics,section8.4.4Liou:TenYearsinthemakingAUSMfamily,NASATM-2001-2109778方法2：VanLeer分裂（不單獨處理壓力）驗證例：92.通量差分分裂（FDS）利用Riemann解，計算通量1）精確Riemann解（Godunov方法）v,w按照被動標(biāo)量處理

滿足：物理意義為平均增長率102）

Roe近似Riemann解uf(u)uLuRuRoeRoe平均Riemann問題：近似：用平均增長率替代瞬時增長率常系數(shù)線性方程組，求解簡單（相似變換解耦求解）113）

HLL近似Riemann解（Harten,Lax&vanLeer）Ref.：E.F.Toro:RiemannSolversandNumericalMethodsforFluidDynamics,Springer,2009(ThirdEdition)基本原理：雙激波近似t=0t=t0激波1，速度Z1激波2，速度Z2假設(shè)間斷面產(chǎn)生兩道激波，速度分別為ZL,ZR根據(jù)質(zhì)量、動量、能量守恒，容易計算出圖中控制體積內(nèi)的總質(zhì)量、總動量、總能量t0

時刻激波才傳到控制體邊界，因此0到t0時刻，控制體邊界處物理量保持0時刻的值。利用總量，求出圖中控制體內(nèi)的平均值，作為該區(qū)域物理量的近似值124）

HLLC近似Riemann解（Toro）發(fā)展了HLL近似解，用三波模型來近似（如圖）三波近似，左、右波的速T時刻的流動狀態(tài)激波激波接觸間斷模型：左右兩道激波，中間有接觸間斷激波速度已知為:ZL,ZR未知數(shù)（4個）：

方程（6個）：兩道激波的RH關(guān)系式方程多了兩個？（因為假設(shè)激波速度已知）常用方法：去掉兩個方程去掉兩個能量方程，4個未知數(shù)，4個方程，求解求解過程簡單，輕易可給出表達式13最終，HLLC公式為：三波近似，左、右波的速T時刻的流動狀態(tài)激波激波接觸間斷接觸間斷移動速度Ref.：E.F.Toro:RiemannSolversandNumericalMethodsforFluidDynamics,Springer,2009(ThirdEdition)14HLL及HLLC均假設(shè)ZL,ZR已知，實際上它們?nèi)孕枰浪銣?zhǔn)確計算ZL,ZR，實際是計算Riemann精確解，計算量大方法1：直接估算1a:假設(shè)以聲速傳播（Davis）

竟然假設(shè)激波以聲速傳播，太OUT了小常識：激波的傳播速度

激波相對于波前介質(zhì)以超聲速傳播，相對于波后介質(zhì)以亞聲速傳播；

弱激波（Ma趨近于1）以聲速傳播。

1b:左、右兩種狀態(tài)聲速的平均(Davis,Einfeldt)要平均嗎？用Roe平均，激波速度介于波后（相對）聲速與波前（相對）聲速之間，平均是個好思路左、右波速ZL,ZR的計算15Roe平均：1c.Roe平均的修正（Einfeldt）

方法2：基于壓力的波速估算法（Toro）已知中心區(qū)壓力，容易計算波傳播速度中心區(qū)的估算：16§8.2加速收斂技術(shù)1.當(dāng)?shù)貢r間步長法時間步長：受制于最小空間步長邊界層近壁空間網(wǎng)格時間步長小，計算速度慢當(dāng)?shù)貢r間步長：每個點采用不同的時間步長推進不同的點采用不同時間步長對于定常問題，收斂后不影響計算精度；可大幅加速收斂；SeeJ.Blazek:6.1.42.隱格式

i,ji-1,ji+1,ji,j+1i,j+118LU-ADIi,ji-1,ji+1,ji,j+1i,j+1Step1i,ji-1,ji+1,ji,j+1Step2三對角方程，追趕法19Step1:求解Step2:求解LU-SGSi,ji-1,ji+1,ji,j+1i,ji-1,ji+1,ji,j+1i,j+1i,j+1203.多重網(wǎng)格細(xì)網(wǎng)格：精度高粗網(wǎng)格：收斂快利用粗網(wǎng)格加速收斂需克服粗網(wǎng)格降低精度的缺點，不能讓粗網(wǎng)格的結(jié)果“污染”細(xì)網(wǎng)格1）不能讓粗網(wǎng)格獨立計算，細(xì)網(wǎng)格的結(jié)果必須控制粗網(wǎng)格；2）不能將粗網(wǎng)格的計算直接插值到細(xì)網(wǎng)格，而是把粗網(wǎng)格的差量插值到細(xì)網(wǎng)格算法步驟守恒變量殘差下標(biāo)“h”細(xì)網(wǎng)格；

“2h”粗網(wǎng)格；

須注意：21Step1:

細(xì)網(wǎng)格上推進1步，并計算出殘差Step2:

將細(xì)網(wǎng)格上的物理量插值到粗網(wǎng)格（作為粗網(wǎng)格的初值）Step3:計算強迫函數(shù)細(xì)網(wǎng)格的殘差（插值到粗網(wǎng)格）粗網(wǎng)格的殘差細(xì)網(wǎng)格對粗網(wǎng)格的修正1步Euler3步RKLU-SGSStep4:粗網(wǎng)格上進行時間推進使用多步法，推進1步即可；使用單步法，推進2-3步，強迫函數(shù)保持不變Step5:計算粗網(wǎng)格上的差量Step6:將差量插值到細(xì)網(wǎng)格，對細(xì)網(wǎng)格進行修正，得到細(xì)網(wǎng)格上的物理量兩重網(wǎng)格的計算步驟細(xì)網(wǎng)格粗網(wǎng)格細(xì)網(wǎng)格22粗、細(xì)網(wǎng)格之間的插值1）細(xì)網(wǎng)格->粗網(wǎng)格物理量：子網(wǎng)格物理量的體積加權(quán)平均殘差：子網(wǎng)格殘差的簡單相加網(wǎng)格內(nèi)的凈流量=子網(wǎng)格凈流量的相加總質(zhì)量(動量、能量）=子網(wǎng)格質(zhì)量（動量、能量）之和232）粗網(wǎng)格->細(xì)網(wǎng)格二維情況：最近的網(wǎng)格權(quán)重9/16;次近3/16;最遠1/16;三維情況：最近網(wǎng)格權(quán)重27/64;次近9/64;再次近3/64;最遠1/64三重網(wǎng)格上的計算步驟細(xì)網(wǎng)格粗網(wǎng)格最粗網(wǎng)格最粗網(wǎng)格上的強迫函數(shù)：修正后粗網(wǎng)格上的殘差準(zhǔn)備初值過程：粗網(wǎng)格->中等網(wǎng)格->細(xì)網(wǎng)格

244.雙時間步法目的：定常問題的加速收斂技術(shù)

非定常問題1）構(gòu)造隱格式例：2階精度隱格式2）構(gòu)造發(fā)展方程定常問題，非定常處理，推進到收斂3）對t*（偽時間）

進行時間推進，直到收斂推進到收斂可使用定常問題的加速收斂手段：局部時間步長法、隱格式、多重網(wǎng)格法……25有限體積法程序編制1.數(shù)據(jù)結(jié)構(gòu)全局可見的數(shù)據(jù)結(jié)構(gòu)

moduleGlobal_Varimplicitnonereal*8,save::Ma,Re,gamma,Pr,dt0,tt!Mach數(shù)，Reynolds數(shù)

integer,save::Num_Block!塊數(shù)（單塊網(wǎng)格為1）

TYPEBlock_TYPE!localvariables(foreachblock)integer::Block_no,nx,ny,subfacereal*8,pointer,dimension(:,:)::x,y,x1,y1！節(jié)點及網(wǎng)格中心坐標(biāo)

real*8,pointer,dimension(:,:)::s0!控制體（網(wǎng)格）面積

real*8,pointer,dimension(:,:,:)::U,Un!守恒變量（當(dāng)前步及前一步的值）

real*8,pointer,dimension(:,:)::dt!時間步長（局部時間步長法需計算該值）

real*8,pointer,dimesnion(:,:,:)::Res!殘差

EndTYPEBlock_TYPE

TYPE(Block_TYPE),save,dimension(:),allocatable,target::Blockreal*8,save,pointer,dimension(:,:)::d,uu,v,T,p,cc,Amu

real*8,save,pointer,dimension(:,:,:)::Fluxi,Fluxj

real*8,save::Res_max(4),Res_rms(4)endmoduleGoem_Var本程序雖只有1塊；為多塊結(jié)構(gòu)網(wǎng)格預(yù)留；可方便地升級為多塊網(wǎng)格參考文獻：J.Blazek,ComputationalFluiddynamics:PrinciplesandApplications.Elsevier(電子版可在“流體中文網(wǎng)”下載）262.程序結(jié)構(gòu)讀入控制信息；

初始化：讀入并計算幾何信息；初始化流場（包括邊界條件）

進行時間推進c1.計算殘差c2.時間推進c3.邊界條件多步法(如R-K)須重復(fù)多次

d.監(jiān)視（最大或均方根）殘差，輸出流場信息（如氣動力、力矩等）網(wǎng)格中心坐標(biāo)、控制體體積通常用來流初始化四邊形面積：中心坐標(biāo)：273.核心模塊：計算殘差無粘通量+粘性通量無粘通量的計算：1）利用差分格式，計算I+1/2點的左、右值：2）坐標(biāo)旋轉(zhuǎn)，得到坐標(biāo)系的物理量3）利用通量技術(shù)，得到坐標(biāo)系下的無粘通量4）坐標(biāo)旋轉(zhuǎn)，得到原坐標(biāo)系下的無粘通量間斷有限元方法(DiscontinuousGalerkin,DG)一、DG基本原理一維雙曲守恒律方程：a.差分法：直接離散導(dǎo)數(shù)（差分化）j-1jj+1j-1/2j+1/2b.有限體積法：求解（控制體）積分方程§8.3間斷有限元方法簡介已知f函數(shù)的點值，計算“重構(gòu)對函數(shù)”h在j+1/2值的過程已知函數(shù)u的均值，計算該函數(shù)

在j+1/2值的過程c.有限元（Galerkin

）方法1）選擇一組基函數(shù)：含義：在該基函數(shù)空間內(nèi)，方程殘差最小全局定義的連續(xù)函數(shù)j-1jj+1x=ax=b……在全局定義的基函數(shù)空間中，選擇最優(yōu)解的過程通常，選取為正交多項式函數(shù)，例如勒讓德多項式2）設(shè)待解變量，“自由度”3）求解方程：如果

是單位正交基函數(shù)特點：理論基礎(chǔ)好：指定空間內(nèi)的“最優(yōu)解”

解足夠光滑時，有較好的收斂性

不足:不易處理間斷問題；基函數(shù)全局定義,不易推廣到多維、復(fù)雜網(wǎng)格；j-1jj+1x=ax=b……分部積分解出在局部區(qū)域?qū)ふ易顑?yōu)解d.間斷有限元（DiscontinuousGalerkin，DG）方法j-1jj+1j-1/2j+1/21）

選擇局部定義的光滑的基函數(shù)：注：與Galerkin法不同，不是全局積分，而是局部積分；

易于推廣到復(fù)雜網(wǎng)格；

如果測試函數(shù)，則為有限體積法通常為正交多項式函數(shù)2）令3）在局部區(qū)域?qū)ふ易顑?yōu)解32ABC（全局）Galerkin方法，基函數(shù)全局定義基函數(shù)選取復(fù)雜；積分復(fù)雜；基函數(shù)全局光滑，無法處理間斷間斷Galerkin方法，基函數(shù)僅在每個單元內(nèi)定義等參變換：

三角形單元

正三角形

四邊形單元

正方形基函數(shù)選取簡單，積分簡單

僅在單元內(nèi)光滑即可，易于處理間斷

j-1jj+1j-1/2j+1/2j-1jj+1j-1/2j+1/2問題：各單元獨立求解，沒有關(guān)聯(lián)！非物理DG具體求解方法單位正交基分部積分單元內(nèi)部積分邊界值？j-1jj+1j-1/2j+1/2j+1/2點處存在兩個不同的通量值：需替換成統(tǒng)一通量！通量分裂或Riemann解；Steger-Warming,VanLeer,LF,Roe,HLL,HLLC,AUSM,…最終離散形式：優(yōu)勢：高度緊致性；

易于推廣到高精度統(tǒng)一通量有限體積法：每個單元僅存儲一個信息（均值）DG:每個單元存儲多個信息（各階矩）步驟：

Step1:在區(qū)間選取基函數(shù)

，區(qū)間內(nèi)物理量分布為：j-1jj+1j-1/2j+1/2例如：（1階精度）：假設(shè)物理量在區(qū)間內(nèi)均勻分布：

（k=1,每個點僅儲存1個信息）

（2階精度）：線性分布：

（k=2,每個點儲存2個信息）

Step2:計算界面通量：Step3:

進行時間積分得到n+1時刻的值：DG的兩個關(guān)鍵步驟：

1）重構(gòu)（選擇基函數(shù)）

2）計算通量（Riemann解）二、DG在非結(jié)構(gòu)網(wǎng)格中的應(yīng)用DG算法具有很好的緊致性，單個單元信息即可完成重構(gòu)設(shè)單元內(nèi)，物理量分布為：其中，為基函數(shù)，通常可取為多項式函數(shù)等價于：例如：分部積分可計算廣義通量必須保證通量的一致性:單元公共邊處的通量必須相同采用通量技術(shù)（例如Riemann解）具體方法：Step1.選取基函數(shù)常用基函數(shù)：正交多項式滿足：由此確定各系數(shù)注：積分方法1）精確積分；2）數(shù)值積分（精度足夠高，常用Gauss積分）……Step2.計算體積分a.計算b.計算

本時間步值（已知）c.計算積分精確或數(shù)值積分Step3.計算面積分（通量）注意：由于交界面為兩相鄰單元共享，因而在交界面處的通量必須統(tǒng)一；采用通量技術(shù)實現(xiàn)（Riemann解或通量分裂方法）Step4:正交基函數(shù)推進求解（Runge-Kutta）時間推進三、DG方法中的限制器（激波捕捉）1.為何要使用限制器？間斷（大梯度）處，高階重構(gòu)預(yù)測值與真實值相差甚遠需要進行限制（切換到低階重構(gòu)）重構(gòu)示意圖j1階2階3階2.DG限制器的困難僅利用本單元信息，很難識別間斷差分法：利用多個網(wǎng)格點信息，易于識別間斷j-1jj+1jDG:僅利用本網(wǎng)格單元信息，不易識別間斷j點上各階導(dǎo)數(shù)信息；3.常用方法2）混合方法例如DG+WENO先利用周圍點