Problem konzole

Greda učvršćena na jednom kraju, a na drugom postoji primjena sile prema dolje. Problem je pronaći optimalnu geometriju koja minimizira deformaciju pod zadanom silom. Račun derivacije funkcionala elastične podložnosti (compliance) možete pronaći na poveznici.

load "iovtk"

ofstream data("data.txt");

real E=15;				// Youngov modul
real nu=0.35;			// Poisson koeficijent
func g1=0;
func g2=-0.5*(y-0.4)*(0.6-y);
//////////////////////////////////////////
// Računanje Lameovih koeficijenta 
//////////////////////////////////////////
real lambda=E*nu/((1.+nu)*(1.-2.*nu));
real mu=E/(2.*(1.+nu));
real lagr=3.8*1e-6;

func truncate=(( (abs(x)<0.002 ) || (abs(x-4)<0.005 && abs(y-0.5)<0.2 )  )? 0:1);
//////////////////////////////////////////
// Mesh
//////////////////////////////////////////

border shapeDirichlet(t=0,1){ x=0; y=1-t; label=1; }
border shapeDown(t=0,4){ x=t; y=0.1*t; label=3; }
border shapeNeumann(t=0.4,0.6){ x=4; y=t; label=2; }
border shapeUp(t=0,4){ x=4-t; y=0.6+0.1*t; label=3; }

border hole01(t=0,2*pi){ x=0.5+0.1*cos(t); y=0.5+0.1*sin(t); label=3;}
border hole02(t=0,2*pi){ x=1+0.1*cos(t); y=0.5+0.1*sin(t); label=3;}
border hole03(t=0,2*pi){ x=1.5+0.1*cos(t); y=0.5+0.1*sin(t); label=3;}
border hole04(t=0,2*pi){ x=0.75+0.1*cos(t); y=0.7+0.1*sin(t); label=3;}
border hole05(t=0,2*pi){ x=0.75+0.1*cos(t); y=0.3+0.1*sin(t); label=3;}
border hole06(t=0,2*pi){ x=1.25+0.1*cos(t); y=0.65+0.1*sin(t); label=3;}
border hole07(t=0,2*pi){ x=1.25+0.1*cos(t); y=0.35+0.1*sin(t); label=3;}
border hole08(t=0,2*pi){ x=2.0+0.1*cos(t); y=0.5+0.1*sin(t); label=3;}
border hole09(t=0,2*pi){ x=1.75+0.1*cos(t); y=0.62+0.1*sin(t); label=3;}
border hole10(t=0,2*pi){ x=1.75+0.1*cos(t); y=0.38+0.1*sin(t); label=3;}
border hole11(t=0,2*pi){ x=2.5+0.1*cos(t); y=0.5+0.1*sin(t); label=3;}
border hole12(t=0,2*pi){ x=2.25+0.1*cos(t); y=0.60+0.05*sin(t); label=3;}
border hole13(t=0,2*pi){ x=2.25+0.1*cos(t); y=0.40+0.05*sin(t); label=3;}
border hole14(t=0,2*pi){ x=0.25+0.1*cos(t); y=0.7+0.05*sin(t); label=3;}
border hole15(t=0,2*pi){ x=0.25+0.1*cos(t); y=0.3+0.05*sin(t); label=3;}



int nn=5;
mesh Omega=buildmesh(shapeDirichlet(nn)+shapeDown(4*nn)+shapeNeumann(nn/2)+shapeUp(4*nn)
				+hole01(-3*nn)+hole02(-3*nn)+hole03(-3*nn)+hole04(-3*nn)+hole05(-3*nn)
				+hole06(-3*nn)+hole07(-3*nn)+hole08(-3*nn)+hole09(-3*nn)+hole10(-3*nn)
				+hole11(-3*nn)+hole12(-3*nn)+hole13(-3*nn)+hole14(-3*nn)+hole15(-3*nn)
				); 
real sizeMesh=0.03;
Omega = adaptmesh (Omega,sizeMesh,IsMetric=1,nbvx=100000);



