An “all Mach” flow solver

We wish to solve the generic momentum equation tq+(qu)=p+(μu)+ρa with q=ρu the momentum, u the velocity, ρ the density, μ the dynamic viscosity, p the pressure and a an acceleration. The pressure is defined through an equation of state and verifies the evolution equation tp+up=ρc2u with c the speed of sound. By default the solver sets c=, ρ=1 and the pressure equation reduces to u=0

The advection of momentum is not performed by this solver (so that different schemes can be used) i.e. in the end, by default, we solve the incompressible (linearised) Euler equations with a projection method.

We build the solver using the generic time loop and the implicit viscous solver (which includes the multigrid Poisson–Helmholtz solver).

#include "run.h"
#include "timestep.h"
#include "viscosity.h"

The primitive variables are the momentum q, pressure p, density ρ, (face) specific volume α=1/ρ, (face) dynamic viscosity μ (which is zero by default) and (face) velocity field uf.

vector q[];
scalar p[];
face vector uf[];
(const) face vector α = unityf, μ = zerof;
(const) scalar ρ = unity;

The equation of state is defined by the pressure field ps and ρc2. In the incompressible limit ρc2. Rather than trying to compute this limit, we set both fields to zero and check for this particular case when computing the pressure (see below). This means that by default the fluid is incompressible.

scalar ps[];
(const) scalar rhoc2 = zeroc;
(const) face vector a = zerof;

We store the combined pressure gradient and acceleration field in g.

vector g[];

event defaults (i = 0) {

The default density field is set to unity (times the metric).

  if (α.x.i == unityf.x.i)
    α = fm;

We apply boundary conditions and define the face velocity field after user initialisation.

event init (i = 0) {
  boundary ({q,p,ρ});

The face velocity field is obtained by simple linear interpolation from the momentum field. We make sure that the specific volume α is defined by calling the “properties” event (see below).

  event ("properties");
    uf.x[] = α.x[]*(q.x[] + q.x[-1])/2.;
  boundary ((scalar *){uf});

The timestep is computed using the CFL condition on the face velocity field.

double dtmax;

event set_dtmax (i++,last) dtmax = DT;

event stability (i++,last) {
  dt = dtnext (timestep (uf, dtmax));

Tracers (including momentum q) are advected by these events.

event vof (i++,last);
event tracer_advection (i++,last);

The equation of state (i.e. fields α, ρ, ρc2 and ps) is defined by this event.

event properties (i++,last)
  boundary ({ρ, rhoc2, ps, α, μ});

If the acceleration is not constant, we reset it to zero.

  if (!is_constant(a.x)) {
    face vector af = a;
      af.x[] = 0.;

This event can be overloaded to add acceleration terms.

event acceleration (i++, last)
  boundary ((scalar *){a});

The equation for the pressure is a Poisson–Helmoltz problem which we will solve with the multigrid solver. The statistics for the solver will be stored in mgp (resp. mgu for the viscosity solver).

mgstats mgp, mgu;

event pressure (i++, last)

If the viscosity is not zero, we use the implicit viscosity solver to obtain the velocity field at time t+Δt. The pressure term is taken into account using the pressure gradient at time t. A provisionary momentum (without the pressure gradient) is then built from the velocity field.

  if (constant(μ.x) != 0.) {
        q.x[] = (q.x[] + dt*g.x[])/ρ[];
    boundary ((scalar *){q});
    mgu = viscosity (q, μ, ρ, dt, mgu.nrelax);
        q.x[] = q.x[]*ρ[] - dt*g.x[];
    boundary ((scalar *){q});

We first define a temporary face velocity field u using simple averaging from q, αn+1 and the acceleration term.

    uf.x[] = α.x[]*(q.x[] + q.x[-1])/2. + dt*fm.x[]*a.x[];
  boundary ((scalar *){uf});

The evolution equation for the pressure is tp+up=ρc2u with ρ the density and c the speed of sound. Following the classical projection method for incompressible flows, we set un+1=uΔt(αp)n+1 The evolution equation for the pressure can then be discretised as pn+1pnΔt+unpn=ρcn+12un+1 which gives, after some manipulations, the Poisson–Helmholtz equation λn+1pn+1+(αp)n+1=λn+1p+1Δtu with p=pnΔtunpn and λ=1Δt2ρc2

  scalar λ = rhoc2, rhs = ps;
  foreach() {

We compute λ and the corresponding term in the right-hand-side of the Poisson–Helmholtz equation.

    if (constant(λ) == 0.)
      rhs[] = 0.;
    else {
      λ[] = - cm[]/(sq(dt)*rhoc2[]);
      rhs[] = λ[]*ps[];

We add the divergence of the velocity field to the right-hand-side.

    double div = 0.;
      div += uf.x[1] - uf.x[];
    rhs[] += div/(dt*Δ);
  boundary ({λ, rhs});

The Poisson–Helmholtz solver is called with a definition of the tolerance identical to that used for incompressible flows.

  mgp = poisson (p, rhs, α, λ, tolerance = TOLERANCE/sq(dt));

The pressure gradient is applied to u to obtain the face velocity field at time n+1.

We also compute the combined face pressure gradient and acceleration field gf.

  face vector gf[];
  foreach_face() {
    double dp = α.x[]*(p[] - p[-1])/Δ;
    uf.x[] -= dt*dp;
    gf.x[] = a.x[] - dp/fm.x[];
  boundary_flux ({gf});

And finally we apply the pressure gradient/acceleration term to the flux/momentum. We also store the centered, combined pressure gradient and acceleration field g.

    foreach_dimension() {
      g.x[] = ρ[]*(gf.x[] + gf.x[1])/2.;
      q.x[] += dt*g.x[];
  boundary ((scalar *){q,g,uf});

After mesh adaptation fluid properties need to be updated.

#if TREE
event adapt (i++,last) {
  event ("properties");