by Mauricio J. Del Razo S. 2014
We consider the one dimensional Riemann problem for the Euler equations with the Tammann EOS. We want to solve 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$ velocity, $E$ the internal energy and $p$ the pressure and the subcripts $x,t$ denote partial derivatives with respect $x$ and $t$. The Tamann EOS is given by $$ p=\rho e(\gamma_k -1) - \gamma_k p_{\infty k}, $$
where $e$ the specific internal energy and $k =L,R$ determines which coefficients to use for the EOS. The initial conditions are given by the left and right constant states $\mathbf{Q}_L=[\rho_L, \rho_L u_L, E_L]$ and $\mathbf{Q}_R=[\rho_R, \rho_R u_R, E_R]$. Note that using the equation of state, the state of the system can also be written on the primitive variables $[\rho, u, p]$.
In order to solve the Riemann problem, we need to find the zero of $\phi(p)=0$. The $\phi(p)$ function is defined below and takes as argument the left and right states in the primitive variables and the left and right parameters of the Tammann equation of state. It returns the value of the function $\phi$ and the corresponding density and velocity in the middle state.
import numpy as np
def phi_exact(p, *initial_data):
rho_l,rho_r,ul,ur,pl,pr,gammal,gammar,pinfl,pinfr = initial_data
global ustar, rhos_l, rhos_r
# Bar pressure (SEE IVINGS paper)
pbl = pl + pinfl
pbr = pr + pinfr
# Star bar pressure
pbsl = p + pinfl
pbsr = p + pinfr
# Useful parameters
gl1 = gammal - 1.0
gr1 = gammar - 1.0
bl = (gammal + 1.0)/(gammal - 1.0)
br = (gammar + 1.0)/(gammar - 1.0)
betal = pbl/bl
betar = pbr/br
al = 2.0/((gammal + 1.0)*rho_l)
ar = 2.0/((gammar + 1.0)*rho_r)
# Calculate velocities
cl = np.sqrt(gammal*(pl + pinfl)/rho_l)
cr = np.sqrt(gammar*(pr + pinfr)/rho_r)
# Calculate the rarefaction phi's
phil_m = ul + 2*cl/gl1*(1 - (pbsl/pbl)**(gl1/(2.0*gammal)))
phir_m = ur - 2*cr/gr1*(1 - (pbsr/pbr)**(gr1/(2.0*gammar)))
# Use shockwave phi's
phil_p = ul - (p - pl)*np.sqrt(al/(pbsl + betal))
phir_p = ur + (p - pr)*np.sqrt(ar/(pbsr + betar))
# Assign value to phi function depending on the case
# Rarefaction - Rarefaction (use both intergal curves)
if p <= pl and p <= pr:
phi = phir_m - phil_m
ustar = 0.5*(phir_m + phil_m)
rhos_l = rho_l*(pbsl/pbl)**(1.0/gammal)
rhos_r = rho_r*(pbsr/pbr)**(1.0/gammar)
#print('Rarefaction-Rarefaction')
# Shockwave - Rarefaction (use Hugoniot locus and integral curve)
elif pl <= p and p <= pr:
phi = phir_m - phil_p
ustar = 0.5*(phir_m + phil_p)
rhos_l = rho_l*(pbsl/pbl + 1.0/bl)/(pbsl/(pbl*bl) + 1.0)
rhos_r = rho_r*(pbsr/pbr)**(1.0/gammar)
# print('Shockwave-Rarefaction')
# Rarefaction - Shockwave (use integral curve and hugoniot locus)
elif pr <= p and p <= pl:
phi = phir_p - phil_m
ustar = 0.5*(phir_p + phil_m)
rhos_l = rho_l*(pbsl/pbl)**(1.0/gammal)
rhos_r = rho_r*(pbsr/pbr + 1.0/br)/(pbsr/(pbr*br) + 1.0)
#print('Rarefaction-Shockwave')
# Shockwave - Shockwave (use both hugoniot locus)
elif p >= pr and p >= pl:
phi = phir_p - phil_p
ustar = 0.5*(phir_p + phil_p)
rhos_l = rho_l*(pbsl/pbl + 1.0/bl)/(pbsl/(pbl*bl) + 1.0)
rhos_r = rho_r*(pbsr/pbr + 1.0/br)/(pbsr/(pbr*br) + 1.0)
# print('Shockwave-Shockwave')
return phi
Now we choose some left and right states and parameters and find a zero fo the phi function,
from scipy.optimize import fsolve
# Define a left state
rho_l = 1.0
ul = 0.0
pl = 3.0
# Define a right state
rho_r = 0.5
ur = 0.0
pr = 1.0
# Define left and right parameters for the EOS
gammal = 1.4
gammar = 1.4
pinfl = 0.0
pinfr = 0.0
# Find the pressure that make phi = 0
initial_data = (rho_l,rho_r,ul,ur,pl,pr,gammal,gammar,pinfl,pinfr)
p0 = (pl + pr)/2.0
pstar = fsolve(phi_exact, p0, args=initial_data)
# Calculate the left and right wave speeds
betal = (pl + pinfl)*(gammal - 1.0)/(gammal + 1.0)
betar = (pr + pinfr)*(gammar - 1.0)/(gammar + 1.0)
alphal = 2.0/(rho_l*(gammal + 1.0))
alphar = 2.0/(rho_r*(gammar + 1.0))
Sl = ul - np.sqrt((pstar + pinfl + betal)/alphal)/rho_l
Sr = ur + np.sqrt((pstar + pinfr + betar)/alphar)/rho_r
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 create a plotting function for our Riemann solution.
# NOT YET SET TO PLOT RAREFACTIONS CORRECTLY **
%matplotlib inline
from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget
import matplotlib.pyplot as plt
from pylab import *
# Create plotting function
def plot_riemann(time):
# Create figure
fig = plt.figure(figsize=(10, 8))
x = np.linspace(-10,10,100)
# Create subplot for (x,t) plane
tmax = 5
ax = fig.add_subplot(2,2,1)
ax.set_xlabel('x')
ax.set_ylabel('time')
ax.axis([min(x),max(x),0,tmax])
# Plot characteristic lines
#
ax.plot(x,x/Sl, '-r', linewidth=2)
ax.plot(x,x/Sr, '-r', linewidth=2)
# Contact discontinuity
ax.plot(x,x/ustar+0.00000001, '--b', linewidth=2)
# Plot time-line in (x,t) plane
ax.plot(x, 0*x+time, 'k', linewidth=3)
# Create Riemann solution vector
sol = np.zeros((3,len(x)))
for i in range(len(x)):
if x[i] < Sl*time:
sol[0,i] = rho_l
sol[1,i] = ul
sol[2,i] = pl
elif x[i] < ustar*time:
sol[0,i] = rhos_l
sol[1,i] = ustar
sol[2,i] = pstar
elif x[i] < Sr*time:
sol[0,i] = rhos_r
sol[1,i] = ustar
sol[2,i] = pstar
else:
sol[0,i] = rho_r
sol[1,i] = ur
sol[2,i] = pr
# Create subplot for density
ax2 = fig.add_subplot(2,2,2)
ax2.set_xlabel('x')
ax2.set_ylabel('density')
ax2.axis([min(x),max(x),-0.5,1.5])
ax2.plot(x,sol[0,:], 'k', linewidth = 3)
ax2.fill_between(x,-20, sol[0,:], facecolor='blue', alpha=0.2)
# Create subplot for velocity
ax3 = fig.add_subplot(2,2,3)
ax3.set_xlabel('x')
ax3.set_ylabel('velocity')
ax3.axis([min(x),max(x),-0.5,1.5])
ax3.plot(x,sol[1,:], 'k', linewidth = 3)
ax3.fill_between(x,-20, sol[1,:], facecolor='blue', alpha=0.2)
# Create subplot for pressure
ax4 = fig.add_subplot(2,2,4)
ax4.set_xlabel('x')
ax4.set_ylabel('pressure')
ax4.axis([min(x),max(x),-5,5])
ax4.plot(x,sol[2,:], 'k', linewidth = 3)
ax4.fill_between(x,-20, sol[2,:], facecolor='blue', alpha=0.2)
return fig
#middle_states = (rhos_l, rhos_r, ustar, pstar)
#all_states = initial_data + middle_states
#plot_riemann(2,Sl,Sr,*all_states)
# Create interactive static widget to visualize solution
StaticInteract(plot_riemann, time=RangeWidget(0,5,0.25))