from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
In this lesson we will look at 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.
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 we discussed in Lesson 1:
$$\rho_t + (\rho u)_x = 0.$$The momentum is given by the product of density and velocity, i.e. $\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; i.e. $\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.$$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))_x = 0.$$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 an ideal gas, for which
$$p = \rho e (\gamma-1).$$Here $\gamma$ is the ratio of specific heats, which for air is approximately 1.4.
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}These equations can be solved in a manner similar to what we used for advection and traffic flow. As you might guess, computing the flux gets significantly more complicated since we now have 3 (or 5) equations and more complicated flux expressions.
Implementing a solver for the Euler equations from scratch would be a lot of fun, but to save some time we'll use a package called PyClaw, which is part of the Clawpack software (Clawpack stands for Conservation LAWs PACKage). PyClaw allows us to quickly and easily set up and solve problems modeled by hyperbolic PDEs.
Now let's get started. First, import the parts of Clawpack that we'll use:
from clawpack import pyclaw
from clawpack import riemann
To solve a problem, we'll need to create the following:
This might sound complicated at first, but stick with me.
Let's start by creating a controller and specifying the simulation end time:
claw = pyclaw.Controller()
claw.tfinal = 0.1 # Simulation end time
claw.keep_copy = True # Keep solution data in memory for plotting
claw.output_format = None # Don't write solution data to file
claw.num_output_times = 50 # Write 50 output frames
The method used to compute the flux between each pair of cells is referred to as a Riemann solver. By specifying a Riemann solver, we will specify the system of PDEs that we want to solve. So far we have only used very simple approximate Riemann solvers. Clawpack includes much more sophisticated Riemann solvers for many hyperbolic systems.
Place your cursor at the end of the line in the box below and hit 'Tab' (for autocompletion). You'll see a dropdown list of all the Riemann solvers currently available in PyClaw. The ones with 'py' at the end of the name are written in pure Python; the others are written in Fortran and wrapped with f2py.
riemann.
We'll start with a simple 1D problem, using the Riemann solver riemann.euler_with_efix_1D
:
riemann_solver = riemann.euler_with_efix_1D
claw.solver = pyclaw.ClawSolver1D(riemann_solver)
We also need to specify boundary conditions. We'll use extrapolation boundary conditions, so that waves simply pass out of the domain:
claw.solver.all_bcs = pyclaw.BC.extrap
Next we need to specify the domain and the grid. We'll solve on the unit line $[0,1]$ using 100 grid cells. Note that each argument to the Domain constructor must be a tuple:
domain = pyclaw.Domain( (0.,), (1.,), (100,))
Next we create a solution object that belongs to the controller and extends over the domain we specified:
claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)
The initial data is specified in an array named solution.q
. The density is contained in q[0,:]
, the momentum in q[1,:]
, and the energy in q[2,:]
.
x=domain.grid.x.centers # grid cell centers
gam = 1.4 # ratio of specific heats
rho_left = 1.0; rho_right = 0.125
p_left = 1.0; p_right = 0.1
claw.solution.q[0,:] = (x<0.5)*rho_left + (x>=0.5)*rho_right
claw.solution.q[1,:] = 0.
claw.solution.q[2,:] = ((x<0.5)*p_left + (x>=0.5)*p_right)/(gam-1.0)
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x, claw.solution.q[0,:],'-o')
This problem is known as the Sod shock-tube. It amounts to setting up a tube with a thin separator between a high-pressure, high-density region and a low-pressure, low-density region, then suddenly removing the separator.
Next we need to specify the value of $\gamma$, the ratio of specific heats.
problem_data = claw.solution.problem_data
problem_data['gamma'] = 1.4
problem_data['gamma1'] = 0.4
Finally, let's run the simulation.
claw.run()
Now we'll plot the results, which are contained in a list called claw.frames
. It's simple to plot a single frame with matplotlib:
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
import numpy as np
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-0.2, 1.2))
frame = claw.frames[0]
pressure = frame.q[0,:]
line, = ax.plot([], [], 'o-', lw=2)
def fplot(frame_number):
frame = claw.frames[frame_number]
pressure = frame.q[0,:]
line.set_data(x,pressure)
return line,
animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=30)
In the solution, 3 waves are visible:
In fact, the solution of any Riemann problem consists of some combination of these three types of waves. In the Euler equations, one of the waves is always a contact discontinuity, but each of the other two waves may be a shock or a rarefaction, depending on the left and right states.
For convenience, all of the code from the cells above to set up and run the shocktube problem is pasted together below. Play around with the code. You might:
%matplotlib inline
import matplotlib.pyplot as plt
from clawpack import pyclaw
from clawpack import riemann
claw = pyclaw.Controller()
claw.tfinal = 0.1
claw.keep_copy = True # Keep solution data in memory for plotting
claw.output_format = None # Don't write solution data to file
claw.num_output_times = 50 # Write 50 output frames
riemann_solver = riemann.euler_with_efix_1D
claw.solver = pyclaw.ClawSolver1D(riemann_solver)
claw.solver.all_bcs = pyclaw.BC.extrap
domain = pyclaw.Domain( (0.,), (1.,), (100,))
x=domain.grid.x.centers # grid cell centers
claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)
gam = 1.4 # ratio of specific heats
claw.solution.problem_data['gamma'] = gam
claw.solution.problem_data['gamma1'] = gam-1.0
rho_left = 1.0; rho_right = 0.125
p_left = 1.0; p_right = 0.1
claw.solution.q[0,:] = (x<0.5)*rho_left + (x>=0.5)*rho_right
claw.solution.q[1,:] = 0.
claw.solution.q[2,:] = ((x<0.5)*p_left + (x>=0.5)*p_right)/(gam-1.0)
status = claw.run()