3d Elasticity problem in beam
![]() |
Initial and distorted mesh |
download example: beam-3d.edp and cube.idp or return to 3D examples
load "medit" include "cube.idp" int[int] Nxyz=[20,5,5]; real [int,int] Bxyz=[[0.,5.],[0.,1.],[0.,1.]]; int [int,int] Lxyz=[[1,2],[2,2],[2,2]]; mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); real E = 21.5e4; real sigma = 0.29; real mu = E/(2*(1+sigma)); real lambda = E*sigma/((1+sigma)*(1-2*sigma)); real gravity = -0.05; fespace Vh(Th,[P1,P1,P1]); Vh [u1,u2,u3], [v1,v2,v3]; cout << "lambda,mu,gravity ="<<lambda<< " " << mu << " " << gravity << endl; real sqrt2=sqrt(2.); macro epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM solve Lame([u1,u2,u3],[v1,v2,v3])= int3d(Th)( lambda*div(u1,u2,u3)*div(v1,v2,v3) +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //') ) - int3d(Th) (gravity*v3) + on(1,u1=0,u2=0,u3=0) ; real dmax= u1[].max; cout << " max deplacement = " << dmax << endl; real coef= 0.1/dmax; int[int] ref2=[1,0,2,0]; mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2); Thm=change(Thm,label=ref2); plot(Th,Thm, wait=1,cmm="coef amplification = "+coef );