-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // Location: 2 : // 3 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/poiseuille/poiseuille.edp 4 : // 5 : // Discussion: 6 : // 7 : // Solve the Navier Stokes equations defining Poiseuille flow in a channel. 8 : // 9 : // This example was adapted from code at 10 : // https://modules.freefem.org/modules/stokes/ 11 : // 12 : // Licensing: 13 : // 14 : // This code is distributed under the GNU LGPL license. 15 : // 16 : // Modified: 17 : // 18 : // 12 February 2021 19 : // 20 : // Reference: 21 : // 22 : cout << "\n"; 23 : cout << "poiseuille:\n"; 24 : cout << " FreeFem++ version.\n"; 25 : cout << " Poiseuille flow through a channel.\n"; 26 : 27 : // Parameters 28 : real uMax = 10.; 29 : 30 : real Mu = 1.; 31 : 32 : func fx = 0.; 33 : func fy = 0.; 34 : // 35 : // Geometry 36 : // 37 : int nn = 4; //Mesh quality 38 : real L = 5.0; //Pipe length 39 : real D = 1.0; //Pipe height 40 : int Wall = 1; //Pipe wall label 41 : int Inlet = 2; //Pipe inlet label 42 : int Outlet = 3; //Pipe outlet label 43 : // 44 : // Wall (y=D) 45 : // ------------------------------ 46 : // :---> : 47 : // Inlet :-----> : Outlet 48 : // (x=0) :---> : (x=L) 49 : // ------------------------------ 50 : // Wall (y=0) 51 : // 52 : border b1 ( t = 0.0, 1.0 ) { x=L*t; y=0.0; label=Wall; }; 53 : border b2 ( t = 0.0, 1.0 ) { x=L; y=D*t; label=Outlet; }; 54 : border b3 ( t = 0.0, 1.0 ) { x=L-L*t; y=D; label=Wall; }; 55 : border b4 ( t = 0.0, 1.0 ) { x=0.0; y=D-D*t; label=Inlet; }; 56 : 57 : int nnL = max ( 2.0, L * nn ); 58 : int nnD = max ( 2.0, D * nn ); 59 : // 60 : // Build the mesh. 61 : // 62 : mesh Th = buildmesh ( b1 ( nnL ) 63 : + b2 ( nnD ) 64 : + b3 ( nnL ) 65 : + b4 ( nnD ) ); 66 : // 67 : // Fespace 68 : // Set the finite element velocity space to piecewise quadratics. 69 : // 70 : fespace Uh ( Th, [ P2, P2 ] ); 71 : Uh [ux, uy]; 72 : Uh [vx, vy]; 73 : // 74 : // Set the finite element pressure space to piecewise linears. 75 : // 76 : fespace Ph(Th, P1); 77 : Ph p; 78 : Ph q; 79 : // 80 : // Macros 81 : // 82 : macro Gradient(u) [dx(u), dy(u)] ) // 83 : macro Divergence(ux, uy) (dx(ux) + dy(uy)) ) // 84 : // 85 : // Define a function for the parabolic inflow velocity. 86 : // 87 : func uIn = uMax * (1.-(y-D/2.)^2 / (D/2.)^2); 88 : // 89 : // Problem: 90 : // Solve for velocity u=(ux,uy) and pressure p, 91 : // using test functions v=(vx,vy) and q, 92 : // which zero out the following function. 93 : // 94 : problem S ( [ux, uy, p],[vx, vy, q]) 95 : = int2d(Th)( 96 : Mu * ( 97 : Gradient(ux) [dx(ux), dy(ux)] ' * Gradient(vx) [dx(vx), dy(vx)] 98 : + Gradient(uy) [dx(uy), dy(uy)] ' * Gradient(vy) [dx(vy), dy(vy)] 99 : ) 100 : - p * Divergence(vx, vy) (dx(vx) + dy( vy)) 101 : - Divergence(ux, uy) (dx(ux) + dy( uy)) * q 102 : ) 103 : - int2d(Th)( 104 : fx*vx + fy*vy 105 : ) 106 : + on ( Inlet, ux = uIn, uy = 0.0 ) 107 : + on ( Wall, ux = 0.0, uy = 0.0 ) 108 : ; 109 : // 110 : // Request the solution ( ux, uy, p ) 111 : // 112 : S; 113 : // 114 : // Plot the velocity vector field. 115 : // 116 : plot ( [ ux, uy ], nbarrow = 1, ps = "poiseuille_velocity.ps", 117 : cmm = "Poiseuille Velocity", wait = 1 ); 118 : // 119 : // Plot the pressure as color contours. 120 : // 121 : plot ( p, nbiso = 10, fill = 1, ps = "poiseuille_pressure.ps", 122 : cmm = "Poiseuille Pressure" ); 123 : // 124 : // Terminate. 125 : // 126 : cout << "\n"; 127 : cout << "poiseuille:\n"; 128 : cout << " Normal end of execution.\n"; 129 : 130 : exit ( 0 ); 131 : sizestack + 1024 =2064 ( 1040 ) poiseuille: FreeFem++ version. Poiseuille flow through a channel. -- mesh: Nb of Triangles = 168, Nb of Vertices 109 -- Solve : min -9.07482e-15 max 10 min -1.37644e-14 max 400 poiseuille: Normal end of execution. current line = 130 exit(0) err code 0 , mpirank 0 CodeAlloc : nb ptr 3728, size :483256 mpirank: 0 Ok: Normal End