Recap

The Euler equations are a set of non-linear PDE's, which comprise an eigenvalue problem. The eigenvalues for the equations are functions of the solution to the equations themselves. This means that the waves, which propagate with speeds = to their eigenvalues, distort the solution, and the solution distorts them over space and time. Thus the solution to the Riemann problem for the non-linear system does not consist of a closed form expression for the values of pstar and ustar like it does for a linear system of equations. To solve the Euler equations exactly then, we have developed the method of characteristics that describe the propagation of waves outside of the star region. To solve for values of the q-array inside of the star region, we used an iterative scheme and then sampled the solution in the different wave regions. We now are concerned with approximations to this exact solution. The method to be discussed here considers a 'linearized' version of the Euler equations, so analytical methods used for linear, constant the averages (tilde expressions in Toro, eqn. 11.118) can be used.

The ROE Solver

The Roe solver is an approximation means for the numerical flux of the Godunov method, which is derived through linearizing a hyperbolic system of equations. For instance, the Euler equations in conservative form are written

\vec{U}_t + \vec{F}(\vec{U})_x = \vec{0}

which using the chain rule is identical to

\vec{U}_t + \frac{\partial \vec{F}} {\partial \vec{U}} ~\vec{U}_x = \vec{0}

If we assume that

\hat{A} = \frac{\partial \vec{F}}{\partial \vec{U}} where, A-hat is a Jacobian matrix of averaged/constant values, we can derive an expression for the numerical flux in terms of 1) wave strengths (alpha), and the 2) eigenvalues (lambda) and 3) right eigenvectors (K) of the 'averaged' Jacobian:

\vec{F}_{i+1/2}= \frac{1}{2}(\vec{F}_L + \vec{F}_R) - \frac{1}{2}\Sigma_{i=1,m} \tilde{\alpha}_i \tilde{\lambda}_i \tilde{\vec{K^{(i)}}} where

\vec{F}_k = \vec{F}(\vec{U}_k) So the goal is to compute the wave speeds and associated eigenvalues and eigenvectors of the Jacobian matrix. There are 2 methods by which we can do this: 1) The 'Roe' approach, which (non-trivially) constructs an averaged Jacobian directly that must satisfy conditions such as hyperbolicity and conservation, and 2), the newer 'Roe-Pike' approach, which avoids solving for the Jacobian and instead develops algebraic expressions for the sought quantities based on averages of the initial data. It is the 2nd, more widely used approach, that we will explore here.

The Roe-Pike Approach

This approach avoids construction of an averaged Jacobian directly. Instead the method involves 1) assuming the Jacobian consists of reference states for which the left and right states are close to, 2) solving for the eigenvalues and vectors for this 'reference' Jacobian, and the wave speeds, 3) evaluating these expressions using averaged versions of the original scalar quantities. In other words, the Roe Pike method consists of first, deriving expressions for the various quantities needed in the flux, then second, solving for the averaged primitive variables, and third, evaluating the expressions at the reference state equal to the average state.

From the equation above of the Euler equations in Jacobian form, we can derive the equations for wave speeds:

\triangle \vec{U} = \vec{U}_R - \vec{U}_L = \Sigma _{i=1,m} \hat{\alpha}_i(\vec{W}) \hat{\vec{K}}^{(i)}(\vec{W})

where the W-vector denotes the vector primitive variables rho, u, and p, as usual. The eigenvalues and vectors are derived straightforwardly from the Jacobian, using typical methods.

Now, we assume the independent variables u, rho, etc. of these functions for alpha, lambda, and K are reference states (associated with hat notation) which we use to solve for any general, averaged versions of the variables (associated with tilde notation), which we use in the numerical flux. Toro goes through the algebraic analysis for the general case for the Euler equations. The results are as follows.

Roe-Pike Approach for Euler Equations

For the x-split, 3-dimensional Euler equations, we have the following eigenvalues and vectors:

\lambda_1 = u-a, ~\lambda_2 = \lambda_3 = \lambda_4 = u, ~\lambda_5 = u+a and