//////////////////////////////////////////
// Elastični sustav 
//////////////////////////////////////////
// Ae(u):e(v)=2mu(e(u):e(v))+lambda(div(u)*div(v))
fespace Vh(Omega,P1);

Vh ux,uy,vx,vy;

problem elasticity([ux,uy],[vx,vy]) =
    int2d(Omega)( 2*mu*( dx(ux)*dx(vx)+dy(uy)*dy(vy)+((dx(uy)+dy(ux))*(dx(vy)+dy(vx)))/2.)
              +lambda*(dx(ux)+dy(uy))*(dx(vx)+dy(vy)) )
    -int1d(Omega,2)(g1*vx+g2*vy)	           
    +on(1, ux=0.0, uy=0.0);

//////////////////////////////////////////
// Konstrukcija vektora silaska 
//////////////////////////////////////////
Vh theta1,theta2;

problem shapeDerivative([theta1,theta2],[vx,vy])=
	int2d(Omega)(sizeMesh*(dx(theta1)*dx(vx)+dy(theta1)*dy(vx)+dx(theta2)*dx(vy)+dy(theta2)*dy(vy))+theta1*vx+theta2*vy)
	+int2d(Omega)(4*mu*(2*(dx(uy)/2 + dy(ux)/2)*((dx(uy)*dx(vx))/2 + (dx(ux)*dy(vx))/2 + (dx(vy)*dy(uy))/2 
					+ (dy(ux)*dy(vy))/2) + dx(ux)*(dx(ux)*dx(vx) + dx(vy)*dy(ux)) + dy(uy)*(dx(uy)*dy(vx) + dy(uy)*dy(vy))) 
				+ 2*lambda*(dx(ux) + dy(uy))*(dx(ux)*dx(vx) + dx(uy)*dy(vx) + dx(vy)*dy(ux) + dy(uy)*dy(vy))
				-(dx(vx)+dy(vy))*(2*mu*( dx(ux)*dx(ux)+dy(uy)*dy(uy)+((dx(uy)+dy(ux))*(dx(uy)+dy(ux)))/2.)
						+lambda*(dx(ux)+dy(uy))*(dx(ux)+dy(uy)))
				)
	+int1d(Omega,3) (lagr*(vx*N.x+vy*N.y))
	+on(1,2,theta1=0.0,theta2=0.0);

//////////////////////////////////////////
// Iteracije 
//////////////////////////////////////////

int MaxIter=650;

real vol0=Omega.measure;
real vol;
for (int i=0; i<MaxIter; ++i)
{
	elasticity;
	shapeDerivative;

	vol=Omega.measure;
	//lagr=lagr*(1+0.05*(vol-vol0)/vol0); // držimo konstantni
	real J=int2d(Omega)( 2*mu*( dx(ux)*dx(ux)+dy(uy)*dy(uy)+((dx(uy)+dy(ux))*(dx(uy)+dy(ux)))/2.)
              +lambda*(dx(ux)+dy(uy))*(dx(ux)+dy(uy)) ) ;
	
	if (i%10==0) 
	{
		savevtk("Omega"+i/10+".vtk",Omega,[theta1,theta2],[ux,uy]);
	}
	
	data<<J<<" "<<vol<<" "<<J+lagr*vol<<" "<<lagr<<" "<<vol-vol0<<endl;
	data.flush;
	
    //////////////////////////////////////////
	//  Update mesh
	//////////////////////////////////////////
    real coef=100;	
	real minT0 = checkmovemesh(Omega, [x, y]);
    while(1){
        real minT = checkmovemesh(Omega, [x+coef*truncate*theta1,y+coef*truncate*theta2]);
        coef/=1.2;
        if (minT > minT0/5.0) break;
    }
    Omega = movemesh(Omega,[x+coef*truncate*theta1,y+coef*truncate*theta2]);
    Omega = adaptmesh (Omega,sizeMesh,IsMetric=1,nbvx=100000,hmin=sizeMesh*3.0);
}

Rezultati