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 -convect([up1,up2],-dt,up1)*v1 -convect([up1,up2],-dt,up2)*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
)
+ 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]);
} ;