This is the first of two notebooks on the Euler equations. In this notebook, we discuss the equations and the structure of the exact solution to the Riemann problem. In the next notebook, we investigate approximate Riemann solvers.

Table of contents

• Fluid dynamics
• The Euler equations
• Hyperbolic structure
• Exact solution of the Riemann problem
• Plot particle trajectories
• Plot Riemann solution with advected colors
• Interactive Riemann solver
• Riemann problems with vacuum

Fluid dynamics

In this chapter we study the system of hyperbolic PDEs that governs the motions of fluids in the absence of viscosity. These consist of conservation laws for mass, momentum, and energy. Together, they are referred to as the compressible Euler equations, or simply the Euler equations. Our discussion here is fairly brief; for much more detail see (Toro, 2013).

Mass conservation

We will use $\rho(x,t)$ to denote the fluid density and $u(x,t)$ for its velocity. Then the equation for conservation of mass is just the continuity equation:

$$\rho_t + (\rho u)_x = 0.$$

Momentum conservation

The momentum is given by the product of density and velocity, $\rho u$. The momentum flux has two components. First, the momentum is transported in the same way that the density is; this flux is given by the momentum times the density: $\rho u^2$.

To understand the second term in the momentum flux, we must realize that a fluid is made up of many tiny molecules. The density and velocity we are modeling are average values over some small region of space. The individual molecules in that region are not all moving with exactly velocity $u$; that's just their average. Each molecule also has some additional random velocity component. These random velocities are what accounts for the pressure of the fluid, which we'll denote by $p$. These velocity components also lead to a net flux of momentum. Thus the momentum conservation equation is

$$(\rho u)_t + (\rho u^2 + p)_x = 0.$$

Energy conservation

The energy has two components: internal energy $\rho e$ and kinetic energy $\rho u^2/2$:

$$E = \rho e + \frac{1}{2}\rho u^2.$$

Like the momentum flux, the energy flux involves both bulk transport ($Eu$) and transport due to pressure ($pu$):

$$E_t + (u(E+p)) = 0.$$

Equation of state

You may have noticed that we have 4 unknowns (density, momentum, energy, and pressure) but only 3 conservation laws. We need one more relation to close the system. That relation, known as the equation of state, expresses how the pressure is related to the other quantities. We'll focus on the case of a polytropic ideal gas, for which

$$p = \rho e (\gamma-1).$$

Here $\gamma$ is the ratio of specific heats, which for air is approximately 1.4.

The Euler equations

We can write the three conservation laws as a single system $q_t + f(q)_x = 0$ by defining

$$\begin{align} q & = \begin{pmatrix} \rho \\ \rho u \\ E\end{pmatrix}, & f(q) & = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E+p)\end{pmatrix}. \end{align}$$

In three dimensions, the equations are similar. We have two additional velocity components $v, w$, and their corresponding fluxes. Additionally, we have to account for fluxes in the $y$ and $z$ directions. We can write the full system as

$$q_t + f(q)_x + g(q)_y + h(q)_z = 0$$

with

$$\begin{align} q & = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ E\end{pmatrix}, & f(q) & = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho u v \\ \rho u w \\ u(E+p)\end{pmatrix} & g(q) & = \begin{pmatrix} \rho v \\ \rho uv \\ \rho v^2 + p \\ \rho v w \\ v(E+p)\end{pmatrix} & h(q) & = \begin{pmatrix} \rho w \\ \rho uw \\ \rho vw \\ \rho w^2 + p \\ w(E+p)\end{pmatrix}. \end{align}$$

Hyperbolic structure of the 1D Euler equations

In our discussion of the structure of these equations, it is convenient to work with the primitive variables $(\rho, u, p)$ rather than the conserved variables. In quasilinear form, we have

$$\begin{align} \begin{bmatrix} \rho \\ u \\ p \end{bmatrix}_t + \begin{bmatrix} u & \rho & 0 \\ 0 & u & 1/\rho \\ 0 & \gamma \rho & u \end{bmatrix} \begin{bmatrix} \rho \\ u \\ p \end{bmatrix}_x & = 0. \end{align}$$

Characteristic velocities

In primitive variables, the eigenvalues of the flux Jacobian for the 1D Euler equations are:

$$\begin{align} \lambda_1 & = u-c & \lambda_2 & = u & \lambda_3 & = u+c \end{align}$$

