discontinuity-advection.c

Advection of two fluids at different pressures
- References

Advection of two fluids at different pressures

This problem, initially proposed by Johnsen and Colonius, 2006 tests the solver for the advection of an interface between two different ideal gases at uniform velocity and pressure quantifying the amplitude of spurious pressure and velocity oscillations induced by the method when advecting an interface with different material properties. See also Figure 4 in Fuster and Popinet, 2018.

\displaystyle (\rho,u,p,\gamma)^T_L = (1, 0.5, 1/1.4, 1.2)^T \displaystyle (\rho,u,p,\gamma)^T_R = (10, 0.5, 1/1.4, 1.4)^T

By construction the method should keep the pressure and velocity uniform. During the advection step the energy is advected avoiding any diffusion and therefore both the provisional pressure and velocity are uniform and do not need to be corrected during the projection step.

#include "grid/multigrid1D.h"
#include "compressible/two-phase.h"
#include "compressible/Mie-Gruneisen.h"

We set the problem parameters, size domain and boundary conditions.

double rhoL = 1., rhoR = 10.;
double pL, pR;
double tend = 8;

int main()
{
  L0 = 2. [1];
  X0 = - L0/2.;
  
  periodic (right);
  
  pL = pR = 1./1.4;
  gamma1 = 1.2;
  gamma2 = 1.4;
 
  N = 128;
  run();
}

The initial conditions are:

event init (i = 0)
{  
  double u0 = 0.5;
  foreach() {
    f[] = (x < 0.);
    p[] = f[]*pL + (1. - f[])*pR;
    frho1[] = f[]*rhoL;
    frho2[] = (1. - f[])*rhoR;
    fE1[] = f[]*pL/(gamma1 - 1.) +  0.5*frho1[]*sq(u0);
    fE2[] = (1. - f[])*pR/(gamma2 - 1.) +  0.5*frho2[]*sq(u0);
    q.x[] = (frho1[] + frho2[])*u0;
  }
}

We output the field variables at the end of the simulation.

event outputdata (t = tend)
{  
  scalar perr[], uerr[];
  
  foreach () {
    perr[] = fabs(p[] - 1./1.4);
    uerr[] = fabs(q.x[]/rho[] - 0.5);

    double Ek = 0.;
    foreach_dimension()
      Ek += sq(q.x[]);

    fprintf (stderr, "%g %g %g %g %g %g %g \n", x, t, p[], rho[], q.x[]/rho[], f[],
	     (fE1[] + fE2[] - 0.5*Ek/rho[])/(f[]/(gamma1 - 1.) + (1. - f[])/(gamma2 - 1.)));
  }

  stats sp = statsf(perr), su = statsf(uerr);
  fprintf (stdout, "%g %g\n", sp.sum/sp.volume, su.sum/su.volume);
  fflush (stdout);
  assert (sp.sum/sp.volume < 2e-9 && su.sum/su.volume < 2e-9);
}

The results below are those displayed in Figure 4 of Fuster and Popinet, 2018.

set xlabel 'x'
set ylabel 'p'
unset key
p "log" u 1:3 w lp

Pressure profile (script)

set xlabel 'x'
set ylabel '{/Symbol r}'
p "log" u 1:4 w lp

Density profile (script)

set xlabel 'x'
set ylabel 'u'
p "log" u 1:5 w lp

Velocity profile (script)

set xlabel 'x'
set ylabel 'VOF function'
p "log" u 1:6 w lp

VOF profile (script)

References

[fuster2018]	Daniel Fuster and Stéphane Popinet. An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension. Journal of Computational Physics, 374:752–768, December 2018. [ DOI \| http \| .pdf ]
[johnsen2006]	Eric Johnsen and Tim Colonius. Implementation of WENO schemes in compressible multicomponent flow problems. Journal of Computational Physics, 219(2):715–732, 2006. [ DOI \| http ]