fluidstruct.edp

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

assert(version>=1.24); //  for big bug is non symetric matrix see HISTORY, and sign in int1d 
include "beam.edp"
// Stokes on square  b,e,f,g  driven cavite on left side g 
border e(t=0,10) { x=t; y=-10; label= 1; };      //  bottom
border f(t=0,10) { x=10; y=-10+t ; label= 1; };   //  right
border g(t=0,10) { x=0; y=-t ;label= 2;};       //  left
border h(t=0,10) { x=t; y=vv(t,0)*( t>=0.001 )*(t <= 9.999); label=3;};   //  top of cavity deforme 


mesh sh = buildmesh(h(-20)+f(10)+e(10)+g(10));
plot(sh,wait=1);

fespace Xh(sh,P2),Mh(sh,P1);
fespace V1h(sh,P2);

Xh u1,u2,v1,v2;
Mh p,q,ppp;


varf bx(u1,q) = int2d(sh)( -(dx(u1)*q)); 

varf by(u1,q) = int2d(sh)( -(dy(u1)*q));

varf Lap(u1,u2)= int2d(sh)(  dx(u1)*dx(u2) + dy(u1)*dy(u2) )
                   +  on(2,u1=1) +  on(1,3,u1=0)   ;

varf Mass(p,q)=int2d(sh)(p*q);

Xh bc1; 
Xh brhs;
                   
matrix A= Lap(Xh,Xh,solver=CG); 
matrix M= Mass(Mh,Mh,solver=CG); 
matrix Bx= bx(Xh,Mh);
matrix By= by(Xh,Mh);

Xh bcx,bcy;
func real[int] divup(real[int] & pp)
{ //  
  int verb=verbosity;
   verbosity=1;
   brhs[]  = Bx'*pp; brhs[] += bc1[] .*bcx[];
   u1[] = A^-1*brhs[];
   brhs[]  = By'*pp; brhs[] += bc1[] .*bcy[];
   u2[] = A^-1*brhs[];
   ppp[] = M^-1*pp;
   ppp[] = 1.e-6*ppp[];
   ppp[] =   Bx*u1[];
   ppp[] +=  By*u2[];
   verbosity=verb;
   return ppp[] ;
};
p=0;q=0;u1=0;v1=0;

int i;
for( i=0;i<2;i++)
 {    
bcx=0;
bcy= (-y)*(10.+y)/25.;

 cout << " loop " << i << endl; 
 bc1[] = Lap(0,Xh);
 q=0;   
 verbosity=50;
 LinearCG(divup,p[],eps=1.e-3,nbiter=50);
 divup(p[]);
 verbosity=1;
 plot([u1,u2],p,wait=1,value=true,coef=0.1,cmm="[u1,u2],p");

 real coef=0.2;

V1h sigma11,sigma22,sigma12;
Vh [uu1,vv1];
  [uu1,vv1]=[uu,vv];

  sigma11([x+uu,y+vv]) = (2*dx(u1)-p);
  sigma22([x+uu,y+vv]) = (2*dy(u2)-p);
  sigma12([x+uu,y+vv]) = (dx(u1)+dy(u2));


solve  bbst([uu,vv],[w,s],init=i)  = 
    int2d(th)(
		  lambda*div(w,s)*div(uu,vv)	
                  +2.*mu*( epsilon(w,s)'*epsilon(uu,vv) ) 
             )
  + int2d(th) (-gravity*s)
  + int1d(th,bottombeam)( coef*(sigma11*N.x*w + sigma22*N.y*s + sigma12*(N.y*w+N.x*s) )  )
  + on(1,uu=0,vv=0) 
  ;
cout << " max deplacement " << uu[].linfty << endl;
 plot([uu,vv],wait=1);
 real  err = sqrt(int2d(th)( (uu-uu1)^2 + (vv-vv1)^2 )); 
 cout << " ----- Iteration = " << i <<  " Erreur L2 = " << err << "----------\n";
 th1 = movemesh(th, [x+uu, y+vv]);
 plot(th1,wait=1);
 sh = buildmesh(h(-20)+f(10)+e(10)+g(10));
 plot(sh);
 p=p;q=q;u1=u1;u2=u2;brhs=brhs;ppp=ppp;
 A= Lap(Xh,Xh,solver=CG); 
 M= Mass(Mh,Mh,solver=CG); 
 Bx= bx(Xh,Mh);
 By= by(Xh,Mh);
 bc1=0;  // for resize
 bc1[] = Lap(0,Xh);
}
cout << " max deplacement " << uu[].linfty << endl;