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

  programmes Maple  -  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

restart;with(linalg):

# 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(%):
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(%):
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(%):
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(%):
k:=jacobian([-Tzi,Tzj],[wi,wj]);
f:=jacobian([Tzi,-Tzj],[pzi,pzj,myi,myj]);
rotyi;rotyj;

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

ber_int

# calcul des fonctions d'interpolation

restart:with(linalg):

# 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(x=0,w):eq2:=wj=subs(x=L,w):
eq3:=rotyi=subs(x=0,roty):eq4:=rotyj=subs(x=L,roty):
solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(%):

# 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

plot([subs(L=1,Nw[1]),subs(L=1,Nw[2]),subs(L=1,Nw[3]),subs(L=1,Nw[4])],
x=0..1,legend=[N1,N2,N3,N4],color=[red,blue,green,cyan],thickness=2,
title="Poutre de Bernoulli : fonctions d'interpolation");

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

ber_int_par

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

restart:with(linalg):

# 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(xi=-1,w):eq2:=wj=subs(xi=1,w):
eq3:=rotyi=subs(xi=-1,roty):eq4:=rotyj=subs(xi=1,roty):
solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(%):

# 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(factor,Nw));
Nroty:=simplify(map(factor,Nroty));
dNroty:=simplify(map(factor,dNroty));

# représentation graphique pour L=1

plot([Nw[1],subs(L=1,Nw[2]),Nw[3],subs(L=1,Nw[4])],xi=-1..1,
legend=[N1,N2,N3,N4],color=[red,blue,green,cyan],thickness=2,
title="Poutre de Bernoulli : fonctions d'interpolation");

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

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

restart:with(linalg):

# 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),shape=symmetric);

# matrice de masse

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

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

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

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

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