Compressible Richtmyer-Meshkov instability with surface tension

We use the compressible solver with VOF interface tracking and surface tension developed in Fuster and Popinet (2018).

#include "grid/multigrid.h"
#include "two-phase-compressible.h"
#include "compressible-tension.h"

Set length scale and resolution for the problem.

#define D 1.0
#define LEVEL 9

Define Mach number and Weber numbers.

#define M 1.2
#define We 6

Some useful reference quantities.

double Ms = M;
double rhoR = 0.1;
double rhoW = 1;
double pR = 1.;
double tend = 4;
double gam = 1.4;

Calculate the postshock conditions based on the Mach number.

double prat(double Ms, double gam){
return 2.*gam*sq(Ms)/(gam+1.) - (gam-1.)/(gam+1.);
}

double rhorat(double Ms, double gam){
return (gam + 1.)*sq(Ms)/((gam - 1.)*sq(Ms) + 2.);
}

double urat(double Ms, double gam){
return 1./rhorat(Ms, gam);
}

double us(double Ms, double gam, double rhoR, double pR){
return Ms*sqrt(gam*pR/rhoR);
}

Some useful fields.

FILE * fp = NULL;

double rho1, rho2, A, dv;
double p0R = 1.;

double rhoL, uL, pL, ushock, qL, EL, eta0, l, x0;
scalar denslap[];
scalar schl[];
scalar totdens[];
scalar vort[];

CFL number needs to be reduced to 0.25 for obtaining physical results at low Mach numbers.

double CFLac = 0.25;

Make sure we meet the CFL condition by correctly estimating maximum signal speeds (|u|+|c|).

event stability (i++) {
stats pstats = statsf(p);
stats fr1stats = statsf(frho1);
stats qstats = statsf(q.x);

double cson = fmax(sqrt(gamma1*(pstats.max)/rhoL), sqrt(gamma2*(pstats.max + PI2)/rhoW));
double u_c = fmax(fabs(qstats.max/rhoL), fabs(qstats.max/rhoW));

foreach () {
DT    = L0*CFLac/(cson + u_c)/pow(2,LEVEL);
dtmax = L0*CFLac/(cson + u_c)/pow(2,LEVEL);
}
}

Boundary conditions.

p[left] = dirichlet(pL);
frho1[left] = dirichlet(rhoL);
q.n[left] = dirichlet(qL);
q.t[left] = dirichlet(0.);
fE1[left] = dirichlet(EL);
q.n[bottom] = neumann(0.);
q.n[top] = neumann(0.);
p[top] = neumann(0.);
p[bottom] = neumann(0.);
frho1[top] = neumann(0.);
frho1[bottom] = neumann(0.);

Main routine, set all the relevant parameters and domain size, etc.

int main() {
gamma1 = gam;
gamma2 = gam;
PI2 = 0.;
l=D;

eta0= 0.01*l;
rhoL = rhorat(Ms, gam)*rhoR;
pL = pR*prat(Ms, gam);
ushock = us(Ms, gam, rhoR, pR);
uL = ushock*(1.-urat(Ms, gam));
qL = uL*rhoL;
EL = pL/(gam-1.) + 0.5*sq(qL)/rhoL;
fprintf(ferr, "%g %g %g %g\n", pL, rhoL, qL, EL);
periodic(top);

Size of the domain:

  size (11.*D);
x0=-eta0-L0/324; origin (x0, -L0/11.);

The density is variable.

  rho1 = rhoR;
rho2 = rhoW;
A=(rho2-rho1)/(rho2+rho1);

Set viscosity (mu1, mu2) and surface tension (f.sigma) via pre-shock Weber numbers.

  dv = Ms*sqrt(gam*pR/rhoR)-sqrt(gam*pR*prat(Ms,gam)/(rhoR*rhorat(Ms,gam)))*sqrt(((gam-1)*sq(Ms)+2)/(2*gam*sq(Ms)-(gam-1)));
fprintf(ferr, "dv=%g\n", dv);
f.sigma = (rho1+rho2)*sq(A*dv)*eta0/We;

mu1=0; mu2=0;

N = (1 << LEVEL)*11;
fprintf(ferr, "uL=%g %g\n", uL, ushock);

We open a file indexed by the level to store the time evolution of the kinetic energy.

  char name[80];
sprintf (name, "k-%d", LEVEL);
fp = fopen (name, "w");
run();
fclose (fp);

We use grep to filter the lines generated by gnuplot containing the results of the fits (see below).

  system ("grep ^fit out >> log");
}