Here $c$ is the sound speed:

$$c = \sqrt{\frac{\gamma p}{\rho}}.$$

The eigenvectors of the flux Jacobian are (again in primitive variables):

$$\begin{align} r_1 & = \begin{bmatrix} -\rho/c \\ 1 \\ - \rho c \end{bmatrix} & r_2 & = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} & r_3 & = \begin{bmatrix} \rho/c \\ 1 \\ \rho c \end{bmatrix}. \end{align}$$

Notice that the second characteristic speed, $\lambda_2$, depends only on $u$ and that $u$ does not change as we move in the direction of $r_2$. In other words, the 2-characteristic velocity is constant on 2-integral curves. We say this characteristic field is linearly degenerate; it admits neither shocks nor rarefactions. In a simple 2-wave, all characteristics are parallel. A jump in this family carries a change only in the density, and is referred to as a contact discontinuity.

The other two fields have characteristic velocities that do vary along the corresponding integral curves; thus the 1-wave and the 3-wave in any Riemann solution will be either a shock or a rarefaction. We say these characteristic fields are genuinely nonlinear.

Mathematically, the $p$th field is linearly degenerate if

$$\nabla \lambda_p(q) \cdot r_p(q) = 0$$

and genuinely nonlinear if

$$\nabla \lambda_p(q) \cdot r_p(q) \ne 0.$$

Riemann invariants

Since the Euler equations have three components, we expect each integral curve (a 1D set in 3D space) to be defined by two Riemann invariants. These are:

$$\begin{align} 1 & : s, u+\frac{2c}{\gamma-1} \\ 2 & : u, p \\ 3 & : s, u-\frac{2c}{\gamma-1}. \end{align}$$

Here $s$ is the specific entropy:

$$s = c_v \log(p/\rho^\gamma) + C.$$

The level sets of these Riemann invariants are two-dimensional surfaces; the intersection of two appropriate level sets defines an integral curve.

Integral curves

The 2-integral curves, of course, are simply lines of constant pressure and velocity (with varying density). Since the field is linearly degenerate, these coincide with the Hugoniot loci. We can determine the form of the 1- and 3-integral curves using the Riemann invariants above. For a curve passing through $(\rho_0,u_0,p_0)$, we find

$$\begin{align} u(p) & = u_0 \pm \frac{2c_0}{\gamma-1}\left(1-(p/p_0)^{(\gamma-1)/(2\gamma)}\right). \end{align}$$

Here the plus sign is for 1-waves and the minus sign is for 3-waves.

Below we plot the projection of some integral curves on the $p-u$ plane.

To do: Discuss how $\rho$ fits into this, or plot 3D integral curves.

%matplotlib inline
from exact_solvers import Euler
import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import widgets
from clawpack import riemann
from utils import riemann_tools
import matplotlib
matplotlib.rcParams.update({'font.size': 12})
from collections import namedtuple
Primitive_State = namedtuple('State', Euler.primitive_variables)
gamma = 1.4

from utils.snapshot_widgets import interact                   # for interactive widgets
#from utils.snapshot_widgets import interact      # for static figure that can be viewed online

Will create static figures with single value of parameters

def plot_integral_curves(plot_1=True,plot_3=False,gamma=1.4,rho_0=1.):
    N = 400
    p = np.linspace(0.,5,N)
    p_0 = 1.
    uu = np.linspace(-3,3,15)
    c_0 = np.sqrt(gamma*p_0/rho_0)
    if plot_1:
        for u_0 in uu:
            u = u_0 + (2*c_0)/(gamma-1.)* \
                (1.-(p/p_0)**((gamma-1)/(2*gamma)))
            plt.plot(p,u,color='coral')
    if plot_3:
        for u_0 in uu:
            u = u_0 - (2*c_0)/(gamma-1.)* \
                (1.-(p/p_0)**((gamma-1)/(2*gamma)))
            plt.plot(p,u,color='maroon')
    plt.xlabel('p'); plt.ylabel('u')
    plt.show()
interact(plot_integral_curves,
         gamma=widgets.FloatSlider(min=1.1,max=3,value=1.4));

The structure of centered rarefaction waves

Rankine-Hugoniot jump conditions

