## Coordonnées Cartésiennes

### Fonctions utiles

U : le déplacement (vecteur) ; X : les coordonnées (vecteur) ;
lambda, mu : les coefficients de Lamé ;
N : le vecteur normal à une surface ;

div(U,X) : divergence de U par rapport à X (scalaire) ;
Grad(U,X) : gradient de U par rapport à X (matrice) ;
grad(phi,X) : gradient de phi par rapport à X (vecteur) ;
Lap(U,X) : laplacien vectoriel de U par rapport à X (vecteur) ;

Navier(U,X,lambda,mu) : equation de Navier (vecteur) ;
Def(U,X) : tenseur des déformation (matrice) ;
tr(t1) : trace de la matrice t1 ;
Hooke(epsilon,lambda,mu) : expression du tenseur des contraintes en fonction du tenseur des déformations epsilon ;
dF(sigma,N) : produit du tenseur des contraintes sigma par le vecteur normal N;

div(U,X):=apply("+",map(lambda([%u,%v],diff(%u,%v)),U,X))$
Grad(U,X):=apply(matrix,map(lambda([%u],map(lambda([%x],diff(%u,%x)),X)),U))$
grad(phi,X):=map(lambda([%u],diff(phi,%u)),X)$
Lap(U,X):=map(lambda([%u],div(grad(%u,X),X)),U)$

Navier(U,X,lambda,mu):=(lambda+mu)*grad(div(U,X),X)+mu*Lap(U,X)$
Def(U,X):=block([eps:Grad(U,X)],eps:(eps+transpose(eps))/2)$
tr(t1):=block([i:0],apply("+",map(lambda([u],i:i+1,u[i]),substpart("[",t1,0))))$
Hooke(epsilon,lambda,mu):=block([unite:matrix([1,0,0],[0,1,0],[0,0,1])],
 lambda*tr(epsilon)*unite+2*mu*epsilon)$
dF(sigma,N):=sigma.N$

LvsYP:[lambda=nu*E/(1-2*nu)/(1+nu),mu=E/(2*(1+nu))]$
YPvsL:[E=mu*(3*lambda+2*mu)/(lambda+mu),nu=lambda/(2*(lambda+mu))]$

DitCeQueCest([s]):=apply(print,append(s,[%]))$DCQ:DitCeQueCest$

### Traction d'une éprouvette   

La forme de U paramétrée par alpha,beta est devinée ; 
ces coefficients sont calculés pour que les équations d'équilibre soient satisfaites exactement.

U:[-alpha*x,-alpha*y,beta*z]$DCQ("Forme devinée pour U : ")$

X:[x,y,z]$
eq:Navier(U,X,lambda,mu)$DCQ("Eq. de Navier pour ce U : ")$

epsilon:Def(U,X)$DCQ("epsilon pour ce U : ")$

sigma:Hooke(epsilon,lambda,mu)$DCQ("sigma pour U : ")$

cll:ratsimp(dF(sigma,[nx,ny,0]))$DCQ("(sigma n) sur Gammal pour U : ")$

clb:ratsimp(dF(sigma,[0,0,1]))$DCQ("(sigma n) sur Gamma1 pour U : ")$

sol:solve([cll[1,1]=0,clb[3,1]=F/S],[alpha,beta])[1]$
DCQ("Résolution : ")$
print("de cll=0 et clb = F/S par rapport à alpha et beta")$

ratsimp(ev(sigma,sol))$DCQ("sigma =")$

factor(ev(ev(epsilon,sol),LvsYP))$DCQ("epsilon =")$

### Torsion d'une éprouvette

La forme de U paramétrée par alpha, phi(x,y) est devinée ; 
Les équations d'équilibre soient satisfaites exactement.

U:[-alpha*z*y,alpha*z*x,alpha*phi(x,y)]$DCQ("U deviné : ")$

X:[x,y,z]$
eq:Navier(U,X,lambda,mu)$DCQ("Eq. Navier pour ce U : ")$

epsilon:Def(U,X)$DCQ("epsilon pour ce U : ")$

sigma:Hooke(epsilon,lambda,mu)$DCQ("sigma pour U : ")$

cll:ratsimp(dF(sigma,[nx,ny,0]))$DCQ("(sigma n) sur Gammal pour U : ")$

clb:ratsimp(dF(sigma,[0,0,1]))$DCQ("(sigma n) sur Gamma1 pour U : ")$

Mb:clb[2,1]*x-clb[1,1]*y$DCQ("Densité de couple à appliquer=")$

### Flexion d'une éprouvette

U:[M/(2*E*I)*(z^2+nu*(x^2-y^2)),nu*M/(E*I)*x*y,-M/(E*I)*x*z];
X:[x,y,z]$
[lambda,mu]:ev([lambda,mu],LvsYP)$
eq:ratsimp(Navier(U,X,lambda,mu));
epsilon:Def(U,X);
sigma:ratsimp(Hooke(epsilon,lambda,mu));
cll:ratsimp(dF(sigma,[nx,ny,0]));
clb:ratsimp(dF(sigma,[0,0,1]));

## Coordonnées cylindriques

Idem cartésien

div(U):=diff(U[1],r)+U[1]/r+diff(U[2],theta)/r+diff(U[3],z)$
Grad(U):=matrix(
 [diff(U[1],r),diff(U[1],theta)/r-U[2]/r,diff(U[1],z)],
 [diff(U[2],r),diff(U[2],theta)/r+U[1]/r,diff(U[2],z)],
 [diff(U[3],r),diff(U[3],theta)/r       ,diff(U[3],z)])$
grad(phi):=[diff(phi,r),diff(phi,theta)/r,diff(phi,z)]$
lap(phi):=div(grad(phi))$
Lap(U):=[lap(U[1])-2*diff(U[2],theta)/r^2-U[1]/r^2,
         lap(U[2])+2*diff(U[1],theta)/r^2-U[2]/r^2,
         lap(U[3])]$
Def(U):=block([eps:Grad(U)],eps:(eps+transpose(eps))/2)$
Navier(U,lambda,mu):=(lambda+mu)*grad(div(U))+mu*Lap(U)$
tr(t1):=block([i:0],apply("+",map(lambda([u],i:i+1,u[i]),substpart("[",t1,0))))$
Hooke(U,lambda,mu):=block([eps:Def(U),unite:matrix([1,0,0],[0,1,0],[0,0,1])],
 lambda*tr(eps)*unite+2*mu*eps)$
dF(sigma,N):=sigma.N$

LvsYP:[lambda=nu*E/(1-2*nu)/(1+nu),mu=E/(2*(1+nu))]$
YPvsL:[E=mu*(3*lambda+2*mu)/(lambda+mu),nu=lambda/(2*(lambda+mu))]$

DitCeQueCest([s]):=apply(print,append(s,[%]))$DCQ:DitCeQueCest$

### Tube sous pression

U:[u,0,0];
depends(u,r);
eq:ratsimp(Navier(U,lambda,mu));
sigma:ratsimp(Hooke(U,lambda,mu));
cll:ratsimp(dF(sigma,[1,0,0]));
eql:ratsimp(ev(cll[1,1],u=alpha*r+beta/r,nouns));
sol:solve([ev(eql=-dP,r=R1),ev(eql=0,r=R2)],[alpha,beta]);
ev(alpha*r+beta/r,sol);
ev(e/2*r+A/r,A=dP*R1^2*R2^2/(R2^2-R1^2)/(2*mu),e=dP*R1^2/(R2^2-R1^2)/(lambda+mu));