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

  programmes MuPAD  -  24/03/06

  Méthode des éléments finis :
  
  Poutre de Bernoulli : flexion dans le plan xz

     Yves DEBARD

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

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

bre_mat

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

reset();export(linalg):

assume(L>0):assume(EIy>0):

// charges

pz:=x->pzi+(pzj-pzi)*x/L;
my:=x->myi+(myj-myi)*x/L;

// effort tranchant et moment fléchissant

Tz:=x->Tzi-int(pz(s),s=0..x);
Mfy:=x->Mfyi+int(Tz(s),s=0..x)-int(my(s),s=0..x);

// pente et flèche

roty:=x->rotyi+int(Mfy(s)/EIy,s=0..x);
w:=x->wi-int(roty(s),s=0..x);

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

solve({rotyj=roty(L),wj=w(L)},{Tzi,Mfyi}):assign(op(%)):
Tzj:=Tz(L):Mfyj:=Mfy(L):

// matrice de rigidité et vecteur force

unod:=[wi,rotyi,wj,rotyj]:
fnod:=[-Tzi,-Mfyi,Tzj,Mfyj]:

k:=jacobian(fnod,unod);
f:=-jacobian(fnod,[pzi,pzj,myi,myj]);

// remarque : fonctions d'interpolation 

Nw:=grad(w(x),unod);
Nroty:=grad(roty(x),unod);

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

ber_rot_rig : élément rotule-rigide

Ajouter les lignes suivantes au programme ber_mat :

// rotule-rigide

solve({Mfyi},{rotyi}):assign(op(%)):
k:=jacobian([-Tzi,Tzj,Mfyj],[wi,wj,rotyj]);
f:=jacobian([Tzi,-Tzj,-Mfyj],[pzi,pzj,myi,myj]);
simplify(rotyi);

ber_rig_rot : élément rigide-rotule

Ajouter les lignes suivantes au programme ber_mat :

// rigide-rotule

solve({Mfyj},{rotyj}):assign(op(%)):
k:=jacobian([-Tzi,-Mfyi,Tzj],[wi,rotyi,wj]);
f:=jacobian([Tzi,Mfyi,-Tzj],[pzi,pzj,myi,myj]);
rotyj;

ber_rot_rot : élément rotule-rotule

Ajouter les lignes suivantes au programme bre_mat :

// rotule-rotule

solve({Mfyi,Mfyj},{rotyi,rotyj}):assign(op(%)):
k:=jacobian([-Tzi,Tzj],[wi,wj]);
f:=jacobian([Tzi,-Tzj],[pzi,pzj,myi,myj]);
rotyi;rotyj;

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

ber_int

// calcul des fonctions d'interpolation

reset():export(linalg):

assume(L>0):

// champ de déplacements

w:=a0+a1*x+a2*x^2+a3*x^3;
roty:=diff(-w,x);
droty:=diff(roty,x):

// calcul des coefficients du polynôme d'interpolation

eq1:=wi=subs(w,x=0):eq2:=wj=subs(w,x=L):
eq3:=rotyi=subs(roty,x=0):eq4:=rotyj=subs(roty,x=L):
solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(op(%)):

// fonctions d'interpolation et dérivées

unod:=[wi,rotyi,wj,rotyj]:

Nw:=grad(w,unod);
Nroty:=grad(roty,unod);
dNroty:=grad(droty,unod);

// représentation graphique pour L=1

plotfunc2d(subs(Nw[1],L=1),subs(Nw[2],L=1),subs(Nw[3],L=1),subs(Nw[4],L=1),x=0..1);

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

ber_int_par

// calcul des fonctions d'interpolation sous forme paramétrique

reset():export(linalg):

assume(L>0):

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

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

// représentation du champ de déplacements

w:=a0+a1*xi+a2*xi^2+a3*xi^3;
roty:=diff(-w,xi)/J;
droty:=diff(roty,xi)/J;

// calcul des coefficients du polynôme d'interpolation

eq1:=wi=subs(w,xi=-1):eq2:=wj=subs(w,xi=1):
eq3:=rotyi=subs(roty,xi=-1):eq4:=rotyj=subs(roty,xi=1):
solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(op(%)):

// fonctions d'interpolation et dérivées

unod:=[wi,rotyi,wj,rotyj]:

Nw:=grad(w,unod):
Nroty:=grad(roty,unod):
dNroty:=grad(droty,unod):

Nw:=simplify(map(Nw,factor));
Nroty:=simplify(map(Nroty,factor));
dNroty:=simplify(map(dNroty,factor));

// représentation graphique pour L=1

plotfunc2d(Nw[1],subs(Nw[2],L=1),Nw[3],subs(Nw[4],L=1),xi=-1..1);

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

ber_interpolation

// fonctions d'interpolation sous forme paramétrique

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

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

// fonctions d'interpolation

Nw:=[(2+xi)*(xi-1)^2/4,-L*(1+xi)*(xi-1)^2/8,
        -(xi-2)*(1+xi)^2/4,-L*(xi-1)*(1+xi)^2/8]:
Nroty:=[-3/2*(-1+xi^2)/L,(3*xi+1)*(xi-1)/4,
        3/2*(-1+xi^2)/L,(1+xi)*(3*xi-1)/4]:
dNroty:=[-6*xi/L^2,(3*xi-1)/L,6*xi/L^2,(3*xi+1)/L]:

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

ber_int_mat

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

reset():export(linalg):

assume(L>0):assume(EIy>0):

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

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

// fonctions d'interpolation

Nw:=[(2+xi)*(xi-1)^2/4,-L*(1+xi)*(xi-1)^2/8,
        -(xi-2)*(1+xi)^2/4,-L*(xi-1)*(1+xi)^2/8]:
Nroty:=[-3/2*(-1+xi^2)/L,(3*xi+1)*(xi-1)/4,
        3/2*(-1+xi^2)/L,(1+xi)*(3*xi-1)/4]:
dNroty:=[-6*xi/L^2,(3*xi-1)/L,6*xi/L^2,(3*xi+1)/L]:

// matrice de rigidité

k:=matrix(4,4,(i,j)->int(dNroty[i]*dNroty[j]*EIy*J,xi=-1..1));

// matrice de masse

mw:=matrix(4,4,(i,j)->int(Nw[i]*Nw[j]*rho*A*J,xi=-1..1));
mroty:=matrix(4,4,(i,j)->int(Nroty[i]*Nroty[j]*rho*Iy*J,xi=-1..1));

pz:=pzi+(pzj-pzi)*x/L:

fpz:=matrix(4,1,i->int(Nw[i]*pz*J,xi=-1..1)):
fpz:=simplify(jacobian(fpz,[pzi,pzj]));

my:=myi+(myj-myi)*x/L:

fmy:=matrix(4,1,i->int(Nroty[i]*my*J,xi=-1..1)):
fmy:=simplify(jacobian(fmy,[myi,myj]));