nsi3d.edp

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

load "msh3"
real nu=0.01,dt=0.3;
real alpha=1./dt,alpha2=sqrt(alpha);

int nn=5;

mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);
Vh2 ux,uz,p2;
int[int] rup=[0,2],  rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;

mesh3 Th=buildlayers(Th2,nn,
  zbound=[zmin,zmax],
  // reftet=r1, 
  reffacemid=rmid, 
  reffaceup = rup,
  reffacelow = rdown);

fespace VVh(Th,[P23d,P23d,P23d,P13d]);
fespace Vh(Th,P23d);
fespace Ph(Th,P13d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM
  
  varf vStokes([u1,u2,u3,p],[v1,v2,v3,q]) = 
  int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3)
             - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
 +  on(2,u1=1.,u2=0,u3= 0)
 + on(1,u1=0,u2=0,u3=0)
 ;

cout << "b  mat " << endl;

matrix A=vStokes(VVh,VVh);
cout << "e  mat " << endl;
set(A,solver=UMFPACK);
cout << "e fac  mat " << endl;
real[int] b= vStokes(0,VVh);

VVh [u1,u2,u3,p];
VVh [X1,X2,X3,Xp];
VVh [x1,x2,x3,xp]=[x,y,z,0];



u1[]= A^-1 * b;

ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5",wait=1);
macro XX1() (x-u1*dt)//
macro XX2() (y-u2*dt)//
macro XX3() (z-u3*dt)//

  varf vNS([uu1,uu2,uu3,p],[v1,v2,v3,q]) = 
  int3d(Th)( alpha*(uu1*v1+uu2*v2+uu3*v3) + nu*(Grad(uu1)'*Grad(v1) +  Grad(uu2)'*Grad(v2) +  Grad(uu3)'*Grad(v3))
  - div(uu1,uu2,uu3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
  + on(2,uu1=1,uu2=0,uu3=0)
  + on(1,uu1=0,uu2=0,uu3=0)
   +  int3d(Th,optimize=1)(   alpha*(  u1(X1,X2,X3)*v1  +  u2(X1,X2,X3)*v2  +  u3(X1,X2,X3)*v3 )  ) ;
//  +  int3d(Th,optimize=1)(   alpha*(  u1(XX1,XX2,XX3)*v1  +  u2(XX1,XX2,XX3)*v2  +  u3(XX1,XX2,XX3)*v3 )  ) ;
//+  int3d(Th,optimize=1)(   alpha*(  u1(x,y,z)*v1  +  u2(x,y,z)*v2  +  u3(x,y,z)*v3 )  ) ;
//+  int3d(Th,optimize=1)(   alpha*(  u1*v1  +  u2*v2  +  u3*v3 )  ) ;

cout << " build  A" << endl;
A = vNS(VVh,VVh);
cout << " fac A" << endl;
set(A,solver=UMFPACK);
real t=0;
for(int i=0;i<10;++i)
  {
    t += dt;
    cout << " iteration " << i << " t = " << t << endl;
    X1[]=x1[]+u1[]*(-dt);
    //    verbosity=200;
    b=vNS(0,VVh);
    verbosity=2;
    u1[]= A^-1 * b;
    ux= u1(x,0.5,y);
    uz= u3(x,0.5,y);
    p2= p(x,0.5,y);
    plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=0);
    {
      exec("mkdir dd");
      string prefu="dd/u-"+(100+i);
      string prefp="dd/p-"+(100+i);
      savemesh(Th,prefu+".mesh");
      savemesh(Th,prefp+".mesh");
     
      ofstream file(prefu+".bb"); 
      ofstream filep(prefp+".bb"); 
      Ph up1=u1,up2=u2,up3=u3,pp=p;
      file << "3 1 3 "<< up1[].n << " 2 \n";
      filep << "3 1 1 "<< pp[].n << " 2 \n";
      for (int j=0;j<up1[].n ; j++)  
	{
	  file << up1[][j] <<" " <<up2[][j] <<" "<< up3[][j] <<"\n";
	  filep << pp[][j] <<  endl; 
	}  
    }
  }
plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=1);