Numerical results error

General Discussion
1

Hello everyone
I modified the code of the elasticity equation in 3D to evaluate the evolution of energy as a function of time…but I don’t know why it does not reflect what I obtained theoretically.

Normally when I put the beta=0 parameter, the energy must be constant and if beta is strictly greater than 0 the energy decreases.

load “msh3”
load “medit”
int nx =10, ny =1, nz = 1;

int[int] rup = [0,6]; // Haut du cube étiqueté par zéro
int[int] rdown = [0,5]; // Bas du cube étiqueté par zéro
int[int] rmid = [1,2,3,4]; // Parois latérales étiquetées

mesh3 Th = buildlayers(square(nx, ny, [0.5x, 0.4y]), nz, zbound=[0., 0.1],
labelmid=rmid, labelup = rup, labeldown = rdown);

real E = 21.5e4,rho=1600;
real sigma = 0.29;
real gravity = -0.05;
// Fespace
fespace Vh(Th, [P1, P1, P1]);
Vh [u1, u2, u3]=[10sin(2pix)sin(2piy)sin(2piz),10 sin(2pix)sin(2piy)sin(2piz), 10]
,[u1old, u2old, u3old]=[10sin(2pix)sin(2piy)sin(2piz),10 sin(2pix)sin(2piy)sin(2piz), 10],[u1oldd, u2oldd, u3oldd], [v1, v2, v3],[up1,up2,up3];

// 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] //
macro div(u1, u2, u3) (dx(u1) + dy(u2) + dz(u3)) //

// Problem
int kk=0;
real mu = E/(2*(1+sigma));
real lambda = Esigma/((1+sigma)(1-2*sigma));
real dt=0.1, beta=0,T=1;
real[int] instT(floor(T/dt) + 1);
real[int] NORML2(floor(T/dt) + 1);

for(real t=0; t <=T; t+=dt )
{

[u1oldd, u2oldd, u3oldd]=[u1old, u2old, u3old];
[u1old, u2old, u3old]=[u1, u2, u3];
solve Lame ([u1, u2, u3], [v1, v2, v3])
= int3d(Th)((rho/dt^2)[u1, u2, u3]'[v1, v2, v3]
+ lambdadiv(u1, u2, u3)div(v1, v2, v3)
+ 2.mu( epsilon(u1, u2, u3)'epsilon(v1, v2, v3) )
)
-int3d(Th)(2(rho/dt^2)[u1old, u2old, u3old]'[v1, v2, v3]
-(rho/dt^2)[u1oldd, u2oldd, u3oldd]'[v1, v2, v3])

     +int2d(Th,1,2,3,4,6)((beta/dt)*[u1, u2, u3]'*[v1, v2, v3]) 
     -int2d(Th,1,2,3,4,6)( (beta/dt)*[u1old, u2old, u3old]'*[v1, v2, v3])
     + on(5, u1=0, u2=0, u3=0) ;

[up1,up2,up3]=[u1-u1old, u2-u2old, u3-u3old];

NORML2[kk]=int3d(Th)(rho*([up1,up2,up3]‘[up1, up2, up3]))
+ int3d(Th)( lambdadiv(u1, u2, u3)*div(u1, u2, u3)
+ 2.mu(epsilon(u1, u2, u3)’*epsilon(u1, u2, u3) ) )
;

cout << " " << NORML2[kk] <<endl;
}
But I don’t know what’s wrong with the code

Thanks for help

2

can you upload the file? copy/paste always has errors.
I guess generally I’d try to breakdown the contributions to total energy
and see what they are doing. It may be something simple like a sign error.
Also you could have hidden sources or sinks
of energy at boundaries etc.

3

Ok this is the attached the source code

(Attachment code_freefem1 is missing)

4

the post says attachment is miessing. Can you put it on github or something?

5

ok i resend

code_freefem.edp.edp (2.23 KB)

6

this is the tex

(Attachment code.txt is missing)

7

code_freefem.edp.edp (2.23 KB)

8

Est ce que vous avez reçu le fichier.
Je n’arrive pas a envoyé sur github.

9

I didn’t get to the math as your value of concern seemed to drop
once I increased the numbers ny and nz and moved the variable kk
init and increment,

code_freefem.edp.edp (2.3 KB)