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

  programmes Maple  -  24/03/06

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

     Yves DEBARD

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

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

ber_mat

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

restart;with(linalg):

# charges

py:=x->pyi+(pyj-pyi)*x/L;
mz:=x->mzi+(mzj-mzi)*x/L;

# effort tranchant et moment fléchissant

Ty:=x->Tyi-int(py(s),s=0..x);
Mfz:=x->Mfzi-int(Ty(s),s=0..x)-int(mz(s),s=0..x);

# rotation des sections autour de l'axe z et déplacement suivant y

rotz:=x->rotzi+int(Mfz(s)/EIz,s=0..x);rotz(x);
v:=x->vi+int(rotz(s),s=0..x);v(x);

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

solve({rotzj=rotz(L),vj=v(L)},{Tyi,Mfzi}):assign(%):
Tyj:=Ty(L):Mfzj:=Mfz(L):

# matrice de rigidité et vecteur force

unod:=[vi,rotzi,vj,rotzj]:
fnod:=[-Tyi,-Mfzi,Tyj,Mfzj]:

k:=jacobian(fnod,unod);
f:=jacobian(-fnod,[pyi,pyj,mzi,mzj]);

# remarque : fonctions d'interpolation 

Nv:=grad(v(x),unod);
Nrotz:=grad(rotz(x),unod);

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

ber_rot_rig : élément rotule-rigide

Ajouter les lignes suivantes au programme BER_MAT :

# rotule-rigide

solve({Mfzi},{rotzi}):assign(%):
k:=jacobian([-Tyi,Tyj,Mfzj],[vi,vj,rotzj]);
f:=jacobian([Tyi,-Tyj,-Mfzj],[pyi,pyj,mzi,mzj]);
simplify(rotzi);

ber_rig_rot : élément rigide-rotule

Ajouter les lignes suivantes au programme BER_MAT :

# rigide-rotule

solve({Mfzj},{rotzj}):assign(%):
k:=jacobian([-Tyi,-Mfzi,Tyj],[vi,rotzi,vj]);
f:=jacobian([Tyi,Mfzi,-Tyj],[pyi,pyj,mzi,mzj]);
rotzj;

ber_rot_rot : élément rotule-rotule

Ajouter les lignes suivantes au programme BER_MAT :

# rotule-rotule

solve({Mfzi,Mfzj},{rotzi,rotzj}):assign(%):
k:=jacobian([-Tyi,Tyj],[vi,vj]);
f:=jacobian([Tyi,-Tyj],[pyi,pyj,mzi,mzj]);
rotzi;rotzj;

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

ber_int

# calcul des fonctions d'interpolation

restart:with(linalg):

# champ de déplacements

v:=a0+a1*x+a2*x^2+a3*x^3;
rotz:=diff(v,x);
drotz:=diff(rotz,x):

# calcul des coefficients du polynôme d'interpolation

eq1:=vi=subs(x=0,v):eq2:=vj=subs(x=L,v):
eq3:=rotzi=subs(x=0,rotz):eq4:=rotzj=subs(x=L,rotz):
solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(%):

# fonctions d'interpolation et dérivées

unod:=[vi,rotzi,vj,rotzj]:

Nv:=grad(v,unod);
Nrotz:=grad(rotz,unod);
dNrotz:=grad(drotz,unod);

# représentation graphique pour L=1

