Hello, guys
Regarding how to use the lame formula (elasiticity) to simulate a 1m x 0.5m x 0.002m plate, the condition of the plate is a simple support condition, and the internal grid of the plate is a tetrahedron, and external force is applied to the center of the plate, and then as time increases, record the midpoint of the plate displacement.
The following is the program I wrote, but now I am very troubled about how to write the simple support conditions of the plate and the relationship between vibration and time
// Parameters
real E = 71.6e9, nu = 0.23;
real f = -100;//gravity
// Mesh
load “msh3”
load “tetgen”
load “medit”
real x0,x1,y0,y1,x2,x3,y2,y3;
x0=0.; x1=1.; y0=0.; y1=0.5;
int[int] lab1 = [20,21,22,23];
mesh Thsq1 = square(100,100,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
x2=0.; x3=1.; y2=0.; y3=0.5;
int[int] lab2 = [24,25,26,27];
mesh Th1 = square(100,100,[x2+(x3-x2)*x,y2+(y3-y2)*y]);
func ZZ1min = 0;
func ZZ1max = 0.002;
func XX1 = x;
func YY1 = y;
int[int] ref31h = [0,11];
int[int] ref31b = [0,12];
meshS Th31h = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1max],region=ref31h,orientation=1);
savemesh(Th31h,“Th31h.mesh”);
meshS Th31b = movemesh23(Th1,transfo=[XX1,YY1,ZZ1min],region=ref31b,orientation=-1);
Th31h = change(Th31h, label=lab1);
Th31b = change(Th31b, label=lab2);
/////////////////////////////////
x0=0.; x1=1.; y0=0.; y1=0.002;
int[int] lab3 = [26,28,22,29];
mesh Thsq2 = square(100,1,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
x2=0.; x3=1.; y2=0.; y3=0.002;
int[int] lab4 = [24,31,20,30];
mesh Thsq2q = square(100,1,[x2+(x3-x2)*x,y2+(y3-y2)*y]);
func ZZ2 = y;
func XX2 = x;
func YY2min = 0.;
func YY2max = 0.5;
int[int] ref32h = [0,13];
int[int] ref32b = [0,14];
meshS Th32h = movemesh23(Thsq2,transfo=[XX2,YY2max,ZZ2],region=ref32h,orientation=-1);
savemesh(Th32h,“Th32h.mesh”);
meshS Th32b = movemesh23(Thsq2q,transfo=[XX2,YY2min,ZZ2],region=ref32b,orientation=1);
Th32h = change(Th32h, label=lab3);
Th32b = change(Th32b, label=lab4);
/////////////////////////////////
x0=0.; x1=0.5; y0=0.; y1=0.002;
int[int] lab5 = [25,31,21,28];
mesh Thsq3 = square(100,1,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
x2=0.; x3=0.5; y2=0.; y3=0.002;
int[int] lab6 = [27,30,23,29];
mesh Thsq3q = square(100,1,[x2+(x3-x2)*x,y2+(y3-y2)*y]);
func XX3min = 0;
func XX3max = 1.;
func YY3 = x;
func ZZ3 = y;
int[int] ref33h = [0,15];
int[int] ref33b = [0,16];
meshS Th33h = movemesh23(Thsq3,transfo=[XX3max,YY3,ZZ3],region=ref33h,orientation=1);
meshS Th33b = movemesh23(Thsq3q,transfo=[XX3min,YY3,ZZ3],region=ref33b,orientation=-1);
Th33h = change(Th33h, label=lab5);
Th33b = change(Th33b, label=lab6);
////////////////////////////////
meshS Th33 = Th31h+Th31b+Th32h+Th32b+Th33h+Th33b; // “gluing” surface meshs to obtain the surface of cube
medit(“glumesh”,Th33);
Th33=change(Th33, rmInternalEdges=true);
savemesh(Th33,“Th33.mesh”);
// build a mesh of a axis parallel box with TetGen
real[int] domain = [1,0.5,0.002,145,0.001];
mesh3 Th44 = tetg(Th33,switch=“pqaAAYYQ”,nbofregions=1,regionlist=domain); // Tetrahelize the interior of the cube with tetgen
medit(“tetg”,Th44);
savemesh(Th44,“Th44.mesh”);
// Problem
real mu = E/(2*(1 + nu));
real lambda = Enu/((1 + nu)(1 - 2*nu));
//Fespace
fespace Vh(Th44,[P1,P1,P1]);
Vh [u1,u2,u3], [v1,v2,v3] ;
varf m([u1,u2,u3],[v1,v2,v3]) = int3d(Th44) ();
matrix mass=m(Vh,Vh);
cout << “lambda,mu,f =”<<lambda<< " " << mu << " " << f << endl;
//Macro
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] // epsilon(u1,u2,u3)’*epsilon(v1,v2,v3) = epsilon(u):epsilon(v)
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM
solve Lame([u1,u2,u3],[v1,v2,v3])=
int3d(Th44)(lambda*div(u1,u2,u3)*div(v1,v2,v3)+2.mu( epsilon(u1,u2,u3)'epsilon(v1,v2,v3) ) )-int3d(Th44) (fv3) + on(ref32h,u1=0,u2=0,u3=0) + on(ref32b,u1=0,u2=0,u3=0)+ on(ref33b,u1=0,u2=0,u3=0)+ on(ref33h,u1=0,u2=0,u3=0) ;
//Display
real dmax= u1[].max;
cout << " max deplacement = " << dmax << endl;
// movemesh
real coef= 0.1/dmax;
int[int] ref2=[100,101,102,103];
mesh3 Thm=movemesh3(Th44,transfo=[x+u1coef,y+u2coef,z+u3*coef],label=ref2);
Thm=change(Thm,label=ref2);
//Plot
plot(Th44,Thm, wait=true,cmm="coef amplification = "+coef );