Waves

The solution to the Riemann problem is built out of waves. These are regions where the variables change, separating constant states. There are two basic types of wave: those across which the variables change continuously (which, for fluid problems, are known as rarefactions), and those across which the variables change discontinuously (which break down into more subclasses).

We will assume that the state is known on one side of the wave, which will be denoted \({\bf U}_k\). The state on the other side will be unknown. Given the value of one variable, typically the pressure \(p\), on the unknown side, the solution across the wave can be found.

The equations here are adapted from the Living Review of Martí and Müller. The approach to reactive waves follows the paper of Zhang and Zheng, extended to the Relativistic case.

Wave sections

In general any wave can be built from multiple sections, pieced together. A wave may “split” into two discontinuous pieces moving at different speeds. A wave may contain a discontinuous section attached to a continuous section. In the most complex cases, the wave may have an arbitrarily large number of sections.

In the inert case, and when the equation of state is convex, each wave contains a single section. If the pressure decreases across the wave, \(p_k > p\), then the wave is a continuous rarefaction. If the pressure increases across the wave, \(p_k < p\), then the wave is a discontinuous shock. If the pressure does not change, the wave is either trivial (the state does not change) or a contact discontinuity (where quantities except \(p\) and \(v_x\) may jump).

In the reactive case, all the reactions that we can model will take place “instantly”. This means that reactions must take place across discontinuous wave sections. The reactive waves may now contain multiple wave sections. If the equation of state is convex, and the reaction exothermic, then across a reactive discontinuous wave the pressure increases. Therefore if the pressure is to decrease across the wave as a whole, there must be an inert precursor shock before the reactive discontinuity.

Finally, we note that not all reactive discontinuous waves are stable. This can be checked by looking at the wave speed on the unknown side of the wave. If the wave speed leads to characteristics coming out of the discontinuity, then the wave is unstable. In this case we need to add another wave section. We find the fastest stable reactive discontinuity (a Chapman-Jouget discontinuity), where the wave speed matches the speed of the discontinuity, and attach to it an inert rarefaction.

Rarefactions

Across a rarefaction the normal velocity satisfies

\[\frac{\text{d} v_x}{\text{d} p} = \pm \frac{1}{\rho_0 h W^2 c_s} \frac{1}{\sqrt{1 + g \left( \xi_{\pm}, v_x, v_t \right)}}.\]

The sign corresponds to the wavenumber: plus for the right acoustic wave and minus for the left. In the code this corresponds to lr_sign = wavenumber - 1.

The function \(g\) quantifies the effect of the tangential velocity, and is

\[g = \frac{\left( v_t \right)^2 \left( \xi_{\pm} - 1 \right)}{\left( 1 - \xi_{\pm} v_x \right)^2}.\]

The wavespeed itself is \(\xi_{\pm}\).

We also solve for the rest mass density and specific internal energy across the wave using

\[\begin{split}\frac{\text{d}}{\text{d} p} \begin{pmatrix} \rho_0 \\ \epsilon \end{pmatrix} = \begin{pmatrix} \frac{1}{h c_s^2} \\ \frac{p}{\rho_0^2 h c_s^2} \end{pmatrix}.\end{split}\]

Having solved across the wave, the tangential velocity follows as

\[v_t = h_k W_k \left( v_t \right)_k \frac{1 - v_x^2}{h^2 + \left( h_k W_k \left( v_t \right)_k \right)^2}.\]

The relation for the tangential velocity holds across both continuous and discontinuous waves.

Shocks

Across an inert shock we first find the post shock density and enthalpy by solving the nonlinear equation

\[h^2 - h_k^2 - \left( \frac{h}{\rho_0} + \frac{h_k}{\left( \rho_0 \right)_k} \right) \left( p - p_k \right) = 0.\]

This assumes the post shock pressure \(p_k\) is known, and that the equation of state gives the enthalpy \(h_k = h \left( \left( \rho_0 \right)_k, p_k \right)\).

From this we can compute the mass flux \(j\) across the shock, as

\[j = \sqrt{ \frac{p_k - p}{ \frac{h^2 - h_k^2}{p - p_k} - \frac{2 h_k}{p_k} } }.\]

This gives the shock velocity as

\[v_S = \pm \frac{ \left( \rho_0 \right)_k^2 W_k^2 \left( v_x \right)_k \pm j^2 \sqrt{ 1 + \frac{\left( \rho_0 \right)_k^2 W_k^2 \left( 1 - \left( v_x \right)^2 \right)}{j^2}}}{\left( \rho_0 \right)_k^2 W_k^2 + j^2}.\]

Again, the sign corresponds to the wave number.

Given the shock velocity we can compute the shock “Lorentz factor” \(W_S = (1 - v_S^2)^{-1/2}\), from which the post shock normal velocity is

\[v_x = \left( h_k W_k v_S \pm \frac{p - p_k}{j \sqrt{1 - v_S^2}} \right) \left[ h_k W_k + \left( p - p_k \right) \left( \frac{1}{\left( \rho_0 \right)_k W_k} \pm \frac{\left( v_x \right)_k W_S}{j} \right) \right]^{-1}.\]

The tangential velocity follows as in the rarefaction case.

Detonation