plot([subs(L=1,Nv[1]),subs(L=1,Nv[2]),subs(L=1,Nv[3]),subs(L=1,Nv[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

v:=a0+a1*xi+a2*xi^2+a3*xi^3;
rotz:=diff(v,xi)/J;
drotz:=diff(rotz,xi)/J:

# calcul des coefficients du polynôme d'interpolation

eq1:=vi=subs(xi=-1,v):eq2:=vj=subs(xi=1,v):
eq3:=rotzi=subs(xi=-1,rotz):eq4:=rotzj=subs(xi=1,rotz):
solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(%):

# fonctions d'interpolation et dérivées

unod:=[vi,rotzi,vj,rotzj]:

Nv:=grad(v,unod):
Nrotz:=grad(rotz,unod):
dNrotz:=grad(drotz,unod):

Nv:=simplify(map(factor,Nv));
Nrotz:=simplify(map(factor,Nrotz));
dNrotz:=simplify(map(factor,dNrotz));

# représentation graphique pour L=1

plot([Nv[1],subs(L=1,Nv[2]),Nv[3],subs(L=1,Nv[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

Nv:=[(xi+2)*(xi-1)^2/4,L*(xi+1)*(xi-1)^2/8,
      -(xi-2)*(xi+1)^2/4,L*(xi-1)*(xi+1)^2/8]:
Nrotz:=[3/2*(xi^2-1)/L,(3*xi+1)*(xi-1)/4,
       -3/2*(xi^2-1)/L,(xi+1)*(3*xi-1)/4]:
dNrotz:=[6/L^2*xi,1/L*(3*xi-1),-6/L^2*xi,1/L*(3*xi+1)]:

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

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

Nv:=[(xi+2)*(xi-1)^2/4,L*(xi+1)*(xi-1)^2/8,
      -(xi-2)*(xi+1)^2/4,L*(xi-1)*(xi+1)^2/8]:
Nrotz:=[3/2*(xi^2-1)/L,(3*xi+1)*(xi-1)/4,
       -3/2*(xi^2-1)/L,(xi+1)*(3*xi-1)/4]:
dNrotz:=[6/L^2*xi,1/L*(3*xi-1),-6/L^2*xi,1/L*(3*xi+1)]:

# matrice de rigidité

k:=Matrix(4,4,(i,j)->int(dNrotz[i]*dNrotz[j]*EIz*J,xi=-1..1),shape=symmetric);

# matrice de masse

mv:=Matrix(4,4,(i,j)->int(Nv[i]*Nv[j]*rho*A*J,xi=-1..1),shape=symmetric);
mrotz:=Matrix(4,4,(i,j)->int(Nrotz[i]*Nrotz[j]*rho*Iz*J,xi=-1..1),shape=symmetric);

# vecteur force

py:=pyi+(pyj-pyi)*x/L:

fpy:=vector(4,i->int(Nv[i]*py*J,xi=-1..1)):
fpy:=simplify(jacobian(fpy,[pyi,pyj]));

mz:=mzi+(mzj-mzi)*x/L:

fmz:=vector(4,i->int(Nrotz[i]*mz*J,xi=-1..1)):
fmz:=simplify(jacobian(fmz,[mzi,mzj]));


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

ber_ex1

restart:        # initialisation

with(linalg):   # calcul matriciel

# calcul des déplacements inconnus

KLL:=matrix([[8*L^2,-6*L,2*L^2],[-6*L,12,-6*L],[2*L^2,-6*L,4*L^2]]):
c:=E*Iz/L^3:
KLL:=scalarmul(KLL,c);
FL:=vector([0,-P,0]);
UL:=linsolve(KLL,FL);

# efforts et déplacements élémentaires

element:=1;

# vecteur déplacement

if element=1
  then unod:=vector([0,0,0,UL[1]]);  # élément 1-2
  else unod:=vector([0,UL[1],UL[2],UL[3]]);  # élément 2-3
fi:

# matrice de rigidité

k:=matrix([[12,6*L,-12,6*L],[6*L,4*L^2,-6*L,2*L^2],
           [-12,-6*L,12,-6*L],[6*L,2*L^2,-6*L,4*L^2]]):
k:=scalarmul(k,E*Iz/L^3);

# vecteur force nodal

fnod:=multiply(k,unod);

# effort tranchant

Ty:=x->-fnod[1];
Ty(x);Ty(L);

# moment fléchissant

Mfz:=x->-fnod[2]-int(Ty(s),s=0..x);
Mfz(x);Mfz(L);

# rotation des sections droites : pente

rotz:=x->unod[2]+int(Mfz(s)/E/Iz,s=0..x);simplify(rotz(x));rotz(L);

# déplacement suivant y : flèche

v:=x->unod[1]+int(rotz(s),s=0..x);simplify(v(x));v(L);

# représentations graphiques

assign(E=1,Iz=1,L=1,P=1);

plot(Ty(x),x=0..L,title="Effort tranchant",
          thickness=2,color=blue,labels=["x","Ty"]);
plot(Mfz(x),x=0..L,title="Moment fléchissant",
          thickness=2,color=blue,labels=["x","Mfz"]);
plot(rotz(x),x=0..L,title="Rotation des sections droites : pente",
          thickness=2,color=blue,labels=["x","rotz"]);
plot(v(x),x=0..L,title="Déplacement suivant y : flèche",
          thickness=2,color=blue,labels=["x","v"]);


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

ber_ex2

restart:        # initialisation

with(linalg):   # calcul matriciel

# calcul des déplacements inconnus

c:=E*Iz/L:
KLL:=matrix([[8*c,2*c],[2*c,4*c]]);
c:=p*L^2/12:
FL:=vector([c,c]);
UL:=linsolve(KLL,FL);

# efforts et déplacements élémentaires

element:=1;

if element=1 then py:=-2*p; else py:=-p;fi:

# vecteur déplacement

if element=1
  then unod:=vector([0,0,0,UL[1]]);  # élément 1-2
  else unod:=vector([0,UL[1],0,UL[2]]);  # élément 2-3
fi:

# matrice de rigidité

k:=matrix([[12,6*L,-12,6*L],[6*L,4*L^2,-6*L,2*L^2],
           [-12,-6*L,12,-6*L],[6*L,2*L^2,-6*L,4*L^2]]):
k:=scalarmul(k,E*Iz/L^3);

# vecteur force dû à la force répartie

c:=py*L/12:f:=vector([6*c,L*c,6*c,-L*c]);

# vecteur force nodal

fnod:=matadd(multiply(k,unod),f,1,-1);

# effort tranchant

Ty:=x->-fnod[1]-int(py,s=0..x);simplify(Ty(x));Ty(L);

# moment fléchissant

Mfz:=x->-fnod[2]-int(Ty(s),s=0..x);simplify(Mfz(x));Mfz(L);
xmax:=-fnod[1]/py;
if (xmax/L>0) and (xmax/L<1) then Mfmax:=Mfz(xmax):fi;
 
# rotation des sections droites : pente

rotz:=x->unod[2]+int(Mfz(s)/E/Iz,s=0..x);simplify(rotz(x));rotz(L);

# déplacement suivant y : flèche

v:=x->unod[1]+int(rotz(s),s=0..x);simplify(v(x));v(L);

# représentations graphiques

assign(E=1,Iz=1,L=1,p=1);

plot(Ty(x),x=0..L,title="Effort tranchant");
plot(Mfz(x),x=0..L,title="Moment fléchissant");

plot([Ty(x),Mfz(x)],x=0..L,title="Effort tranchant et moment fléchissant",
          color=[red,blue],legend=["Ty","Mfz"]);

plot(rotz(x),x=0..L,title="Rotation des sections droites : pente");
plot(v(x),x=0..L,title="Déplacement suivant y : flèche");