==================================================================== programmes MuPAD - 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 reset();export(linalg): assume(L>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); // 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(op(%)): 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(op(%)): 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(op(%)): 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(op(%)): k:=jacobian([-Tyi,Tyj],[vi,vj]); f:=jacobian([Tyi,-Tyj],[pyi,pyj,mzi,mzj]); rotzi;rotzj; --------------------------------------------------------- ber_int // calcul des fonctions d'interpolation reset():export(linalg): assume(L>0): // 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(v,x=0):eq2:=vj=subs(v,x=L): eq3:=rotzi=subs(rotz,x=0):eq4:=rotzj=subs(rotz,x=L): solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(op(%)): // 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 plotfunc2d(subs(Nv[1],L=1),subs(Nv[2],L=1),subs(Nv[3],L=1),subs(Nv[4],L=1),x=0..1); --------------------------------------------------------- ber_int_par // calcul des fonctions d'interpolation sous forme paramétrique reset():export(linalg): assume(L>0): // 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(v,xi=-1):eq2:=vj=subs(v,xi=1): eq3:=rotzi=subs(rotz,xi=-1):eq4:=rotzj=subs(rotz,xi=1): solve({eq1,eq2,eq3,eq4},{a0,a1,a2,a3}):assign(op(%)): // 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(Nv,factor)); Nrotz:=simplify(map(Nrotz,factor)); dNrotz:=simplify(map(dNrotz,factor)); // représentation graphique pour L=1 plotfunc2d(Nv[1],subs(Nv[2],L=1),Nv[3],subs(Nv[4],L=1),xi=-1..1); ---------------------------------------------------------------------- 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 reset():export(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)); // matrice de masse mv:=matrix(4,4,(i,j)->int(Nv[i]*Nv[j]*rho*A*J,xi=-1..1)); mrotz:=matrix(4,4,(i,j)->int(Nrotz[i]*Nrotz[j]*rho*Iz*J,xi=-1..1)); // vecteur force py:=pyi+(pyj-pyi)*x/L: fpy:=matrix(4,1,i->int(Nv[i]*py*J,xi=-1..1)): fpy:=simplify(jacobian(fpy,[pyi,pyj])); mz:=mzi+(mzj-mzi)*x/L: fmz:=matrix(4,1,i->int(Nrotz[i]*mz*J,xi=-1..1)): fmz:=simplify(jacobian(fmz,[mzi,mzj])); --------------------------------------------------------------------- ber_ex1 reset(): // initialisation export(linalg): // calcul matriciel // calcul des déplacements inconnus KLL:=matrix(3,3,[[8*L^2,-6*L,2*L^2],[-6*L,12,-6*L],[2*L^2,-6*L,4*L^2]]): c:=EIz/L^3: KLL:=KLL*c; FL:=matrix(3,1,[0,-P,0]); UL:=matlinsolve(KLL,FL); // efforts et déplacements élémentaires element:=1; // vecteur déplacement if element=1 then unod:=matrix(4,1,[0,0,0,UL[1]]); // élément 1-2 else unod:=matrix(4,1,[0,UL[1],UL[2],UL[3]]); // élément 2-3 end_if: // matrice de rigidité k:=matrix(4,4,[[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:=k*EIz/L^3; // vecteur force nodal fnod:=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)/EIz,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([EIz=1,L=1,P=1]); plotfunc2d(Ty(x),x=0..L); plotfunc2d(Mfz(x),x=0..L); plotfunc2d(rotz(x),x=0..L); plotfunc2d(v(x),x=0..L); -------------------------------------- ber_ex2 reset(): // initialisation export(linalg): // calcul matriciel // calcul des déplacements inconnus c:=EIz/L: KLL:=matrix(2,2,[[8*c,2*c],[2*c,4*c]]); c:=p*L^2/12: FL:=matrix(2,1,[c,c]); UL:=matlinsolve(KLL,FL); // efforts et déplacements élémentaires element:=1; if element=1 then py:=-2*p; else py:=-p;end_if: // vecteur déplacement if element=1 then unod:=matrix(4,1,[0,0,0,UL[1]]); // élément 1-2 else unod:=matrix(4,1,[0,UL[1],0,UL[2]]); // élément 2-3 end_if: // matrice de rigidité k:=matrix(4,4,[[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:=k*EIz/L^3; // vecteur force dû à la force répartie c:=py*L/12:f:=matrix(4,1,[6*c,L*c,6*c,-L*c]); // vecteur force nodal fnod:=k*unod-f; // 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):end_if; // rotation des sections droites : pente rotz:=x->unod[2]+int(Mfz(s)/EIz,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([EIz=1,L=1,p=1]); plotfunc2d(Ty(x),x=0..L); plotfunc2d(Mfz(x),x=0..L); plotfunc2d(rotz(x),x=0..L); plotfunc2d(v(x),x=0..L);