高等計(jì)算流體力學(xué)-06

1、第六講Riemann問題的近似解1十八、Riemann問題的近似解Approx. Riemann Solver2為什么要研究Riemann問題的近似解? 精確求解Riemann問題計(jì)算量大；某些雙曲型守恒律的Riemann問題無精確解！ 由于在時(shí)間積分、通量積分、重構(gòu)等步驟已經(jīng)引入了各種近似，精確求解Riemann問題并不會導(dǎo)致整體精度的提高（但對一維問題而言，基于Riemann問題精確解的Godunov格式的穩(wěn)健性和精度的平衡可能較好。）3Riemann 問題4擴(kuò)展一維系統(tǒng)擴(kuò)展一維系統(tǒng)初始條件初始條件 的解對應(yīng)有限體積格的解對應(yīng)有限體積格式的數(shù)值通量式的數(shù)值通量Riemann問題的近似解法1

2、：HLL方法 Harten-Lax-Van Leer， Einfeldt（1988） 基本思路：雙波近似（不考慮接觸間斷）5tRxLxHLLULURULdxSdtRdxSdt為了避免奇異性，不對左右波的類型作假定。左右波只表示間斷初始值的影為了避免奇異性，不對左右波的類型作假定。左右波只表示間斷初始值的影響范圍；響范圍；優(yōu)點(diǎn)：簡單；優(yōu)點(diǎn)：簡單；已知問題：對接觸間斷分辨率低。已知問題：對接觸間斷分辨率低。假定SL，SR已知積分關(guān)系 取 積分區(qū)域 方程6相容條件相容條件7 把前述積分區(qū)域分為 左右兩部分，分別積分8相容條件相容條件HLL Riemann Solver9一般情況一般情況通量計(jì)算時(shí)考慮

3、通量計(jì)算時(shí)考慮x/t=0波速的確定10Roe平均平均不推薦！不推薦！推薦！推薦！還有其他方法，見還有其他方法，見Toro的書。的書。HLL的變形11Rusanov（Lax）LaxRiemann問題的近似解法2：HLLC方法 Toro 基本思路：三波近似（在HLL基礎(chǔ)上加接觸間斷）1213相容條件相容條件14超定還是欠定？超定還是欠定？假定假定SL，SR已知，關(guān)鍵問題是求已知，關(guān)鍵問題是求S*利用上二式的第一、二兩利用上二式的第一、二兩項(xiàng)項(xiàng)中間波是接觸間斷中間波是接觸間斷接觸間斷接觸間斷滿足相容條件滿足相容條件15*22*()()()()()()LLLLLLLLLLLLSuSuSuupSuup*

4、22*()()()() ()() LLLLLLLLLLLLLuuSSppSuuuu22* ()() ()()()()()()()()()()()()()LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLppS SuupSuSuSuSupSuSuSuSupSuSuSupSuuu16變形變形117變形變形2kS18變形變形3假定假定*狀態(tài)有單一壓力狀態(tài)有單一壓力 ，不不直接用直接用Riemann問題的近似解法3：Roe方法19考慮擴(kuò)張一維守恒型考慮擴(kuò)張一維守恒型Euler 方程方程 ：0UFtx 對應(yīng)的擬線性形式為：對應(yīng)的擬線性形式為：0UUAtx Roe 把

5、上式中的矩陣把上式中的矩陣A 用一個(gè)常矩陣用一個(gè)常矩陣 代替代替 ，即，即 (,)LRAA U U,:LRU URiemann問題左右狀態(tài)0UUAtx則則0UFtx守恒形式守恒形式FAU基本思路基本思路:在控制體界面處把方程線性化在控制體界面處把方程線性化通過求解線性通過求解線性Riemann問題計(jì)算通量問題計(jì)算通量00( ,0)0txLRUAUUifxU xUifx線性線性Riemann問題問題關(guān)鍵問題關(guān)鍵問題: 如何線性化如何線性化?稱為稱為Roe Jacobian Matrix, Roe要求它滿足下列條件：要求它滿足下列條件：( ,)A U UA()()()RLRLF UF UA UU（

