高等计算流体力学-06

1、第六讲Riemann问题的近似解1十八、Riemann问题的近似解Approx. Riemann Solver2为什么要研究Riemann问题的近似解? 精确求解Riemann问题计算量大；某些双曲型守恒律的Riemann问题无精确解！ 由于在时间积分、通量积分、重构等步骤已经引入了各种近似，精确求解Riemann问题并不会导致整体精度的提高（但对一维问题而言，基于Riemann问题精确解的Godunov格式的稳健性和精度的平衡可能较好。）3Riemann 问题4扩展一维系统扩展一维系统初始条件初始条件 的解对应有限体积格的解对应有限体积格式的数值通量式的数值通量Riemann问题的近似解法1

2、：HLL方法 Harten-Lax-Van Leer， Einfeldt（1988） 基本思路：双波近似（不考虑接触间断）5tRxLxHLLULURULdxSdtRdxSdt为了避免奇异性，不对左右波的类型作假定。左右波只表示间断初始值的影为了避免奇异性，不对左右波的类型作假定。左右波只表示间断初始值的影响范围；响范围；优点：简单；优点：简单；已知问题：对接触间断分辨率低。已知问题：对接触间断分辨率低。假定SL，SR已知积分关系 取 积分区域 方程6相容条件相容条件7 把前述积分区域分为 左右两部分，分别积分8相容条件相容条件HLL Riemann Solver9一般情况一般情况通量计算时考虑

3、通量计算时考虑x/t=0波速的确定10Roe平均平均不推荐！不推荐！推荐！推荐！还有其他方法，见还有其他方法，见Toro的书。的书。HLL的变形11Rusanov（Lax）LaxRiemann问题的近似解法2：HLLC方法 Toro 基本思路：三波近似（在HLL基础上加接触间断）1213相容条件相容条件14超定还是欠定？超定还是欠定？假定假定SL，SR已知，关键问题是求已知，关键问题是求S*利用上二式的第一、二两利用上二式的第一、二两项项中间波是接触间断中间波是接触间断接触间断接触间断满足相容条件满足相容条件15*22*()()()()()()LLLLLLLLLLLLSuSuSuupSuup*

4、22*()()()() ()() LLLLLLLLLLLLLuuSSppSuuuu22* ()() ()()()()()()()()()()()()()LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLppS SuupSuSuSuSupSuSuSuSupSuSuSupSuuu16变形变形117变形变形2kS18变形变形3假定假定*状态有单一压力状态有单一压力 ，不不直接用直接用Riemann问题的近似解法3：Roe方法19考虑扩张一维守恒型考虑扩张一维守恒型Euler 方程方程 ：0UFtx 对应的拟线性形式为：对应的拟线性形式为：0UUAtx Roe 把

5、上式中的矩阵把上式中的矩阵A 用一个常矩阵用一个常矩阵 代替代替 ，即，即 (,)LRAA U U,:LRU URiemann问题左右状态0UUAtx则则0UFtx守恒形式守恒形式FAU基本思路基本思路:在控制体界面处把方程线性化在控制体界面处把方程线性化通过求解线性通过求解线性Riemann问题计算通量问题计算通量00( ,0)0txLRUAUUifxU xUifx线性线性Riemann问题问题关键问题关键问题: 如何线性化如何线性化?称为称为Roe Jacobian Matrix, Roe要求它满足下列条件：要求它满足下列条件：( ,)A U UA()()()RLRLF UF UA UU（

6、A） 双曲性双曲性 即即 的特征值均为实数，且存在完备的左右特征向量；的特征值均为实数，且存在完备的左右特征向量；（B） 相容性相容性 （C） 守恒性守恒性 如果左右状态满足激波关系，则如果左右状态满足激波关系，则 ()()()()RLRLRLF UF UD UUA UU D是是 的特征值，此时线性的特征值，此时线性Riemann问题的解与问题的解与Riemann问题精确解相同！问题精确解相同！ 根据这些条件可根据这些条件可确定确定 （先假定（先假定为为已知）已知）(,)( )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 当x/t 当x/t22为了计算通量，只需为了计算通量，只需在在

8、x=0处处的解：的解：( )( )( )( )( )( , )jjLjjjRlUvx tlU 当x/t 当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线化方程线化方程的数值通量的数值通量11/2011()22iLRRLFFAUAULL UU如何计算原始如何计算原始Euler方程的数值通量？方程的数值通量？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方程（无线性化，但假定与线化方程有相同的解的结构）方程（无线性化，但假定与线化方程有相同的解的结构）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编程实施方案：编程实施方案：()1-miRLmlUUUr(1)(2)(1)(2)()()(1)(2)(1)(2)()()1(),mmmimmlmllURL UrrrUllUlUrrrrlU( ).iilU( ) ( )( )1( )( )mjjjrl 事实上，对于任事实上，对于任意的意的m维列向量维列向量11( )( )1( )( )1()()RLmjjRLjm

