wiki:u/erica/ExactRiemannSolver

Version 35 (modified by Erica Kaminski, 12 years ago) ( diff )

Intro

The Riemann problem is an IVP for the Euler Equations, which consists of 2 constant initial data states separated by a discontinuity between them, say at x=0. For x<x=0, we say the data state is XL, for which there are initial variables WL=(rho_L, p_L, u_L). Similarly for the initial right data state, XR=x>x=0, the initial variables are WR=(rho_R, p_R, u_R). The solution of this IVP consists of 3 nonlinear waves, a left wave that is either a shock or a rarefaction, a center contact discontinuity, and a right wave that is either a shock or a rarefaction wave. Depending on which type of L- or R-wave is present, different expressions exist that describe the change in variables across them. The contact discontinuity is special in that the pressure (p) and velocity (u) are constant across it.

Program Outline for Exact Riemann Solver

Given the above properties of the Riemann problem, an algebraic expression can be derived which gives p in the central "star-region" (denoted by '*'). The overall structure of the Riemann problem then is to solve this algebraic equation for p*. Once p* is known, u* follows immediately. The remaining rho*_L and rho*_R follow from expressions valid for the specific L- or R- wave present. Specifically, this means that the code will determine at every point (x,t), that point's relative position to the different waves present. Once the position is determined, the solution is given by analytic expressions for the following five possibilities: pre- or post- shock, ahead of rarefaction head, behind rarefaction tail, or within the rarefaction fan. The position of each sampling point (x,t) is determined by a characteristic speed in the grid given by s = dx/t, where dx is the distance from the initial discontinuity to the sampling point (x,t), and t is the simulation end time. This s for every point on the grid can be compared to the known present waves, as their speeds are known exactly. In this way, the relative position of the sampling point with respect to the waves on the grid is determined, and hence the fluid variables are given at that point.

Exact Riemann Solver

Attached here is the Exact RS program I wrote. It reproduces the results for the Toro tests, as can be seen in the next section. The program is broken up into the following routines:

  1. Estimate an initial starting P*. This subroutine was taken from Toro almost exactly, as I have not yet covered the details for approximation methods in my studies. Once this P* is found, it is fed into the FindPstar routine, whereby a Newton Rhapson iteration method corrects this P* until it converges.
  1. Using this p*, u* is found in the star region given a simple analytical expression.
  1. Using p* and u*, the solution is sampled at every point on the grid. The position is determined by comparing speeds in the grid. First it is determined whether the point is on the left or right side of the contact discontinuity. (The speed of the contact disc. is given by u*). If s<u*, the point is on left of the contact. If s>u*, it is on right.
  1. Next, given the initial L- and R- pressure compared to P*, the type of wave is determined. If p_L < p* - left wave is shock, else it is a rarefaction. If p_R < p*, wave is shock, else it is a rarefaction. Note, this assumes waves always travel away from the contact.
  1. Once the waves are identified, then the different speeds in the grid are compared to the speeds of the different left and right waves. The appropriate analytic expression is then chosen, based on relative position to wave.

Results

These tests were taken from Toro, chapter 4. They are on a domain of length = 1 computational unit, with 1000 zones, and a gamma = 1.4. The final time in each of them varies, and were chosen to match those in Toro, as were the initial values for the left and right data states. These results match Toro's exactly, the slight mismatch in solutions at the rarefaction-tail, contact-discontinuity boundary was fixed, but the plots haven't been updated to show this change.

Note, the region of the contact discontinuity is identified by constant p and u between the two flanking outgoing waves.

Questions

Is the pressure condition for determining types of waves sufficient? How would the algorithm handle the colliding shocks problem (Test 5)?

Attachments (5)

Download all attachments as: .zip

Note: See TracWiki for help on using the wiki.