6、A） 雙曲性雙曲性 即即 的特征值均為實(shí)數(shù)，且存在完備的左右特征向量；的特征值均為實(shí)數(shù)，且存在完備的左右特征向量；（B） 相容性相容性 （C） 守恒性守恒性 如果左右狀態(tài)滿足激波關(guān)系，則如果左右狀態(tài)滿足激波關(guān)系，則 ()()()()RLRLRLF UF UD UUA UU D是是 的特征值，此時(shí)線性的特征值，此時(shí)線性Riemann問題的解與問題的解與Riemann問題精確解相同！問題精確解相同！ 根據(jù)這些條件可根據(jù)這些條件可確定確定 （先假定（先假定為為已知）已知）(,)( )LRA U UA U20A()m(1)(2)()m(),()mmlr的特征值為的特征值為 左、右特征值向量為左、右特征

7、值向量為 1ALL由由雙曲性，可知：雙曲性，可知： 0UULLtx由于由于A是常矩陣，所以是常矩陣，所以L, L-1也是常矩陣也是常矩陣 定義特征變量定義特征變量V=LU0VVtx(1)(1)(2)(2)()()mmlUvlUvVvlU( )( )( )0jjjvvtx2100( ,0)0txLRUAUUifxU xUifx( )( )( )( )( )( )( )( )0(1,2,)0( ,0)0jjjjLjLjjRjRvvtxjmvlUifxvxvlUifx( )( )( )( )( )( , )jjLjjjRlUvx tlU 當(dāng)x/t 當(dāng)x/t22為了計(jì)算通量，只需為了計(jì)算通量，只需在在

8、x=0處處的解：的解：( )( )( )( )( )( , )jjLjjjRlUvx tlU 當(dāng)x/t 當(dāng)x/t)( )( )1/2)( )0(0)0jjLjijjRlUvlU（若若( )( )1/2( )( )11( , )()(/ )()22jjijLRjRLvx tlUUsignx t lUU1/211( / )()( )( / ) ()22iLRRLVx tL UUsignx t I L UU 111( , )()( / ) ) ()22LRRLU x tUUL signx t I L UU11/211(0, )()() ()22iLRRLUUtUULsignL UU23FAU注意到注

9、意到11/211(0, )()() ()22iLRRLUUtUULsignL UU線化方程線化方程的數(shù)值通量的數(shù)值通量11/2011()22iLRRLFFAUAULL UU如何計(jì)算原始如何計(jì)算原始Euler方程的數(shù)值通量？方程的數(shù)值通量？2425txRxLxRTSLTST未擾動未擾動的值的值, 0, LRxxT( , )()RLTSRRLLLRTSU x T dxT S US UFF ,00,LxT 0,0,RxT001( , )RTSRRRRFFS UU x T dxTRoe格式（格式（Euler方程）：方程）：(1)()(2)(1),mLRmSS001( , )LLLLLTSFFS UU

10、x T dxTEuler方程（無線性化，但假定與線化方程有相同的解的結(jié)構(gòu)）方程（無線性化，但假定與線化方程有相同的解的結(jié)構(gòu)）Euler方程（線性化）方程（線性化）001( , )LLLLLTSFAUS UU x T dxT001( , )RTSRRRRFAUS UU x T dxT000LRFFF1/2000iLRFFFF001( , )RTSRRRRFFS UU x T dxT001( , )LLLLLTSFFS UU x T dxT001( , )LLLLLTSFAUS UU x T dxT001( , )RTSRRRRFAUS UU x T dxT00LLLLFFFAU00RRRRFFF

11、AU1/2011()()22iLRLRFFFFAUAU11/2011()22iLRRLFFAUAULL UU11/211()()()22iLRRLFF UF ULL UU2611/211()()()22iLRRLFF UF ULL UU編程實(shí)施方案：編程實(shí)施方案：()1-miRLmlUUUr(1)(2)(1)(2)()()(1)(2)(1)(2)()()1(),mmmimmlmllURL UrrrUllUlUrrrrlU( ).iilU( ) ( )( )1( )( )mjjjrl 事實(shí)上，對于任事實(shí)上，對于任意的意的m維列向量維列向量11( )( )1( )( )1()()RLmjjRLjm

