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

  programmes MuPAD  -  24/03/06

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

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

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(fy>0):assume(EIz>0):

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

// formules de Bresse

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

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

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

// matrice de rigidité

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

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

// vecteur force

fpy:=-jacobian(fnod,[pyi,pyj]):
simplify(fpy*(1+fy)*120/L);
fmz:=-jacobian(fnod,[mzi,mzj]):
simplify(fmz*(1+fy)*12);

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

tim_int

// fonctions d'interpolation en fonction de x

reset():export(linalg):

assume(L>0):assume(fy>0):

// champ de déplacements

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

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

eq1:=subs(v,x=0)=vi;
eq2:=subs(v,x=L)=vj;
eq3:=subs(rotz,x=0)=rotzi;
eq4:=subs(rotz,x=L)=rotzj;
eq5:=diff(rotz,x$2)=12*c/fy/L^2;

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

// fonctions d'interpolation en fonction de x

unod:=[vi,rotzi,vj,rotzj]:
Nv:=simplify(grad(v,unod));
Nrotz:=simplify(grad(rotz,unod));
dv:=diff(v,x):
dNv:=simplify(grad(dv,unod));
drotz:=diff(rotz,x):
dNrotz:=simplify(grad(drotz,unod));
Bc:=simplify(dNv-Nrotz);

// factorisation

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

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

tim_int_par

// fonctions d'interpolation sous forme paramétrique

reset():export(linalg):

assume(L>0):assume(fy>0):

// jacobien de la transformation géométrique

J:=L/2;

// champ de déplacements

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

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

eq1:=subs(v,xi=-1)=vi;
eq2:=subs(v,xi=1)=vj;
eq3:=subs(rotz,xi=-1)=rotzi;
eq4:=subs(rotz,xi=1)=rotzj;
eq5:=diff(rotz,xi$2)/J^2=12*c/fy/L^2;

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

// fonctions d'interpolation

unod:=[vi,rotzi,vj,rotzj]:
Nv:=simplify(grad(v,unod)):
Nrotz:=simplify(grad(rotz,unod)):
dv:=diff(v,xi)/J:
dNv:=simplify(grad(dv,unod)):
drotz:=diff(rotz,xi)/J:
dNrotz:=simplify(grad(drotz,unod)):

Bc:=simplify(dNv-Nrotz);

// factorisation

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

simplify(Nv*(1+fy)*8);
simplify(dNv*(1+fy)*L*4);
simplify(Nrotz*(1+fy)*L*4);
simplify(dNrotz*(1+fy)*L^2);

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

tim_int_rem

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

reset():export(linalg):

assume(L>0):assume(fy>0):

//  formules de Bresse

rotz:=rotzi+int(Mfzi-Tyi*s,s=0..x):      
v:=vi+rotzi*x+int((Mfzi-Tyi*s)*(x-s),s=0..x)+Tyi*x*fy*L^2/12:

solve({rotzj=subs(rotz,x=L),vj=subs(v,x=L)},{Tyi,Mfzi}):assign(op(%)):

// fonctions d'interpolation en fonction de x

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

unod:=[vi,rotzi,vj,rotzj]:
Nv:=simplify(grad(v,unod));
Nrotz:=simplify(grad(rotz,unod));
dv:=diff(v,x):
dNv:=simplify(grad(dv,unod));
drotz:=diff(rotz,x):
dNrotz:=simplify(grad(drotz,unod));
Bc:=simplify(dNv-Nrotz);

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

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/8/(1+fy):
Nv:=[2*(xi*xi+xi-2-2*fy)*(-1+xi)*c,(-1+xi*xi)*(xi-1-fy)*L*c,
     -2*(-2-xi+xi*xi-2*fy)*(1+xi)*c,(-1+xi*xi)*(xi+1+fy)*L*c]:
c:=1/4/L/(1+fy):
dNv:=[(6*xi*xi-6-4*fy)*c,-L*(1-3*xi*xi+2*xi+2*fy*xi)*c,
      (-6*xi*xi+6+4*fy)*c,L*(-1+3*xi*xi+2*xi+2*fy*xi)*c]:
