CFD Online Discussion Forums - Sod Shock Tube problem

Hi all,

I'm trying to model the Sod Shock Tube Riemann Problem using a Riemann decomposition method. The theory is outlined in EF Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, section 3.1.

The problem consists of an area of density, pressure of 1 and initial velocity of 0.75 on the left of a membrane, and an area with density of 0.125, and pressure of 0.1 with a speed of 0 on the right.

At time zero the membrane is removed.

I'm using a Riemann method to solve the problem, first for the Riemann variables, which is the matrix of right eigenvectors (of the Jacobian matrix of the flux) premultiplied by the variables (rho, rho.u, E)^T where the right eigenvectors are

| 1 , 1 , 1 |

| u-a , u , u+a |

| H -ua , 0.5u^2 , H+ua |

a = soundspeed, H = (E + p)/rho, p = pressure, E = total energy = rho.(0.5u^2 + specific internal energy).

The eigenvalues of this system are (u-a, u, u+a).

For the next timestep.

If the eigenvalue is greater than zero, we take INITOUTPUT(i,k) = INITOUTPUT(i,k)*(1 - EIG(i,k)*(DT/DX)) + EIG(i,k)*(DT/DX)*(INITOUTPUT(i-1,k).

If the eigenvalue is less than zero, we take

INITOUTPUT(i,k) = INITOUTPUT(i,k)*(1 + (DT/DX)*(EIG(i,k)) - EIG(i,k)*(DT/DX)*INITOUTPUT(i+1,k).

Where INITOUTPUT(i,k) are the Riemann variables, i is the x variable.

Then we use a newton solver to solve the non-linear relationship between the Riemann variables and the original variables, (rho, rho.u, E).

Regardless of the timestep or grid spacing, I seem to get the density decreasing smoothly towards the discontinuity over time, then sharply reducing to much lower than the original 0.125 on the right hand side just after the discontinuity and only then going back to the original 0.125. The solution to this problem is well known and not this.

I'm using reflective boundary conditions for the near wall rho(0) = rho(1), u(0) = -u(1), E(0)=E(1) and the same on the far wall.

Can anyone help? I would really appreciate it.