12、jjjjjLL UUrlLL UUr ( )1/2( )11122mjiLRjjjFFFr 27( )1/2( )11122mjiLRjjjFFFr ( )( ),jjjr只需計(jì)算：Roe-Pike方法方法: 無需計(jì)算無需計(jì)算 ,直接計(jì)算上述三個(gè)量。直接計(jì)算上述三個(gè)量。考慮三維問題對應(yīng)的擴(kuò)展一維問題，有：考慮三維問題對應(yīng)的擴(kuò)展一維問題，有：A右特征向量右特征向量（注意記號變化）（注意記號變化）(,)()(,)LRLRA UUA UUUU28利用：利用：29關(guān)鍵問題：如何確定關(guān)鍵問題：如何確定Roe-Pike方法：求平均狀方法：求平均狀態(tài)態(tài)，使下列關(guān)系成立：使下列關(guān)系成立：(,)LRUU U其中

13、：其中：11( )( )1( )( )1()()RLmjjRLjmjjjjjLL UUrlLL UUr 11( )( )1( )( )1()()()RLRLmjjRLjmjjjjA UULL UUrlLL UUr ()()()RLRLF UF UA UU3031,ii j iAii當(dāng)當(dāng)中所有的對角元素中所有的對角元素都充分大于零時(shí)，格式粘性項(xiàng)都充分大于零時(shí)，格式粘性項(xiàng)是正定的，從而保證了計(jì)算格式有足夠的數(shù)值耗散，是正定的，從而保證了計(jì)算格式有足夠的數(shù)值耗散，變?yōu)榱?。此時(shí)，格式粘性項(xiàng)的系數(shù)矩陣是半正定的，變?yōu)榱?。此時(shí)，格式粘性項(xiàng)的系數(shù)矩陣是半正定的，某些特征波的數(shù)值粘性過某些特征波的數(shù)值粘性過小

14、，從而有可能得到非物理解。所以，一般情況下，非物理解出現(xiàn)在某些小，從而有可能得到非物理解。所以，一般情況下，非物理解出現(xiàn)在某些非常接近零的膨脹流動區(qū)域，表現(xiàn)為類似于非常接近零的膨脹流動區(qū)域，表現(xiàn)為類似于“膨脹激波膨脹激波”的系數(shù)矩陣的系數(shù)矩陣可以得到物理解。但是，在某些流動區(qū)域，如音速點(diǎn)附近，一些可以得到物理解。但是，在某些流動區(qū)域，如音速點(diǎn)附近，一些（過激波熵減少）的數(shù)值解的跳躍。（過激波熵減少）的數(shù)值解的跳躍。 1/21/21/21/21/211 ()()()22LRRLjjjjjAFF UF UUU,ii jARL 熵修正熵修正32進(jìn)行限制，使其不致過小。例如，可以指定一個(gè)小正數(shù)進(jìn)行限制

15、，使其不致過小。例如，可以指定一個(gè)小正數(shù)“熵修正（熵修正（entropy fixentropy fix）”是避免非物理解出現(xiàn)的一個(gè)有效手段。一種簡單的熵修是避免非物理解出現(xiàn)的一個(gè)有效手段。一種簡單的熵修正方法正方法是對是對ii*222iiiiifotherwise*ii使限制后的使限制后的*i（記為（記為）為）為用用代替代替進(jìn)行數(shù)值通量的計(jì)算，對于避免進(jìn)行數(shù)值通量的計(jì)算，對于避免FVSFVS格式和格式和RoeRoe格式出現(xiàn)格式出現(xiàn)非物理解都是有效的。非物理解都是有效的。3334評述及擴(kuò)展 基于Riemann問題精確或者近似解計(jì)算數(shù)值通量的方法稱為Godunov類型格式或者通量差分裂（Flux

