-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // two_fluid.edp 2 : // 3 : // Discussion: 4 : // 5 : // Two−fluid incompressible flow 6 : // 7 : // div u = 0 in Omega 8 : // D_t rho = 0 in Omega 9 : // D_t(rho u) + div (mu grad u) + grad p = rho g. 10 : // 11 : // MU is the dynamic viscosity, 12 : // MU = RHO * NU, where NU is the kinematic viscosity. 13 : // 14 : // D_t is the material derivative: 15 : // 16 : // D_t rho = d/dt rho + d/dx rho * v 17 : // 18 : // where v is the convection velocity. 19 : // 20 : // The fluids have different densities, 1 and 10, so the local 21 : // density, between 1 and 10, reflects the relative proportion of 22 : // the two fluids. 23 : // 24 : // Location: 25 : // 26 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/two_fluid/two_fluid.edp 27 : // 28 : // Modified: 29 : // 30 : // 25 June 2015 31 : // 32 : // Author: 33 : // 34 : // Florian De Vuyst 35 : // 36 : // Reference: 37 : // 38 : // Florian De Vuyst, 39 : // Numerical modeling of transport problems using freefem++ software - 40 : // with examples in biology, CFD, traffic flow and energy transfer, 41 : // HAL id: cel-00842234 42 : // https://cel.archives-ouvertes.fr/cel-00842234 43 : // 44 : cout << "\n"; 45 : cout << "two_fluid\n"; 46 : cout << " FreeFem++ version\n"; 47 : cout << " Flow of two incompressible fluids of different ... : densitites.\n"; 48 : 49 : real gravity = 9.81; 50 : real rhom = 1.0; 51 : real rhop = 10.0; 52 : real num = 0.05; 53 : real nup = 0.01; 54 : real dt = 0.2; 55 : // 56 : // Define a fine mesh. 57 : // 58 : real Lx = 2; 59 : real Ly = 1; 60 : mesh Th = square ( 120, 60, [ x * Lx, y * Ly ] ); 61 : plot ( Th, cmm = "two_fluid Mesh", wait = 1, ps = "two_fluid_mesh.ps" ); 62 : // 63 : // Define the finite element space. 64 : // 65 : fespace Uh ( Th, P1b ); 66 : fespace Vh ( Th, P1 ); 67 : fespace Wh ( Th, P2 ); 68 : // 69 : // Define test and trial variables. 70 : // 71 : Uh u1; 72 : Uh u2; 73 : Uh u1old; 74 : Uh u2old; 75 : Uh u1h; 76 : Uh u2h; 77 : Vh p; 78 : Vh ph; 79 : Wh phi; 80 : Wh phih; 81 : Wh phiold; 82 : Wh rho; 83 : Wh rhoold; 84 : Wh mu; 85 : // 86 : // Initialization 87 : // 88 : phi = 0.75 * Ly / Lx * x - y; 89 : phiold = phi; 90 : u1 = 0.0; 91 : u1old = u1; 92 : u2 = 0.0; 93 : u2old = u2; 94 : // 95 : // Define the weak form of the convection equation. 96 : // 97 : problem ConvectLevel ( phi, phih ) = 98 : int2d(Th) ( phi * phih / dt ) 99 : - int2d(Th) ( convect ( [u1old,u2old], -dt, phiold ) * phih / dt ) 100 : + int2d(Th) ( 0.001 * dx(phi)*dx(phih) + 0.001 * dy(phi)*dy(phih) ); 101 : // 102 : // Define the weak form of the Navier Stokes equations. 103 : // 104 : problem NavierStokes ( [u1, u2, p], [u1h, u2h, ph] ) = 105 : int2d(Th) (rho*u1*u1h / dt) 106 : - int2d(Th) (rho*convect([u1old,u2old], -dt, u1old)*u1h/dt) 107 : + int2d(Th) (mu*dx(u1)*dx(u1h)+mu*dy(u1)*dy(u1h)) 108 : - int2d(Th) (p*dx(u1h)) 109 : + int2d(Th) (rho*u2*u2h/dt) 110 : - int2d(Th) (rho*convect([u1old,u2old], -dt, u2old)*u2h/dt) 111 : + int2d(Th) (mu*dx(u2)*dx(u2h)+mu*dy(u2)*dy(u2h)) 112 : - int2d(Th) (p*dy(u2h)) 113 : + int2d(Th) (rho*gravity*u2h) 114 : + int2d(Th) ((dx(u1)+dy(u2))*ph+0.000001*p*ph) 115 : + on ( 1, 2, 3, 4, u1=0.0, u2=0.0 ); 116 : // 117 : for ( int it = 0; it < 30; it++ ) 118 : { 119 : for ( int innerloop = 0; innerloop < 1; innerloop++ ) 120 : { 121 : // 122 : // Convect PHI using the current velocity field. 123 : // 124 : phiold = phi; 125 : ConvectLevel; 126 : // 127 : // Update RHO and MU based on the value of PHI. 128 : // 129 : rhoold = rho; 130 : rho = rhom * (phi<=0) + rhop * (phi>0); 131 : 132 : mu = rhom * num * (phi<=0) + rhop * nup * (phi>0); 133 : // 134 : // Solve the Navier Stokes equations for the new velocity field. 135 : // 136 : NavierStokes; 137 : u1old = u1; 138 : u2old = u2; 139 : } 140 : // 141 : // Plot the density field. 142 : // 143 : plot ( rho, nbiso = 40, fill = 1, value = 1, 144 : cmm = "two_fluid iteration " + it, 145 : ps = "two_fluid_rho_"+it+".ps" ); 146 : 147 : cout << " IT = " << it 148 : << ", pmin = " << p[].min 149 : << ", pmax = " << p[].max << endl; 150 : } 151 : // 152 : // Terminate. 153 : // 154 : cout << "\n"; 155 : cout << "two_fluid:\n"; 156 : cout << " Normal end of execution.\n"; 157 : 158 : sizestack + 1024 =8672 ( 7648 ) two_fluid FreeFem++ version Flow of two incompressible fluids of different densitites. -- Square mesh : nb vertices =7381 , nb triangles = 14400 , nb boundary edges 360 -- Solve : min -0.980507 max 0.730507 -- Solve : min -0.56857 max 0.620852 min -0.581025 max 0.251749 min -14.9432 max 46.346 IT = 0, pmin = -14.9432, pmax = 46.346 -- Solve : min -0.970798 max 0.720726 -- Solve : min -1.0224 max 0.942156 min -0.752543 max 0.593848 min -14.1722 max 44.71 IT = 1, pmin = -14.1722, pmax = 44.71 -- Solve : min -0.963484 max 0.713251 -- Solve : min -1.28041 max 1.01542 min -0.731331 max 1.26126 min -12.8428 max 40.4151 IT = 2, pmin = -12.8428, pmax = 40.4151 -- Solve : min -0.957333 max 0.706875 -- Solve : min -1.13663 max 0.930816 min -0.631655 max 1.7337 min -12.1561 max 39.1485 IT = 3, pmin = -12.1561, pmax = 39.1485 -- Solve : min -0.951877 max 0.701161 -- Solve : min -0.752739 max 0.652685 min -0.492338 max 1.06626 min -11.8058 max 44.6694 IT = 4, pmin = -11.8058, pmax = 44.6694 -- Solve : min -0.946884 max 0.695876 -- Solve : min -0.40872 max 0.358527 min -0.357062 max 0.375007 min -12.2491 max 42.7361 IT = 5, pmin = -12.2491, pmax = 42.7361 -- Solve : min -0.942249 max 0.69089 -- Solve : min -0.409772 max 0.415168 min -0.494369 max 0.461502 min -12.6371 max 41.8575 IT = 6, pmin = -12.6371, pmax = 41.8575 -- Solve : min -0.937921 max 0.686117 -- Solve : min -0.540564 max 0.4933 min -0.456019 max 0.398145 min -12.5306 max 39.8124 IT = 7, pmin = -12.5306, pmax = 39.8124 -- Solve : min -0.933853 max 0.681512 -- Solve : min -0.620099 max 0.61253 min -0.475364 max 0.386954 min -12.1249 max 36.2822 IT = 8, pmin = -12.1249, pmax = 36.2822 -- Solve : min -0.930003 max 0.677074 -- Solve : min -0.64159 max 0.697629 min -0.491826 max 0.440389 min -11.8113 max 34.4881 IT = 9, pmin = -11.8113, pmax = 34.4881 -- Solve : min -0.926338 max 0.672856 -- Solve : min -0.57503 max 0.634355 min -0.415272 max 0.697954 min -11.6773 max 34.4433 IT = 10, pmin = -11.6773, pmax = 34.4433 -- Solve : min -0.922833 max 0.668878 -- Solve : min -0.30102 max 0.280396 min -0.276714 max 0.699287 min -11.8323 max 39.9171 IT = 11, pmin = -11.8323, pmax = 39.9171 -- Solve : min -0.919467 max 0.665067 -- Solve : min -0.353316 max 0.278436 min -0.465512 max 0.304944 min -11.9919 max 38.632 IT = 12, pmin = -11.9919, pmax = 38.632 -- Solve : min -0.916222 max 0.661341 -- Solve : min -0.42525 max 0.357269 min -0.649286 max 0.185181 min -11.9832 max 37.444 IT = 13, pmin = -11.9832, pmax = 37.444 -- Solve : min -0.913085 max 0.657654 -- Solve : min -0.358198 max 0.322015 min -0.461266 max 0.198668 min -11.8372 max 35.2838 IT = 14, pmin = -11.8372, pmax = 35.2838 -- Solve : min -0.910041 max 0.654006 -- Solve : min -0.374861 max 0.342957 min -0.257683 max 0.193499 min -11.6532 max 34.1217 IT = 15, pmin = -11.6532, pmax = 34.1217 -- Solve : min -0.907076 max 0.650419 -- Solve : min -0.38153 max 0.337527 min -0.194218 max 0.305093 min -11.5765 max 33.5531 IT = 16, pmin = -11.5765, pmax = 33.5531 -- Solve : min -0.904175 max 0.646907 -- Solve : min -0.286575 max 0.256146 min -0.137271 max 0.37323 min -11.5445 max 34.9396 IT = 17, pmin = -11.5445, pmax = 34.9396 -- Solve : min -0.901327 max 0.643468 -- Solve : min -0.107192 max 0.107782 min -0.0816166 max 0.28786 min -11.7014 max 36.794 IT = 18, pmin = -11.7014, pmax = 36.794 -- Solve : min -0.898525 max 0.640085 -- Solve : min -0.150252 max 0.177965 min -0.205963 max 0.102182 min -11.7435 max 36.4057 IT = 19, pmin = -11.7435, pmax = 36.4057 -- Solve : min -0.895771 max 0.636749 -- Solve : min -0.196678 max 0.218063 min -0.269354 max 0.166882 min -11.6895 max 35.0127 IT = 20, pmin = -11.6895, pmax = 35.0127 -- Solve : min -0.893063 max 0.633463 -- Solve : min -0.17886 max 0.194859 min -0.208939 max 0.100307 min -11.5686 max 33.2956 IT = 21, pmin = -11.5686, pmax = 33.2956 -- Solve : min -0.890401 max 0.630234 -- Solve : min -0.175634 max 0.204576 min -0.106326 max 0.136523 min -11.4949 max 33.1995 IT = 22, pmin = -11.4949, pmax = 33.1995 -- Solve : min -0.88778 max 0.627069 -- Solve : min -0.161557 max 0.183793 min -0.137927 max 0.149665 min -11.4877 max 33.1512 IT = 23, pmin = -11.4877, pmax = 33.1512 -- Solve : min -0.885197 max 0.623971 -- Solve : min -0.0877914 max 0.0963357 min -0.0853345 max 0.136682 min -11.5306 max 34.4556 IT = 24, pmin = -11.5306, pmax = 34.4556 -- Solve : min -0.882647 max 0.620933 -- Solve : min -0.0616607 max 0.0536106 min -0.0913303 max 0.0701756 min -11.6395 max 35.0162 IT = 25, pmin = -11.6395, pmax = 35.0162 -- Solve : min -0.880129 max 0.617939 -- Solve : min -0.122801 max 0.102138 min -0.122628 max 0.085611 min -11.6292 max 34.5549 IT = 26, pmin = -11.6292, pmax = 34.5549 -- Solve : min -0.87764 max 0.614976 -- Solve : min -0.159911 max 0.131984 min -0.15099 max 0.124965 min -11.5632 max 33.6015 IT = 27, pmin = -11.5632, pmax = 33.6015 -- Solve : min -0.875178 max 0.612032 -- Solve : min -0.145127 max 0.127161 min -0.105908 max 0.0795402 min -11.5053 max 32.6946 IT = 28, pmin = -11.5053, pmax = 32.6946 -- Solve : min -0.87274 max 0.609107 -- Solve : min -0.127463 max 0.103554 min -0.0864743 max 0.111821 min -11.4821 max 32.8929 IT = 29, pmin = -11.4821, pmax = 32.8929 two_fluid: Normal end of execution. times: compile 0.005919s, execution 233.612s, mpirank:0 CodeAlloc : nb ptr 3993, size :499944 mpirank: 0 Ok: Normal End