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

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

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

# 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(%):
Tyj:=Ty(L):Mfzj:=Mfz(L):

# matrice de rigidité

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

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

# vecteur force

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

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

tim_int

# fonctions d'interpolation en fonction de x

restart:with(linalg):

# 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(x=0,v)=vi;
eq2:=subs(x=L,v)=vj;
eq3:=subs(x=0,rotz)=rotzi;
eq4:=subs(x=L,rotz)=rotzj;
eq5:=diff(rotz,x$2)=12*c/fy/L^2;

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

# 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(matadd(dNv,Nrotz,1,-1));

# factorisation

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

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

tim_int_par

# fonctions d'interpolation sous forme paramétrique

restart:with(linalg):

# jacobien de la transformation

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(xi=-1,v)=vi;
eq2:=subs(xi=1,v)=vj;
eq3:=subs(xi=-1,rotz)=rotzi;
eq4:=subs(xi=1,rotz)=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(%):

# 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(matadd(dNv,Nrotz,1,-1));

# factorisation

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

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

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

tim_int_rem

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

restart:with(linalg):

#  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(x=L,rotz),vj=subs(x=L,v)},{Tyi,Mfzi}):assign(%):

# 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(matadd(dNv,Nrotz,1,-1));

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

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

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/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),shape=symmetric):
GAky:=12*EIz/L^2/fy:
kc:=Matrix(4,4,(i,j)->int(Bc[i]*Bc[j]*GAky*J,xi=-1..1),shape=symmetric):
k:=matadd(kf,kc):k:=simplify(k);
simplify(scalarmul(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),shape=symmetric):
mv:=simplify(mv);
simplify(scalarmul(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),shape=symmetric):
mrotz:=simplify(mrotz);
simplify(scalarmul(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:=vector(4,i->int(Nv[i]*py*J,xi=-1..1)):
fpy:=jacobian(fpy,[pyi,pyj]):fpy:=simplify(fpy);
simplify(scalarmul(fpy,(1+fy)*120));

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

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

fmz:=vector(4,i->int(Nrotz[i]*mz*J,xi=-1..1)):
fmz:=jacobian(fmz,[mzi,mzj]):fmz:=simplify(fmz);
simplify(scalarmul(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

restart:with(linalg):

# 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),shape=symmetric):
GAky:=12*EIz/L^2/fy:
kc:=Matrix(2,2,(i,j)->int(Bc[i]*Bc[j]*GAky*J,xi=-1..1),shape=symmetric):
kD:=matadd(kf,kc):kD:=simplify(kD);

k12:=multiply(ht,kD):
k11:=multiply(k12,h):
k:=blockmatrix(2,2,[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),shape=symmetric):
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),shape=symmetric):
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),shape=symmetric):
mRD:=simplify(mRD);
mDR:=transpose(mRD):

m11:=evalm(mR+mRD&*h+ht&*mDR+ht&*mD&*h):
m12:=evalm(mRD+mD&*h):
m:=blockmatrix(2,2,[m11,m12,transpose(m12),mD]):m:=simplify(m);

# vecteur force

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

fR:=vector(2,i->int(NRv[i]*py*J+NRrotz[i]*mz*J,xi=-1..1)):
fD:=vector(2,i->int(NDv[i]*py*J+NDrotz[i]*mz*J,xi=-1..1)):
f1:=evalm(fR+ht&*fD):
f:=vector([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

restart:with(linalg):

# 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),shape=symmetric):
GAky:=12*EIz/L^2/fy:
kc:=Matrix(2,2,(i,j)->int(Bc[i]*Bc[j]*GAky*J,xi=-1..1),shape=symmetric):
kD:=matadd(kf,kc):kD:=simplify(kD);

k:=multiply(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),shape=symmetric):
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),shape=symmetric):
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),shape=symmetric):
mRD:=simplify(mRD);
mDR:=transpose(mRD):

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

# vecteur force

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

fR:=vector(2,i->int(NRv[i]*py*J+NRrotz[i]*mz*J,xi=-1..1)):
fD:=vector(2,i->int(NDv[i]*py*J+NDrotz[i]*mz*J,xi=-1..1)):
f:=evalm(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:=inverse(c);

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

tim_ex1

# influence de l'effort tranchant

restart:with(linalg):

# calcul des déplacements inconnus

c:=E*Iz/(L^3*(1+fy));
KLL:=matrix([[12*c,-6*L*c],[-6*L*c,L^2*(4+fy)*c]]);
FL:=vector([p*L/2,-p*L^2/12]);
UL:=linsolve(KLL,FL);

f:=1-subs(fy=0,UL[1])/UL[1];

# la section droite est rectangulaire

A:=b*h;Iz:=b*h^3/12;
ky:=5/6;
#ky:=10*(1+nu)/(12+11*nu);

# influence de l'effort tranchant

L:=x*h;
fy:=24*(1+nu)*Iz/L^2/ky/A;
f:=simplify(f);

f1:=subs(nu=0,f*100):
f2:=subs(nu=0.2,f*100):
f3:=subs(nu=0.4,f*100):
f4:=subs(nu=0.5,f*100):

plot([f1,f2,f3,f4],x=3..10,
title="Influence en % de l'effort tranchant",
labels=["L/h",""],
color=[green,red,blue,gray],thickness=2,

legend=["nu=0","nu=0.2","nu=0.4","nu=0.5"]);

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

tim_ex2

restart:

Ty:=F;
Mfz:=M+F*(L-x);

dEdef:=Ty^2/(2*G*A*ky)+Mfz^2/(2*E*Iz);

v2:=diff(dEdef,F);
rotz2:=diff(dEdef,M);
v2:=subs(M=0,v2);
rotz2:=subs(M=0,rotz2);

A:=c^2*(2-x/L)^2;Iz:=c^4*(2-x/L)^4/12;G:=E/(1+nu)/2;
ky:=5/6;
#ky:=10*(1+nu)/(12+11*nu);

v2:=int(v2,x=0..L,continuous);
rotz2:=int(rotz2,x=0..L,continuous);

restart:

Ty:=F+p*(L-x);
Mfz:=M+F*(L-x)+p*(L-x)^2/2;

dEdef:=Ty^2/(2*G*A*ky)+Mfz^2/(2*E*Iz);

v2:=diff(dEdef,F);
rotz2:=diff(dEdef,M);
v2:=subs(F=0,M=0,v2);
rotz2:=subs(F=0,M=0,rotz2);

A:=c^2*(2-x/L)^2;Iz:=c^4*(2-x/L)^4/12;G:=E/(1+nu)/2;
#ky:=5/6;
ky:=10*(1+nu)/(12+11*nu);

v2:=int(v2,x=0..L,continuous);
rotz2:=int(rotz2,x=0..L,continuous);