16、Difference Splitting，F(xiàn)DS）格式。 優(yōu)點(diǎn)： Godunov、Roe、HLLC對激波和接觸間斷都有很好的分辨能力；對粘性問題的邊界層和剪切層計(jì)算效果較好。 缺點(diǎn)：有可能出現(xiàn)奇偶失連或者Carbuncle現(xiàn)象 35 其他迎風(fēng)型的通量計(jì)算方法還包括： Flux Vector Splitting：Steger-Warming， Van Leer，Lax（Rusanov） 不出現(xiàn)奇偶失連或者Carbuncle現(xiàn)象，可計(jì)算強(qiáng)激波 對接觸間斷、邊界層和剪切層的分辨率低 壓力-對流分別分裂的方法：AUSM類（Liou等），CUSP（Jameson）， LDFSS（Edwards）等等，這

17、些方法均試圖結(jié)合FDS和FVS的某些優(yōu)點(diǎn)而克服其缺點(diǎn)。36 解決FDS方法奇偶失連或者Carbuncle現(xiàn)象的一種有效方法是采用旋轉(zhuǎn)Riemann求解器 任玉新（2003）， Nishikawa等（2008），Wu（2010） Godunov、Lax、HLL等方法滿足離散的熵條件，但分辨率較低。 沒有一種在任何情況下都最優(yōu)的方法！3738 FVS方法方法 （流通矢量分裂（流通矢量分裂 逐點(diǎn)分裂）逐點(diǎn)分裂） fff 具體方法：具體方法： Steger-Warming 分裂分裂 Lax-Friedrichs分裂分裂 Van Leer分裂：分裂：2kkkwcucuucucuu232221321321

18、)(2)(2) 1()()() 1(2) 1(22)(f2/ )(*Uff根據(jù)當(dāng)?shù)馗鶕?jù)當(dāng)?shù)豈ach數(shù)分裂數(shù)分裂保證保證 的的Jocabian陣特征值為正，陣特征值為正， 的為負(fù)的為負(fù)ffUAf)2/12/12/1RjLjnj(Uf(Uff一個(gè)參數(shù)，反映全部特征一個(gè)參數(shù)，反映全部特征39AUSM:Liou-Steffen分裂分裂)()(2200)()(pcFFpupEuuupEpuuf(U)對流項(xiàng)壓力項(xiàng)思路：思路： 決定特征的關(guān)鍵參數(shù)決定特征的關(guān)鍵參數(shù) 當(dāng)?shù)禺?dāng)?shù)豈ach數(shù)數(shù)1 1 , 00 , 11cuMa超音速，超音速，x-方向方向超音速，超音速，x+方向方向0, 0, 0321cucuu32