The Hugoniot loci for 1- and 3-shocks are $$\begin{align} u(p) & = u_0 \pm \frac{2c_0}{\sqrt{2\gamma(\gamma-1)}} \left(\frac{1-p/p_0}{\sqrt{1+\beta p/p_0}}\right), \ \end{align}$$ where $\beta = (\gamma+1)/(\gamma-1)$. Here the plus sign is for 1-shocks and the minus sign is for 3-shocks.

To do: Discuss how $\rho$ varies, and maybe plot 3D integral curves.

def plot_hugoniot_loci(plot_1=True,plot_3=False,gamma=1.4,rho_0=1.):
    N = 400
    p = np.linspace(1.e-3,5,N)
    p_0 = 1.
    uu = np.linspace(-3,3,15)
    c_0 = np.sqrt(gamma*p_0/rho_0)
    beta = (gamma+1.)/(gamma-1.)
    if plot_1:
        for u_0 in uu:
            u_1 = u_0 + (2*c_0)/np.sqrt(2*gamma*(gamma-1.))* \
                (1.-p/p_0)/(np.sqrt(1+beta*p/p_0))
            plt.plot(p,u_1,color='coral')
    if plot_3:
        for u_0 in uu:
            u_1 = u_0 - (2*c_0)/np.sqrt(2*gamma*(gamma-1.))* \
                (1.-p/p_0)/(np.sqrt(1+beta*p/p_0))
            plt.plot(p,u_1,color='maroon')
    plt.xlabel('p'); plt.ylabel('u')
    plt.show()
interact(plot_hugoniot_loci,
         gamma=widgets.FloatSlider(min=1.1,max=3,value=1.4));

Entropy condition

Exact solution of the Riemann problem

Executing the cell below loads some subroutines that find the exact solution of the Riemann problem. In brief, the Riemann solution is found as follows:

1. Define a piecewise function giving the middle state velocity $u_m$ that can be connected to the left state by an entropy-satisfying shock or rarefaction, as a function of the middle-state pressure $p_m$.
2. Define a piecewise function giving the middle state velocity $u_m$ that can be connected to the right state by an entropy-satisfying shock or rarefaction, as a function of the middle-state pressure $p_m$.
3. Use an iterative solver to find the intersection of the two functions defined above.
4. Use the Riemann invariants to find the intermediate state densities and the solution structure inside any rarefaction waves.

Execute the cell below (after removing #) to bring the Euler solver into the notebook with syntax highlighting, or you can examine it by looking at this file: exact_solvers/Euler.py

#%load exact_solvers/Euler.py

Examples of Riemann solutions

Problem 1: Sod shock tube

First we consider the classic shock tube problem, with high density and pressure on the left, low density and pressure on the right. Both sides are initially at rest. The solution includes a rarefaction, a contact, and a shock.

def riemann_solution(left_state, right_state):
    q_left  = Euler.primitive_to_conservative(*left_state)
    q_right = Euler.primitive_to_conservative(*right_state)

    ex_states, ex_speeds, reval, wave_types = Euler.exact_riemann_solution(q_left ,q_right, gamma)
    
    plot_function = riemann_tools.make_plot_function(ex_states, ex_speeds, reval, wave_types,
                                                     layout='vertical', 
                                                     variable_names=Euler.primitive_variables,
                                                     plot_chars=[Euler.lambda1,Euler.lambda2,Euler.lambda3],
                                                     derived_variables=Euler.cons_to_prim)

    interact(plot_function, t=widgets.FloatSlider(value=0.1,min=0,max=.9),
             which_char=widgets.Dropdown(options=[None,1,2,3],description='Show characteristics'))

left_state  = Primitive_State(Density = 3.,
                              Velocity = 0.,
                              Pressure = 3.)
right_state = Primitive_State(Density = 1.,
                              Velocity = 0.,
                              Pressure = 1.)

riemann_solution(left_state,right_state)

Here is a plot of the solution in the phase plane, showing the integral curve connecting the left and middle states, and the Hugoniot locus connecting the middle and right states.

Euler.phase_plane_plot(left_state, right_state)

Problem 2: Symmetric expansion

Next we consider the case of 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).

left_state  = Primitive_State(Density = 1.,
                              Velocity = -3.,
                              Pressure = 1.)
right_state = Primitive_State(Density = 1.,
                              Velocity = 3.,
                              Pressure = 1.)

riemann_solution(left_state,right_state);

