Solving 2D Euler Euqations Using 2D MUSCL Scheme

1 Abstract

Euler equations have been used as a benchmark for testing the performances of various numercial schemes. Here, the 2D MUSCL scheme, short for Monotone Upstream-centered Schemes for Conservation Laws, will be utilized to solve 2D Euler equations. A 2D Riemann problem(6.1 ), a 2D Lax problem(6.2 ), and a 2D Shu-Osher problem(6.3 ) will be test examples.

2 Introduction to the finite volume method

The idea of the finite volume method is: if the domain is divided into small cells, the change in cell average of a unit of measure which obeys the law of conservation, say mass, equals the net flux change on all the boundaries of the cell, as described in the following equation:

             ∫
d¯ui-    --1-      ⃗ ⃗
dt  = − |Di|  ∂D F (U ) ⋅ ˆnds                            (1)
                i

where ui =  1
|Di|Diud⃗A is the cell average of u in Di, and ˆn is the unit normal vector of the boundary ∂Di.
Equation (1 ) can be applied to cells of any shape. However, since the initial conditions (6.1 )(6.2 )(6.3 ) are defined in rectangular domains, it is reasonable for one to choose rectangular cell partitions in order to exactly fit the domain boundaries. Hence equation (1 ) can be rewritten as

              ∫                  ∫                   ∫                    ∫
d¯ui-    --1-        ⃗ ⃗                ⃗  ⃗                ⃗  ⃗                ⃗  ⃗
 dt = − |Di |{      F (U ) ⋅ ˆneds +     F (U) ⋅ ˆnnds +      F (U) ⋅ ˆnwds +      F (U ) ⋅ ˆnsds} (2)
               ∂Die                ∂Din                 ∂Diw                 ∂Dis

where the subscripts e,n,w,s represent the east, north, west, and south boundary of a cell, respectively. In the Cartesian coordinate system, ˆne = [  ]
  1
  0 , ˆnn = [  ]
  0
  1 , nˆw = [   ]
 − 1
  0 , and ˆns = [   ]
  0
 − 1 .

3 2D Euler equations

2D Euler equations can be written in the following way

∂-⃗U-      ⃗  ⃗
 ∂t + ∇  ⋅F (U) = 0                                  (3)

where
                      ⌊    ⌋
                         ρ
                 ⃗    || ρu ||
                 U =  ⌈ ρv ⌉

                        ρE
                     [⃗ ⃗  ]
            ⃗F (⃗U) =   f(U )
                      ⃗g(⃗U )
                ⌊          ⌋
                     ρu
                |  ρu2 + p |
        f⃗(U⃗) = |⌈          |⌉
                     ρuv
                  u(E +  p)
                 ⌊         ⌋
                     ρu
        ⃗g(⃗U) =  ||   ρuv   ||
                 ⌈ ρv2 + p ⌉
                  v(E +  p)

                1-   2    2
p = (ν − 1)(E − 2 ρ(u  + v ))   ν =  1.4

There are five unknowns and five equations, so 2D Euler equations are closed. After applying the concept of equation (2 ) to equations above, we get:

               ∫ uby               ∫ ubx               ∫ uby               ∫  ubx
d¯ui-= − --1-- {     F⃗(⃗U ) ⋅ ˆnedy +     ⃗F (⃗U) ⋅ ˆnndx +      ⃗F (U⃗) ⋅ ˆnwdy +      ⃗F (U⃗) ⋅ ˆnsdx} (4)
 dt     dxdy    lby                 lbx                  lby                  lbx

For example, if we write the second equation of the Euler equations in the form of equation (4 ), we get

 

dρ¯u        1    ∫ uby              ∫ ubx         ∫  uby               ∫ ubx
---- = − -----{      (ρu2 + p)dy +      ρuvdx  −       (ρu2 + p)dy −      ρuvdx }   (5)
 dt      dxdy    lby                 lbx            lby                 lbx

and the equation for 1D conservative schemes (can be found on most tutorials)

d-¯u     -1--                                 ˜         ˜
 dt = − Δx {Fj+ 12 − Fj− 12}  where Fj+ 12 = F (U(xj+ 12− ),U (xj+12+))         (6)

is just a special case of equation (4 ).
Numerically solving the Euler equations requires a trick: solving each equation independently by using the variables solved at the previous time step. Take the following equation as an example. In theory, at tn one needs to know u, v, and p, which are solutions (solved at tn as well) from other equations in order to solve for ρ

d(ρu )        2
------=  − (ρu + p )x − (ρuv )y
  dt

