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

  programmes MuPAD  -  24/03/06

  Méthode des éléments finis :
  
  Poutre de Timoshenko : 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

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

tim_mat

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

reset();export(linalg):

assume(L>0):assume(fz>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);

// formules de Nawier-Bresse

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

// 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:=Mfz(L):

// matrice de rigidité

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

k:=jacobian(fnod,unod):
simplify(k*L^3*(1+fz)/EIy);

// vecteur force

fpz:=-jacobian(fnod,[pzi,pzj]):
simplify(fpz*(1+fz)*120/L);
fmy:=-jacobian(fnod,[myi,myj]):
simplify(fmy*(1+fz)*12);

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

tim_int

// fonctions d'interpolation en fonction de x

reset():export(linalg):

assume(L>0):assume(fz>0):

// champ de déplacements

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

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

eq1:=subs(w,x=0)=wi;
eq2:=subs(w,x=L)=wj;
eq3:=subs(roty,x=0)=rotyi;
eq4:=subs(roty,x=L)=rotyj;
eq5:=diff(roty,x$2)=12*c/fz/L^2;

solve({eq1,eq2,eq3,eq4,eq5},{a0,a1,a2,a3,c}):assign(op(%));

// fonctions d'interpolation en fonction de x

unod:=[wi,rotyi,wj,rotyj]:
Nw:=simplify(grad(w,unod));
Nroty:=simplify(grad(roty,unod));
dw:=diff(w,x):
dNw:=simplify(grad(dw,unod));
droty:=diff(roty,x):
dNroty:=simplify(grad(droty,unod));
Bc:=simplify(dNw+Nroty);

// factorisation

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

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

tim_int_rem

// calcul des fonctions d'interpolation et de leurs dérivées

reset():export(linalg):

assume(L>0):assume(fz>0):

Mfy:=x->Mfyi+Tzi*x:

//  formules de Navier-Bresse avec EIy=1

roty:=x->rotyi+int(Mfy(s),s=0..x):      
w:=x->wi-rotyi*x-int(Mfy(s)*(x-s),s=0..x)+Tzi*x*fz*L^2/12:

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

// fonctions d'interpolation en fonction de x

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

Nw:=grad(w(x),unod):Nw:=simplify(Nw);
Nroty:=grad(roty(x),unod):Nroty:=simplify(Nroty);
dw:=diff(w(x),x):
dNw:=simplify(grad(dw,unod));
droty:=diff(roty(x),x):
dNroty:=simplify(grad(droty,unod));
Bc:=simplify(dNw+Nroty);


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

tim_int_par

// fonctions d'interpolation sous forme paramétrique

reset():export(linalg):

assume(L>0):assume(fz>0):

// jacobien de la transformation géométrique

J:=L/2;

// champ de déplacements

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

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

eq1:=subs(w,xi=-1)=wi;
eq2:=subs(w,xi=1)=wj;
eq3:=subs(roty,xi=-1)=rotyi;
eq4:=subs(roty,xi=1)=rotyj;
eq5:=diff(roty,xi$2)/J^2=12*c/fz/L^2;

solve({eq1,eq2,eq3,eq4,eq5},{a0,a1,a2,a3,c}):assign(op(%)):

// fonctions d'interpolation

unod:=[wi,rotyi,wj,rotyj]:
Nw:=simplify(grad(w,unod)):
Nroty:=simplify(grad(roty,unod)):
dw:=diff(w,xi)/J:
dNw:=simplify(grad(dw,unod)):
droty:=diff(roty,xi)/J:
dNroty:=simplify(grad(droty,unod)):

Bc:=simplify(dNw+Nroty);

// factorisation

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

simplify(Nw*(1+fz)*8);
simplify(dNw*(1+fz)*L*4);
simplify(Nroty*(1+fz)*L*4);
simplify(dNroty*(1+fz)*L^2);

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

tim_interpolation

// fonctions d'interpolation et dérivées sous forme paramétrique

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

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

// fonctions d'interpolation

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

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

tim_int_mat

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

reset():export(linalg):

assume(L>0):assume(fz>0):

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

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

// fonctions d'interpolation

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

// matrice de rigidité

Bf:=dNroty:
kf:=matrix(4,4,(i,j)->int(Bf[i]*Bf[j]*EIy*J,xi=-1..1)):
GAkz:=12*EIy/L^2/fz:
kc:=matrix(4,4,(i,j)->int(Bc[i]*Bc[j]*GAkz*J,xi=-1..1)):
k:=kf+kc:k:=simplify(k);
simplify(k*L^3*(1+fz)/EIy);

// matrice de masse

mv:=matrix(4,4,(i,j)->int(Nw[i]*Nw[j]*rho*A*J,xi=-1..1)):
mv:=simplify(mv);
simplify(mv*420*(1+fz)^2/rho/A/L);
mroty:=matrix(4,4,(i,j)->int(Nroty[i]*Nroty[j]*rho*Iy*J,xi=-1..1)):
mroty:=simplify(mroty);
simplify(mroty*30*(1+fz)^2*L/rho/Iy);

// vecteur force

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

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

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

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