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
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
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}
Table of content