Riemann problem for Euler equations with the Tammann EOS¶
In this chapter, we consider the one-dimensional Riemann problem for the Euler equations with the Tammann equation of state (EOS). The Tammann EOS, also known as stiffened gas EOS, is used to model almost-incompressible fluids. It can be physically understood as modeling an ideal gas under very high ambient pressures. The Riemann problem we want to solve consists of the one-dimensional Euler equations,
\begin{align*} \left[\begin{array}{c} \rho \\ \rho u \\ E \end{array} \right]_t + \left[\begin{array}{c} \rho u \\ \rho u^2 + p \\ u(E+p) \end{array} \right]_x = 0, \end{align*}where $\rho$ is density, $u$ the velocity, $E$ the internal energy, and $p$ the pressure; the subcripts $x,t$ denote partial derivatives. The initial conditions are given by the left and right constant states $q_l=[\rho_l, \rho_l u_l, E_l]$ and $q_r=[\rho_r, \rho_r u_r, E_r]$. The Tamann EOS is given by
\begin{align*} p (\rho,e)=\rho e(\gamma -1) - \gamma p_{\infty}, \end{align*}where $e$ is the specific internal energy. For a general equation of state $p (\rho,e)$, the speed of sound is given by the isentropic derivative (keeping entropy constant) of the pressure,
\begin{align*} c = \sqrt{\frac{\partial p(\rho,e)}{\partial \rho} \bigg|_s} = \sqrt{\frac{\partial p(\rho,e)}{\partial \rho} + \frac{p(\rho,e)}{\rho^2}\frac{\partial p(\rho,e)}{\partial e}}. \end{align*}In the case of the Tammann EOS, we obtain
\begin{align*} c= \sqrt{\gamma\frac{p+p_{\infty }}{\rho}}. \end{align*}Note this is almost the same as the sound speed in an ideal gas except for the additional $p_{\infty}$ term. This means that the fluid behaves similarly to an ideal gas that is already at pressure $p_{\infty }$. When $p_{\infty } \gg 1$, the Tammann EOS emulates an ideal gas at a very high pressure which can be used to model almost-incompressible fluids like water. Note the sound speed (as well as the density and any other physical quantity) is well defined for negative values of $p$, with the only limitation that $p+p_{\infty}>0$. We will revisit this issue in the last section of this chapter.
Using the equation of state, the state of the system can also be written in the primitive variables $[\rho, u, p]$.
Exact solution¶
The exact solution to the Riemann problem with a Tammann EOS can be found in (Ivings & Toro 1998). In the supplemental material of the work (del Razo, 2017), the reader can find an extension of the analytical solution to the case where the parameters $\gamma$ and $p_\infty$ in the Tammann EOS have a discontinuity. In this chapter we will focus in the latter, since it is a more general case. We will reproduce the derivation of the exact solution while at the same time showing some sricpts and plots to illustrate the structure of the solution.
As we know from the previous examples, the solution of the Euler equations will consist of three waves: the 1-wave and 3-wave are genuinely nonlinear (rarefactions or shocks), while the 2-wave is a linearly degenerate contact discontinuity. The solution will have four different solution states, $q_l,q_{*l},q_{*r},q_r$ separated by the three waves. Just as in the case of the ideal gas, the velocity and pressure do not change across the contact wave, so we denote them simply by $u_*$ and $p_*$. In order to determine whether the 1-wave and 3-wave are rarefactions or shocks, we will need to create a function of the middle state pressure $p_*$ that ensures that the velocity $u_*$ across the contact discontinuity is consistent. Since we know that the left state velocity $u_l$ should be connected by a rarefaction or shock to $u_*$, we write $u_*=u_l + [u]_1$, where $[u]_1$ is the jump of the velocity across the 1-wave. In a similar manner, we also know the 3-wave should be a shock or rarefaction, so we write $u_*= u_r - [u]_3$. Therefore, it is useful to define
\begin{align} \phi_l(p_*) = u_* = u_l - \mathcal{F}_l(p_*), \label{CDvel-1} \\ \phi_r(p_*) = u_* = u_r + \mathcal{F}_r(p_*). \label{CDvel-2} \end{align}where the value of $\mathcal{F}_{l,r}(p_*) = -[u]_{1,3}$ depends on whether the wave in question is a shock or a rarefaction (signs were chosen for notation consistency). As we expect these two equations yield the same contact discontinuity velocity $u_*$, we have
\begin{align} \Phi(p_*)= \phi_r(p_*) - \phi_l(p_*) = 0. \label{Phi} \end{align}This nonlinear equation will yield the pressure $p_*$ that provides consistency between the type of waves (rarefactions or shocks), their speeds and the contact discontinuity velocity $u_*$. As we mentioned before, the shape of $\phi_k(p_*)$ will depend on whether the states are connected by a shock wave or rarefaction. Once the $p_*$ has been found, the contact discontinuity velocity can be found from the expressions just derived. However, it is not yet clear how to calculate the density or the speeds of the 1-wave and 3-wave. It will become obvious how to obtain these quantities once we write the explicit equations for our system further below.
Before coding the exact solution, we should first note that having a rarefaction or shock in the 1-wave and 3-wave will depend on the pressure $p_*$, in the same way it did for an ideal gas. If the pressure is higher on the side toward which the wave is propagating (i.e., in front of the wave), it must be a rarefaction. If the pressure is lower in front of the wave, it must be a shock. In the Euler equations, this yields four possible cases for the value $\Phi(p_*)$, just as in the solution with an ideal gas EOS:
- 1-rarefaction, 3-rarefaction: $p_*< p_l$ and $p_*
- 1-shock, 3-rarefaction $p_l \le p_* \le p_r$
$\Phi(p_*)= \phi_r^R(p_*) - \phi_l^S(p_*)$, - 1-rarefaction, 3-shock $p_r \le p_* \le p_l$
$\Phi(p_*)= \phi_r^S(p_*) - \phi_l^R(p_*)$, - 1-shock, 3-shock: $p_* > p_l$ and $p_*>p_r$
$\Phi(p_*)= \phi_r^S(p_*) - \phi_l^S(p_*)$,
where the index $S,R$ indicates if the $\phi$ was obtained by using the Rankine-Hugoniot equations to connect states by shocks or the Riemann invariants to connect them by rarefactions respectively.
We will now derive the functions $\phi_{k}^{\mu}$ for all four cases with $k=l,r$ (left and right) and $\mu=R,S$ (rarefaction and shock). We will also show how to obtain the density and the missing wave speeds. We denote the speed of the 1-wave, $s_1$, the 2-wave, $s_2$, and the 3-wave $s_3$.
Rankine-Hugoniot conditions for shock waves¶
As we know the 1-wave and the 3-wave could each be a shock. In that case, the velocity of the wave, i.e. the shock, will be given by the Rankine-Hugoniot conditions. These conditions will also allow us to calculate $\phi_k^S(p_*)$ with $k=l,r$. The following analysis is written to be applicable to both the 1-wave velocity $s_1$ and the 3-wave velocity $s_3$, by writing $s_n$, with $n=1,3$.
The Rankine-Hugoniot conditions are in general given by $s_n\left(q_k - q_{*k}\right) = f(q_k) - f(q_{*k})$, where $q$ is the vector state variable, $f(q)$ the vector state flux and $s_n$ the shock velocity. For the Euler equations this can be easily rewritten as ((Ivings & Toro 1998)),
\begin{align} \rho_k \omega_k = \rho_{*k} \omega_*, \label{RH-mass} \\ \rho_k \omega_k^2 + p_k = \rho_{*k} \omega_*^2 + p_{*k}, \label{RH-mom} \\ \frac{1}{2} \omega_k^2 + h_k = \frac{1}{2} \omega_*^2 + h_{*k}, \label{RH-ene} \end{align}where $k=l,r$, $\omega_k = u_k - s_n$ , $\omega_* = u_* - s_n$, with $n=1$ when $k=l$ and $n=3$ when $k=r$, and the specific enthalpy is given by $h = e+(p+p_\infty)/\rho$ with $e$ the specific internal energy that relates to the internal energy of our original variables by $E = \rho e + \rho u^2/2$. We will use these relations to find $\phi_K^S(p_*)$ and the wave speeds $s_n$.
Finding $\phi_l^S(p_*)$ and $\phi_r^S(p_*)$ and $s_1$ and $s_3$:¶
We can start by defining the mass fluxes $F_k$ for $k=l,r$ as
\begin{align} F_l = \rho_l \omega_l = \rho_{*l} \omega_* \label{Ql} \\ F_r = -\rho_r \omega_r = -\rho_{*r} \omega_*. \label{Qr} \end{align}As $\omega_k = u_k - s_n$, from these two equations we can obtain the wave speeds in terms of $F_l$ and $F_r$,
\begin{align} s_1 = u_l - \frac{F_l}{\rho_l}, \ \ \ s_3 = u_r + \frac{F_r}{\rho_r}. \label{RH-speeds} \end{align}We still need to find $F_l$ and $F_r$, so we substitute equation (\ref{Ql}) and (\ref{Qr}) into (\ref{RH-mom}) to immediately obtain
\begin{align} F_l &= \frac{\tilde{p}_{*l} -\tilde{p}_l}{\omega_l - \omega_*} = \frac{p_* -p_l}{u_l - u_*} \label{Ql-1} \\ F_r &= -\frac{\tilde{p}_{*r} -\tilde{p}_r}{\omega_r - \omega_*} = -\frac{p_* -p_r}{u_r - u_*}, \label{Qr-1} \end{align}where $\tilde{p}_\kappa = p_\kappa + p_{\infty \kappa}$ is defined to simplify future notation with $\kappa = l, *l, *r$ and $r$. Note that $\tilde{p}_{*l} \neq \tilde{p}_{*r}$ and that $\tilde{p}_{*k} -\tilde{p}_k=p_* -p_k$ since $p_*=p_{*l}=p_{*r}$ and $p_{\infty k} = p_{\infty *k}$ ($k=l,r)$. Solving for $u_*$ we obtain the equations,
\begin{align} u_* = u_l - \frac{\tilde{p}_{*l} - \tilde{p}_l}{F_l} = \phi_l^S(p_*) \\ u_* = u_r + \frac{\tilde{p}_{*r} - \tilde{p}_r}{F_r} = \phi_r^S(p_*) \label{RH-phis} \end{align}Comparing to equations (\ref{CDvel-1}, \ref{CDvel-2}), we notice $\mathcal{F}_k(p_*) = \frac{p_* - p_k}{F_k}$. We also notice we have almost obtained the $\phi$ functions we are looking for, though we still need to find $F_l$ and $F_r$ in terms of known variables.
Finding $F_l$ and $F_r$:¶
From equations (\ref{Ql}) and (\ref{Ql-1}), we know that
\begin{align*} \frac{\tilde{p}_{*l} -\tilde{p}_l}{\omega_l - \omega_*} = \rho_l \omega_l. \end{align*}Solving for $\omega_*$, substituting the solution into (\ref{Ql}) and subtituting the $w_l$ for $F_l/\rho_l$, we obtain a new equation that we can solve for $F_l$ that yields,
\begin{align} F_k = \sqrt{\rho_k \rho_{*k} \frac{\tilde{p}_{*k} -\tilde{p}_k}{\rho_{*k} - \rho_k}}. \label{Qk} \end{align}with $k=l,r$, since we repeated the same process for $F_r$ and obtained exactly the same equation. However, we still don't know $\rho_{*k}$, for this we will need our third Rankine-Hugoniot condition (\ref{RH-ene}).
Finding $\rho_{*l}$ and $\rho_{*r}$:¶
From equation (\ref{RH-ene}), we can obtain
\begin{align} h_{*k} - h_k &= \frac{1}{2}\left(w_k^2 - w_*^2\right) \nonumber\\ &= \frac{1}{2}\left(\pm \frac{F_k^2}{\rho_k^2} \mp \frac{F_k^2}{\rho_{*k}^2}\right) \nonumber \\ &= \frac{1}{2}\left(\frac{1}{\rho_k} + \frac{1}{\rho_{*k}}\right)(\tilde{p}_{*k}-\tilde{p}_k), \label{delta-h} \end{align}where the sign above is used for $k=l$ and the one below for $k=r$, and we used equations (\ref{Ql}) and (\ref{Qr}) for the second line and (\ref{Qk}) for the third line. We can now substitute the specific enthalpy $h=\gamma \tilde{p}/(\rho(\gamma -1))$ in equation (\ref{delta-h}) to obtain,
\begin{align*} \frac{\gamma_{*k}}{\gamma_{*k} - 1}\frac{\tilde{p}_{*k}}{\rho_{*k}} - \frac{\gamma_{k}}{\gamma_{k} - 1}\frac{\tilde{p}_k}{\rho_{k}} =\frac{1}{2}\left(\frac{1}{\rho_k} + \frac{1}{\rho_{*k}}\right)(\tilde{p}_{*k}-\tilde{p}_k). \end{align*}As the interface is the contact discontinuity, the jump in the parameters is only across the contact discontinuity, so $\gamma_{*k}=\gamma_k$. Now we can solve for the unknown density,
\begin{align} \rho_{*k} = \rho_k \left(\frac{\frac{\tilde{p_*}}{\tilde{p_k}} + \frac{\gamma_k-1}{\gamma_k+1}}{\frac{\tilde{p_*}}{\tilde{p_k}}\frac{\gamma_k-1}{\gamma_k+1}+1}\right). \label{rho-k} \end{align}Replacing this result into equation (\ref{Qk}), we obtain $F_k$ in terms of $p_*$ and known variables,'
\begin{align} F_k = \sqrt{\rho_k \frac{\tilde{p}_{*k} + \tilde{p}_{k}\frac{\gamma_k - 1}{\gamma_k +1}}{\frac{2}{\gamma_k+1}}}. \label{Qk-f} \end{align}Subsitituting (\ref{Qk-f}) into (\ref{RH-phis}), we can calculate the $\phi^S_{l,r}$ nonlinear functions of $p_*$ in terms of known variables, which yields
\begin{align} \phi_l^S(p_*) = u_l - (p - pl)\sqrt{\frac{\alpha_l}{p_* + p_{\infty l} + \beta_l}} \\ \phi_r^S(p_*) = u_r + (p - pr)\sqrt{\frac{\alpha_r}{p_* + p_{\infty r} + \beta_r}} \end{align}with $\alpha_k = 2.0/((\gamma_k + 1.0)\rho_l)$ and $\beta_k = (p_l + p_{\infty l})(\gamma_l - 1.0) / (\gamma_l + 1.0)$. The $\phi_k^S$ functions are called the Hugoniot Loci, and they allow us to construct equation (\ref{Phi}) and solve it using a Newton method or other root finder in order to obtain the value of $p_*$.
Riemann invariants for rarefaction waves¶
Riemann invariants are variables that remain constant through simple waves such as rarefactions. The Riemann invariants across the 2-wave are the pressure $p_*$ and the normal velocity $u_*$. The Riemann invariants for the 1-wave and 3-wave are the entropy and the quantities,
\begin{align} u_l+\frac{2c_l}{\gamma_l-1} = u_* + \frac{2c_{*l}}{\gamma_l-1}, \label{RI1} \\ u_r-\frac{2c_r}{\gamma_r-1} = u_* - \frac{2c_{*r}}{\gamma_r-1}, \label{RI2} \end{align}correspondingly. The speed of sound $c_k$ with the Tammann EOS is given by
\begin{align} c_k= \sqrt{\gamma_k\frac{p_k+p_{\infty k}}{\rho_k}}. \label{sndsp} \end{align}As the entropy is invariant, we can use the Tammann EOS isentropic relation to obtain the density in the middle states,
\begin{align} \rho_{*k} = \rho_k \left(\frac{\tilde{p}_{*k}}{\tilde{p}_k}\right)^{1/\gamma}. \label{isentrop} \end{align}Solving (\ref{RI1}) and (\ref{RI2}) for $u_*$ and using equations (\ref{sndsp}) and (\ref{isentrop}), we immediately obtain
\begin{align*} u_* = u_l + \frac{2c_l}{\gamma_l-1} \left[ 1-\left(\frac{\tilde{p}_{*l}}{\tilde{p}_l}\right)^{\frac{\gamma_l-1}{2\gamma_l}} \right] = \phi_l^R(p_*), \\ u_* = u_r - \frac{2c_r}{\gamma_r-1} \left[ 1-\left(\frac{\tilde{p}_{*r}}{\tilde{p}_r}\right)^{\frac{\gamma_r-1}{2\gamma_r}} \right] = \phi_r^R(p_*). \end{align*}These functions are called the integral curves, and they tell you which values are acceptable to connect two states with a rarefaction curve. Note that as $\tilde{p}_\kappa = p_\kappa + p_{\infty \kappa}$, when we compare to equations (\ref{CDvel-1} and \ref{CDvel-2}), we obtain the $\phi_{l,r}^R$ functions. The rarefaction head velocities will be given by $u_l-c_l$ and $u_r+c_r$; the tail velocities will be $u_*-c_{*l}$ and $u_*+c_{*r}$. For numerical purposes, a simple approximate velocity is provided for $s_1$ and $s_3$ as the average between the head and tail velocity.
In order to compute the complete structure of the rarefaction wave (Ivings & Toro 1998, LeVeque 2002), we can use the Riemann invariants from Eqs. \ref{RI1} and \ref{RI2}, along with Eq. \ref{sndsp} and the isentropic relation from Eq. \ref{isentrop}. The solution for the 1-rarefaction wave along the rays $x/t=\xi= u_{rar1} - c_{rar1}$ is then
\begin{align*} u_{rar1}(\xi) &= \frac{u_l(\gamma_l - 1) + 2(\xi + c_l)}{\gamma_l + 1}, \\ \rho_{rar1}(\xi) &= \rho_l \left[\frac{u_{rar1}(\xi) - \xi}{c_l}\right]^{\frac{2}{\gamma_l -1}}, \\ p_{rar1}(\xi) &= \tilde{p}_l \left[\frac{u_{rar1}(\xi) - \xi}{c_l}\right]^{\frac{2\gamma_l}{\gamma_l -1}} - p_{\infty l}, \end{align*}and for a 3-rarefaction wave along the rays $x/t=\xi= u_{rar3} + c_{rar3}$ is,
\begin{align*} u_{rar3}(\xi) &= \frac{u_r(\gamma_r - 1) + 2(\xi - c_r)}{\gamma_r + 1}, \\ \rho_{rar3}(\xi) &= \rho_r \left[\frac{\xi - u_{rar3}(\xi)}{c_r}\right]^{\frac{2}{\gamma_r -1}}, \\ p_{rar3}(\xi) &= \tilde{p}_r \left[\frac{\xi - u_{rar3}(\xi)}{c_r}\right]^{\frac{2\gamma_r}{\gamma_r -1}} - p_{\infty r}. \end{align*}Now that we know the integral curves given by $\phi^{s,r}_{l,r}$ (for the rarefactions), we can construct the function $\Phi(p_*)$ from (\ref{Phi}) for any of the 4 possible scenarios. The value of $p_*$ will be found by numerically finding the roots of $\Phi(p_*)=0$. Note which case to employ to calculate $\Phi(p_*)$ might change in each iteration of the root finder. Once $p_*$ is found, $u_*$, $\rho_{*l}$, $\rho_{*r}$, $s_1$ and $s_3$ can be found using the relations above. As we know the three wave speeds $s_1$, $s_2$ and $s_3$ and the values of the primitive variables $[\rho,u,p]$ in all the 4 states, we have solved the Riemann problem.
Solution in the pressure-velocity plane¶
The easiest way to visualize the solution is in the phase plane whose axes are pressure and velocity. We have a left and a right state such that each one is connected to some middle state by a shock or a rarefaction, depending on the entropy conditions. The connection is established by curves in the phase plane: hugoniot loci in the case of shocks and integral curves in the case of rarefactions. The intersection between these curves yields the middle state, which is the most important element of the solution to the Riemann problem. Lets look at an interactive plot illustrating the solution in the phase plane. The initial states and the middle state are indicated as black points; the Hugoniot loci are plotted in red and the integral curves in blue. The physical solution satisfying the entropy conditions is plotted as a continous line, whereas the unphysical solution is shown as a dashed line. Note that, in this case, the phase plane is three dimensional, but it is easierly understood when projecting into a two dimensional phase space. Nonetheless, the variable $\rho$ in the flattened dimension can still be modified as an additional paramater of the solution.
%matplotlib inline
from utils.snapshot_widgets import interact
from ipywidgets import widgets
import matplotlib.pyplot as plt
from exact_solvers import euler_tammann, interactive_pplanes
from utils import riemann_tools
# Initial states [density, velocity, pressure]
ql = [600.0, 10.0, 50000.0]
qr = [50.0, -10.0, 25000.0]
# Tammann EOS parameters [gamma, pinf]
paramsl = [1.4, 0.0]
paramsr = [7.0, 100.0]
interactive_pplanes.euler_tammann_interactive_phase_plane(ql,qr,paramsl,paramsr)
Full solution¶
In order to get the full solution, we still need a couple more quantities. Equations (\ref{RH-phis}) will yield the contact discontinuity speed $s_2=u_*$ in terms of $p_*$. We can also calculate $F_l$ and $F_r$ from (\ref{Qk-f}), and we can substitute in (\ref{RH-speeds}) to obtain the corresponding wave speeds. With this, we can provide the full solution of the Riemann problem. In the next section, we will show the full solution for different examples. It should be noted that this model is to be handled carefully because it could yield unphysical results (negative pressures) for some initial conditions. The cases where this happens are in the limit of validity for modeling almost incompressible media with the Tammann EOS. We will explore these limits in Problem 4. Here we define a function to plot the exact solution of the Riemann problem and parameters that we will use for the different examples.
def plot_riemann_solution(q_l, q_r, gamma, pinf):
ex_states, ex_speeds, reval, wave_types, varsout = \
euler_tammann.exact_riemann_solution(q_l ,q_r, gamma, pinf,
varin='primitive', varout='primitive')
plot_function = riemann_tools.make_plot_function(ex_states, ex_speeds,
reval, wave_types,
layout='vertical',
variable_names=euler_tammann.primitive_variables,
plot_chars=[euler_tammann.lambda1,
euler_tammann.lambda2,
euler_tammann.lambda3])
return interact(plot_function, t=widgets.FloatSlider(value=0.0,min=0,max=1.0),
which_char=widgets.Dropdown(options=[None,1,2,3],description='Show characteristics'));
Examples of Riemann Solutions¶
Problem 1: Shock tube in water¶
We first consider a shock tube filled with water, and a shock traveling from left two right. The propagation of waves in water can be accurately modeled using $\gamma = 7.15$, $p_{\infty} = 3\times 10^8$ and $\rho\approx 1000 kg/m^3$, see (del Razo 2016). Atmospheric pressure is modeled at $p_{atm}=101325 Pa$ and speed in $m/s$.
gamma_water = 7.15
gamma_air = 1.4
pinf_water = 3.e8
pinf_air = 0.
patm = 101325.0
q_l = [1010.0, 0.0, 3*patm]
q_r = [1000.0, 0.0, patm]
gamma = [gamma_water, gamma_water]
pinf = [pinf_water, pinf_water]
plot_riemann_solution(q_l, q_r, gamma, pinf);
The density difference between the states is very small because water is almost incompressible. Note the small waves generated in the density plot and the very compressed rarefaction fan, which are consequence of the low compressibility in water.
Problem 2: Air-water interface in shock tube¶
We can use the Euler equations with the Tammann EOS to model wave propagation in almost-incompressible media, such as water. As we provided the solution with different EOS parameters on the left and right states, it can couple different materials across the interface. For instance, we can use this model to study the propagation of waves from air into water. Air is modeled by the ideal gas EOS, which corresponds to $\gamma=1.4$ and $p_{\infty} = 0$, and water can be modeled using the same parameters as the example before. Consider a shock traveling from left to right; the left state has only air while the right state has only water. Although in a realistic two or three dimensional setting, the water would spill into the left state, we can consider the water is being held in place by a container with a thin wall. In this setting, the current model has been shown to be a valid approximation (del Razo 2016). The solution has the following form:
q_l = [1.0, 350.0, 2*patm] # Initial condition for shock in air
q_r = [1000.0, 0.0, patm] # Initial condition in still water
gamma = [gamma_air, gamma_water]
pinf = [pinf_air, pinf_water]
plot_riemann_solution(q_l, q_r, gamma, pinf);
Problem 3: Water-air interface in shock tube¶
The same setup as in Problem 2, but now the shock moves from water (on the left) into air (on the right).
q_l = [1000.0, 350.0, 2*patm] # Initial condition for shock in air
q_r = [1.0, 0.0, patm] # Initial condition in still water
gamma = [gamma_water, gamma_air]
pinf = [pinf_water, pinf_air]
plot_riemann_solution(q_l, q_r, gamma, pinf);
Negative pressure and vacuum states¶
As we mentioned when we introduced the Tammann EOS, the relevant physical quantities remain well-defined for negative pressures, as long as $p+p_{\infty}>0$. For instance, the densities \ref{rho-k}, \ref{isentrop} and the sound speed \ref{sndsp} are all well defined in this case. This is consistent with the interpretation of the Tammann EOS as a shift of the zero pressure point of an ideal gas.
As an example of this scenario, we consider a strong symmetric expansion wave in water. Analogous to the Euler equations using an ideal gas EOS, we consider equal densities and pressures, and equal and opposite velocities, with the initial states moving away from each other. The result is two rarefaction waves (the contact has zero strength). However, notice the value in the pressure plot is negative; it reaches a value of about $-2.8e8$, so $p+p_{\infty}$ is still positive. Of course, it is questionable whether the model remains accurate when $p+p_{\infty}$ is nearly zero, as it is in this example.
velocity = 350
q_l = [1000.0, -velocity, 2*patm] # Initial condition for shock in air
q_r = [1000.0, velocity, 2*patm] # Initial condition in still water
gamma = [gamma_water, gamma_water]
pinf = [pinf_water, pinf_water]
plot_riemann_solution(q_l, q_r, gamma, pinf);
It may be possible to extend the solution to cases where the Riemann solver yields $p+p_{\infty} \le 0$. As in the Euler equations chapter, the vacuum states arise in regions with zero density. As vacuum states can only connected through rarefaction waves, the middle state densities can only be given by \ref{isentrop}. The middle state density will only yield zero when $\tilde{p_*} = p_*+p_{\infty}=0$ or $p_*=-p_{\infty}$. The convergence failure of the Riemann solution for this scenario can be observed if the velocity in the example above is increased one order of magnitude. How to deal with vacuum states when using a Tammann EOS with a jump in the parameters is beyond the scope of this book.
Full solution¶
In order to get the full solution, we still need a couple more quantities. Equations (\ref{RH-phis}) will yield the contact discontinuity speed $s_2=u_*$ in terms of $p_*$. We can also calculate $F_l$ and $F_r$ from (\ref{Qk-f}), and we can substitute in (\ref{RH-speeds}) to obtain the corresponding wave speeds. With this, we now provide the full solution of the Riemann problem.
# %load exact_solvers/euler_tammann.py
ql = [1.0, -3.0, 1.0]
qr = [1.0, 3.0, 1.0]
gamma = [1.4, 1.4]
pinf = [0.0, 0.0]
ex_states, ex_speeds, reval, wave_types, varsout = euler_tammann.exact_riemann_solution(ql ,qr, gamma, pinf,
varin = 'primitive', varout = 'conservative')
plot_function = riemann_tools.make_plot_function(ex_states, ex_speeds, reval, wave_types,
layout='vertical', variable_names=varsout)
interact(plot_function, t=widgets.FloatSlider(value=0.0,min=0,max=1.0));
to be modified¶
Now we need to iterate this function using a Newton method, to find which value of $p_*$ yields $\Phi(p_*)=0$ We just obtained the pressure and velocity middle states $p_*$ and $u_*$ that are continuous across the contact disconitnuity and the middle states for the density $\rho_{*l}$ and $\rho_{*r}$ have also been obtained from the function and saved as global variables. We also know the speed the left and right waves are traveling. We can now plot the solution to the Riemann problem.
*Test: Interactive pplane for ideal gas Euler eqs.¶
# Initial states [density, velocity, pressure]
ql = [1.0, -10.0, 100.0]
qr = [1.0, 10.0, 100.0]
# EOS parameter
gamma = 1.4
interactive_pplanes.euler_interactive_phase_plane(ql,qr,gamma)
