Request help with a 3D model

[EChase (Ed)]

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

[fb77 (François Bouchut)]

You just miss a closing parenthese for the term with lambda in factor.

[Monirul25 (Monirul)]

Here, is your Code with several corrections.
FEM_elasticity-3D.edp (2.0 KB)
(Now, error is coming higher as i have not included the source term in the problem. If you give that i hope error will come nice as exact solution and approximate solution are matching).

If you give me your equations or provide your source, i will check order of convergence.

Thanks in advance!.

[EChase (Ed)]

Thank you both very much!

My equations are all written down on paper, and I need to formalize them. I will post the math and ask you to look at it. Thank you for making the offer.

Thank you

[EChase (Ed)]

Math Notes:

Grad u= ∇u= ∂u/∂x i+ ∂u/∂y j+ ∂u/∂z k
div u= ∇∙u= ∂u/∂x+ ∂u/∂y+ ∂u/∂z
curl u= ∇×u=((∂u_3)/∂y- (∂u_2)/∂z,(∂u_3)/∂z- (∂u_1)/∂x,(∂u_2)/∂x- (∂u_1)/∂y)
curl u= ∇×u= (∂/∂x,∂/∂y,∂/∂z)×(u_1,u_2,u_3 )
∇×u= ((∂u_3)/∂y- (∂u_2)/∂z)i-((∂u_3)/∂z- (∂u_1)/∂x)j+((∂u_2)/∂x- (∂u_1)/∂y)k

Kronecker Delta=δ_ij= ■(1 if i=j@0 otherwise)

Relationship between stress tensor and strain tensor

stress tensor= σ_ij= λδ_ij ∇∙u+2μϵ_ij (u)
strain tensor= ϵ_ij= 1/2 (〖∂u〗_i/〖∂x〗_j +〖∂u〗_j/〖∂x〗_i )

Lame constants
λ= Eν/(1+ν)(1-2ν)
μ=E/2(1+ν)

Where:
E = Young’s Modulus
 = poisson’s ratio

Displacement formula : Navier-Lamé equations

(λ+μ) u_(k,ki)+μu_(i,kk)+F_i=ρ (∂^2 u_i)/〖∂t〗^2
(λ+μ)∇∇∙u+μ∇^2×u+F_i=ρ (∂^2 u_i)/〖∂t〗^2
(λ+μ)grad div u+μcurl curl u+F_i=ρ (∂^2 u_i)/〖∂t〗^2

For static conditions

ρ (∂^2 u_i)/〖∂t〗^2 =0

FreeFEM does not use the Navier-Lame equation for a structural problem, because the associated variational form does not give the right boundary conditions. Instead:

-div(σ)=f in Ω

As FreeFEM solves the div value to begin with, the matrix definition of  is required.

σ_ij= λδ_ij (∂u/∂x+ ∂u/∂y+ ∂u/∂z)+2/2 μ(〖∂u〗_i/〖∂x〗_j +〖∂u〗_j/〖∂x〗_i )

σ_11= λ(〖∂u〗_x/∂x)+2μ(〖∂u〗_x/∂x)
σ_12= μ(〖∂u〗_x/∂y+〖∂u〗_y/∂x)
σ_13= μ(〖∂u〗_x/∂z+〖∂u〗_z/∂x)
σ_22= λ( 〖∂u〗_y/∂y)+2μ(〖∂u〗_y/∂y)
σ_23= μ(〖∂u〗_y/∂z+〖∂u〗_z/∂y)
σ_33= λ( 〖∂u〗_z/∂z)+2μ(〖∂u〗_z/∂z)

These values are used in the div and epsilon matrices which are fed to the div and epsilon matrices which are fed to the solve function

epsilon=[■(〖∂u〗_x/∂x&〖∂u〗_x/∂y+〖∂u〗_y/∂x&〖∂u〗_x/∂z+〖∂u〗_z/∂x@〖∂u〗_x/∂y+〖∂u〗_y/∂x&〖∂u〗_y/∂y&〖∂u〗_y/∂z+〖∂u〗_z/∂y@〖∂u〗_x/∂z+〖∂u〗_z/∂x&〖∂u〗_y/∂z+〖∂u〗_z/∂y& 〖∂u〗_z/∂z)]

div=[■(〖∂u〗_x/∂x&0&0@0&〖∂u〗_y/∂y&0@0&0& 〖∂u〗_z/∂z)]
Note: The coefficients ,2 and are handles in the solve sequence. These matrices just separate the terms.

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]

// epsilon(u1,u2)’*epsilon(v1,v2) = epsilon(u): epsilon(v)

macro div(u,v,w) (dx(u)+dy(v)+dz(w))

and the solve function looks like this:
// Problem
solve lame([u, v, w], [uu, vv, ww])
= int3d(Sp)(lambda*(div(u, v, w) * div(uu, vv, ww)))
+ int3d(Sp)(2.*mu * ( epsilon(u,v,w)’ * epsilon(uu, vv, ww)))

…
Minus forces
…
On boundary conditions
…
;

lame; // solve

These are the notes that go with the script. It is a simple cantilevered beam structural statics model. I was wondering if the definitions are correct, especially that for epsilon. Does the order of the terms matter?