However, in other equations, one needs to know ρ value (solved at tn) before u or v or p is solved. Numerically solving the coupling equations can be done byThis cou instead of doing so, we obtain their values (u, v, and p) from the preivous time step at tn−1. This way, the coupled Euler equations become four separate equations and we can easily solve for the four elements in ⃗
U by vectorizing Euler equations. In MATLAB, this will significantly accelerate the computation process.

4 MUSCL scheme

MUSCL scheme is a Riemann-free-solver scheme. In the 2D MUSCL scheme, one needs to reconstruct a 2D piece-wise polynomial for each cell. This polynomial is determined by the cell averages of the current cell and its neighbors. Therefore, the degree of the piece-wise polynomial, its covering range, and the cell configuration determine the complexity of the reconstruction of a solution. For the 2D MUSCL scheme, a linear polynomial is adopted, so it can be written in the form of P(x,y) = a0 + a1(x − xi) + a2(y − yj). The three coefficients a0,a1,a2 are determined by different stencil configurations (see Appendix(7 )). A MINMOD function (the slope limiter) is a reasonable choice to determine the coefficients. After the coefficients are determined, the reconstructed solution for the specific cell is determined.
When the reconstructed piece-wise solutions are available for each cell, one then uses the reconstructed solution on the boundaries to determine the flux for each cell. There are several types of flux functions to choose, and they must satisfy

  • Consistency,
  • Lipschitz continuity, and
  • Monotonicity.

Lax-Friedrichs flux is a good choice here.
Lax-Friedriches flux has the following form:

 ⃗ ⃗        1- ⃗ ˜ in       ⃗  ˜out        ˜a-˜ in   ˜ out
F (U ) ⋅ ˆn = 2{F (U ) ⋅ ˆn + F (U  ) ⋅ ˆn } + 2(U − U   )                (7)

where Ũ is the reconstructed solution at the boundaries and superscripts in or out represents the values in or out of the current cell, respectively.
ã is a factor which determines the flux. It determines the size of time-marching step. In order to keep this scheme stable (See equation (9 )), it should be the largest characteristic speed. The characteristic speed can be represented as follows:

                ∘ ---
    √ --2---2-      p
˜a =   u  + v  ±   γ ρ-                                (8)

For all the elements in ⃗U and all the cells in the domain, at moment tn, the time step Δt for marching to the next step should obey the following rule:

           dx-  dy-
Δt < min (c ˜a ,c ˜a )                                 (9)

where c is the CFL number.

5 Time Discretization

Since the MUSCL scheme is a second order scheme, using the improved Euler (predictor-corrector) method for discretizing time will be effective and efficient (using Runge-Kutta 3rd order is overkill for this problem). The improved Euler method contains three steps:

GivendUdt- = L (U)

  1. U∗n+1 = Un + ΔtL (Un)
  2. U∗n+2 = U∗n+1 + ΔtL (U∗n+1)
  3. Un+1 = 1
2(Un + U∗n+2)

where Δt satisfies equation (9 )

6 Results

The three test examples all have discontinuities in the initial conditions. Overcoming discontinuities (Riemann problems) in numerical schemes has been a challenging task since an exact Riemann solver requires Raphson-Newton iterations and the information of wave structure and wave speeds. On the other hand, the MUSCL scheme uses predictor-corrector strategy, which saves computation.
The MUSCL is a second order conservative scheme, so it doesn’t catch as many fine details as the high order scheme in [1]. However, this scheme is relatively computationally inexpensive (compared to the shemes in [1]).

6.1 2D Riemann problem

2D Riemann problem is defined as follows: the domain is [0, 1] × [0, 1]. The initial states are constant within each of the 4 quadrants. Counterclockwise from the upper right quadrant, they are initialized with (ρi,ui,vi,pi), i = 1,2,3,4. Initially, ρ1 = 1.1,u1 = 0,v1 = 0,p1 = 1.1; ρ2 = 0.5065,u2 = 0.8939,v2 = 0,p2 = 0.35; ρ3 = 1.1,u3 = 0.8939,v3 = 0.8939,p3 = 1.1; ρ4 = 0.5065,u4 = 0,v4 = 0.8939,p4 = 0.35.


Figure 1: Riemann problem: the initial condition of density


PIC

Figure 2: Riemann problem: the dentisty contours at Tn = 0.1


PIC

Figure 3: Riemann problem: the dentisty in 3D at Tn = 0.1


PIC

Figure 4: Riemann problem: the dentisty in 3D at Tn = 0.25, Δx = Δy = -1-
400, CFL# =0.4


PIC

