====================================================================

  programmes Maple  

  24 mars 2006 - 7 novembre 2006

  Méthode des éléments finis : élasticité à une dimension

     Yves DEBARD

     Institut Universitaire de Technologie du MANS
     Département Génie Mécanique et Productique
     Avenue Olivier Messiaen
     72085 LE MANS CEDEX 9

====================================================================

3n_int

# élément à 3 noeuds 

# fonctions d'interpolation 

restart:with(linalg):

x:=xi->a0+a1*xi+a2*xi*xi;
solve({x1=x(-1),x2=x(0),x3=x(1)},{a0,a1,a2}):assign(%):
N:=grad(x(xi),[x1,x2,x3]):
N:=map(factor,N);
dN:=map(diff,N,xi);

plot([N[1],N[2],N[3]],xi=-1..1,legend=[N1,N2,N3],
color=[red,blue,green],thickness=2,
title="Elément à 3 noeuds : fonctions d'interpolation");

-------------------------------------------------------------

3n_mat

# élément à 3 noeuds 

# calcul des matrices élémentaires

restart:with(linalg):

# représentation de la géométrie et jacobien

x:=(1+xi)*L/2;J:=L/2;

# fonctions d'interpolation

N:=vector([xi*(-1+xi)/2,-(-1+xi)*(xi+1),xi*(xi+1)/2]);
dN:=vector([-1/2+xi,-2*xi,xi+1/2]);

# matrice de rigidité

B:=scalarmul(dN,1/J);
k:=Matrix(3,3,(i,j)->int(B[i]*B[j]*E*A*J,xi=-1..1),shape=symmetric);

# matrice de masse

m:=Matrix(3,3,(i,j)->int(N[i]*N[j]*rho*A*J,xi=-1..1),shape=symmetric);

# vecteur force

px:=pxi*(1-xi)/2+pxj*(1+xi)/2;
f:=vector(3,i->int(N[i]*px*J,xi=-1..1)):f:=simplify(f);
f:=jacobian(f,[pxi,pxj]);

===========================================================

3n_milieu

# influence de la position du noeud milieu
# sur la performance d'un élément à 3 noeuds

restart:with(linalg):with(plots):

# fonctions d'interpolation

N:=[-xi*(1-xi)/2,1-xi^2,xi*(1+xi)/2];
dN:=[(2*xi-1)/2,-2*xi,(2*xi+1)/2];

# transformation géométrique et jacobien

x:=-N[1]*L/2+N[2]*x2+N[3]*L/2:x:=simplify(%);
J:=diff(x,xi);

# matrice [B]

B:=simplify([seq(dN[i]/J,i=1..3)]);

# calcul de la matrice de rigidité 
# par intégration numérique avec 2 points de Gauss

G1:=-sqrt(3)/3;G2:=-G1;
Kij:=(i,j)->subs(xi=G1,E*A*B[i]*B[j]*J)+subs(xi=G2,E*A*B[i]*B[j]*J);
K:=simplify(matrix(3,3,Kij));
Fi:=i->subs(xi=G1,p*N[i]*J)+subs(xi=G2,p*N[i]*J);
F:=simplify(vector(3,Fi));

# calcul du déplacement du noeud 2

U2:=simplify(F[2]/K[2,2]);

# champ de déplacements

u:=N[2]*U2;

# calcul des efforts normaux aux noeuds par la méthode 2

N1:=-simplify(K[1,2]*U2-F[1]);
N3:=simplify(K[3,2]*U2-F[3]);

# calcul des efforts normaux par la méthode 1

Nx:=A*E*B[2]*U2;

# représentations graphiques

L:=1;A:=1;E:=1;p:=1;

eq:=x2=0.15*L;

# champ de déplacements

plotuexact:=plot([subs(x2=0,x),subs(x2=0,u),xi=-1..1],
  title="Champ de déplacements u(x)",color=green,thickness=2):
plotu:=plot([subs(eq,x),subs(eq,u),xi=-1..1],color=blue,thickness=2):
plotx2:=plot([subs(eq,xi=0,x),xi,xi=0..p*L^2/8/E/A],color=red):
display(plotuexact,plotu,plotx2);

# efforts normaux : méthode 1 et méthode 2

plotNexact:=plot([subs(x2=0,x),subs(x2=0,Nx),xi=-1..1],
  title="Effort normal N(x)",color=green,thickness=2):
plotN:=plot([subs(eq,x),subs(eq,Nx),xi=-1..1],color=blue,thickness=2):
plotG1:=plot([subs(eq,xi=G1,x),xi,xi=0..0.5],color=red):
plotG2:=plot([subs(eq,xi=G2,x),xi,xi=0..-0.5],color=red):
display(plotNexact,plotN,plotG1,plotG2);