Good morning,
Please be patient with this request. I am trying to learn Freefem, and decided to try to convert the 2D example in the documentation (under the system of elasticity) into a 3d problem. I am stuck at the point where the boundary conditions are applied.
Here is my attempt at the model:
// Three dimensional plate model
// a trial model
//
// metric model kg,m,s unit system
load “msh3”
// r3 is the surface label vector
int[int] r3=[11,12,13,14,15,16];
// face y=0 is 11
// face x=1 is 12
// face y=1 is 13
// face x=0 is 14
// face z=0 is 15
// face z=1 is 16
int r5 = 5;
// r5 is the volume label
// Material Properties
real Yngs = 210e9; // Pa
real nu = 0.3;
real rho = 7750; // kg/m^3
real mu = Yngs/(2*(1+nu)); // lame coefficient mu
real lambda = Yngsnu/((1+nu)(1-2*nu)); // lame coefficient lambda
// Geometry
real Wdth=0.42; // width
real Lngth=0.42; // length
real Th=0.01; // Thickness
int NW = 20; // Grid resolution width
int NL = 20; // Grid resolution length
int NT = 5; // grid resolution thickness
// Meshing
mesh3 Sp = cube(NW,NL,NT, [Wdthx, Lngthy, Th*z], label=r3, flags=3, region=r5);
plot(Sp);
// gravity type load
real f = -1;
// Define Fespace
fespace Vh(Sp, P2);
Vh u, v, w;
Vh uu, vv, ww;
// define the differential operators
real sqrt2=(2.)^0.5;
macro epsilon(u1,u2,u3) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2, dz(u3), (dz(u2)+dy(u3))/sqrt2, (dz(u1)+dx(u2))/sqrt2] // diff op
macro div(u,v,w) (dx(u)+dy(v)+dz(w))
// Problem
solve lame([u, v, w], [uu, vv, ww])
= int3d(Sp)(
lambda*(div(u, v, w) * div(uu, vv, ww)
+ 2.mu * ( epsilon(u,v,w)’ * epsilon(uu, vv, ww) )
)
- int3d(Sp)(
fww
)
+ on(14, u=0)
+ on(14, v=0)
+ on(14, w=0)
;
plot([u, v, w],wait 1, coef=coef);
// move mesh
mesh SP2 = movemesh(Sp, [x+ucoef, y+vcoef, z+w*coef]);
Freefem crashes at the on command, and I am stuck.
I would appreciate some assistance in getting this example model to work, and would greatly appreciate any advice in the model in general.
Thank you in advance for your time and assistance.
Ed