// Resolution d'un probleme d'elasticite // Schema de Newmark /////////////////////////// real E=15; //Module de Young (toujours positif) real nu=0.35; //Coefficient de Poisson (toujours compris entre -1 et 1/2) real lambda=E*nu/((1.+nu)*(1.-2.*nu));//coef de Lame real mu=E/(2.*(1.+nu)); func g1=-50.; //Forces appliquées func g2=-1.; real T =50.; real dt=.2; real theta=0.8; real delta=0.5; //////////////////////////// real exa=0.05; mesh Th=square(20,5,[4.*x,y]); mesh Sh; fespace Vh(Th,[P1,P1]); //definition de l'espace Vh [ux,uy],[vx,vy],[ux0,uy0],[ux1,uy1];//definition de l'inconnu et des fonctions tests problem elasticite([ux0,uy0],[vx,vy],solver=LU) = int2d(Th)( 2.*mu*(dx(ux0)*dx(vx)+dy(uy0)*dy(vy)+((dx(uy0)+dy(ux0))*(dx(vy)+dy(vx)))/2.) +lambda*(dx(ux0)+dy(uy0))*(dx(vx)+dy(vy)) ) -int1d(Th,2)(g1*vx+g2*vy) +on(4,ux0=0,uy0=0) ; ////////////////////////////////////////////////////////////////////////////: problem elasticiteevolv([ux,uy],[vx,vy],solver=LU) = int2d(Th)( (ux*vx+uy*vy) +(dt*dt)*theta*(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))) ) +int2d(Th)( (dt*dt)*(0.5+delta-2*theta)*(2.*mu*(dx(ux1)*dx(vx)+dy(uy1)*dy(vy)+((dx(uy1)+dy(ux1))*(dx(vy)+dy(vx)))/2.) +lambda*(dx(ux1)+dy(uy1))*(dx(vx)+dy(vy))) +(dt*dt)*(0.5-delta+theta)*(2.*mu*(dx(ux0)*dx(vx)+dy(uy0)*dy(vy)+((dx(uy0)+dy(ux0))*(dx(vy)+dy(vx)))/2.) +lambda*(dx(ux0)+dy(uy0))*(dx(vx)+dy(vy))) -2.*(ux1*vx+uy1*vy) +(ux0*vx+uy0*vy) ) +on(4,ux=0,uy=0) ; elasticite; [ux1,uy1]=[ux0,uy0]; Sh=movemesh(Th,[x+exa*ux1,y+exa*uy1]); plot(Sh,wait=1,bb=[[0,-2.5],[5.,2.5]]); int Niter=T/dt; for(int i=1;i