event init (i = 0) {
fprintf(ferr, "start\n");
fprintf(ferr, "l=%g\n", l); fprintf(ferr, "L0=%g\n", L0); fprintf(ferr, "sigma=%g\n", f.sigma);
fprintf(ferr, "mu=%g\n", mu1);

We initialise the shape of the interface.

  fraction (f, -x+eta0*cos(2*pi*y/l) );
boundary({f});

foreach() {
double pLap = p0R;
double m = (x < -7.);
p[] = f[]*(m*pL + (1. - m)*pR) + (1.-f[])*pLap;
frho1[] = f[]*(m*rhoL + (1. - m)*rhoR);
frho2[] = (1. - f[])*rhoW;
q.x[] = m*frho1[]*uL;
q.y[] = 0.;
fE1[] = f[]*(p[]/(gam - 1.) + 0.5*pow(q.x[],2)/rhoL);
fE2[] = (1.-f[])*(pLap + PI2*gamma2)/(gamma2-1.);
}
boundary ((scalar *){q.x, q.y, frho1, frho2, p, fE1, fE2});
}

At each timestep we output the kinetic energy.

//event logfile (i++; t <= 3.0) {
event logfile (i++; t <= 12.0) {
stats pstats = statsf(p);
double ke = 0.;
foreach (reduction(+:ke))
ke += sq(Delta)*(sq(q.x[]) + sq(q.y[]))/(frho1[]+frho2[]);
fprintf (fp, "%g %g %g %d\n", t, ke, mgp.i, pstats.max);
fflush (fp);
}

event outputInterface(t += 0.02) {
char resultname[100];
sprintf( resultname, "results%4.2f_%d.dat", t, pid() );
FILE * fp = fopen(resultname, "w");
scalar xpos[];
scalar ypos[];
position (f, xpos, {1, 0});
position (f, ypos, {0, 1});
foreach()
if (xpos[] != nodata)
fprintf (fp, "%g %g\n", xpos[], ypos[]);
fclose (fp);
}

event output (t = 0.00; t += 0.01){
srand(time(NULL));
scalar vort[];
vector uu[];
dump("dump");
fprintf(ferr, "t=%g\n", t);

foreach()
{
totdens[] = frho1[] + frho2[];
vort[] = (q.y[1] - q.y[-1])/(2.*Delta) - (q.x[0,1]-q.x[0,-1])/(2.*Delta);
foreach_dimension() {
}
denslap[] = -fabs(totdens[1] + totdens[0,1] + totdens[-1] + totdens[0,-1] - 4.*totdens[])/sq(Delta);
}

foreach()
{
double kscale = 40.;
}

static FILE * fp = fopen("p.ppm", "w");
static FILE * fuf = fopen("f.ppm", "w");
static FILE * fr1 = fopen("frho1.ppm", "w");
static FILE * fr2 = fopen("frho2.ppm", "w");
static FILE * fqx = fopen("fqx.ppm", "w");
static FILE * fe1 = fopen("fe1.ppm", "w");
static FILE * fe2 = fopen("fe2.ppm", "w");
static FILE * fdl = fopen("fdl.ppm", "w");
static FILE * fsc = fopen("fsc.ppm", "w");
static FILE * fom = fopen("fom.ppm", "w");
output_ppm (p, fp, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (f, fuf, n=512*11, min=0., max=1., box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (frho1, fr1, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (frho2, fr2, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (q.x, fqx, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (fE1, fe1, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (fE2, fe2, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (denslap, fdl, n=512*11, map=gray, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (schl, fsc, n=512*11, map=gray, box = {{x0,-L0/11},{x0+L0,0}});
output_ppm (vort, fom, n=512*11, box = {{x0,-L0/11},{x0+L0,0}});
}

event end (t==tend) {
}

#if TREE
#endif //TREE