<jupyter_start><jupyter_text>Nearly-incompressible flow: the Tammann equation of stateIn 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; as usual, 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 that this is almost the same as the [sound speed in an ideal gas](Euler_equations.ipynbCharacteristic-velocities) 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](Euler_equations_TammannEOS.ipynbNegative-pressure-and-vacuum-states).Using the equation of state, the state of the system can also be written in the primitive variables $[\rho, u, p]$.   Exact solutionThe exact solution to the Riemann problem with the Tammann EOS when the fluid parameters $\gamma and p_\infty$ do not vary can be found in (Ivings & Toro 1998). For the more general case when the left and right states may have different values of $\gamma$ and $p_\infty$, the solution can be found in the supplemental material of the work (del Razo, 2017). In this chapter we will focus on the general case. The presentation follows that of (del Razo, 2017).Just as in [the earlier chapter on ideal gas flow](Euler_equations.ipynb), the solution of the Riemann problem considered here consists 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 thus includes four constant regions with solutions $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](Euler_equations.ipynbEntropy-condition).  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_*<p_r$   $\Phi(p_*)= \phi_r^R(p_*) - \phi_l^R(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 whether the function $\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 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 also show how to obtain the density and the missing wave speeds. We write $s_n$ for the speed of the wave belonging to family $n$ ($n=1,2,3$).  Rankine-Hugoniot conditions for shock wavesLet us first consider the case in which both the 1-wave and the 3-wave are shocks. The velocity of eachshock 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 considering $s_n$ and $q_k$, with $n=1, k=l$ for the 1-wave and $s=3, k=r$ for the 3-wave.The Rankine-Hugoniot conditions for the Euler equations can be written as (see e.g. (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 $\omega_k = u_k - s_n$ , $\omega_* = u_* - s_n$, 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 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 equations (\ref{Ql}) and (\ref{Qr}) into (\ref{RH-mom}) to 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 almostobtained 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 knowthat\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}Substituting 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 nonlinear functions $\phi^S_{l,r}(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 functions $\phi_k^S$ 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 wavesRiemann 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}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 along the rarefaction wave, 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 indicate which pairs of states can be connected by a rarefaction. Note that since $\tilde{p}_\kappa = p_\kappa + p_{\infty \kappa}$, when we compare to equations (\ref{CDvel-1}) and (\ref{CDvel-2}), we obtain the functions $\phi_{l,r}^R$. The rarefaction head velocities are given by $u_l-c_l$ and $u_r+c_r$; the tail velocities are $u_*-c_{*l}$ and $u_*+c_{*r}$.In order to compute the complete structure of the rarefaction wave (Ivings & Toro 1998,  LeVeque 2002), we can use the Riemann invariants from Eqs. \eqref{RI1} and \eqref{RI2}, along with Eq. \eqref{sndsp} and the isentropic relation from Eq. \eqref{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. \begin{align*}\Phi(p_*) & = \begin{cases} \phi_r^R(p_*) - \phi_l^R(p_*) & \text{if } p_* < p_l \text{ and } p_*<p_r \\ \phi_r^R(p_*) - \phi_l^S(p_*) & \text{if } p_l \le p_* \le p_r \\ \phi_r^S(p_*) - \phi_l^R(p_*) & \text{if } p_r \le p_* \le p_l \\ \phi_r^S(p_*) - \phi_l^S(p_*) & \text{if } p_* > p_l \text{ and } p_* > p_r              \end{cases}\end{align*}The value of $p_*$ will be found by numerically finding the roots of $\Phi(p_*)=0$. 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 planeAs we did in [the chapter on ideal gas flow](Euler_equations.ipynb), here we visualize the solution in the 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.<jupyter_code>%matplotlib inline
from ipywidgets 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)<jupyter_output><empty_output><jupyter_text>Full solutionIn order to get the full solution, we need some additional 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.  Here we define a function to plot the exact solution of the Riemann problem and parameters that we will use for the different examples.<jupyter_code>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')
    varnames = euler_tammann.primitive_variables
    plot_function = \
            riemann_tools.make_plot_function(ex_states, ex_speeds, 
                                             reval, wave_types, 
                                             layout='vertical', 
                                             variable_names=varnames,
                                             plot_chars=[euler_tammann.lambda1,
                                                         euler_tammann.lambda2,
                                                         euler_tammann.lambda3])
        
    dropdown = widgets.Dropdown(options=[None,1,2,3],
                                description='Show characteristics')
    return interact(plot_function, 
                    t=widgets.FloatSlider(value=0.0,min=0,max=1.0), 
                    which_char=dropdown);<jupyter_output><empty_output><jupyter_text>Examples of Riemann Solutions Problem 1: Shock tube in waterWe first consider a shock tube filled with water, and a shock traveling from left to 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). A value of $p_{atm}=101325 Pa$ is used for Atmospheric pressure.<jupyter_code>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);<jupyter_output><empty_output><jupyter_text>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 consequences of the low compressibility of water.  Problem 2: Air-water interface in shock tubeWe can use the Euler equations with the Tammann EOS to model wave propagation in almost-incompressible media, such as water. Since we also have the solution for different EOS parameters in the left and right states, we can model interaction of different fluids 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 above. 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 suppose that 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 2017). The solution has the following form:<jupyter_code>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);<jupyter_output><empty_output><jupyter_text>Problem 3: Water-air interface in shock tubeHere we consider a setup similar to that of Problem 2, but now the shock moves from water (on the left) into air (on the right).<jupyter_code>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);<jupyter_output><empty_output><jupyter_text>Strong rarefaction wavesAs 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 \eqref{rho-k}, \eqref{isentrop} and the sound speed \eqref{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). Notice that the value in the pressure plot is negative; it reaches a value of about $-2.8\times 10^8$, 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.<jupyter_code>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);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13778758808048222, 'fraction_numerical': 0.015026742507852959, 'mean_word_length': 4.01361991913173, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '468481', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 8028, 'n_tokens_neox': 7579, 'n_words': 2907, 'path': 'Euler_equations_TammannEOS.ipynb', 'repo_name': 'srinivasvl81/riemann_book', 'size': '23558'}