The Riemann problem in special relativistic hydrodynamics

In this project we study the evolution of a relativistic fluid with an ultra-relativistic equation of state in slab symmetry. We will use one of the so called HRSC (High Resolution Shock Capturing) methods which are suited for the evolution of discontinuities. These methods have been proven to be useful in both the special relativistic case [1] and the general relativistic one [2]. The equations of motion can be derived from the following conservation laws:

$$\left(\rho_o u^\mu \right)_{;\mu} = 0,$$ (48)

$$T^{\mu \nu}{}_{;\mu} = 0,$$ (49)

where here $\mu, \nu$ take values $\{0,1,2,3\}$, $\rho_o$ is the proper rest mass density on a local inertial frame, $u^\mu$ stands for the four velocity of the fluid and $T^{\mu \nu}$ is the energy momentum tensor. Equation (49) represents the conservation of baryons in our system and (50) the conservation of the energy momentum.

For a perfect fluid (we will only consider an adiabatic fluid, i.e. we do not consider heat exchange or viscous terms) the stress energy tensor can be written as [3]: :

$$T^{\mu \nu} = \left( \rho +P \right) u^\mu u^\nu + P g^{\mu \nu},$$ (50)

where $\rho=\rho_o\left(\epsilon+1 \right)$ is the energy density of the fluid, $g^{\mu \nu}$ is the inverse of the metric of the space-time and $\epsilon$ is the specific internal energy. In our case we consider the fluid on Minkowski space-time, therefore $g^{\mu \nu}=\gamma^{\mu \nu}=diag\{-1,1,1,1\}$. In order to completely describe the fluid we need to support these equations of motion with and equation of state that relates the pressure with the rest of the variables that describe the fluid. the equation of state for the relativistic ideal fluid is

$$P = (\Gamma-1)\rho_o\epsilon$$ (51)

with $1 < \Gamma < 2$ as the adiabatic constant, in our case we used $1.5$, which is recommended so to avoid severe jump situations where the code will breakdown and the solution turns extremely sensitive to initial conditions⁴. We will apply our Roe type methods for this problem.

The special relativistic hydrodynamic (SRHD) equations are hyperbolic for causal equations of state⁵. In conservative form the $1D$ SRHD equations are:

$$\frac{\partial {\bf U}}{\partial t}+\frac{{\bf F}}{\partial x} =0$$ (53)

with the conservative variables

$${\bf U} = \begin{array}{cc} \left(\begin{array}{c} D\\ ... ... v \\ (\rho+P)W^2-P-D \\ \end{array}\right) \\ \end{array}$$ (54)

and the fluxes,

$${\bf F}=\left(\begin{array}{c} Dv\\ Sv+P \\ S-Dv\\ \end{array}\right)$$ (55)

The system is balanced by the equation of state (53). The transform from conserved variables to primitive variables $(\rho,v,P)$ is not straightforward and a non-linear equation must be solved numerically.

The first modification to the code for Burgers equation was just to take care now of 3 equations instead of one. Then we have to figure out the calculation of both primitive and conservative variables. The primitive variables are reconstructed at the boundaries of the cells using a slope limiter. However, in order to compute the physical fluxes , we need to compute the conservative variables as well using equations (55).

Conversely, at every half and full step in our update procedure we need to calculate the primitive variables after the conservative variables have been evolved. So we need to invert the equations that define the conservative variables in order to get the primitive ones.

The primitive variables are needed for the spectral decomposition (exact Riemann solver). To calculate the primitive variables from the conservative ones is provided in David Neilsen's thesis [10], page 85, as pointed out by Martin. This is basically a Newton iterative procedure:

$$\fbox{ \begin{minipage}{12cm} \begin{center} Conservative $... ...{new} = \omega - \frac{f(\omega)}{f'(\omega)}$ \end{minipage}}$$ (56)

The auxiliary variable $\mu$ helps calculating the derivative of the function $f$. The inequality had the wrong orientation for the while statement, which did not make sense at the beginning.

The spectral decomposition was followed from Ian Hawke's thesis [12], pag. 94-95 and feeded into spec_decom.f The characteristic structure of the system (eigenvalues and eigenvectors) are:

eigenvalues,

$$\lambda_+ = \frac{v+c}{1+vc}$$ (57)

$$\lambda_o = v$$ (58)

$$\nonumber \lambda_- = \frac{v-c}{1-vc}$$

The right and left eigenvectors are written as follows,

$$r = \left(\begin{array}{ccc} \left(\begin{array}{c} 1 \\ ... ...- \\ hWA_--1 \\ \end{array}\right) \\ \end{array}\right)$$ (59)

$$l = \left(\begin{array}{c} -\frac{h^2}{\Delta}\left(\begin... ...A_+\lambda_+-v \\ \end{array}\right) \\ \end{array}\right)$$ (60)

we have for the Jacobian matrix (see appendix on Godunov methods),

$$A_\pm = \frac{1-v^2}{1-v\lambda_\pm}$$ (61)

and $\Delta$, c y h are

$$\Delta = h^3W(h-1)(1-v^2)(\lambda_+A_+-\lambda_-A_-)$$ (62)

$$c = \sqrt{\frac{P\Gamma(\Gamma-1)}{(\Gamma-1)\rho_o + \Gamma P}}$$ (63)

$$h = 1+\frac{P\Gamma}{(\Gamma-1)\rho_o}$$ (64)

The SRHD system of equations (4.1) together with the initial conditions:

$$(P,\rho,v)= \bigg\{ \begin{array}{lr} (P_L,\rho_L,v_L)& if x<0\\ (P_L,\rho_L,v_L)& if x>0 \end{array}$$ (65)

constitute the special relativistic Riemann problem. As in the classical Riemann problem it contains three types of waves: shock, rarefraction and contact discontinuity waves. The solution is schematically represented in fig (8).

Figure 8: The $L,L^*,R^*,R$ represent constant state solutions, the dashed line is a contact discontinuity and the double lines represent a shock - or a rarefraction wave. The states L and R are known. also, plotted in the figure are the labels $a,b,$ to indicate a state (a)head or (b)ehind a shock/rarefraction. Across the contact discontinuity $P_L^* =P_R^*$ and $v_L^*=v_R^*$; only the density makes a jump across the contact discontinuity. A rarefraction is formed if $P_b\le p_A$, else a shock is formed.

So we have the ambiguity about which variables to reconstruct, and which to average for use in the spectral decomposition. We decided to implement the four combinations. We test the code using the models provided in the assignment, which we summarize in the next table (4.1) We will provide some comments about the overall procedure in each code and we refer the reader to the web site to observe the different evolutions we encountered. We continue using a reconstruction based on the slope limiter.

Model	A	B	C
$\rho_{oL}$	10	10	10
$\rho_{oR}$	15	15	15
$v_L$	0.6	0	-0.6
$v_R$	-0.6	0	0.6
$P_L$	50	50	50
$P_R$	10	10	10

We will provide some comments about the overall procedure in each code and we refer the reader to the web site to observe the different evolutions we encountered. We continue using a reconstruction based on the slope limiter. All reconstructions in each step are carried out in the local domain for each step and on each variable, of course.