\vec{K}^1 = <1, u-a, v, w, H-ua>, \vec{K}^2 = <1,u,v,w,V^2/2>, \vec{K}^3 = <0,0,1,0,v>, {\vec{K}^4 = <0,0,0,1,w>, {\vec{K}^5 = <1, u+a, v, w, H+ua>

where u, v, w, are velocities in x, y, z directions, respectively, a is the local sound speed, given by

a^2 = \frac{\gamma P}{\rho}

H is the enthalpy given by

H = \frac{E+p}{\rho} E is the total energy per unit volume,

E = \frac{1}{2}\rho \vec{V}^2 + \rho e with \vec{V}^2 = u^2 + v^2 + w^2 The Euler equations are a set of non-linear PDE's, which comprise an eigenvalue problem. The eigenvalues for the equations are functions of the solution to the equations themselves. This means that the waves, which propagate with speeds = to their eigenvalues, distort the solution, and the solution distorts them over space and time. Thus the solution to the Riemann

and e being the specific internal energy, which for ideal gases is,

e = \frac{p}{(\gamma - 1)\rho} Recall that the numerical flux is given by,

\vec{F}_{i+1/2}= \frac{1}{2}(\vec{F}_L + \vec{F}_R) - \frac{1}{2}\Sigma_{i=1,m} \tilde{\alpha}_i \tilde{\lambda}_i \tilde{\vec{K^{(i)}}} where the tilde's denote derived expressions evaluated at the average primitive variables:

\tilde{\alpha}_i = \hat{\alpha}_i(\tilde{\vec{W}}), ~\tilde{\lambda}_i = \lambda_i(\tilde{\vec{W}}), ~\tilde{\vec{K}^{(i)}} = \vec{K}^{(i)}(\tilde{\vec{W}})

Note that the vector equation for the numerical flux makes up 3 equations for the 3D Euler equations, one for each component of the following vectors:

the vector of numerical fluxes in conserved form, \vec{F}_{i+1/2} = <\rho_{i+1/2}, ~p_{i+1/2}, ~E_{i+1/2}> vector of left-fluxes in conserved form, \vec{F}_L = <\rho_L, ~\rho_L u_L, ~0.5\rho_L u_L^2 + \frac{p_L}{\gamma - 1}>

vector of right-fluxes in conserved form, \vec{F}_R = <\rho_R, ~\rho_R u_R, ~0.5\rho_R u_R^2 + \frac{p_R}{\gamma - 1}> and eigenvectors,

\vec{K}^{(1)} = <1, ~u-a, ~H-ua>

\vec{K}^{(2)} = <1, ~u, ~0.5 u^2>

\vec{K}^{(3)} = <1, ~u+a, ~H+ua>

Finally, the derived expressions for the wave-speeds are as follows,

\hat{\alpha}_1 = \frac{\triangle p - \hat{\rho} \hat{a} \triangle u}{2 \hat{a} ^2}

\hat{\alpha}_2 = \triangle \rho - \triangle p / \hat{a}^2

\hat{\alpha}_3 = \hat{\rho} \triangle v

\hat{\alpha}_4 = \hat{\rho} \triangle w

\hat{\alpha}_5 = \frac{\triangle p + \hat{\rho} \hat{a} \triangle u}{2 \hat{a} ^2}

and the corresponding expressions for rho, u, v, w, H, and a that we would plug into these wave speeds and the above expressions for the eigenvalues and eigenvectors are listed in Toro, eqn. 11.118. These in turn would be used in the numerical flux for the Godunov method.

Note in the above formalism, the following definition is used,

\triangle r = r_R - r_L

A Sample Algorithm

1. Compute the averages (tilde expressions in Toro, eqn. 11.118) for the various quantities
2. Compute the eigenvalues, eigenvectors, and wave speeds as functions of these averages
3. Compute the numerical flux given in above formula

Short-comings of the Roe Solver

As with all methods of linearizing the Euler equations, discontinuities are resolved, but continuous fan-like structures are not. Thus the Roe solver has trouble solving the RP inside of a sonic (aka transonic) rarefaction wave (see below figure), but less of a problem with contacts and shocks. This problem is refered to as the 'entropy problem', because the method treats the rarefaction wave as a 'rarefaction shock', which is entropy violating. Recall that the entropy condition that must be satisfied for a discontinuous wave jump is that the characteristics run into each other. That is, Sb>S>Sa, where Sb is the speed BEHIND the wave, S is the speed of the wave, and Sa is the speed ahead of the wave. There was a later 'entropy fix' developed for the Roe solver that corrected this issue.

Here is a diagram of the different types of waves:

The waves are as follows: (a) left supersonic shock, (b) right supersonic shock, © left supersonic rarefaction, (d) right supersonic rarefaction, (e) transonic or sonic rarefaction.

Here is an example of the Roe solver poorly computing the flow inside of a left sonic rarefaction. As you can see, it computes the solution as if there was a discontinuity inside of the wave, not a continuous fan. Note the sonic point inside of the wave is where the eigenvalue for the left wave, Lambda1 = u-a = 0. Throughout the wave, the eigenvalue should change sign through this sonic point.

The code - Entropy Fix

My Roe solver + entropy fix essentially does the following: 1) sets up left and right local Riemann problems as usual, 2) for each intercell boundary determines if the contact is left or right traveling, 3) based on the position of the contact, chooses which wave would cross the intercell boundary (the left wave associated with the (u-a) Euler eigenvalue for a right traveling contact, and vice versa), 4) decides if this wave is a sonic rarefaction, and if so, applies the Harten- Hyman entropy fix (described below), 5) if the wave is not a sonic rarefaction it determines if it crossed the intercell boundary, and if so, calculates the Roe flux, 6) if the wave has not crossed the intercell boundary, then it can be thought of as not having disturbed the adjacent cell's state, and so we just take the numerical flux to be the physical flux of that adjacent cell.

As stated above, the Rho flux is given by:

\vec{F}_{i+1/2} = \vec{F}_L + \lambda*alpha*\vec{K}

Where \vec{K}_{i+1/2} is the vector of numerical fluxes, \vec{F}_L is the vector of physical fluxes for the left state, lambda and vec{K are the eigenvalues and eigenvectors, and alpha is the wave speed, all with corresponding expressions given in detail above.

For the case of sonic rarefactions, the lambda in the above equation for Roe's flux must be adjusted, and following the Harten-Hyman approach, is taken to be:

\lambda ' = (u_L-a_L)*\frac{ustar-astarL-\lambda}{ustar-astarL-u_L+a_L}

for a left sonic rarefaction. For a right sonic rarefaction, the coefficients are slightly adjusted. Note that the entropy fix requires estimates for the velocity and sound speed of the star region, ustar and astarL. We get estimates of these using the Roe approach for linearizing the Euler equations. They are listed in Toro, section 11.4.3.

Results, Entropy Fix Code

Test 1 - Sod Test

Test 2

Fails with the Roe Solver.

Test 3

Test 4

Test 5

Discussion