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

  programmes MuPAD  -  24 mars 2006 - 9 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 

reset():export(linalg):

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

plotfunc2d(N[1],N[2],N[3],xi=-1..1);

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

3n_mat

// élément à 3 noeuds 

// calcul des matrices élémentaires

reset():export(linalg):

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

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

// fonctions d'interpolation

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

// matrice de rigidité

B:=dN/J;
k:=matrix(3,3,(i,j)->int(B[i]*B[j]*EA*J,xi=-1..1));

// matrice de masse

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

// vecteur force

px:=pxi*(1-xi)/2+pxj*(1+xi)/2;
f:=matrix(3,1,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

reset():export(linalg):

// 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([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(EA*B[i]*B[j]*J,xi=G1)+subs(EA*B[i]*B[j]*J,xi=G2);
K:=simplify(matrix(3,3,Kij));
Fi:=i->subs(p*N[i]*J,xi=G1)+subs(p*N[i]*J,xi=G2);
F:=simplify(matrix(3,1,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:=EA*B[2]*U2;

// représentations graphiques

L:=1;EA:=1;p:=1;

eq:=x2=0.15*L;

// champ de déplacements

plotuexact:=plot::Curve2d([subs(x,x2=0),subs(u,x2=0)],xi=-1..1):
plotu:=plot::Curve2d([subs(x,eq),subs(u,eq)],xi=-1..1):
plotx2:=plot::Curve2d([subs(x,eq,xi=0),xi],xi=0..p*L^2/8/EA):
plot(plotuexact,plotu,plotx2);

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

plotNexact:=plot::Curve2d([subs(x,x2=0),subs(Nx,x2=0)],xi=-1..1):
plotN:=plot::Curve2d([subs(x,eq),subs(Nx,eq)],xi=-1..1):
plotG1:=plot::Curve2d([subs(x,eq,xi=G1),xi],xi=0..0.5):
plotG2:=plot::Curve2d([subs(x,eq,xi=G2),xi],xi=0..-0.5):
plot(plotNexact,plotN,plotG1,plotG2);