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

  programmes MuPAD  -  24/03/06

  Méthode des éléments finis : poutre soumise à un effort normal

     Yves DEBARD

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

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

mat_elem1

// calcul de la matrice de rigidité et du vecteur force 
// d'un élément de poutre à section constante

reset():export(linalg):

assume(L>0):

// charges

px:=x->pxi+(pxj-pxi)*x/L;

// effort normal

N:=x->Ni-int(px(s),s=0..x);N(x);

// champ de déplacements

u:=x->ui+int(N(s)/EA+alpha*DT,s=0..x);u(x);

// calcul des efforts nodaux en fonction des déplacements nodaux

solve({uj=u(L)},{Ni}):assign(op(%)):
Nj:=N(L):

// matrice de rigidité

k:=jacobian([-Ni,Nj],[ui,uj]);

// vecteur force

f:=jacobian([Ni,-Nj],[pxi,pxj,DT]);

// remarque : fonctions d'interpolation

Nu:=grad(u(x),[ui,uj]);


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

mat_elem2

// calcul des matrices élémentaires
// d'un élément de poutre à section constante
// à l'aide des fonctions d'interpolation

reset():export(linalg):

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

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

// fonctions d'interpolation

Nu:=[(1-xi)/2,(1+xi)/2];

// matrice de rigidité

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

// matrice de masse

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

// vecteur force

px:=pxi+(pxj-pxi)*x/L:
f:=matrix(2,1,i->int(Nu[i]*px*J,xi=-1..1)):simplify(f);
fth:=matrix(2,1,i->int(B[i]*EA*alpha*DT*J,xi=-1..1));

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

mat_var1

// calcul de la matrice de rigidité et du vecteur force 
// d'un élément de poutre à section variable

reset():export(linalg):

assume(L>0):assume(c>0):assume(alpha>0):

// charges

px:=x->pxi+(pxj-pxi)*x/L;

// effort normal

N:=x->Ni-int(px(s),s=0..x);

// champ de déplacements

A:=x->c^2*(1+x/L)^2;
u:=x->ui+int(N(s)/E_/A(s)+alpha*DT,s=0..x,Continuous);

// calcul des efforts nodaux en fonction des déplacements nodaux

solve({uj=u(L)},{Ni}):assign(op(%)):
Nj:=N(L):

// matrice de rigidité

k:=jacobian([-Ni,Nj],[ui,uj]);

// vecteur force

f:=jacobian([Ni,-Nj],[pxi,pxj,DT]):simplify(f);


mat_var2

// calcul des matrices élémentaires
// d'un élément de poutre à section variable
// à l'aide des fonctions d'interpolation

reset():export(linalg):

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

x:=(1+xi)*L/2;J:=L/2;
Nu:=[(1-xi)/2,(1+xi)/2];

A:=c^2*(1+x/L)^2:

// matrice de rigidité

B:=[-1/L,1/L];
k:=matrix(2,2,(i,j)->int(B[i]*B[j]*E_*A*J,xi=-1..1));

// matrice de masse

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

// vecteur force

px:=pxi+(pxj-pxi)*x/L:
f:=matrix(2,1,i->int(Nu[i]*px*J,xi=-1..1)):simplify(f);
fth:=matrix(2,1,i->int(B[i]*E_*A*alpha*DT*J,xi=-1..1));

mat_var3

// calcul de la matrice de rigidité et du vecteur force 
// d'un élément de poutre à section variable
// à l'aide du théorème de Castigliano

reset():export(linalg):

A:=c^2*(1+x/L)^2;

// matrice de rigidité

C:=int(1/E_/A,x=0..L,Continuous):
kD:=1/C:
k:=matrix(2,2,[[kD,-kD],[-kD,kD]]);

// vecteur force

px:=x->pxi+(pxj-pxi)*x/L:
Fpx:=x->int(px(s),s=x..L):
fR:=simplify(Fpx(0));
upD:=int(Fpx(x)/E_/A,x=0..L,Continuous):

fD:=kD*upD;
f:=matrix(2,1,[fR-fD,fD]):
f:=jacobian(f,[pxi,pxj]):simplify(f);

fthD:=kD*int(alpha*DT,x=0..L);