FreeFEM Documentation on GitHub

stars - forks

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 \(\displaystyle \bar f=\frac12\left(f^++f^-\right)\) are used at shocks in the scheme, and finally mesh adaptation.

(6)\[\begin{split}\begin{array}{rcl} \partial_t\rho+\bar u\nabla\rho + \bar\rho\nabla\cdot u &=& 0\nonumber\\ \bar\rho( \partial_t u+\frac{\overline{\rho u}}{\bar\rho}\nabla u +\nabla p &=& 0\nonumber\\ \partial_t p + \bar u\nabla p +(\gamma-1)\bar p\nabla\cdot u &=& 0\\ \end{array}\end{split}\]

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

(7)\[\begin{split}\begin{array}{rcl} \frac 1{(\gamma-1)\delta t\bar p^m} (p^{m+1}-p^m \circ X^m) + \nabla\cdot u^{m+1} &=& 0\nonumber\\ \frac{\bar\rho^m}{\delta t}(u^{m+1}-u^m \circ {\tilde X}^m ) +\nabla 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{array}\end{split}\]

A numerical result is given on Fig. 42 and the FreeFEM script is

  1 // Parameters
  2 verbosity = 1;
  3 int anew = 1;
  4 int m = 5;
  5 real x0 = 0.5;
  6 real y0 = 0.;
  7 real rr = 0.2;
  8 real dt = 0.01;
  9 real u0 = 2.;
 10 real err0 = 0.00625;
 11 real pena = 2.;
 12 
 13 // Mesh
 14 border ccc(t=0, 2){x=2-t; y=1;};
 15 border ddd(t=0, 1){x=0; y=1-t;};
 16 border aaa1(t=0, x0-rr){x=t; y=0;};
 17 border cercle(t=pi, 0){x=x0+rr*cos(t); y=y0+rr*sin(t);}
 18 border aaa2(t=x0+rr, 2){x=t; y=0;};
 19 border bbb(t=0, 1){x=2; y=t;};
 20 
 21  mesh Th;
 22 if(anew)
 23    Th = buildmesh (ccc(5*m) + ddd(3*m) + aaa1(2*m) + cercle(5*m) + aaa2(5*m) + bbb(2*m));
 24 else
 25    Th = readmesh("Th_circle.mesh"); plot(Th);
 26 
 27 // fespace
 28 fespace Wh(Th, P1);
 29 Wh u, v;
 30 Wh u1, v1;
 31 Wh uh, vh;
 32 
 33 fespace Vh(Th, P1);
 34 Vh r, rh, r1;
 35 
 36 // Macro
 37 macro dn(u) (N.x*dx(u)+N.y*dy(u)) //
 38 
 39 // Initialization
 40 if(anew){
 41    u1 = u0;
 42    v1 = 0;
 43    r1 = 1;
 44 }
 45 else{
 46    ifstream g("u.txt"); g >> u1[];
 47    ifstream gg("v.txt"); gg >> v1[];
 48    ifstream ggg("r.txt"); ggg >> r1[];
 49    plot(u1, ps="eta.eps", value=1, wait=1);
 50    err0 = err0/10;
 51    dt = dt/10;
 52 }
 53 
 54 // Problem
 55 problem euler(u, v, r, uh, vh, rh)
 56    = int2d(Th)(
 57         (u*uh + v*vh + r*rh)/dt
 58       + ((dx(r)*uh + dy(r)*vh) - (dx(rh)*u + dy(rh)*v))
 59    )
 60    + int2d(Th)(
 61       - (
 62            rh*convect([u1,v1],-dt,r1)
 63          + uh*convect([u1,v1],-dt,u1)
 64          + vh*convect([u1,v1],-dt,v1)
 65       )/dt
 66    )
 67    +int1d(Th, 6)(
 68         rh*u
 69    )
 70    + on(2, r=0)
 71    + on(2, u=u0)
 72    + on(2, v=0)
 73    ;
 74 
 75 // Iterations
 76 int j = 80;
 77 for(int k = 0; k < 3; k++){
 78    if(k==20){
 79       err0 = err0/10;
 80       dt = dt/10;
 81       j = 5;
 82    }
 83 
 84    // Solve
 85    for(int i = 0; i < j; i++){
 86       euler;
 87       u1=u;
 88       v1=v;
 89       r1=abs(r);
 90       cout << "k = " << k << " E = " << int2d(Th)(u^2+v^2+r) << endl;
 91       plot(r, value=1);
 92    }
 93 
 94    // Mesh adaptation
 95    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);
 96    plot(Th);
 97    u = u;
 98    v = v;
 99    r = r;
100 
101    // Save
102    savemesh(Th, "Th_circle.mesh");
103    ofstream f("u.txt"); f << u[];
104    ofstream ff("v.txt"); ff << v[];
105    ofstream fff("r.txt"); fff << r[];
106    r1 = sqrt(u*u+v*v);
107    plot(r1, ps="mach.eps", value=1);
108    r1 = r;
109 }
../_images/mach_2r.png

Fig. 42 Pressure for a Euler flow around a disk at Mach 2 computed by (7)

Table of content