beam.edp

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

//  Fluid-structure interaction for a weighting beam sitting on a
// square cavity filled with a fluid.

int bottombeam = 2;
border a(t=2,0)  { x=0; y=t ;label=1;};        //  left beam
border b(t=0,10) { x=t; y=0 ;label=bottombeam;};        //  bottom of beam
border c(t=0,2)  { x=10; y=t ;label=1;};       //  rigth beam
border d(t=0,10) { x=10-t; y=2; label=3;};     //  top beam 
real E = 21.5;
real sigma = 0.29;
real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));
real gravity = -0.05;
mesh th = buildmesh( b(20)+c(5)+d(20)+a(5));
fespace Vh(th,[P1,P1]);
Vh [uu,vv], [w,s];
cout << "lambda,mu,gravity ="<<lambda<< " " << mu << " " << gravity << endl;
// deformation of a beam under its own weight 
// definition of 2 macro :
real sqrt2=sqrt(2.);
macro epsilon(u1,u2)  [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u,v) ( dx(u)+dy(v) ) // EOM

solve bb([uu,vv],[w,s],solver=CG)=
	int2d(th)(  
		  lambda*div(w,s)*div(uu,vv)	
                  +2.*mu*( epsilon(w,s)'*epsilon(uu,vv) ) 
                 )
  - int2d(th) (gravity*s)
  + on(1,uu=0,vv=0)
  ;

//verbosity=100;
plot([uu,vv],wait=1);
plot([uu,vv],wait=1,bb=[[-0.5,2.5],[2.5,-0.5]]);
mesh th1 = movemesh(th, [x+uu, y+vv]);
plot(th1,wait=1);
cout << " max deplacement = " << uu[].linfty << endl;