Figure 5: Riemann problem: the dentisty contours at Tn = 0.25, Δx = Δy = -1-
400, CFL# =0.4


PIC

Figure 6: Riemann problem: the dentisty along x = 0.8 at Tn = 0.25, Δx = Δy = -1-
400, CFL# =0.4


6.2 2D Lax problem

2D Lax problem is defined as follows:

              {
                 (0.445, 0.311, 0,8.928), for 0 ≤ x < 0.5and  0 ≤ y ≤ 1
[ρ ρu ρv  E] =    (0.5,0,0, 1.4275 ), for 0.5 ≤ x ≤ 1and 0 ≤ y ≤ 1              (10)

PIC

Figure 7: Lax problem: the initial condition of density


PIC

Figure 8: Lax problem: the dentisty along x direction at Tn = 0.16, Δx = Δy = -1-
200, CFL# =0.5


6.3 2D Shu-Osher problem

2D Shu-Osher problem is defined as follows:

           {
              (3.857143, 2.629369, 0,10.333333), for − 5 ≤ x < − 4and 0 ≤ y ≤ 2.5
[ρ u v p] =    (1 + 0.2 sin(5x ),0,0,1), for − 4 ≤ x ≤ 5and 0 ≤ y ≤ 2.5               (11)

PIC

Figure 9: Shu-Osher problem: the initial condition of density


PIC

Figure 10: Shu-Osher problem: the dentisty along x direction at Tn = 1.8, Δx = Δy = -1-
160, CFL#= 0.9

7 Appendix

In order to determine the coefficients of the piece-wise polynomial for each cell, one has to consider four stencil configurations: Figure (11 ), (12 ), (13 ), (14 ), and

                   ∫  uby ∫ ubx
        U¯ =  --1--          {a  + a  (x − x ) + a (y − y )}dxdy
          i   dxdy   lby   lbx    0    1      i    2      j
                                2                     2
= --1--{a0x |ubxy|uby+ a1 (x-−--xi)-|ubxy|uby + a2(y-−-yj)-|ubyx|ubx}
  dxdy      lbx  lby         2     lbx   lby         2     lby  lbx


Figure 11: Stencil 1

PIC

(        -1--
{  U¯1 =  dxdy{a0dxdy +  0 + 0} = a0
   U¯2 =  d1xdy{a0dxdy +  a1dx2dy + 0} =  a0 + a1dx
(  U¯ =  -1--{a dxdy +  0 + a dy2dx} =  a + a  dy
     3   dxdy  0             2           0    2



Figure 12: Stencil 2

PIC

(        -1--
{  U¯1 =  dxdy{a0dxdy +  0 + 0} = a0
   U¯2 =  d1xdy{a0dxdy +  0 + a2dy2dx} =  a0 + a2dy
(  U¯ =  -1--{a dxdy −  a dx2dy + 0 } = a −  a dx
     3   dxdy  0         1               0    1



Figure 13: Stencil 3

PIC

(        -1--
{  U¯1 =  dxdy{a0dxdy +  0 + 0} = a0
   U¯2 =  d1xdy{a0dxdy −  a1dx2dy + 0 } = a0 − a1dx
(  U¯ =  -1--{a dxdy +  0 − a dy2dx } = a −  a dy
     3   dxdy  0             2           0    2



Figure 14: Stencil 4

PIC

(        -1--
{  U¯1 =  dxdy{a0dxdy +  0 + 0} = a0
   U¯2 =  d1xdy{a0dxdy +  0 − a2dy2dx + 0} =  a0 − a2dy
(  U¯ =  -1--{a dxdy +  0 + a dx2dy} =  a + a  dx
     3   dxdy  0             1           0    1


Obviously, for each cell, there are two possible values for a1, two possible values for a2, and only one possible value for a0. After all the possible values are available, we use MINMOD to determine one set of {a0,a1,a2}. MINMOD is used to maintain the TVD (total variation diminishing) property and is defined as follows:

                           (
                           { min {x1,x2,...,xn}  if xi > 0, for all i
MINMOD    {x1,x2,...,xn } = ( max {x1,x2,...,xn}  if xi < 0, for all i       (12)
                                     0              otherwise
It is clear that using linear piece-wise polynomials in rectangular cell configuration could save a lot of computation because solving a1 and a2 doesn’t require any Gaussian elimination.

References

[1]   Y. Liu, C.W. Shu, E. Tadmor, M. Zhang, et al. Non-oscillatory hierarchical reconstruction for central and finite volume schemes. Comm. Comput. Phys, 2:933–963, 2007.