Mie-Gruneisen.h

Go to the documentation of this file.

/** @file Mie-Gruneisen.h
 */
/**
# The Mie--Gruneisen Equation of State
 
This EOS is typically used in combination with the
[two-phase compressible solver](/src/compressible/two-phase.h).
 
The general form of the
[Mie--Gruneisen](https://en.wikipedia.org/wiki/Mie%E2%80%93Gr%C3%BCneisen_equation_of_state)
EOS can be written
\f[
\rho_i e_i = \frac{p_i + \Gamma_i \Pi_i}{\Gamma_i - 1}
\f]
with \f$\rho_i\f$, \f$e_i\f$ and \f$p_i\f$ the densities, internal energies and
pressures of each phase.
 
These are the coefficients of the Mie-Gruneisen EOS for each phase. */
 
double gamma1 = 1.4 [0], gamma2 = 1.4 [0], PI1 = 0., PI2 = 0.;
 
/**
## Sound speed
 
In mixture cells, this function returns the maximum between the speeds
in both phases. */
 
double sound_speed (Point point)
{
  double fc = clamp (f[],0.,1.);
  double c2speed1 = 0., c2speed2 = 0.;
 
  double Ek = 0.;
  for (int _d = 0; _d < dimension; _d++)
    Ek += sq(q.x[]);
  Ek /= 2.*(frho1[] + frho2[]);
  
  if (fc > 0.00001) {
    double fe1 = fE1[] - fc*Ek;
    double p  = fe1/fc*(gamma1 - 1.) - gamma1*PI1;
    c2speed1 = fc*gamma1*(p + PI1)/frho1[];
  }
  
  if (fc < 0.99999) {
    double fe2 = fE2[] - (1. - fc)*Ek;
    double p  = fe2/(1. - fc)*(gamma2 - 1.) - gamma2*PI2;
    c2speed2 = (1. - fc)*gamma2*(p + PI2)/frho2[];
  }
 
  return sqrt (max (c2speed1, c2speed2));
}
 
/**
## Average pressure */
 
#define PIGAMMA  double invgammaavg = fc/(gamma1 - 1.) + (1. - fc)/(gamma2 - 1.), \
    PIGAMMAavg = fc*PI1*gamma1/(gamma1 - 1.) + (1. - fc)*PI2*gamma2/(gamma2 - 1.)
 
double average_pressure (Point point)
{
  double fc = clamp (f[],0.,1.);
  PIGAMMA;
  double Ek = 0.;
  for (int _d = 0; _d < dimension; _d++)
    Ek += sq(q.x[]);
  Ek /= 2.*(frho1[] + frho2[]);
  return (fE1[] + fE2[] - Ek - PIGAMMAavg)/invgammaavg;
}
 
/**
## Bulk compressibility of the mixture
 
i.e. \f$\rho c^2\f$. */
 
double bulk_compressibility (Point point)
{
  double fc = clamp (f[],0.,1.);
  PIGAMMA;
  return (p[]*(invgammaavg + 1.) + PIGAMMAavg)/invgammaavg;
}
 
/**
## Internal energy */
 
double internal_energy (Point point, double fc)
{
  PIGAMMA;
  return p[]*invgammaavg + PIGAMMAavg;
}