19、1,0, 0, 0321因此，對因此，對Mach數(shù)進(jìn)行分裂更為簡潔！數(shù)進(jìn)行分裂更為簡潔！1當(dāng)01當(dāng)4/) 1(1當(dāng)2MMMMMMaHauaMFc)(114/) 1(102MMMMMM1012/ )1 (1MMMpMpp112/ )1 (10MpMMpMp顯然：顯然： pppMMMfff010paHauaMf參考文獻(xiàn)：參考文獻(xiàn)：Toro: Riemann Solvers and Numerical Methods for Fluid Dynamics, section 8.4.4Liou: Ten Years in the making AUSM family, NASA TM-2001-210

20、977類似類似 Van Leer分裂，但壓力單獨(dú)處理分裂，但壓力單獨(dú)處理MM保證光滑過渡保證光滑過渡M=1Journal of Computational Physics 208 (2005) 527569, Kim&Kim4041三維問題的旋轉(zhuǎn)不變性三維問題的旋轉(zhuǎn)不變性100000coscoscossinsin00sincos000sincossinsincos000001fT,fxfyfzfnnnncoscoscossinsinxfyfzfnnn1xfxfyfyfzfzfffnnnFGHTF T Q十九、十九、A Shock Instability Free FDS Scheme

21、Based on Rotated Riemann Solvers1.Introduction Flux-Difference Splitting Schemes: Advantage and Disadvantage Resolve shock waves and contact discontinuities in higher resolution. Admit non-physical numerical solution such as the rarefaction shock for some Riemann solvers. Entropy fix is needed in th

22、is case. May produce the shock instability in multi-dimensional problems using the grid-aligned method: the “oddeven decoupling” and the “carbuncle” phenomena The first order grid-aligned Roe scheme The first order rotated Roe scheme with the direction of upwinding determined by the velocity-differe

23、nce vector. The Cure for the Shock Instability Detect the cell faces deemed as susceptible to the shock instabilities; Modify the flux functions with some special entropy fix procedures, or Replace the flux functions with more dissipative flux functions. The Motivations of the Present Work Shock ins

24、tability: a particular problem for multi-dimensional computations. It is helpful to take the multi-dimensional effects into consideration. The Rotated Riemann Solver is the simplest Riemann solver take multi-dimensional effect into account. What We Have Done Reformulation and generalization of the R

25、otated Riemann Solvers; Determination of the direction of rotation: crucial for the Rotated Riemann Solver to be shock instability free; Application to Roes Riemann Solver and the construction of the corresponding first and second order finite volume schemes; Numerical experiments.2. The Generalized

26、 Rotated Flux Function 2.1 The Gas Dynamic Euler Equation0txyUFG 2.2 The Finite Volume Scheme,0i ji jdxdydltUH nHFiGj For 2D Problems,41,44,11,11i jki ji ji jkIki ji ji ji jkkkIkki ji jdldldxdydxdydxdydllttH nH nUUUUH nH n 2.3 The Grid Aligned Riemann Solvers Introduce the rotation matrix11000100000

27、00000000010001xkykxkykkkykxkykxknnnnnnnnTT Then we have1()kkkxkkykkknnH nFGT F U4,11,1()i jkkkki jlt UT F Ukkkcossincossin()kkkTkkuvEuuvvuv UT UUT UFF U The numerical flux is computed by solving the following Riemann problem in terms of the augmented one-dimensional system ()()00( ,0)0kktxLkkRkif xx

28、if xUFUUU()kkkkkUT UFF U()kF U()()()kktxkyUFG0 2.4 Generalized Rotated Riemann Solvers4,1,11,(2): unit vectori jkkkki jMmmkkkmmkltM UHnnnn Decomposing the outward unit vector normal to face Ik into the summation of two or more vectors Introduce the rotation matrixes corresponding to . Then we have m

29、kn1114,111,()()()1()(),MMmmmmmmkkkkxkkykkkkkmmMi jmmmkkkkkmi jmmmmkkkkknnlt HnFGTF T UUTFFF UUT UmkT In each direction, is evaluated by solving the following Riemann problem()()00(,0)01,2,.,mmmkktxLmkmkmRkmif xxif xmMUFUUU()mmkkFF U We call it the Generalized Rotated Riemann Solvers. If it is called

30、 the Rotated Riemann Solvers. 122,0,kkM nn11kkn22kkn In what follows, we consider the Rotated Riemann Solvers only. The determination of will be discussed later.12,kkn n3. The Rotated Form of Roe Riemann Solver 3.1 The Flux Function of Rotated Roe Riemann Solver24111/2,1/2,1/2,1/2,1/2,2211()()LRmmmm

31、ijijijijlllijmlR21/2,1/2,1/21/2,1,(),(),()mmmmmLLRRkkkkkkijijiijmsignnnnn U U1234,mmmmmmmmmxymyxmmmqunvnrunvnqaqaq T1T2T322T1421,1,0,1, , , ()mmmxymmmmxymmmmyxmmuanvanHaquanvanHaqananaru vuvRRRR 22221()21()211()mmmmpa qapa qaraapa 3.2 The Choice of Upwinding Direction (1,2)mkm n21112221,kkkkkkkkmmkk

32、km nknkijn nn nnn If is known1kn How to choose ?1kn1/2,1/2,11/2,221/2,1/2,1/2,11/2,1/2,1/2,1/2,11/2,1/2,1/2,221/2,1/2,1/2,11/2,11.()()2.3.()()4.:ijijijijijijijijijijijijijijijijijiuvuvppuvuvifuvOur Choiceu ijnnijnijnn/2,1/2,221/2,1/2,()()jijijijvotherwiseuv ijCant not remove the shock instabilityThe

33、 shock instability remains but is confined in small regions near the shock wave Why we make such choice? Physical significance: If there is a shock or a shear-wave at the cell interface, the shock will propagate in the direction of and the shear-wave will move in the direction of . A simple analysis

34、 of the dissipation properties of the Rotated Riemann Solver indicates that at locations where shock instabilities may happen, the use of the rotated flux function effectively increases the numerical dissipations for linear degenerate field of the Euler equations. Numerical experiments show that the

35、 shock instabilities can be eliminated completely using this procedure.1kn2kn 3.3 Time Stepping Scheme and MUSCL Interpolation for 2nd Order Scheme(0),(1)(0)(0),(2)(0)(1)(1)11,221(2),()().ni ji ji ji ji ji ji ji ji jni ji jtt UUUULUUUULUUU,1/2,1,1,1,1/2,1,1,1,(1)(1)433(1)(1)433iiii ji ji jLiiiji ji

36、ji jiiiijijijRiiijijijijssssssUUUU 3.4 Efficiency: About 1.5 times more costly than the grid-aligned method. A proposal to implement the Rotated Riemann Solver more efficiently: 221/2,1/2,1/2,11/2,1/2,1/2,221/2,1/2,()()()()ijijijijijijijijifuvuvotherwiseuv nnij4. Numerical Results 4.1 Odd-Even Grid

37、Perturbation Problem We consider a plane shock wave that propagates downstream in a straight duct at the speed of . This problem is computed on a grid of 80020 cells. Each cell is a square with unit side, except those on the centerline where the grid is perturbed in the following manner: 3,3,1010i j

38、midi jmidi jmidyforievenyyforiodd6Mxy50055060065070051015xy50055060065070051015xy50055060065070051015xy50055060065070051015xy50055060065070051015xy3500355036003650370051015 the rotated Roe scheme with the direction of upwinding determined by velocity vector, t=100 the grid-aligned Roe scheme t=100 t

39、he rotated Roe scheme with the direction of upwinding determined by pressure gradient, t=100 the rotated Roe scheme with the direction of upwinding determined by velocity-magnitude gradient, t=100 the rotated Roe scheme with the direction of upwinding determined by velocity-difference vector, t=100

40、the rotated Roe scheme with the direction of upwinding determined by velocity-difference vector, t=600 4.2 Inviscid Supersonic Flow (M=20) around a Circular Cylinder The free stream Mach number is 20 , and this test problem is computed on a grid with 20 cells in the radial direction and 720 cells in

41、 the circumferential direction with the first order grid-aligned Roe scheme and the first order rotated Roe scheme. The first order grid-aligned Roe scheme The first order rotated Roe scheme with the direction of upwinding determined by the velocity-difference vector. 4.3 The Diffraction of a Supers

42、onic Shock Moving over a Corner The shock Mach number is 5.09. The computational domain is a unit square 0, 1 0, 1 that is discretized into a 400400 uniform cells. The corner is at (x, y)=(0.05, 0.625). Initially, the shock is at x=0.05. To the right of the shock, the flow field is initialized to To

43、 the left of the shock, the flow variables are computed using moving shock relations ( , , , )(1.4,0,0,1)u v pxy05xy05 The first order grid-aligned Roe scheme The first order rotated Roe scheme with the direction of upwind

44、ing determined by the velocity-difference vector. 4.4 Mach 3 Wind Tunnel with a Step The computational domain is 0, 3 0, 1. The corner of the step is located at (x,y)=(0.6,0.2). The initial conditions are( , , , )(1.4,3,0,1)u v pxyxy The second order grid-aligned Roe sche