-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae) file : thermal_convection.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // thermal_convection.edp 2 : // 3 : // Discussion: 4 : // 5 : // Forced + Natural heat convection in a pipe 6 : // Navier−Stokes equations + convection−diffusion on temperature 7 : // 8 : // Licensing: 9 : // 10 : // This code is distributed under the MIT license. 11 : // 12 : // Modified: 13 : // 14 : // 26 June 2015 15 : // 16 : // Author: 17 : // 18 : // Florian De Vuyst 19 : // 20 : // Reference: 21 : // 22 : // Florian De Vuyst, 23 : // Numerical modeling of transport problems using freefem++ software - 24 : // with examples in biology, CFD, traffic flow and energy transfer, 25 : // HAL id: cel-00842234 26 : // https://cel.archives-ouvertes.fr/cel-00842234 27 : // 28 : cout << "\n"; 29 : cout << "thermal_convection():\n"; 30 : cout << " FreeFem++ version\n"; 31 : cout << " Thermal convection and diffusion by Navier-Stok ... : es flow.\n"; 32 : // 33 : // Define the boundary lines. 34 : // 35 : real lx = 0.25; 36 : real Lx = 3.0; 37 : real Ly = 1.0; 38 : 39 : border c1 (t=0.0, lx) {x = t; y = 0.0;} 40 : border c2 (t=lx, Lx-lx) {x = t; y = 0.0;} 41 : border c3 (t=Lx-lx,Lx) {x = t; y = 0.0;} 42 : border c4 (t=0.0,Ly) {x = Lx; y = t;} 43 : border c5 (t=Lx,0.0) {x = t; y = Ly;} 44 : border c6 (t=Ly,0.0) {x = 0.0; y = t;} 45 : // 46 : // Define the mesh. 47 : // 48 : mesh Th = buildmesh ( c1(10) + c2(200) + c3(10) + c4(30) + c5(100) + c6(30) ); 49 : // 50 : // Plot the mesh. 51 : // 52 : plot ( Th, wait = false, ps = "thermal_convection_mesh.ps" ); 53 : // 54 : // Define the finite element spaces: 55 : // UH for velocity fields, 56 : // XH for pressure P, temperature THETA, and specific volume TAU. 57 : // 58 : fespace Uh(Th, P1b ); 59 : Uh u, v, uold, vold, uh, vh; 60 : fespace Xh(Th, P1 ); 61 : Xh p, ph; 62 : Xh theta, thetaold, thetah; 63 : Xh tau; 64 : // 65 : // Initialize the fields. 66 : // 67 : real uin0 = 0.2; 68 : func uin = uin0 * 4.0 * (y/Ly) * (1.0-y/Ly); 69 : u = uin; 70 : uold = u; 71 : v = 0.0; 72 : vold = v; 73 : 74 : real thetain = 20.0; 75 : theta = thetain; 76 : thetaold = theta; 77 : // 78 : // Define the weak form of the heat equation. 79 : // 80 : real dt = 0.05; 81 : real heatflux = 100.0; 82 : real kappa = 0.001; 83 : 84 : problem HeatStep ( theta, thetah ) = 85 : int2d (Th) ( theta * thetah / dt ) 86 : - int2d (Th) ( convect ([u,v], -dt, thetaold)*thetah/dt ) 87 : + int2d (Th) ( kappa * dx(theta) * dx(thetah) + kappa * dy(theta) * dy(thetah) ) 88 : - int1d (Th, c2) ( kappa * heatflux * thetah ) 89 : + on ( c6, theta = thetain ); 90 : // 91 : // Define the weak form of the Navier Stokes equations. 92 : // 93 : real gravity = 9.81; 94 : real nu = 0.001; 95 : 96 : problem NavierStokesStep ( [u, v, p], [uh, vh, ph] ) = 97 : int2d (Th) ( u * uh / dt ) 98 : - int2d (Th) ( convect ( [uold, vold], -dt, uold ) * uh/dt ) 99 : + int2d (Th) ( nu * dx(u) * dx(uh) + nu * dy(u) * dy(uh) ) 100 : + int2d (Th) ( tau * dx(p) * uh ) 101 : + int2d (Th) ( v * vh / dt ) 102 : - int2d (Th) ( convect ([uold, vold], -dt, vold) * vh / dt ) 103 : + int2d (Th) ( nu * dx(v) * dx(vh) + nu * dy(v) * dy(vh) ) 104 : + int2d (Th) ( tau * dy(p) * vh ) 105 : + int2d (Th) ( gravity * vh ) 106 : + int2d (Th) ( dx(u) * ph + dy(v) * ph ) 107 : + on ( c6, u=uin, v=0 ) 108 : + on ( c1, c2, c3, c5, u=0, v=0 ); 109 : // 110 : // Time steps 111 : // 112 : real cst = 1.0; 113 : real dtsnap = 0.5; 114 : real t = 0.0; 115 : real ttosnap = 2.0; 116 : 117 : for ( int it = 1; it <= 100; it++ ) 118 : { 119 : t = it * dt; 120 : cout << it << " T = " << t << "\n"; 121 : 122 : HeatStep; 123 : 124 : thetaold = theta; 125 : tau = cst * theta; 126 : 127 : NavierStokesStep; 128 : 129 : uold = u; 130 : vold = v; 131 : 132 : plot ( [ u, v ], theta, nbiso = 80, value = false ); 133 : 134 : if ( ttosnap <= t ) 135 : { 136 : cout << " Plotting T = " << t << "\n"; 137 : ttosnap = ttosnap + dtsnap; 138 : plot ( [u,v], theta, nbiso = 80, value = false, 139 : ps = "thermal_convection_"+it+".ps", 140 : cmm = "Step "+it ); 141 : } 142 : } 143 : // 144 : // Terminate. 145 : // 146 : cout << "\n"; 147 : cout << "thermal_convection():\n"; 148 : cout << " Normal end of execution.\n"; 149 : 150 : exit ( 0 ); 151 : sizestack + 1024 =7024 ( 6000 ) thermal_convection(): FreeFem++ version Thermal convection and diffusion by Navier-Stokes flow. -- mesh: Nb of Triangles = 15176, Nb of Vertices 7779 1 T = 0.05 -- Solve : min 19.9863 max 20.697 -- Solve : min -2.81046e-35 max 0.225722 min -0.060748 max 0.16685 min -3.67314e+10 max -3.67314e+10 2 T = 0.1 -- Solve : min 19.9962 max 21.0788 -- Solve : min -0.00026222 max 0.214376 min -0.0779192 max 0.0567255 min -2.17108e+10 max -2.17108e+10 3 T = 0.15 -- Solve : min 19.9969 max 21.3532 -- Solve : min -0.00117345 max 0.219653 min -0.0681308 max 0.0679538 min 4.54199e+10 max 4.54199e+10 4 T = 0.2 -- Solve : min 19.9975 max 21.577 -- Solve : min -0.00318042 max 0.218825 min -0.151085 max 0.0803894 min 9.34453e+10 max 9.34453e+10 5 T = 0.25 -- Solve : min 19.9988 max 21.7711 -- Solve : min -1.74673e-35 max 0.211645 min -0.160392 max 0.031055 min -3.17776e+10 max -3.17776e+10 6 T = 0.3 -- Solve : min 19.9986 max 21.945 -- Solve : min -2.23281e-35 max 0.208349 min -0.191443 max 0.0479605 min -4.23814e+10 max -4.23814e+10 7 T = 0.35 -- Solve : min 19.999 max 22.1041 -- Solve : min -1.84797e-35 max 0.216804 min -0.064856 max 0.0561387 min -2.90488e+10 max -2.90488e+10 8 T = 0.4 -- Solve : min 19.9996 max 22.2516 -- Solve : min -3.65856e-05 max 0.250065 min -0.0569217 max 0.1331 min 2.30037e+10 max 2.30037e+10 9 T = 0.45 -- Solve : min 19.9997 max 22.3898 -- Solve : min -0.00118847 max 0.209453 min -0.0380499 max 0.0557881 min -3.24223e+10 max -3.24223e+10 10 T = 0.5 -- Solve : min 19.9997 max 22.5203 -- Solve : min -0.0024451 max 0.205592 min -0.027432 max 0.0189372 min -8.91012e+09 max -8.91012e+09 11 T = 0.55 -- Solve : min 19.9997 max 22.6442 -- Solve : min -0.00261342 max 0.241459 min -0.091613 max 0.0443903 min -2.8758e+10 max -2.8758e+10 12 T = 0.6 -- Solve : min 19.9998 max 22.7625 -- Solve : min -0.00549416 max 0.215497 min -0.0456777 max 0.0632652 min -2.10931e+10 max -2.10931e+10 13 T = 0.65 -- Solve : min 19.9999 max 22.8759 -- Solve : min -0.00239848 max 0.29757 min -0.0355404 max 0.0824904 min 3.11767e+10 max 3.11767e+10 14 T = 0.7 -- Solve : min 19.9999 max 22.985 -- Solve : min -0.00599389 max 0.221352 min -0.0205985 max 0.0391449 min 6.93234e+10 max 6.93234e+10 15 T = 0.75 -- Solve : min 19.9999 max 23.09 -- Solve : min -0.0106369 max 0.230999 min -0.0285129 max 0.102421 min 2.38048e+10 max 2.38048e+10 16 T = 0.8 -- Solve : min 19.9999 max 23.1915 -- Solve : min -0.0122546 max 0.231284 min -0.0381023 max 0.0393722 min 5.81345e+10 max 5.81345e+10 17 T = 0.85 -- Solve : min 19.9999 max 23.2899 -- Solve : min -0.0165869 max 0.261107 min -0.125676 max 0.0791357 min 1.23072e+11 max 1.23072e+11 18 T = 0.9 -- Solve : min 19.9999 max 23.3853 -- Solve : min -0.0134045 max 0.21878 min -0.0531634 max 0.0356195 min -4.30696e+10 max -4.30696e+10 19 T = 0.95 -- Solve : min 20 max 23.4781 -- Solve : min -0.015173 max 0.214214 min -0.042227 max 0.0156637 min 2.38871e+10 max 2.38871e+10 20 T = 1 -- Solve : min 20 max 23.5686 -- Solve : min -0.022381 max 0.227157 min -0.0564961 max 0.107002 min -7.72266e+10 max -7.72266e+10 21 T = 1.05 -- Solve : min 20 max 23.6569 -- Solve : min -0.0223892 max 0.265083 min -0.0541264 max 0.104326 min -1.99683e+10 max -1.99683e+10 22 T = 1.1 -- Solve : min 20 max 23.7432 -- Solve : min -0.199478 max 1.26777 min -1.22154 max 0.399861 min 5.18385e+11 max 5.18385e+11 23 T = 1.15 -- Solve : min 20 max 23.8508 -- Solve : min -0.235246 max 0.460067 min -0.356255 max 0.420641 min -1.17445e+11 max -1.17445e+11 24 T = 1.2 -- Solve : min 20 max 23.9423 -- Solve : min -0.576236 max 1.00824 min -0.808017 max 1.45994 min -7.46166e+11 max -7.46166e+11 25 T = 1.25 -- Solve : min 20 max 24.0233 -- Solve : min -0.138471 max 0.388089 min -0.349485 max 1.1663 min -3.14353e+11 max -3.14353e+11 26 T = 1.3 -- Solve : min 20 max 24.0989 -- Solve : min -0.0292591 max 0.280856 min -0.191438 max 0.170283 min -2.0538e+11 max -2.0538e+11 27 T = 1.35 -- Solve : min 20 max 24.1783 -- Solve : min -0.0484757 max 0.255891 min -0.31318 max 0.257819 min 2.59352e+11 max 2.59352e+11 28 T = 1.4 -- Solve : min 20 max 24.2522 -- Solve : min -0.0549778 max 0.241511 min -0.110925 max 0.0519308 min -3.763e+10 max -3.763e+10 29 T = 1.45 -- Solve : min 20 max 24.3278 -- Solve : min -0.394753 max 0.8588 min -0.264927 max 1.88541 min -2.6981e+11 max -2.6981e+11 30 T = 1.5 -- Solve : min 20 max 24.4065 -- Solve : min -0.131223 max 0.290611 min -0.258218 max 0.352683 min 1.51511e+11 max 1.51511e+11 31 T = 1.55 -- Solve : min 20 max 24.4868 -- Solve : min -0.506159 max 1.41437 min -1.67384 max 1.14255 min -7.18872e+11 max -7.18872e+11 32 T = 1.6 -- Solve : min 20 max 24.5688 -- Solve : min -0.0391601 max 0.365789 min -0.294989 max 0.225591 min -7.46471e+10 max -7.46471e+10 33 T = 1.65 -- Solve : min 20 max 24.651 -- Solve : min -0.0684153 max 0.343546 min -0.171737 max 0.577209 min -2.72421e+11 max -2.72421e+11 34 T = 1.7 -- Solve : min 20 max 24.7329 -- Solve : min -0.0695827 max 0.27591 min -0.287283 max 0.401207 min -1.28561e+11 max -1.28561e+11 35 T = 1.75 -- Solve : min 20 max 24.8145 -- Solve : min -0.0583076 max 0.273951 min -0.168147 max 0.234318 min 3.38126e+11 max 3.38126e+11 36 T = 1.8 -- Solve : min 20 max 24.897 -- Solve : min -0.0602533 max 0.277021 min -0.14158 max 0.336562 min -1.55748e+11 max -1.55748e+11 37 T = 1.85 -- Solve : min 20 max 24.9802 -- Solve : min -0.0590477 max 0.24522 min -0.0576279 max 0.126515 min 8.61054e+10 max 8.61054e+10 38 T = 1.9 -- Solve : min 20 max 25.0639 -- Solve : min -0.062818 max 0.231645 min -0.0509552 max 0.10123 min 2.5229e+10 max 2.5229e+10 39 T = 1.95 -- Solve : min 20 max 25.1481 -- Solve : min -0.0715465 max 0.370924 min -0.369377 max 0.442397 min -8.30661e+10 max -8.30661e+10 40 T = 2 -- Solve : min 20 max 25.2319 -- Solve : min -0.0793033 max 0.254058 min -0.123272 max 0.0971548 min 5.93652e+10 max 5.93652e+10 Plotting T = 2 41 T = 2.05 -- Solve : min 20 max 25.3157 -- Solve : min -0.0792795 max 0.240811 min -0.0459141 max 0.0882246 min 8.15505e+10 max 8.15505e+10 42 T = 2.1 -- Solve : min 20 max 25.3979 -- Solve : min -0.0807606 max 0.215218 min -0.0391735 max 0.0806496 min -1.60648e+10 max -1.60648e+10 43 T = 2.15 -- Solve : min 20 max 25.4814 -- Solve : min -0.0807274 max 0.224211 min -0.0412664 max 0.0735849 min 2.71915e+10 max 2.71915e+10 44 T = 2.2 -- Solve : min 20 max 25.57 -- Solve : min -0.085191 max 0.212195 min -0.0430483 max 0.0706384 min -9.16202e+09 max -9.16202e+09 45 T = 2.25 -- Solve : min 20 max 25.6568 -- Solve : min -0.0903324 max 0.237514 min -0.0450332 max 0.0762124 min 1.47435e+10 max 1.47435e+10 46 T = 2.3 -- Solve : min 20 max 25.7398 -- Solve : min -0.0945768 max 0.24143 min -0.0475524 max 0.0813579 min -2.30028e+10 max -2.30028e+10 47 T = 2.35 -- Solve : min 20 max 25.8278 -- Solve : min -0.096504 max 0.253104 min -0.12797 max 0.0876053 min 8.55339e+10 max 8.55339e+10 48 T = 2.4 -- Solve : min 20 max 25.9194 -- Solve : min -26.6226 max 15.5495 min -81.5365 max 21.9361 min -8.24977e+12 max -8.24977e+12 49 T = 2.45 -- Solve : min 19.9998 max 26.1203 -- Solve : min -48.467 max 58.4933 min -69.5363 max 35.775 min 3.38579e+13 max 3.38579e+13 50 T = 2.5 -- Solve : min 19.996 max 26.2868 -- Solve : min -70.9334 max 87.1923 min -96.2246 max 421.801 min 7.39674e+13 max 7.39674e+13 Plotting T = 2.5 51 T = 2.55 -- Solve : min 19.9954 max 26.5695 -- Solve : min -44.9006 max 51.8261 min -105.974 max 85.8047 min 6.89377e+13 max 6.89377e+13 52 T = 2.6 -- Solve : min 19.9718 max 26.4651 -- Solve : min -21.6164 max 24.4574 min -21.9879 max 96.5548 min -2.57297e+13 max -2.57297e+13 53 T = 2.65 -- Solve : min 19.9952 max 26.4932 -- Solve : min -4.28893 max 5.33583 min -6.11996 max 5.4402 min 3.40838e+12 max 3.40838e+12 54 T = 2.7 -- Solve : min 19.9921 max 26.5615 -- Solve : min -28.4734 max 37.8 min -31.7333 max 81.9512 min -4.48902e+13 max -4.48902e+13 55 T = 2.75 -- Solve : min 19.9787 max 26.5836 -- Solve : min -21.9375 max 18.7849 min -40.0069 max 51.1932 min 1.1497e+13 max 1.1497e+13 56 T = 2.8 -- Solve : min 19.9813 max 26.718 -- Solve : min -25.3904 max 26.0839 min -31.858 max 18.1744 min -1.36011e+13 max -1.36011e+13 57 T = 2.85 -- Solve : min 19.9797 max 26.7777 -- Solve : min -8.01807 max 3.86676 min -4.03997 max 5.21179 min -3.35638e+11 max -3.35638e+11 58 T = 2.9 -- Solve : min 19.9845 max 26.9079 -- Solve : min -29.8331 max 70.8342 min -176.904 max 36.8792 min -4.85694e+13 max -4.85694e+13 59 T = 2.95 -- Solve : min 19.9609 max 27.0889 -- Solve : min -120.123 max 80.8621 min -100.657 max 111.008 min 7.54761e+13 max 7.54761e+13 60 T = 3 -- Solve : min 19.9544 max 27.3598 -- Solve : min -21.9384 max 34.8995 min -69.262 max 77.2632 min 2.51833e+13 max 2.51833e+13 Plotting T = 3 61 T = 3.05 -- Solve : min 19.9847 max 27.412 -- Solve : min -13.6369 max 9.57183 min -10.803 max 11.1987 min -1.93364e+13 max -1.93364e+13 62 T = 3.1 -- Solve : min 19.8979 max 27.496 -- Solve : min -9.42367 max 11.9203 min -14.6393 max 55.0899 min 1.98922e+13 max 1.98922e+13 63 T = 3.15 -- Solve : min 19.9776 max 27.5056 -- Solve : min -69.6881 max 31.1194 min -85.1621 max 231.999 min 1.66303e+14 max 1.66303e+14 64 T = 3.2 -- Solve : min 19.9243 max 27.4914 -- Solve : min -24.9462 max 67.9601 min -116.728 max 32.6775 min -4.22125e+14 max -4.22125e+14 65 T = 3.25 -- Solve : min 19.9374 max 27.1278 -- Solve : min -26.3223 max 66.2509 min -32.9489 max 55.4068 min 4.08856e+13 max 4.08856e+13 66 T = 3.3 -- Solve : min 19.9083 max 27.0136 -- Solve : min -13.1004 max 18.6966 min -13.5292 max 26.6985 min -1.20085e+13 max -1.20085e+13 67 T = 3.35 -- Solve : min 19.9556 max 26.9599 -- Solve : min -31.61 max 41.523 min -61.5432 max 71.335 min -4.40945e+13 max -4.40945e+13 68 T = 3.4 -- Solve : min 19.9605 max 26.8619 -- Solve : min -82.8602 max 54.9428 min -68.4696 max 192.3 min -8.6003e+13 max -8.6003e+13 69 T = 3.45 -- Solve : min 19.9586 max 26.8223 -- Solve : min -381.506 max 263.879 min -706.942 max 677.659 min -2.60027e+14 max -2.60027e+14 70 T = 3.5 -- Solve : min 19.9408 max 26.7412 -- Solve : min -433.556 max 537.393 min -524.196 max 946.879 min 5.44239e+14 max 5.44239e+14 Plotting T = 3.5 71 T = 3.55 -- Solve : min 19.9762 max 25.7175 -- Solve : min -1076.97 max 2489.08 min -823.545 max 5381.8 min 6.85464e+14 max 6.85464e+14 72 T = 3.6 Fatal error in Convect (R2) operator: loop => velocity too high ???? or NaN F. Hecht current line = 122 Assertion fail : (0) line :5177, in file lgfem.cpp Assertion fail : (0) line :5177, in file lgfem.cpp err code 6 , mpirank 0