Elasticity problem

General Discussion
1

Bonjour cher Professeur

Ma compilation affiche je n’ai pas compris

– Square mesh : nb vertices =22 , nb triangles = 20 , nb boundary edges 22
– FESpace: Nb of Nodes 44 Nb of DoF 132
– Solve :
min -1.3237e-12 max 0.345492
current line = 68
Exec error : Try to get unset x,y, …
– number :1
Exec error : Try to get unset x,y, …
– number :1
err code 8 , mpirank 0
try getConsole C:\Users\ADMIN\Desktop\simulations\codee.edp
save log in : ‘C:\Users\ADMIN\Desktop\simulations\codee.log’
wait enter ?

voici mon code:

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],[u1oldd, u2oldd, u3oldd]=[sin(2pix)sin(2piy)sin(2piz), sin(2pix)sin(2piy)sin(2piz), 0]
,[u1old, u2old, u3old]=[sin(2pix)sin(2piy)sin(2piz), sin(2pix)sin(2piy)sin(2piz), 0], [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)) //

int kk=0;

// Problem
real mu = E/(2*(1+sigma));
real lambda = Esigma/((1+sigma)(1-2*sigma));
real dt=0.00000001, beta=0,T=0.000001, alpha=0;

real[int] instT(floor(T/dt) + 1);
real[int] NORML2(floor(T/dt) + 1);

real[int] NORML(floor(T/dt) + 1);

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

solve Lame ([u1, u2, u3], [v1, v2, v3])
= int3d(Th)((rho/dt^2)[u1, u2, u3]'[v1, v2, v3]
+ 0.5lambdadiv(u1, u2, u3)div(v1, v2, v3)
+0.52.mu( epsilon(u1, u2, u3)'*epsilon(v1, v2, v3) )

              )
        +int2d(Th,1,2,3,4,6)((beta/dt)*[u1old, u2old, u3old]'*[v1, v2, v3]) 
    
        -int3d(Th)( (alpha/dt)*([u1old, u2old, u3old]'*[v1, v2, v3])
                    + 2*(rho/dt^2)*[u1old, u2old, u3old]'*[v1, v2, v3]
                    -(rho/dt^2)*[u1oldd, u2oldd, u3oldd]'*[v1, v2, v3]
                    - 0.5* lambda*div(u1oldd, u2oldd, u3oldd)*div(v1, v2, v3)
                    - 0.5*2.*mu*( epsilon(u1oldd, u2oldd, u3oldd)'*epsilon(v1, v2, v3) )
                     +(alpha/dt)*[u1old, u2old, u3old]'*[v1, v2, v3]
                 )

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

cout << " " << u1old<< u2old<< u3old<< u1<< u2 <<endl;

[up1,up2,up3]=[u1old-u1oldd, u2old-u2oldd, u3old-u3oldd];

instT[kk]=t;
NORML2[kk]=int3d(Th)((rho/dt^2)([up1,up2,up3]'[up1, up2, up3]))
+ int3d(Th)( lambda*div(u1, u2, u3)*div(u1, u2, u3)
+ 2.mu(epsilon(u1, u2, u3)'*epsilon(u1, u2, u3) ) )
;

NORML[kk]= int2d(Th,1,2,3,4,6)((beta/dt)[u1, u2, u3]'[u1, u2, u3]) ;

[u1oldd, u2oldd, u3oldd]=[u1old, u2old, u3old];
[u1old, u2old, u3old]=[u1, u2, u3];

cout << " " << NORML2[kk]<< " " << instT[kk] <<endl;

kk++;

}

cout << " " << NORML2 <<endl;

J’aimerais envoyer le fichier.loq mais ça dit que je suis nouveau je ne peux pas envoyer un fichier.
S’il ya lien pour faire l’envoie je serai ravi.

2

did you try making ny and nz bigger like the 10 I suggested?

3

No i will try it now to see.

4

that alone won’t fix your problem see if this helps,

code_freefem.edp.edp (2.3 KB)

5

I increased but I still observe the same problem where the numerical results do not reflect the theoretical results.

6

no the problem still persists

7

If you ran the last code I posted do you get a result that quickly
decreases towards zero? Now its just a matter of getting the numbers
right? Just looking at the results they could be exponentially decreasing
which may be what you expect except for the time constant being wrong.
It may help to post the theory too.

8

\Omega est un ouvert de \mathbb{R}^3
\begin{equation}
\begin{array}{ll}
\rho \frac{\partial^2 u}{\partial t^2}- div (Ce(u))=0\\
u=0 \mbox{sur} \Gamma_D\\
Ce(u).n= -\beta u~~\mbox{sur} ~~\Gamma_N
\end{array}
\end{equation}

\texbf{Formulation variationnelle}

\begin{equation}
\displaystyle \int_\Omega \rho \frac{\partial^2 u}{\partial t^2}\cdot v+ \displaystyle \int_\Omega Ce(u): Ce(v)= -\beta\displaystyle \int_{\Gamma_N} \frac{\partial u}{\partial t} \cdot v
\end{equation}

\textbf{\estimation d’Energie}

On pose
Energie = \displaystyle \int_\Omega \rho| \frac{\partial u}{\partial t}|^2 + \displaystyle \int_\Omega Ce(u): Ce(u)
\begin{equation}
\displaystyle \frac{\partial }{\partial t} \left(
\displaystyle \int_\Omega \rho| \frac{\partial u}{\partial t}|^2 + \displaystyle \int_\Omega Ce(u): Ce(u)
\right)= -\beta\displaystyle \int_{\Gamma_N}| \frac{\partial u}{\partial t}|^2
\end{equation}

Donc ici on a une décroissance simple et non une décroissance exponentielle.

Deplus lorsque \beta=0 la dérivée de l’énergie est nulle par conséquent l’énergie doit être une constante.
C’est ce que dit la theory

9

Oui je suis censée obtenu que l’énergie est constant.

Par besoin d’une décroissance.

10

Can you describe the error or post the working code? I was curious now.
Thanks.

11

I haven’t found a good code yet… I myself cannot detect an error in the code that you published… I don’t know why but the result is not correct… when the parameters beta =0 the energy is theoretically constant.but this code prove that the energy are decrease.

12

Your kinetic term doesn’t seem to have a dt in it. If I reduce the time
step to .0001 or so it still loses some energy. I guess you could check the
time and length scales and see if the mesh and time step are fine enough
or see how the energy loss varies with those parameters.
Apparently your initial conditions are at zero derivative or maximum deflection
with no instaneous motion.