-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // shock.edp 2 : // 3 : // Discussion: 4 : // 5 : // Supersonic perfect gas flow around an ellipse. 6 : // 7 : // Be careful : nonconservative formulation of the compressible 8 : // Euler equations. 9 : // 10 : // A conservative Finite Volume method should be used. 11 : // 12 : // Location: 13 : // 14 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/shock/shock.edp 15 : // 16 : // Modified: 17 : // 18 : // 20 May 2015 19 : // 20 : // Author: 21 : // 22 : // Florian De Vuyst 23 : // 24 : // Reference: 25 : // 26 : // Florian De Vuyst, 27 : // Numerical modeling of transport problems using freefem++ software - 28 : // with examples in biology, CFD, traffic flow and energy transfer, 29 : // HAL id: cel-00842234 30 : // https://cel.archives-ouvertes.fr/cel-00842234 31 : // 32 : // Parameters: 33 : // 34 : // real ALPHA, the angle of attack, in radians. 35 : // 36 : // real CINF, the speed of sound. 37 : // 38 : // real DT, time step. 39 : // 40 : // real RA, major radius of the ellipse. 41 : // 42 : // real RADIUS, the radius of the infinite flow boundary. 43 : // 44 : // real RB, minor radius of the ellipse. 45 : // 46 : // real UINF, the infinite flow Mach number is 2. 47 : // 48 : // real X0, X translation of the body. 49 : // 50 : cout << "\n"; 51 : cout << "shock\n"; 52 : cout << " FreeFem++ version\n"; 53 : cout << " Solve the Euler equations for perfect gas flow\n"; 54 : cout << " around an elliptical obstacle.\n"; 55 : 56 : real gamma = 1.4; 57 : real pinf = 100.0; 58 : real rhoinf = 0.3; 59 : real cinf = sqrt ( gamma * pinf / rhoinf ); 60 : real uinf = 1.5 * cinf; 61 : real radius = 10.0; 62 : real alpha = 0.3; 63 : real x0 = -4.0; 64 : real ra = 2.0; 65 : real rb = 0.3; 66 : real dt = 0.01; 67 : // 68 : // External "infinite" boundary is a circle 69 : // 70 : border binf (t=0,2*pi) {x=radius*cos(t); y=radius*sin(t);} 71 : // 72 : // Obstacle is an ellipse 73 : // 74 : border wall (t=0,2*pi) { x = x0+cos(alpha)*ra*cos(t)+sin(alpha)*rb*sin(t); 75 : y = -sin(alpha)*ra*cos(t)+cos(alpha)*rb*sin(t);} 76 : // 77 : // Define the mesh. 78 : // 79 : mesh Th = buildmesh ( binf ( 100 ) + wall ( -80 ) ); 80 : plot ( Th, wait = 1, ps = "shock_mesh.ps" ); 81 : // 82 : // Define the finite element space. 83 : // 84 : fespace Vh ( Th, P1 ); 85 : Vh ap; 86 : Vh arho; 87 : Vh p; 88 : Vh rho; 89 : Vh T; 90 : Vh u1; 91 : Vh u2; 92 : Vh rhoold, pold, Told, arhoold, apold, u1old, u2old; 93 : // 94 : // Define test functions. 95 : // 96 : Vh v; 97 : Vh v1; 98 : Vh v2; 99 : Vh w; 100 : // 101 : // Field initialization 102 : // 103 : rho = rhoinf; 104 : p = pinf; 105 : u1 = uinf; 106 : u2 = 0; 107 : arho = log ( rho ); 108 : ap = log ( p ); 109 : T = p / rho; // temperature−like 110 : 111 : plot ( T, nbiso=50, fill=1, value=1, wait=0 ); 112 : // 113 : // Define the weak form of the Euler equations. 114 : // 115 : problem euler ( [arho, u1, u2, ap], [v, v1, v2, w] ) = 116 : int2d ( Th ) ( arho*v/dt ) 117 : - int2d ( Th ) ( convect([u1old,u2old], -dt, arhoold)*v/dt ) 118 : + int1d ( Th, binf) ( (uinf*N.x)*v ) 119 : - int2d ( Th ) ( u1*dx(v)+u2*dy(v) ) 120 : + int2d ( Th ) ( u1*v1/Told/dt ) 121 : - int2d ( Th ) ( convect([u1old,u2old], -dt, u1old)*v1/Told/dt ) 122 : + int2d ( Th ) ( u2*v2/Told/dt ) 123 : - int2d ( Th ) ( convect([u1old,u2old], -dt, u2old)*v2/Told/dt ) 124 : + int2d ( Th ) ( dx(ap)*v1+dy(ap)*v2 ) 125 : + int2d ( Th ) ( ap*w/dt ) 126 : - int2d ( Th ) ( convect([u1old,u2old], -dt, apold)*w/dt ) 127 : + int1d ( Th, binf) ( gamma*(uinf*N.x)*w ) 128 : - int2d ( Th ) ( gamma*u1*dx(w)+gamma*u2*dy(w) ); 129 : 130 : for ( int it = 0; it < 80; it++ ) 131 : { 132 : u1old = u1; 133 : u2old = u2; 134 : arhoold = arho; 135 : apold = ap; 136 : Told = T; 137 : euler; 138 : if ( 70 < it ) 139 : { 140 : Th = adaptmesh ( Th, ap ); 141 : u1 = u1; 142 : u2 = u2; 143 : arho = arho; 144 : ap = ap; 145 : } 146 : rho = exp ( arho ); 147 : p = exp ( ap ); 148 : T = p / rho; 149 : plot ( rho, nbiso=50, fill=0, value=1 ); 150 : if ( it < 20 ) 151 : { 152 : plot ( rho, nbiso=50, fill=1, value=1 ); 153 : } 154 : } 155 : // 156 : // Plot final results. 157 : // 158 : plot ( T, nbiso = 50, fill = 1, value = 1, ps = "shock_t.ps" ); 159 : plot ( T, nbiso = 50, fill = 0, value = 1, ps = "shock_tiso.ps" ); 160 : plot ( p, nbiso = 50, fill = 1, value = 1, ps = "shock_p.ps" ); 161 : plot ( rho, nbiso = 50, fill = 1, value = 1, ps = "shock_rho.ps" ); 162 : plot ( [ u1, u2 ], ps = "shock_velocity.ps" ); 163 : plot ( Th, ps = "shock_mesh2.ps" ); 164 : // 165 : // Terminate. 166 : // 167 : cout << "\n"; 168 : cout << "shock:\n"; 169 : cout << " Normal end of execution.\n"; 170 : 171 : sizestack + 1024 =9288 ( 8264 ) shock FreeFem++ version Solve the Euler equations for perfect gas flow around an elliptical obstacle. -- mesh: Nb of Triangles = 7696, Nb of Vertices 3938 kk 4 3 : -- Solve : min -1.94322 max -0.461697 min -1.65075 max 40.4072 min -15.8547 max 25.6546 min 3.57022 max 5.64436 -- Solve : min -1.92602 max -0.489338 min -2.2905 max 41.3387 min -16.9227 max 24.461 min 3.59431 max 5.60566 -- Solve : min -1.92442 max -0.47489 min -3.50496 max 43.7808 min -17.6832 max 23.1432 min 3.59654 max 5.62589 -- Solve : min -1.9363 max -0.474641 min -4.56684 max 45.7317 min -18.4157 max 23.049 min 3.57991 max 5.62623 -- Solve : min -1.95137 max -0.467409 min -5.81349 max 43.861 min -18.6386 max 21.34 min 3.55882 max 5.63636 -- Solve : min -1.95829 max -0.462065 min -14.9283 max 44.0594 min -18.873 max 20.281 min 3.54912 max 5.64384 -- Solve : min -1.95477 max -0.455095 min -8.67772 max 43.3352 min -19.3384 max 19.6745 min 3.55406 max 5.6536 -- Solve : min -1.95748 max -0.449378 min -13.9185 max 44.6726 min -19.5762 max 19.5838 min 3.55027 max 5.6616 -- Solve : min -1.97149 max -0.445905 min -13.5785 max 41.8976 min -19.7899 max 19.3415 min 3.53065 max 5.66646 -- Solve : min -2.0027 max -0.443888 min -12.4393 max 41.4741 min -19.9955 max 19.0481 min 3.48695 max 5.66929 -- Solve : min -2.0231 max -0.441861 min -11.5935 max 41.6094 min -19.8474 max 18.9095 min 3.45839 max 5.67213 -- Solve : min -2.03625 max -0.439793 min -12.3424 max 41.4361 min -19.6252 max 18.9281 min 3.43998 max 5.67502 -- Solve : min -2.0566 max -0.437583 min -12.363 max 40.8117 min -17.8042 max 18.9458 min 3.41149 max 5.67812 -- Solve : min -2.08282 max -0.435059 min -12.8215 max 40.3998 min -17.2088 max 18.9959 min 3.37478 max 5.68165 -- Solve : min -2.11227 max -0.432277 min -10.3419 max 40.4692 min -17.1259 max 19.0594 min 3.33356 max 5.68554 -- Solve : min -2.14377 max -0.429912 min -10.4873 max 40.164 min -17.2645 max 19.1037 min 3.28946 max 5.68886 -- Solve : min -2.17503 max -0.427728 min -12.8473 max 40.0331 min -17.3337 max 19.1584 min 3.24568 max 5.69191 -- Solve : min -2.20156 max -0.425666 min -18.2443 max 40.7546 min -17.2194 max 19.1895 min 3.20855 max 5.6948 -- Solve : min -2.21818 max -0.423854 min -21.7451 max 40.2645 min -21.0628 max 19.1986 min 3.18527 max 5.69734 -- Solve : min -2.22733 max -0.422112 min -15.3222 max 40.7501 min -16.7867 max 19.2111 min 3.17246 max 5.69977 -- Solve : min -2.24349 max -0.42051 min -15.1479 max 39.8697 min -16.8129 max 19.2951 min 3.14984 max 5.70202 -- Solve : min -2.25754 max -0.419062 min -14.2148 max 40.6421 min -19.7009 max 19.2729 min 3.13018 max 5.70404 -- Solve : min -2.26342 max -0.41771 min -12.8417 max 40.8917 min -17.8406 max 19.3049 min 3.12194 max 5.70594 -- Solve : min -2.27495 max -0.416487 min -12.7073 max 40.071 min -18.2113 max 19.3449 min 3.1058 max 5.70765 -- Solve : min -2.26036 max -0.415308 min -17.0586 max 41.7006 min -16.5886 max 19.3983 min 3.12622 max 5.7093 -- Solve : min -2.24659 max -0.414216 min -21.2342 max 40.3591 min -16.3101 max 20.605 min 3.14551 max 5.71083 -- Solve : min -2.24592 max -0.413236 min -16.2356 max 40.4301 min -16.1681 max 25.5759 min 3.14645 max 5.7122 -- Solve : min -2.24608 max -0.412266 min -15.7719 max 40.2183 min -16.1328 max 19.686 min 3.14621 max 5.71356 -- Solve : min -2.24597 max -0.41133 min -13.5929 max 39.9995 min -16.2192 max 19.5779 min 3.14638 max 5.71487 -- Solve : min -2.24827 max -0.41048 min -13.7274 max 41.2645 min -16.3449 max 19.6203 min 3.14315 max 5.71606 -- Solve : min -2.25011 max -0.409694 min -14.2935 max 40.9293 min -16.2302 max 19.6587 min 3.14058 max 5.71716 -- Solve : min -2.24219 max -0.408972 min -18.2288 max 41.5468 min -16.2506 max 19.693 min 3.15167 max 5.71817 -- Solve : min -2.24063 max -0.40829 min -18.0022 max 40.7439 min -16.0133 max 21.2827 min 3.15385 max 5.71913 -- Solve : min -2.23211 max -0.407628 min -16.5626 max 41.3553 min -15.979 max 20.9046 min 3.16578 max 5.72005 -- Solve : min -2.23119 max -0.407 min -15.3793 max 39.9384 min -15.9429 max 19.8021 min 3.16707 max 5.72093 -- Solve : min -2.2341 max -0.406426 min -12.3293 max 41.8239 min -15.9222 max 19.8379 min 3.163 max 5.72174 -- Solve : min -2.23085 max -0.405911 min -15.7794 max 42.0974 min -15.8962 max 19.8731 min 3.16754 max 5.72246 -- Solve : min -2.22976 max -0.405444 min -17.0412 max 40.4573 min -15.9262 max 19.9073 min 3.16907 max 5.72311 -- Solve : min -2.22417 max -0.405055 min -17.4416 max 40.8712 min -15.8529 max 19.9405 min 3.1769 max 5.72365 -- Solve : min -2.22235 max -0.404726 min -20.2203 max 39.4434 min -15.8743 max 21.9119 min 3.17944 max 5.72412 -- Solve : min -2.21267 max -0.404441 min -18.2206 max 40.2678 min -15.7654 max 20.0005 min 3.193 max 5.72451 -- Solve : min -2.20409 max -0.404179 min -16.2865 max 39.9737 min -15.7206 max 20.0304 min 3.20501 max 5.72488 -- Solve : min -2.20381 max -0.403946 min -14.78 max 41.0612 min -15.6746 max 20.0585 min 3.2054 max 5.72521 -- Solve : min -2.1929 max -0.403738 min -10.4277 max 43.1228 min -15.5605 max 20.0837 min 3.22068 max 5.7255 -- Solve : min -2.19164 max -0.403548 min -12.1706 max 42.8034 min -15.5721 max 20.1069 min 3.22243 max 5.72576 -- Solve : min -2.18567 max -0.403365 min -13.9448 max 38.9234 min -15.4175 max 20.1302 min 3.2308 max 5.72602 -- Solve : min -2.18049 max -0.403178 min -17.9497 max 40.8135 min -15.368 max 20.1564 min 3.23805 max 5.72628 -- Solve : min -2.16775 max -0.40299 min -14.4183 max 39.7818 min -15.1917 max 20.1857 min 3.25588 max 5.72655 -- Solve : min -2.16376 max -0.402797 min -12.1551 max 40.2004 min -15.1253 max 20.2089 min 3.26147 max 5.72682 -- Solve : min -2.16273 max -0.40261 min -12.2054 max 41.8645 min -15.1085 max 20.2294 min 3.2629 max 5.72708 -- Solve : min -2.16378 max -0.402435 min -16.1389 max 39.089 min -15.0597 max 20.2449 min 3.26144 max 5.72732 -- Solve : min -2.16264 max -0.402272 min -14.5714 max 42.0716 min -15.135 max 20.2551 min 3.26303 max 5.72755 -- Solve : min -2.15997 max -0.402121 min -12.8016 max 38.4784 min -15.0717 max 20.2619 min 3.26677 max 5.72776 -- Solve : min -2.15537 max -0.402004 min -12.9793 max 41.7159 min -15.0545 max 20.2644 min 3.27321 max 5.72793 -- Solve : min -2.15464 max -0.401872 min -13.3004 max 41.6026 min -14.9185 max 20.2671 min 3.27424 max 5.72811 -- Solve : min -2.15035 max -0.401709 min -12.5968 max 40.603 min -14.8543 max 20.2689 min 3.28024 max 5.72834 -- Solve : min -2.15353 max -0.401537 min -12.9333 max 41.3309 min -14.8265 max 20.2732 min 3.27579 max 5.72858 -- Solve : min -2.15121 max -0.401376 min -17.1239 max 39.0505 min -14.8732 max 20.2793 min 3.27904 max 5.72881 -- Solve : min -2.15275 max -0.401222 min -15.5895 max 38.3039 min -14.8457 max 20.2782 min 3.27688 max 5.72902 -- Solve : min -2.14721 max -0.401085 min -12.4941 max 37.7604 min -14.839 max 20.2823 min 3.28463 max 5.72921 -- Solve : min -2.14393 max -0.400923 min -13.2591 max 37.9146 min -14.7863 max 20.2914 min 3.28923 max 5.72944 -- Solve : min -2.14091 max -0.400761 min -12.2381 max 39.2232 min -14.662 max 20.2942 min 3.29346 max 5.72967 -- Solve : min -2.14039 max -0.400604 min -12.2322 max 40.8397 min -14.6609 max 20.2954 min 3.29418 max 5.72989 -- Solve : min -2.14022 max -0.400437 min -15.4515 max 37.8138 min -14.6612 max 20.2951 min 3.29443 max 5.73012 -- Solve : min -2.13893 max -0.400255 min -12.5074 max 38.4535 min -14.6458 max 20.2945 min 3.29623 max 5.73037 -- Solve : min -2.13764 max -0.400059 min -10.6555 max 40.1142 min -14.6327 max 20.2891 min 3.29804 max 5.73065 -- Solve : min -2.13681 max -0.39986 min -13.3836 max 39.2776 min -14.6289 max 20.2835 min 3.2992 max 5.73093 -- Solve : min -2.13506 max -0.399657 min -11.8302 max 38.7599 min -14.6157 max 20.2814 min 3.30165 max 5.73121 -- Solve : min -2.13308 max -0.399451 min -14.8688 max 37.6022 min -14.5966 max 20.2787 min 3.30442 max 5.7315 -- Solve : min -2.1314 max -0.399251 min -14.9747 max 39.06 min -14.5862 max 20.2778 min 3.30677 max 5.73178 -- Solve : min -2.13013 max -0.39905 min -13.9267 max 39.5337 min -14.577 max 20.2773 min 3.30854 max 5.73206 -- Solve : min -2.12935 max -0.398853 min -16.4209 max 40.9988 min -14.5749 max 20.2775 min 3.30965 max 5.73234 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 6132 , h min 0.0072112 , h max 2.79515 area = 310.608 , M area = 2704.23 , M area/( |Khat| nt) 1.01845 infiny-regularity: min 0.331109 max 2.26989 anisomax 10.3456, beta max = 1.4791 min 0.719112 -- mesh: Nb of Triangles = 6132, Nb of Vertices 3142 kk 4 3 : -- Solve : min -2.12772 max -0.395628 min -14.1404 max 42.4638 min -14.9397 max 19.5627 min 3.31192 max 5.73685 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 5923 , h min 0.00448495 , h max 4.26832 area = 310.591 , M area = 2663.78 , M area/( |Khat| nt) 1.03862 infiny-regularity: min 0.287998 max 2.54602 anisomax 16.3895, beta max = 1.59741 min 0.709148 -- mesh: Nb of Triangles = 5923, Nb of Vertices 3041 kk 4 3 : -- Solve : min -2.12421 max -0.395303 min -16.2121 max 42.1828 min -14.745 max 20.1073 min 3.31684 max 5.73731 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 6024 , h min 0.00303431 , h max 4.87794 area = 310.563 , M area = 2654.55 , M area/( |Khat| nt) 1.01767 infiny-regularity: min 0.27462 max 2.42153 anisomax 14.3988, beta max = 1.3632 min 0.72303 -- mesh: Nb of Triangles = 6024, Nb of Vertices 3089 kk 4 3 : -- Solve : min -2.12986 max -0.394174 min -19.9487 max 40.0885 min -14.6139 max 19.0328 min 3.30893 max 5.73889 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 6579 , h min 0.00197598 , h max 4.33502 area = 310.532 , M area = 2948.36 , M area/( |Khat| nt) 1.03495 infiny-regularity: min 0.333367 max 2.8977 anisomax 13.9447, beta max = 1.36392 min 0.666561 -- mesh: Nb of Triangles = 6579, Nb of Vertices 3368 kk 4 3 : -- Solve : min -2.11862 max -0.396297 min -15.5231 max 42.6733 min -14.6724 max 19.2723 min 3.32466 max 5.73592 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 6412 , h min 0.00161394 , h max 4.23317 area = 310.538 , M area = 2854.45 , M area/( |Khat| nt) 1.02808 infiny-regularity: min 0.304634 max 2.67444 anisomax 13.1702, beta max = 1.51733 min 0.700188 -- mesh: Nb of Triangles = 6412, Nb of Vertices 3282 kk 4 3 : -- Solve : min -2.10826 max -0.396642 min -21.1776 max 43.2749 min -14.6271 max 19.4392 min 3.33917 max 5.73543 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 6710 , h min 0.00179391 , h max 4.19834 area = 310.537 , M area = 2997.99 , M area/( |Khat| nt) 1.03183 infiny-regularity: min 0.317366 max 3.09818 anisomax 14.4026, beta max = 1.37087 min 0.671921 -- mesh: Nb of Triangles = 6710, Nb of Vertices 3431 kk 4 3 : -- Solve : min -2.10694 max -0.398155 min -22.0976 max 42.3222 min -14.4717 max 19.7004 min 3.34102 max 5.73331 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 7093 , h min 0.00186884 , h max 4.15223 area = 310.521 , M area = 3224.66 , M area/( |Khat| nt) 1.04991 infiny-regularity: min 0.311096 max 2.73618 anisomax 14.839, beta max = 1.40097 min 0.695603 -- mesh: Nb of Triangles = 7093, Nb of Vertices 3625 kk 4 3 : -- Solve : min -2.09046 max -0.398288 min -18.7216 max 43.1648 min -14.2344 max 19.6292 min 3.36409 max 5.73313 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 7181 , h min 0.00145496 , h max 4.23492 area = 310.524 , M area = 3282.41 , M area/( |Khat| nt) 1.05562 infiny-regularity: min 0.298818 max 2.76759 anisomax 14.7679, beta max = 1.37989 min 0.664673 -- mesh: Nb of Triangles = 7181, Nb of Vertices 3670 kk 4 3 : -- Solve : min -2.0883 max -0.396031 min -20.6548 max 44.2971 min -14.154 max 18.9869 min 3.36711 max 5.73629 number of required edges : 0 -- adaptmesh Regulary: Nb triangles 7719 , h min 0.00154832 , h max 4.06971 area = 310.516 , M area = 3571.16 , M area/( |Khat| nt) 1.06843 infiny-regularity: min 0.338874 max 3.17906 anisomax 12.8808, beta max = 1.43531 min 0.676817 -- mesh: Nb of Triangles = 7719, Nb of Vertices 3938 shock: Normal end of execution. times: compile 0.005176s, execution 260.812s, mpirank:0 CodeAlloc : nb ptr 4125, size :507488 mpirank: 0 Ok: Normal End