cavity.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

mesh Th=square(8,8);
fespace Xh(Th,P2);
fespace Mh(Th,P1);
Xh u2,v2;
Xh u1,v1; 
Mh p,q;

solve Stokes ([u1,u2,p],[v1,v2,q],solver=Crout) =
    int2d(Th)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            + p*dx(v1)+ p*dy(v2)
            + dx(u1)*q+ dy(u2)*q
           )
  + on(3,u1=1,u2=0) 
  + on(1,2,4,u1=0,u2=0);
  
plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[u1,u2]);

Xh psi,phi;

solve streamlines(psi,phi) = 
      int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi))
   +  int2d(Th)( -phi*(dy(u1)-dx(u2)))
   +  on(1,2,3,4,psi=0);

plot(psi,wait=1);
int i=0;
real  nu=1./100.;
real dt=0.1;
real alpha=1/dt;

Xh up1,up2; 

problem  NS ([u1,u2,p],[v1,v2,q],solver=Crout,init=i) =
    int2d(Th)(
             alpha*( u1*v1 + u2*v2) 
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(0.000001) 
            + p*dx(v1)+ p*dy(v2)
            + dx(u1)*q+ dy(u2)*q
           )
  + int2d(Th) ( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 )
  + on(3,u1=1,u2=0) 
  + on(1,2,4,u1=0,u2=0)
;

for (i=0;i<=20;i++)
 {
   up1=u1;
   up2=u2;
   NS;
   if ( !(i % 10)) 
    plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[u1,u2]);  
    
 } ;
if ( (i % 10))  
  plot(coef=0.2,cmm=" [u1,u2] et p  ",p,[u1,u2]);  
streamlines;
plot(psi,wait=1);