Nrotz:=[(-6+6*xi*xi)*c,L*(-1+xi)*(3*xi+1-2*fy)*c,
        (6-6*xi*xi)*c,(1+xi)*L*(-1+2*fy+3*xi)*c]:
c:=1/L^2/(1+fy):
dNrotz:=[6*xi*c,-L*(-3*xi+1+fy)*c,-6*xi*c,L*(1+fy+3*xi)*c]:
c:=fy/2/L/(1+fy):
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(fy>0):

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

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

// fonctions d'interpolation

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

// matrice de rigidité

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

// matrice de masse

mv:=matrix(4,4,(i,j)->int(Nv[i]*Nv[j]*rho*A*J,xi=-1..1)):
mv:=simplify(mv);
simplify(mv*420*(1+fy)^2/rho/A/L);

mrotz:=matrix(4,4,(i,j)->int(Nrotz[i]*Nrotz[j]*rho*Iz*J,xi=-1..1)):
mrotz:=simplify(mrotz);
simplify(mrotz*30*(1+fy)^2*L/rho/Iz);

// force linéairement répartie sur toute la longueur de l'élément

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

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

// couple linéairement réparti sur toute la longueur de l'élément

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

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

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

tim_partition_tan

// partition du champ de déplacements 
// en mouvement de corps rigide et mouvement de déformation pure
// méthode de la tangente
// calcul des matrices élémentaires à l'aide des fonctions d'interpolation
// section droite constante

reset():export(linalg):

assume(L>0):assume(fy>0):

// matrice cinématique

h:=matrix([[-1,-L],[0,-1]]):
ht:=transpose(h):

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

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

// fonctions d'interpolation

NRv:=[1,(1+xi)*L/2]:
NRrotz:=[0,1]:

c:=1/8/(1+fy):
NDv:=[-2*(-2-xi+xi*xi-2*fy)*(1+xi)*c,(-1+xi*xi)*(xi+1+fy)*L*c]:
c:=1/4/L/(1+fy):
dNDv:=[(-6*xi*xi+6+4*fy)*c,L*(-1+3*xi*xi+2*xi+2*fy*xi)*c]:
NDrotz:=[(6-6*xi*xi)*c,(1+xi)*L*(-1+2*fy+3*xi)*c]:
c:=1/L^2/(1+fy):
dNDrotz:=[-6*xi*c,L*(1+fy+3*xi)*c]:
c:=fy/2/L/(1+fy):
Bc:=[2*c,-L*c]:

// matrice de rigidité

Bf:=dNDrotz:
kf:=matrix(2,2,(i,j)->int(Bf[i]*Bf[j]*EIz*J,xi=-1..1)):
GAky:=12*EIz/L^2/fy:
kc:=matrix(2,2,(i,j)->int(Bc[i]*Bc[j]*GAky*J,xi=-1..1)):
kD:=kf+kc:kD:=simplify(kD);

k12:=ht*kD:
k11:=k12*h:
k:=stackMatrix(k11.k12,transpose(k12).kD):simplify(k);

// matrice de masse

mR:=matrix(2,2,(i,j)->int(NRv[i]*NRv[j]*rho*A*J+NRrotz[i]*NRrotz[j]*rho*Iz*J,xi=-1..1)):
mR:=simplify(mR);

mD:=matrix(2,2,(i,j)->int(NDv[i]*NDv[j]*rho*A*J+NDrotz[i]*NDrotz[j]*rho*Iz*J,xi=-1..1)):
mD:=simplify(mD);

mRD:=matrix(2,2,(i,j)->int(NRv[i]*NDv[j]*rho*A*J+NRrotz[i]*NDrotz[j]*rho*Iz*J,xi=-1..1)):
mRD:=simplify(mRD);
mDR:=transpose(mRD):

m11:=mR+mRD*h+ht*mDR+ht*mD*h:
m12:=mRD+mD*h:
m:=stackMatrix(m11.m12,transpose(m12).mD):m:=simplify(m);

// vecteur force

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

