Skip to content

A Flow with Shocks#

Compressible Euler equations should be discretized with Finite Volumes or FEM with flux up-winding scheme but these are not implemented in FreeFem++. Nevertheless acceptable results can be obtained with the method of characteristics provided that the mean values \bar f=\frac12(f^++f^-) are used at shocks in the scheme, and finally mesh adaptation.

\begin{eqnarray}\label{euler} \p_t\rho+\bar u\n\rho + \bar\rho\n\cdot u &=& 0\nonumber\\ \bar\rho( \p_t u+\frac{\overline{\rho u}}{\bar\rho}\n u +\n p &=& 0\nonumber\\ \p_t p + \bar u\n p +(\gamma-1)\bar p\n\cdot u &=& 0\\ \end{eqnarray}

One possibility is to couple u,p and then update \rho, i.e.

\begin{eqnarray}\label{eulalgo} \frac 1{(\gamma-1)\delta t\bar p^m} (p^{m+1}-p^m \circ X^m) + \n\cdot u^{m+1} &=& 0\nonumber\\ \frac{\bar\rho^m}{\delta t}(u^{m+1}-u^m \circ {\tilde X}^m ) +\n p^{m+1} &=& 0\nonumber\\ \rho^{m+1} = \rho^m \circ X^m + \frac{\bar\rho^m}{(\gamma-1)\bar p^m}(p^{m+1}-p^m &\circ& X^m) \end{eqnarray}

A numerical result is given on figure 1 and the FreeFem++ script is

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
// Parameters
verbosity = 1;
int anew = 1;
int m = 5;
real x0 = 0.5;
real y0 = 0.;
real rr = 0.2;
real dt = 0.01;
real u0 = 2.;
real err0 = 0.00625;
real pena = 2.;

// Mesh
border ccc(t=0, 2){x=2-t; y=1;};
border ddd(t=0, 1){x=0; y=1-t;};
border aaa1(t=0, x0-rr){x=t; y=0;};
border cercle(t=pi, 0){x=x0+rr*cos(t); y=y0+rr*sin(t);}
border aaa2(t=x0+rr, 2){x=t; y=0;};
border bbb(t=0, 1){x=2; y=t;};

 mesh Th;
if(anew)
    Th = buildmesh (ccc(5*m) + ddd(3*m) + aaa1(2*m) + cercle(5*m) + aaa2(5*m) + bbb(2*m));
else
    Th = readmesh("Th_circle.mesh"); plot(Th);

// fespace
fespace Wh(Th, P1);
Wh u, v;
Wh u1, v1;
Wh uh, vh;

fespace Vh(Th, P1);
Vh r, rh, r1;

// Macro
macro dn(u) (N.x*dx(u)+N.y*dy(u)) //

// Initialization
if(anew){
    u1 = u0;
    v1 = 0;
    r1 = 1;
}
else{
    ifstream g("u.txt"); g >> u1[];
    ifstream gg("v.txt"); gg >> v1[];
    ifstream ggg("r.txt"); ggg >> r1[];
    plot(u1, ps="eta.eps", value=1, wait=1);
    err0 = err0/10;
    dt = dt/10;
}

// Problem
problem euler(u, v, r, uh, vh, rh)
    = int2d(Th)(
          (u*uh + v*vh + r*rh)/dt
        + ((dx(r)*uh + dy(r)*vh) - (dx(rh)*u + dy(rh)*v))
    )
    + int2d(Th)(
        - (
              rh*convect([u1,v1],-dt,r1)
            + uh*convect([u1,v1],-dt,u1)
            + vh*convect([u1,v1],-dt,v1)
        )/dt
    )
    +int1d(Th, 6)(
          rh*u
    )
    + on(2, r=0)
    + on(2, u=u0)
    + on(2, v=0)
    ;

// Iterations
int j = 80;
for(int k = 0; k < 3; k++){
    if(k==20){
        err0 = err0/10;
        dt = dt/10;
        j = 5;
    }

    // Solve
    for(int i = 0; i < j; i++){
        euler;
        u1=u;
        v1=v;
        r1=abs(r);
        cout << "k = " << k << " E = " << int2d(Th)(u^2+v^2+r) << endl;
        plot(r, value=1);
    }

    // Mesh adaptation
    Th = adaptmesh (Th, r, nbvx=40000, err=err0, abserror=1, nbjacoby=2, omega=1.8, ratio=1.8, nbsmooth=3, splitpbedge=1, maxsubdiv=5, rescaling=1);
    plot(Th);
    u = u;
    v = v;
    r = r;

    // Save
    savemesh(Th, "Th_circle.mesh");
    ofstream f("u.txt"); f << u[];
    ofstream ff("v.txt"); ff << v[];
    ofstream fff("r.txt"); fff << r[];
    r1 = sqrt(u*u+v*v);
    plot(r1, ps="mach.eps", value=1);
    r1 = r;
}
Fig. 1: Pressure for a Euler flow around a disk at Mach 2 computed by \eqref{eulalgo}
Mach 2