nsstepad.pde

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

border(1,0,1,6,5)  begin x:=0;      y:=1-t   end;
border(2,0,1,15,2) begin x:=2*t;    y:=0     end;
border(2,0,1,10,2) begin x:=2;      y:=-t    end;
border(2,0,1,20,2) begin x:=2+3*t;  y:=-1    end;
border(2,0,1,35,2) begin x:=5+15*t; y:=-1    end;
border(4,0,1,10,2) begin x:=20;     y:=-1+2*t end;
border(5,0,1,35,5) begin x:=5+15*(1-t); y:=1 end;
border(5,0,1,40,5) begin x:=5*(1-t);y:=1     end;

buildmesh(1500);
nu := 0.001; dt := 0.4;

/*initial values */
u = y*(1-y)*one(y>0);
v = 0;
p = 2*nu*x*one(y>0);
un = u; vn = v; 
i:=1; j:=2; k:=3;count:=1;
iter(300)
  begin 
if (count%100 == 0) then { /* adapt every 50 iterations */
   savemesh('stepmesh');
   adaptmesh(u,p);
   i:=1; j:=2;k:=3;
   wait;
   plot(u);plot(v);plot(p);
   save('stepu',u,count);
   nowait;
 };
 /* increment the counter */
 count:=count+1;

 f=convect(un,u,v,dt);  g=convect(vn,u,v,dt);

 solve(u,i){         /*Horizontal velocity*/
          onbdy(1) u = y*(1-y);
          onbdy(2,5) u = 0;
          onbdy(4)dnu(u)=0;
   pde(u) id(u)/dt-laplace(u)*nu = f/dt -dx(p);
 };
 
 solve(v,j){        /* Vertical velocity */
          onbdy(1,2,4,5) v = 0;
   pde(v) id(v)/dt-laplace(v)*nu = g/dt -dy(p);
 };

         solve(p,k) {      /*  Pressure */
          onbdy(1,2,5) dnu(p) = 0;
          onbdy(4) p=0;
   pde(p)  -laplace(p)= -(dx(f) + dy(g))/dt;
  };
 un = u;  vn = v; i:=-1; j:=-2; k:=-3;
 end ;
wait;
plot(u);plot(v);plot(p);