fR:=matrix(2,1,i->int(NRv[i]*py*J+NRrotz[i]*mz*J,xi=-1..1)):
fD:=matrix(2,1,i->int(NDv[i]*py*J+NDrotz[i]*mz*J,xi=-1..1)):
f1:=fR+ht*fD:
f:=matrix([f1[1],f1[2],fD[1],fD[2]]):f:=simplify(f);

jacobian(f,[pyi,pyj,mzi,mzj]);

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

tim_partition_sec

// partition du champ de déplacements 
// en mouvement de corps rigide et mouvement de déformation pure
// méthode de la sécante
// calcul des matrices élémentaires à l'aide des fonctions d'interpolation
// section droite constante

reset():export(linalg):

assume(L>0):assume(fy>0):assume(EIz>0):

// matrice cinématique

aR:=matrix([[1,0,0,0],[0,0,1,0]]):
aD:=matrix([[1/L,1,-1/L,0],[1/L,0,-1/L,1]]):

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

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

// fonctions d'interpolation

NRv:=[(1-xi)/2,(1+xi)/2]:
NRrotz:=[-1/L,1/L]:

c:=L*(xi*xi-1)/8/(1+fy):
NDv:=[(xi-1-fy)*c,(xi+1+fy)*c]:
c:=1/4/(1+fy):
dNDv:=[(3*xi*xi-2*xi-2*fy*xi-1)*c,(3*xi*xi+2*xi+2*fy*xi-1)*c]:
NDrotz:=[(xi-1)*(3*xi+1-2*fy)*c,(xi+1)*(3*xi-1+2*fy)*c]:
c:=1/L/(1+fy):
dNDrotz:=[(3*xi-1-fy)*c,(3*xi+1+fy)*c]:
c:=-fy/2/(1+fy):
Bc:=[c,c]:

// matrice de rigidité

Bf:=dNDrotz:
kf:=matrix(2,2,(i,j)->int(Bf[i]*Bf[j]*EIz*J,xi=-1..1)):
GAky:=12*EIz/L^2/fy:
kc:=matrix(2,2,(i,j)->int(Bc[i]*Bc[j]*GAky*J,xi=-1..1)):
kD:=kf+kc:kD:=simplify(kD);

k:=transpose(aD)*kD*aD:simplify(k);

// matrice de masse

mR:=matrix(2,2,(i,j)->int(NRv[i]*NRv[j]*rho*A*J+NRrotz[i]*NRrotz[j]*rho*Iz*J,xi=-1..1)):
mR:=simplify(mR);

mD:=matrix(2,2,(i,j)->int(NDv[i]*NDv[j]*rho*A*J+NDrotz[i]*NDrotz[j]*rho*Iz*J,xi=-1..1)):
mD:=simplify(mD);

mRD:=matrix(2,2,(i,j)->int(NRv[i]*NDv[j]*rho*A*J+NRrotz[i]*NDrotz[j]*rho*Iz*J,xi=-1..1)):
mRD:=simplify(mRD);
mDR:=transpose(mRD):

m1:=transpose(aR)*mR*aR:
m2:=transpose(aD)*mD*aD:
m3:=transpose(aR)*mRD*aD:
m4:=transpose(m3):
m:=m1+m2+m3+m4:m:=simplify(m);

// vecteur force

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

fR:=matrix(2,1,i->int(NRv[i]*py*J+NRrotz[i]*mz*J,xi=-1..1)):
fD:=matrix(2,1,i->int(NDv[i]*py*J+NDrotz[i]*mz*J,xi=-1..1)):
f:=transpose(aR)*fR+transpose(aD)*fD:f:=simplify(f);

// calcul de kD à l'aide du théorème de Castigliano

cT:=int(1/GAky/L*J,xi=-1..1);
c11:=int((1-x/L)^2/EIz*J,xi=-1..1)+cT;
c22:=int((x/L)^2/EIz*J,xi=-1..1)+cT;
c12:=int((1-x/L)*x/L/EIz*J,xi=-1..1)-cT;
c:=matrix([[c11,c12],[c12,c22]]):
kD:=simplify(c^-1);