Euler.phase_plane_plot(left_state, right_state)

Problem 3: Colliding flows

Next, consider the case in which the left and right states are moving toward eachother. This leads to a pair of shocks, with a high-density, high-pressure state in between.

left_state  = Primitive_State(Density = 1.,
                              Velocity = 3.,
                              Pressure = 1.)
right_state = Primitive_State(Density = 1.,
                              Velocity = -3.,
                              Pressure = 1.)

riemann_solution(left_state,right_state)

Euler.phase_plane_plot(left_state, right_state)

Plot particle trajectories

In the next plot of the Riemann solution in the $x$-$t$ plane, we also plot the trajectories of a set of particles initially distributed along the $x$ axis at $t=0$, with the spacing inversely proportional to the density. The evolution of the distance between particles gives an indication of how the density changes.

left_state  = Primitive_State(Density = 3.,
                              Velocity = 0.,
                              Pressure = 3.)
right_state = Primitive_State(Density = 1.,
                              Velocity = 0.,
                              Pressure = 1.)

q_left  = Euler.primitive_to_conservative(*left_state)
q_right = Euler.primitive_to_conservative(*right_state)

ex_states, ex_speeds, reval, wave_types = Euler.exact_riemann_solution(q_left ,q_right, gamma)

def reval_rho_u(x): 
    q = reval(x)
    rho = q[0]
    u = q[1]/q[0]
    rho_u = np.vstack((rho,u))
    return rho_u

# Specify density of trajectories to left and right:
rho_l = q_left[0] / 10.
rho_r = q_right[0] / 10.
x_traj, t_traj, xmax = riemann_tools.compute_riemann_trajectories(ex_states, ex_speeds, reval_rho_u, wave_types,
                              i_vel=1, rho_left=rho_l, rho_right=rho_r)
riemann_tools.plot_riemann_trajectories(x_traj, t_traj, ex_speeds, wave_types)
                                        

Recall that the evolution of the distance between particles gives an indication of how the density changes. Note that it increases across the shock wave and decreases through the rarefaction wave, and that in general there is a jump in density across the contact discontinuity.

Plot Riemann solution with advected colors

The next cell defines a function to plot the Riemann solution with the density plot also showing an advected color to help visualize the flow better. The fluid to the left of $x=0$ initially is colored red and to the right of $x=0$ is colored blue, with stripes of different shades of these colors to help visualize the motion of the fluids.

def plot_exact_riemann_solution_stripes(rho_l=3.,u_l=0.,p_l=3.,rho_r=1.,u_r=0.,p_r=1.,t=0.4):    
    q_l = Euler.primitive_to_conservative(rho_l,u_l,p_l)
    q_r = Euler.primitive_to_conservative(rho_r,u_r,p_r)
    
    from matplotlib.mlab import find
    
    x = np.linspace(-1.,1.,1000)
    states, speeds, reval, wave_types = Euler.exact_riemann_solution(q_l, q_r, gamma=gamma)
    q = reval(x/t)
    primitive = Euler.conservative_to_primitive(q[0],q[1],q[2])
    
    # compute particle trajectories:
    def reval_rho_u(x): 
        q = reval(x)
        rho = q[0]
        u = q[1]/q[0]
        rho_u = np.vstack((rho,u))
        return rho_u
    
    # Specify density of trajectories to left and right:
    num_left = 10
    num_right = 10
    rho_left = q_l[0] / 10.
    rho_right = q_r[0] / 10.
    x_traj, t_traj, xmax = riemann_tools.compute_riemann_trajectories(states, speeds, reval_rho_u, wave_types,
                                  i_vel=1, xmax=1, rho_left=rho_left, rho_right=rho_right)
                                                                          
    fig = plt.figure(figsize=(18,6))
    names = ['Density','Velocity','Pressure']
    axes = [0]*3
    for i in range(3):
        axes[i] = fig.add_subplot(1,3,i+1)
        q = primitive[i]
        plt.plot(x,q,linewidth=3)
        plt.title(names[i])
        qmax = max(q)
        qmin = min(q)
        qdiff = qmax - qmin
        axes[i].set_ylim((qmin-0.1*qdiff,qmax+0.1*qdiff))
        axes[i].set_xlim(-xmax,xmax)
                
        if i==0:
            # plot stripes only on density plot
            n = find(t > t_traj)
            if len(n)==0:
                n = 0
            else:
                n = min(n.max(), len(t_traj)-1)

            for i in range(1, x_traj.shape[1]-1):
                j1 = find(x_traj[n,i] > x)
                if len(j1)==0:
                    j1 = 0
                else:
                    j1 = min(j1.max(), len(x)-1)
                j2 = find(x_traj[n,i+1] > x)
                if len(j2)==0:
                    j2 = 0
                else:
                    j2 = min(j2.max(), len(x)-1)

                # set advected color for density plot:
                if x_traj[0,i]<0: 
                    # shades of red for fluid starting from x<0
                    if np.mod(i,2)==0:
                        c = [1,0,0]
                    else:
                        c = [1,0.8,0.8]
                else:
                    # shades of blue for fluid starting from x<0
                    if np.mod(i,2)==0:
                        c = [0,0,1]
                    else:
                        c = [0.8,0.8,1]
                plt.fill_between(x[j1:j2],q[j1:j2],0,color=c)
    plt.show()
 

