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

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

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);

# 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(%):
Tzj:=Tz(L):Mfyj:=Mfy(L):

# matrice de rigidité

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

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

# vecteur force

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

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

tim_int

# fonctions d'interpolation en fonction de x

restart:with(linalg):

# 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(x=0,w)=wi;
eq2:=subs(x=L,w)=wj;
eq3:=subs(x=0,roty)=rotyi;
eq4:=subs(x=L,roty)=rotyj;
eq5:=diff(roty,x$2)=12*c/fz/L^2;

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

# 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(matadd(dNw,Nroty));

# factorisation

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

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

tim_int_rem

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

restart:with(linalg):

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(%):

# 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);
dNw:=map(diff,Nw,x);
dNroty:=map(diff,Nroty,x);
Bc:=simplify(matadd(dNw,Nroty));

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

tim_int_par

# fonctions d'interpolation sous forme paramétrique

restart:with(linalg):

# 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(xi=-1,w)=wi;
eq2:=subs(xi=1,w)=wj;
eq3:=subs(xi=-1,roty)=rotyi;
eq4:=subs(xi=1,roty)=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(%):

# 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(matadd(dNw,Nroty));

# factorisation

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

simplify(scalarmul(Nw,(1+fz)*8));
simplify(scalarmul(dNw,(1+fz)*L*4));
simplify(scalarmul(Nroty,(1+fz)*L*4));
simplify(scalarmul(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

restart:with(linalg):

# 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),shape=symmetric):
GAkz:=12*EIy/L^2/fz:
kc:=Matrix(4,4,(i,j)->int(Bc[i]*Bc[j]*GAkz*J,xi=-1..1),shape=symmetric):
k:=matadd(kf,kc):k:=simplify(k);
simplify(scalarmul(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),shape=symmetric):
mv:=simplify(mv);
simplify(scalarmul(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),shape=symmetric):
mroty:=simplify(mroty);
simplify(scalarmul(mroty,30*(1+fz)^2*L/rho/Iy));

# vecteur force

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

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

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

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