A detonation is a discontinuous reactive wave section across which the pressure increases. The equations are identical to those for the shock case, but the interpretation changes. All “known” (pre shock) variables have the reactive equation of state. All “unknown” (post shock) variables use the inert equation of state - the reaction has taken place across the discontinuity.

As noted above, it is possible for the resulting detonation wave section to be unstable. In general, detonations fall into two classes: weak detonations (which are unstable) and strong detonations (which are stable). For the stable strong detonations the characteristic waves impinge on the discontinuity from both sides. For the unstable weak detonations the characteristics only enter the discontinuity on one side.

If an unstable weak detonation is found, the section should be replaced by the fastest detonation that is stable. This Chapman-Jouget, or CJ detonation, is where the characteristic wave is parallel to the discontinuity. In this case, the post detonation pressure will not match the required post wave pressure, and an additional rarefaction wave section is needed to complete the wave.

Deflagration

A deflagration is a discontinuous wave section across which the pressure decreases. As a reaction cannot take place across a rarefaction wave, this requires a discontinuity. However, this discontinuity will reduce the temperature along with the pressure. This means that, unless the material was already at the right temperature to react, any reaction across this wave would be unphysical.

The solution is to start with an inert precursor shock which raises the temperature of the material to the ignition temperature. This first section follows the exact equations in the shock section above. Next, there is a deflagration wave section, across which the reaction takes place and the pressure drops. Again, this follows the shock equations, but with the same interpretation as in the detonation case.

Again, the deflagration wave section need not be stable. As with detonations, deflagrations fall into two classes: weak deflagrations (which are stable), and strong deflagrations (which are unstable). For the stable weak deflagrations the characteristic waves from one side impinge on the discontinuity, but not the other. For the unstable strong deflagrations, neither set of characteristic waves impinges on the discontinuity.

Again, if an unstable strong deflagration is found it is replaced with a CJ deflagration, and the full wave is completed with a rarefaction wave section.

Code

The Wave class in r3d2 in intended for internal use, but can be used to construct single waves directly. We need a known state with an equation of state first:

from r3d2 import eos_defns, State, wave

eos = eos_defns.eos_gamma_law(5.0/3.0)
U = State(1.0, 0.1, 0.05, 0.2, eos, label="known")
U

\begin{equation}\begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \end{pmatrix}_{known} = \begin{pmatrix} 1.0000 \\ 0.1000 \\ 0.0500 \\ 0.2000 \end{pmatrix}\end{equation}

We then solve for a left going (wavenumber 0) inert acoustic wave with post-wave pressure of \(10\) or \(0.1\):

wave_1 = wave.Wave(U, unknown_value=10, wavenumber=0)
wave_2 = wave.Wave(U, unknown_value=0.1, wavenumber=0)

The first wave is a shock. When output in the notebook, it gives the type, wave number, and wave speed:

wave_1

\begin{equation}{\cal S}_{\leftarrow}: \lambda^{(0)}= -0.9614\end{equation}

The second wave is a rarefaction. In this case, as it’s a continuous wave, its wave speed is a range:

wave_2

\begin{equation}{\cal R}_{\leftarrow}: \lambda^{(0)}\in [-0.3209, -0.2383]\end{equation}

The right going (wavenumber 2) acoustic waves follow the same structure:

wave_right = wave.Wave(U, unknown_value=100, wavenumber=2)
wave_right

\begin{equation}{\cal S}_{\rightarrow}: \lambda^{(2)}= 0.9978\end{equation}

The advective waves are different. The pressure must be constant across the wave, so if anything jumps it must be related to the state: either the thermodynamic variables or the equation of state:

eos2 = eos_defns.eos_gamma_law(4.0/3.0)
U2 = State(2.0, U.v, -0.3, U.p/2.0*3.0, eos2)
wave_contact = wave.Wave(U, unknown_value=U2, wavenumber=1)
wave_contact

\begin{equation}{\cal C}: \lambda^{(1)}= 0.1000\end{equation}

For a reactive wave, we first need a reactive equation of state:

eos_reactive = eos_defns.eos_gamma_law_react(5.0/3.0, 0.1, 1.0, 1.0, eos)
U_reactive = State(5.0, 0.0, 0.0, 2.0, eos_reactive)
U_reactive

\begin{equation}\begin{pmatrix} \rho \\ v_x \\ v_t \\ \epsilon \\ q \end{pmatrix}= \begin{pmatrix} 5.0000 \\ 0.0000 \\ 0.0000 \\ 2.0000 \\ 0.1000 \end{pmatrix}\end{equation}

We can now connect this to states with large (\(p=10\)) and small (\(p=0.1\)) pressures again:

wave_1_reactive = wave.Wave(U_reactive, unknown_value=10, wavenumber=0)
wave_2_reactive = wave.Wave(U_reactive, unknown_value=0.1, wavenumber=0)

The analogue of the inert shock is the strong detonation:

wave_1_reactive

\begin{equation}{\cal SDT}_{\leftarrow}: \lambda^{(0)}= -0.7979\end{equation}

The analogue of the rarefaction is the weak deflagration:

wave_2_reactive

\begin{equation}\left({\cal CJDF}_{\leftarrow}{\cal R}_{\leftarrow}\right) : \lambda^{(0)}\in [-0.6097, 0.7607]\end{equation}

In this case we have a CJ deflagration, which is attached to a rarefaction to complete the wave.