Make a plot with only a time slider to illustrate this viewpoint with the Sod shock tube data:

def plot_exact_riemann_solution_stripes_t_slider(t):
    plot_exact_riemann_solution_stripes(rho_l=3.,u_l=0.,p_l=3.,rho_r=1.,u_r=0.,p_r=1.,t=t)
    
interact(plot_exact_riemann_solution_stripes_t_slider, 
         t=widgets.FloatSlider(min=0.1,max=1.,step=0.1,value=0.5));

Note the following in the figure above:

• The edges of each stripe are being advected with the fluid velocity, so you can visualize how the fluid is moving.
• The width of each stripe initially is inversely proportional to the density of the fluid, so that the total mass of gas within each stripe is the same.
• The total mass within each stripe remains constant as the flow evolves, and the width of each stripe remains inversely proportional to the local density.
• The interface between the red and blue gas moves with the contact discontinuity. The velocity and pressure are constant but the density can vary across this wave.

Interactive Riemann solver

Here you can set up your own Riemann problem and immediately see the solution. If you don't want to download and run the notebook, an online interactive version is here.

interact(plot_exact_riemann_solution_stripes,
         rho_l=widgets.FloatSlider(min=1.,max=10.,step=0.1,value=3.,description=r'$\rho_l$'),
         u_l=widgets.FloatSlider(min=-10.,max=10.,step=0.1,value=0.,description=r'$u_l$'),
         p_l=widgets.FloatSlider(min=1.,max=10.,step=0.1,value=3.,description=r'$p_l$'),
         rho_r=widgets.FloatSlider(min=1.,max=10.,step=0.1,value=1.,description=r'$\rho_r$'),
         u_r=widgets.FloatSlider(min=-10.,max=10.,step=0.1,value=0.,description=r'$u_r$'),
         p_r=widgets.FloatSlider(min=1.,max=10.,step=0.1,value=1.,description=r'$p_r$'),
         t=widgets.FloatSlider(min=0.1,max=1.,step=0.1,value=0.5));

Riemann problems with vacuum

A vacuum state (with zero pressure and density) can arise in the solution of the Riemann problem in two ways:

1. An initial left or right vacuum state: in this case the Riemann solution consists of a single rarefaction, connecting the non-vacuum state to vacuum.
2. A problem where the left and right states are not vacuum but middle states are vacuum. Since this means the middle pressure is smaller than that to the left or right, this can occur only if the 1- and 3-waves are both rarefactions. These rarefactions are precisely those required to connect the left and right states to the middle vacuum state.

Initial vacuum state

The velocity plot looks a bit strange, but note that the velocity is undefined in vacuum.

left_state  = Primitive_State(Density =0.,
                              Velocity = 0.,
                              Pressure = 0.)
right_state = Primitive_State(Density = 1.,
                              Velocity = -3.,
                              Pressure = 1.)

riemann_solution(left_state,right_state)

Euler.phase_plane_plot(left_state, right_state)

Middle vacuum state

left_state  = Primitive_State(Density =1.,
                              Velocity = -10.,
                              Pressure = 1.)
right_state = Primitive_State(Density = 1.,
                              Velocity = 10.,
                              Pressure = 1.)

riemann_solution(left_state,right_state)

Euler.phase_plane_plot(left_state, right_state)