src/compressible.h

Compressible gas dynamics

The Euler system of conservation laws for a compressible gas can be written

t(ρEwxwy)+x(wxwxρ(E+p)wx2ρ+pwywxρ)+y(wywyρ(E+p)wywxρwy2ρ+p)=0

with ρ the gas density, E the total energy, w the gas momentum and p the pressure given by the equation of state

p=(γ1)(Eρu2/2)

with γ the polytropic exponent. This system can be solved using the generic solver for systems of conservation laws.

#include "conservation.h"

The conserved scalars are the gas density ρ and the total energy E. The only conserved vector is the momentum w. The constant γ is represented by gammao here, with a default value of 1.4.

scalar ρ[], E[];
vector w[];
scalar * scalars = {ρ, E};
vector * vectors = {w};
double gammao = 1.4 ;

The system is entirely defined by the flux() function called by the generic solver for conservation laws. The parameter passed to the function is the array s which contains the state variables for each conserved field, in the order of their definition above (i.e. scalars then vectors).

void flux (const double * s, double * f, double e[2])
{

We first recover each value (ρ, E, wx and wy) and then compute the corresponding fluxes (f[0], f[1], f[2] and f[3]).

  double ρ = s[0], E = s[1], wn = s[2], w2 = 0.;
  for (int i = 2; i < 2 + dimension; i++)
    w2 += sq(s[i]);
  double un = wn/ρ, p = (gammao - 1.)*(E - 0.5*w2/ρ);

  f[0] = wn;
  f[1] = un*(E + p);
  f[2] = un*wn + p;
  for (int i = 3; i < 2 + dimension; i++)
    f[i] = un*s[i];

The minimum and maximum eigenvalues for the Euler system are the characteristic speeds u±(γp/ρ).

  double c = sqrt(gammao*p/ρ);
  e[0] = un - c; // min
  e[1] = un + c; // max
}

Usage

Tests