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
  2verbosity = 1;
  3int anew = 1;
  4int m = 5;
  5real x0 = 0.5;
  6real y0 = 0.;
  7real rr = 0.2;
  8real dt = 0.01;
  9real u0 = 2.;
 10real err0 = 0.00625;
 11real pena = 2.;
 12
 13// Mesh
 14border ccc(t=0, 2){x=2-t; y=1;};
 15border ddd(t=0, 1){x=0; y=1-t;};
 16border aaa1(t=0, x0-rr){x=t; y=0;};
 17border cercle(t=pi, 0){x=x0+rr*cos(t); y=y0+rr*sin(t);}
 18border aaa2(t=x0+rr, 2){x=t; y=0;};
 19border bbb(t=0, 1){x=2; y=t;};
 20
 21 mesh Th;
 22if(anew)
 23   Th = buildmesh (ccc(5*m) + ddd(3*m) + aaa1(2*m) + cercle(5*m) + aaa2(5*m) + bbb(2*m));
 24else
 25   Th = readmesh("Th_circle.mesh"); plot(Th);
 26
 27// fespace
 28fespace Wh(Th, P1);
 29Wh u, v;
 30Wh u1, v1;
 31Wh uh, vh;
 32
 33fespace Vh(Th, P1);
 34Vh r, rh, r1;
 35
 36// Macro
 37macro dn(u) (N.x*dx(u)+N.y*dy(u)) //
 38
 39// Initialization
 40if(anew){
 41   u1 = u0;
 42   v1 = 0;
 43   r1 = 1;
 44}
 45else{
 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
 55problem 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
 76int j = 80;
 77for(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