12、jjjjjLL UUrlLL UUr ( )1/2( )11122mjiLRjjjFFFr 27( )1/2( )11122mjiLRjjjFFFr ( )( ),jjjr只需计算：Roe-Pike方法方法: 无需计算无需计算 ,直接计算上述三个量。直接计算上述三个量。考虑三维问题对应的扩展一维问题，有：考虑三维问题对应的扩展一维问题，有：A右特征向量右特征向量（注意记号变化）（注意记号变化）(,)()(,)LRLRA UUA UUUU28利用：利用：29关键问题：如何确定关键问题：如何确定Roe-Pike方法：求平均状方法：求平均状态态，使下列关系成立：使下列关系成立：(,)LRUU U其中

13、：其中：11( )( )1( )( )1()()RLmjjRLjmjjjjjLL UUrlLL UUr 11( )( )1( )( )1()()()RLRLmjjRLjmjjjjA UULL UUrlLL UUr ()()()RLRLF UF UA UU3031,ii j iAii当当中所有的对角元素中所有的对角元素都充分大于零时，格式粘性项都充分大于零时，格式粘性项是正定的，从而保证了计算格式有足够的数值耗散，是正定的，从而保证了计算格式有足够的数值耗散，变为零。此时，格式粘性项的系数矩阵是半正定的，变为零。此时，格式粘性项的系数矩阵是半正定的，某些特征波的数值粘性过某些特征波的数值粘性过小

14、，从而有可能得到非物理解。所以，一般情况下，非物理解出现在某些小，从而有可能得到非物理解。所以，一般情况下，非物理解出现在某些非常接近零的膨胀流动区域，表现为类似于非常接近零的膨胀流动区域，表现为类似于“膨胀激波膨胀激波”的系数矩阵的系数矩阵可以得到物理解。但是，在某些流动区域，如音速点附近，一些可以得到物理解。但是，在某些流动区域，如音速点附近，一些（过激波熵减少）的数值解的跳跃。（过激波熵减少）的数值解的跳跃。 1/21/21/21/21/211 ()()()22LRRLjjjjjAFF UF UUU,ii jARL 熵修正熵修正32进行限制，使其不致过小。例如，可以指定一个小正数进行限制

15、，使其不致过小。例如，可以指定一个小正数“熵修正（熵修正（entropy fixentropy fix）”是避免非物理解出现的一个有效手段。一种简单的熵修是避免非物理解出现的一个有效手段。一种简单的熵修正方法正方法是对是对ii*222iiiiifotherwise*ii使限制后的使限制后的*i（记为（记为）为）为用用代替代替进行数值通量的计算，对于避免进行数值通量的计算，对于避免FVSFVS格式和格式和RoeRoe格式出现格式出现非物理解都是有效的。非物理解都是有效的。3334评述及扩展 基于Riemann问题精确或者近似解计算数值通量的方法称为Godunov类型格式或者通量差分裂（Flux

16、Difference Splitting，FDS）格式。 优点： Godunov、Roe、HLLC对激波和接触间断都有很好的分辨能力；对粘性问题的边界层和剪切层计算效果较好。 缺点：有可能出现奇偶失连或者Carbuncle现象 35 其他迎风型的通量计算方法还包括： Flux Vector Splitting：Steger-Warming， Van Leer，Lax（Rusanov） 不出现奇偶失连或者Carbuncle现象，可计算强激波 对接触间断、边界层和剪切层的分辨率低 压力-对流分别分裂的方法：AUSM类（Liou等），CUSP（Jameson）， LDFSS（Edwards）等等，这

17、些方法均试图结合FDS和FVS的某些优点而克服其缺点。36 解决FDS方法奇偶失连或者Carbuncle现象的一种有效方法是采用旋转Riemann求解器 任玉新（2003）， Nishikawa等（2008），Wu（2010） Godunov、Lax、HLL等方法满足离散的熵条件，但分辨率较低。 没有一种在任何情况下都最优的方法！3738 FVS方法方法 （流通矢量分裂（流通矢量分裂 逐点分裂）逐点分裂） fff 具体方法：具体方法： Steger-Warming 分裂分裂 Lax-Friedrichs分裂分裂 Van Leer分裂：分裂：2kkkwcucuucucuu232221321321

18、)(2)(2) 1()()() 1(2) 1(22)(f2/ )(*Uff根据当地根据当地Mach数分裂数分裂保证保证 的的Jocabian阵特征值为正，阵特征值为正， 的为负的为负ffUAf)2/12/12/1RjLjnj(Uf(Uff一个参数，反映全部特征一个参数，反映全部特征39AUSM:Liou-Steffen分裂分裂)()(2200)()(pcFFpupEuuupEpuuf(U)对流项压力项思路：思路： 决定特征的关键参数决定特征的关键参数 当地当地Mach数数1 1 , 00 , 11cuMa超音速，超音速，x-方向方向超音速，超音速，x+方向方向0, 0, 0321cucuu32

19、1,0, 0, 0321因此，对因此，对Mach数进行分裂更为简洁！数进行分裂更为简洁！1当01当4/) 1(1当2MMMMMMaHauaMFc)(114/) 1(102MMMMMM1012/ )1 (1MMMpMpp112/ )1 (10MpMMpMp显然：显然： pppMMMfff010paHauaMf参考文献：参考文献：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分裂，但压力单独处理分裂，但压力单独处理MM保证光滑过渡保证光滑过渡M=1Journal of Computational Physics 208 (2005) 527569, Kim&Kim4041三维问题的旋转不变性三维问题的旋转不变性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