LECCION 6 - EJERCICIO 2 (17.2) v.2005

INICIO

Off[General::"spell1"]
Off[General::"spell"]

DEFINICION  ELEMENTO CUADRILATERO

Cua4=;

Show[Cua4, ImageSize->600]

NO SE DEFINEN DATOS NUMERICOS - SE ELIMINAN LOS PREVIAMENTE DEFINIDOS

ClearAll[Em,ν,a,b,h];

RELACION ENTRE LONGITUD Y ALTURA DEL ELEMENTO

b=γ*a;

COORDENADAS NODOS

ncoor={{0,0},{a,0},{a,b},{0,b}};

MATRIZ DE PROPIEDADES DEL MATERIAL

Emat=Em/(1-ν^2)*{{1,ν,0},{X,1,0},{0,0,(1-ν)/2}};

FORMULACION ELEMENTO CUADRILATERO 4 NODOS - PROBLEMA TENSION PLANA

Quad4IsoPMembraneStiffness[ncoor_,mprop_,fprop_,options_]:=Module[{i,k,p=2,numer=False,Emat,th=1,h,qcoor,c,w,Nf,dNx,dNy,Jdet,B,Ke=Table[0,{8},{8}]},Emat=mprop[[1]];
If[Length[options]==2,{numer,p}=options,{numer}=options];
If[Length[fprop]>0,th=fprop[[1]]];
If[p<1||p>4,Print["p out of range"];Return[Null]];
For[k=1,k≤p*p,k++,{qcoor,w}=QuadGaussRuleInfo[{p,numer},k];
{Nf,dNx,dNy,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
If[Length[th]==0,h=th,h=th.Nf];c=w*Jdet*h;
B={Flatten[Table[{dNx[[i]],0},{i,4}]],Flatten[Table[{0,dNy[[i]]},{i,4}]],Flatten[Table[{dNy[[i]],dNx[[i]]},{i,4}]]};
Ke+=Simplify[c*Transpose[B].(Emat.B)];];Return[Simplify[Ke]]];

Quad4IsoPShapeFunDer[ncoor_,qcoor_]:=Module[{Nf,dNx,dNy,dNξ,dNη,i,J11,J12,J21,J22,Jdet,ξ,η,x,y},{ξ,η}=qcoor;
Nf={(1-ξ)*(1-η),(1+ξ)*(1-η),(1+ξ)*(1+η),(1-ξ)*(1+η)}/4;
dNξ={-(1-η),(1-η),(1+η),-(1+η)}/4;
dNη={-(1-ξ),-(1+ξ),(1+ξ),(1-ξ)}/4;
x=Table[ncoor[[i,1]],{i,4}];y=Table[ncoor[[i,2]],{i,4}];
J11=dNξ.x;J21=dNξ.y;J12=dNη.x;J22=dNη.y;
Jdet=Simplify[J11*J22-J12*J21];
dNx=(J22*dNξ-J21*dNη)/Jdet;dNx=Simplify[dNx];
dNy=(-J12*dNξ+J11*dNη)/Jdet;dNy=Simplify[dNy];
Return[{Nf,dNx,dNy,Jdet}]];

QuadGaussRuleInfo[{rule_,numer_},point_]:=Module[{xi,eta,p1,p2,i1,i2,w1,w2,k,info=Null},If[Length[rule]==2,{p1,p2}=rule,p1=p2=rule];
If[Length[point]==2,{i1,i2}=point,k=point;i2=Floor[(k-1)/p1]+1;i1=k-p1*(i2-1)];
{xi,w1}=LineGaussRuleInfo[{p1,numer},i1];
{eta,w2}=LineGaussRuleInfo[{p2,numer},i2];
info={{xi,eta},w1*w2};
If[numer,Return[N[info]],Return[Simplify[info]]];];

LineGaussRuleInfo[{rule_,numer_},point_]:=Module[{g2={-1,1}/Sqrt[3],w3={5/9,8/9,5/9},g3={-Sqrt[3/5],0,Sqrt[3/5]},w4={(1/2)-Sqrt[5/6]/6,(1/2)+Sqrt[5/6]/6,(1/2)+Sqrt[5/6]/6,(1/2)-Sqrt[5/6]/6},g4={-Sqrt[(3+2*Sqrt[6/5])/7],-Sqrt[(3-2*Sqrt[6/5])/7],Sqrt[(3-2*Sqrt[6/5])/7],Sqrt[(3+2*Sqrt[6/5])/7]},i,info=Null},i=point;
If[rule==1,info={0,2}];
If[rule==2,info={g2[[i]],1}];
If[rule==3,info={g3[[i]],w3[[i]]}];
If[rule==4,info={g4[[i]],w4[[i]]}];
If[numer,Return[N[info]],Return[Simplify[info]]];];

REGLA DE GAUSS p = 2 (2 x 2) -  PROPORCIONA SUFICIENCIA DE RANGO

p=2;

Ke=Quad4IsoPMembraneStiffness[ncoor,{Emat,0,0},{h},{X,p}];

Ke //MatrixForm

ESCALAMOS LA MATRIZ DIVIENDO POR EL FACTOR DE ESCALA

scaledKe=Simplify[Ke*(24*γ*(1-(ν^2))/(Em*h))];

RESULTADO A COMPARAR - (E17.4)

Print["Ke=",Em*h/(24*γ*(1-(ν^2))),"*\n",scaledKe //MatrixForm];