-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae) file : shock.edp 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 : // Licensing: 13 : // 14 : // This code is distributed under the MIT license. 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 = true, 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 = true, value = true, wait = false ); 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 = false, value = true ); 150 : if ( it < 20 ) 151 : { 152 : plot ( rho, nbiso = 50, fill = true, value = true ); 153 : } 154 : } 155 : // 156 : // Plot final results. 157 : // 158 : plot ( T, nbiso = 50, fill = true, value = true, ps = "shock_t.ps" ); 159 : plot ( T, nbiso = 50, fill = false, value = true, ps = "shock_tiso.ps" ); 160 : plot ( p, nbiso = 50, fill = true, value = true, ps = "shock_p.ps" ); 161 : plot ( rho, nbiso = 50, fill = true, value = true, 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 : exit ( 0 ); 172 : 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.9449 max 38.9234 min -15.4175 max 20.1302 min 3.2308 max 5.72602 -- Solve : min -2.18049 max -0.403178 min -17.9496 max 40.8135 min -15.368 max 20.1564 min 3.23805 max 5.72628 -- Solve : min -2.16775 max -0.40299 min -14.4185 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.8642 min -15.1085 max 20.2294 min 3.2629 max 5.72708 -- Solve : min -2.16378 max -0.402435 min -16.1393 max 39.0887 min -15.0597 max 20.2449 min 3.26144 max 5.72732 -- Solve : min -2.16264 max -0.402272 min -14.5697 max 42.0723 min -15.135 max 20.2551 min 3.26303 max 5.72755 -- Solve : min -2.15997 max -0.402121 min -12.8013 max 38.4789 min -15.0717 max 20.2619 min 3.26677 max 5.72776 -- Solve : min -2.15537 max -0.402004 min -12.9789 max 41.715 min -15.0545 max 20.2644 min 3.27321 max 5.72793 -- Solve : min -2.15464 max -0.401872 min -13.3005 max 41.6025 min -14.9185 max 20.2671 min 3.27424 max 5.72811 -- Solve : min -2.15035 max -0.401709 min -12.5969 max 40.6029 min -14.8542 max 20.2689 min 3.28025 max 5.72834 -- Solve : min -2.15353 max -0.401537 min -12.9377 max 41.3305 min -14.8265 max 20.2732 min 3.27579 max 5.72858 -- Solve : min -2.15121 max -0.401376 min -17.1314 max 39.0501 min -14.8732 max 20.2793 min 3.27904 max 5.72881 -- Solve : min -2.15275 max -0.401222 min -15.5698 max 38.3041 min -14.8457 max 20.2782 min 3.27688 max 5.72902 -- Solve : min -2.14721 max -0.401085 min -12.4944 max 37.761 min -14.839 max 20.2823 min 3.28463 max 5.72921 -- Solve : min -2.14393 max -0.400923 min -13.2594 max 37.9127 min -14.7863 max 20.2914 min 3.28923 max 5.72944 -- Solve : min -2.14091 max -0.400761 min -12.2362 max 39.2212 min -14.662 max 20.2942 min 3.29346 max 5.72967 -- Solve : min -2.14039 max -0.400604 min -12.2334 max 40.8251 min -14.6608 max 20.2954 min 3.29419 max 5.72989 -- Solve : min -2.14021 max -0.400437 min -15.4335 max 37.8161 min -14.6612 max 20.2951 min 3.29443 max 5.73012 -- Solve : min -2.13895 max -0.400255 min -12.5035 max 38.41 min -14.646 max 20.2945 min 3.2962 max 5.73037 -- Solve : min -2.13765 max -0.400059 min -10.648 max 40.0823 min -14.6328 max 20.2891 min 3.29803 max 5.73065 -- Solve : min -2.13679 max -0.39986 min -13.4131 max 39.2726 min -14.6287 max 20.2835 min 3.29922 max 5.73093 -- Solve : min -2.13507 max -0.399657 min -11.7879 max 38.7669 min -14.6157 max 20.2814 min 3.30164 max 5.73121 -- Solve : min -2.13308 max -0.399451 min -14.8597 max 37.6409 min -14.5966 max 20.2787 min 3.30442 max 5.7315 -- Solve : min -2.13142 max -0.399251 min -15.0184 max 39.098 min -14.5863 max 20.2778 min 3.30675 max 5.73178 -- Solve : min -2.13013 max -0.39905 min -13.9192 max 39.5624 min -14.5769 max 20.2773 min 3.30855 max 5.73206 -- Solve : min -2.12936 max -0.398853 min -16.3015 max 40.9952 min -14.5753 max 20.2775 min 3.30962 max 5.73234 number of required edges : 0 -- adaptmesh statistics: Nb triangles 6128 , h min 0.00660137 , h max 2.79515 area = 310.608 , M area = 2700.56 , M area/( |Khat| nt) 1.01773 infiny-regularity: min 0.309091 max 2.33784 anisomax 12.7485, beta max = 1.40978 min 0.705754 -- mesh: Nb of Triangles = 6128, Nb of Vertices 3140 kk 4 3 : -- Solve : min -2.12908 max -0.395019 min -13.2917 max 41.7529 min -15.3655 max 19.1957 min 3.31002 max 5.7377 number of required edges : 0 -- adaptmesh statistics: Nb triangles 5633 , h min 0.00522674 , h max 4.69734 area = 310.596 , M area = 2530.78 , M area/( |Khat| nt) 1.03756 infiny-regularity: min 0.264734 max 2.57609 anisomax 21.0718, beta max = 1.63112 min 0.671508 -- mesh: Nb of Triangles = 5633, Nb of Vertices 2891 kk 4 3 : -- Solve : min -2.1224 max -0.392401 min -15.1235 max 43.2889 min -14.4461 max 20.178 min 3.31938 max 5.74137 number of required edges : 0 -- adaptmesh statistics: Nb triangles 6056 , h min 0.00350207 , h max 4.70725 area = 310.569 , M area = 2695.89 , M area/( |Khat| nt) 1.02805 infiny-regularity: min 0.316605 max 2.45241 anisomax 19.4848, beta max = 1.4081 min 0.690603 -- mesh: Nb of Triangles = 6056, Nb of Vertices 3103 kk 4 3 : -- Solve : min -2.12612 max -0.393968 min -18.8864 max 40.4843 min -14.6301 max 19.8193 min 3.31417 max 5.73918 number of required edges : 0 -- adaptmesh statistics: Nb triangles 6366 , h min 0.00344729 , h max 4.79637 area = 310.531 , M area = 2817.91 , M area/( |Khat| nt) 1.02226 infiny-regularity: min 0.325709 max 2.47754 anisomax 19.163, beta max = 1.34129 min 0.721861 -- mesh: Nb of Triangles = 6366, Nb of Vertices 3259 kk 4 3 : -- Solve : min -2.12116 max -0.394795 min -18.7409 max 41.7581 min -14.35 max 19.5358 min 3.32111 max 5.73802 number of required edges : 0 -- adaptmesh statistics: Nb triangles 6248 , h min 0.00226227 , h max 4.53135 area = 310.529 , M area = 2775.26 , M area/( |Khat| nt) 1.0258 infiny-regularity: min 0.310997 max 2.84892 anisomax 18.2171, beta max = 1.41181 min 0.705812 -- mesh: Nb of Triangles = 6248, Nb of Vertices 3202 kk 4 3 : -- Solve : min -2.10447 max -0.395477 min -20.1578 max 40.3441 min -14.3662 max 20.2221 min 3.34447 max 5.73706 number of required edges : 0 -- adaptmesh statistics: Nb triangles 6235 , h min 0.00215549 , h max 4.42957 area = 310.521 , M area = 2759.79 , M area/( |Khat| nt) 1.02221 infiny-regularity: min 0.321439 max 2.56945 anisomax 19.4158, beta max = 1.37586 min 0.726636 -- mesh: Nb of Triangles = 6235, Nb of Vertices 3194 kk 4 3 : -- Solve : min -2.09975 max -0.395245 min -16.7555 max 42.9165 min -14.2681 max 20.6457 min 3.35108 max 5.73739 number of required edges : 0 -- adaptmesh statistics: Nb triangles 6599 , h min 0.00122754 , h max 4.3966 area = 310.527 , M area = 3022.7 , M area/( |Khat| nt) 1.05783 infiny-regularity: min 0.298433 max 2.7827 anisomax 19.4126, beta max = 1.31491 min 0.705859 -- mesh: Nb of Triangles = 6599, Nb of Vertices 3377 kk 4 3 : -- Solve : min -2.09727 max -0.395221 min -18.86 max 43.3604 min -14.1139 max 20.3118 min 3.35455 max 5.73742 number of required edges : 0 -- adaptmesh statistics: Nb triangles 7095 , h min 0.00168344 , h max 4.3589 area = 310.512 , M area = 3261.12 , M area/( |Khat| nt) 1.06149 infiny-regularity: min 0.292642 max 2.57333 anisomax 15.1619, beta max = 1.36726 min 0.705258 -- mesh: Nb of Triangles = 7095, Nb of Vertices 3625 kk 4 3 : -- Solve : min -2.08823 max -0.398927 min -18.2558 max 42.3562 min -14.0679 max 20.3233 min 3.36721 max 5.73223 number of required edges : 0 -- adaptmesh statistics: Nb triangles 7561 , h min 0.00134733 , h max 4.3875 area = 310.51 , M area = 3453.32 , M area/( |Khat| nt) 1.05477 infiny-regularity: min 0.297265 max 2.69928 anisomax 15.6542, beta max = 1.38997 min 0.644608 -- mesh: Nb of Triangles = 7561, Nb of Vertices 3860 shock(): Normal end of execution. current line = 171 exit(0) err code 0 , mpirank 0 CodeAlloc : nb ptr 4465, size :555120 mpirank: 0 WARNING NUMBER bad SearchMethod cas in 2d: 3692 int 3d 0(essai d2: 0 3d: 0 ) Ok: Normal End