ELEMENTO CUADRILATERO REGULAR - 9 NODOS

Implementación en Mathematica

v.2018

1. DATOS INICIALES

INICIO

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

SetDirectory[NotebookDirectory[]]

DEFINICION  ELEMENTO CUADRILATERO REGULAR DE 9 NODOS

DEFINICION GRAFICA

CuaR9=;

CuaR9r=Show[CuaR9, ImageSize->250]

COORDENADAS NATURALES NODOS

Cn={{-1,-1},{1,-1},{1,1},{-1,1},{0,-1},{1,0},{0,1},{-1,0},{0,0}};

NNodos=Dimensions[Cn][[1]]

IMAGEN DEL ELEMENTO - COMPROBACION

FUNCION REPRESENTACION GRAFICA ELEMENTOS Y NODOS

ElementPlot[b_List,options___]:=Module[{asa,color,nr,circles,lines},asa=Select[{options},((!SameQ[#[[1]],NodeColor])&&(!SameQ[#[[1]],NodeSize]))&];
{color,nr}={NodeColor,NodeSize}/.{options}/.{NodeColor->GrayLevel[0],NodeSize->PointSize[0.06]};
circles=Map[Point[#]&,Partition[Flatten[b],2]];
lines=Line[Append[b[[1]],First[b[[1]]]]];
Show[Graphics[{nr,color,circles}],Graphics[lines],Evaluate[Sequence[##]&@@asa]]];

DEFINICION VECTOR DE NODOS

ptsexteriores={Cn[[1]],Cn[[5]],Cn[[2]],Cn[[6]],Cn[[3]],Cn[[7]],Cn[[4]],Cn[[8]]};

ptsinteriores={Cn[[9]]};

IMAGEN DE COMPROBACION

Imagen=ElementPlot[{ptsexteriores,ptsinteriores},AspectRatio->1,PlotRange->{{-1.5,1.5},{-1.5,1.5}},ImageSize-> 150,Frame->False,NodeColor->RGBColor[0,0,1]]

3. FUNCIONES DE FORMA - METODO PRODUCTO DE CURVAS

CURVAS A CONSIDERAR

Cu=Table[0,{i,6}];

LADOS

Cu[[1]]=(η+1);Cu[[2]]=(ξ-1);Cu[[3]]=(η-1);Cu[[4]]=(ξ+1);

MEDIANAS

Cu[[5]]=η;Cu[[6]]=ξ;

DEFINICION PRODUCTOS DE CURVAS EN CADA NODO

Nc=Table[0,{i,NNodos}];

Tipo 1 - ESQUINA

Nc[[4]]=Cu[[1]]*Cu[[2]]*Cu[[5]]*Cu[[6]]

Nc[[3]]=Cu[[1]]*Cu[[4]]*Cu[[5]]*Cu[[6]]

Nc[[2]]=Cu[[3]]*Cu[[4]]*Cu[[5]]*Cu[[6]]

Nc[[1]]=Cu[[2]]*Cu[[3]]*Cu[[5]]*Cu[[6]]

Tipo 2 - LADOS

Nc[[5]]=Cu[[2]]*Cu[[3]]*Cu[[4]]*Cu[[5]]

Nc[[6]]=Cu[[1]]*Cu[[3]]*Cu[[4]]*Cu[[6]]

Nc[[7]]=Cu[[1]]*Cu[[2]]*Cu[[4]]*Cu[[5]]

Nc[[8]]=Cu[[1]]*Cu[[2]]*Cu[[3]]*Cu[[6]]

Tipo 3 - INTERIOR

Nc[[9]]=Cu[[1]]*Cu[[2]]*Cu[[3]]*Cu[[4]]

OBTENCION FUNCIONES DE FORMA

Clear[Nf]

Nfp=Table[0,{i,NNodos}];

Nf=Table[0,{i,NNodos}];

Do[
Nfp[[i]]=a*Nc[[i]];
eq=1==Nfp[[i]]/.{ξ->Cn[[i,1]],η->Cn[[i,2]]};
as=a/.Solve[eq,a][[1]];Print["Nodo ",i];
Nf[[i]]=Simplify[Nfp[[i]]/.{a->as}],
{i,NNodos}
];

MatrixForm[Nf]

COMPROBACION SUMA UNIDAD

Suma=(∑_(i=1)^NNodos)Nf[[i]]

Simplify[%]

OK.

REPRESENTACION GRAFICA

Función Representación Gráfica Funciones de Forma

PlotQuadrilateralShapeFunction[xyquad_,f_,Nsub_,aspect_]:=Module[{Ne,Nev,line3D={},poly3D={},xyf1,xyf2,xyf3,i,j,n,ixi,ieta,xi,eta,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,xc,yc},{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4}}=Take[xyquad,4];xc={x1,x2,x3,x4};yc={y1,y2,y3,y4};Ne[xi_,eta_]:=N[{(1-xi)*(1-eta),(1+xi)*(1-eta),(1+xi)*(1+eta),(1-xi)*(1+eta)}/4];n=Nsub;Do[Do[ixi=(2*i-n-1)/n;ieta=(2*j-n-1)/n;{xi,eta}=N[{ixi-1/n,ieta-1/n}];Nev=Ne[xi,eta];xyf1={xc.Nev,yc.Nev,f[xi,eta]};{xi,eta}=N[{ixi+1/n,ieta-1/n}];Nev=Ne[xi,eta];xyf2={xc.Nev,yc.Nev,f[xi,eta]};{xi,eta}=N[{ixi+1/n,ieta+1/n}];Nev=Ne[xi,eta];xyf3={xc.Nev,yc.Nev,f[xi,eta]};{xi,eta}=N[{ixi-1/n,ieta+1/n}];Nev=Ne[xi,eta];xyf4={xc.Nev,yc.Nev,f[xi,eta]};AppendTo[poly3D,Polygon[{xyf1,xyf2,xyf3,xyf4}]];AppendTo[line3D,Line[{xyf1,xyf2,xyf3,xyf4,xyf1}]],{i,1,Nsub}],{j,1,Nsub}];Show[Graphics3D[RGBColor[1,1,0]],Graphics3D[poly3D],Graphics3D[Thickness[.002]],Graphics3D[line3D],Graphics3D[RGBColor[0,0,0]],Graphics3D[Thickness[.005]],Graphics3D[Line[xyquad]],PlotRange->All,BoxRatios->{1,1,aspect},Boxed->False]];

Representación Gráfica Funciones Forma Elemento.

Ng=Table[0,{i,NNodos}];

xyc1={0,0,0};xyc2={3,0,0};xyc3={3,3,0};xyc4={0,3,0};xyquad=N[{xyc1,xyc2,xyc3,xyc4,xyc1}];

Control de Cuadricula

Nsub=10;

Do[
fi[ξ_,η_]=Nf[[i]];
Ng[[i]]=PlotQuadrilateralShapeFunction[xyquad,fi,Nsub,1/2],
{i,NNodos}
];

4. RESULTADOS INTERACTIVOS - #

Manipulate[{CuaR9r,Ng[[n]],Nf[[n]]},{n,1,Dimensions[Nf][[1]],1},{n,Range[Dimensions[Nf][[1]]]},FrameLabel->{"FUNCION DE FORMA EN NODO n - CUADRILATERO REGULAR 9 NODOS"}, SaveDefinitions->True]

5. DERIVADAS FUNCIONES DE FORMA Y JACOBIANO

FUNCIONES DE FORMA

Nf

DERIVADAS FUNCIONES DE FORMA / COORDENADAS NATURALES

iNICIALIZACIONES

Clear[DNfξ,DNfη]

DNfξ=Table[0,{i,NNodos}];

DNfη=Table[0,{i,NNodos}];

PROCESO DE CALCULO

Do[
Nfe=Expand[Nf[[i]]];
DNfξ[[i]]=Simplify[D[Nf[[i]],ξ]];
DNfη[[i]]=Simplify[D[Nf[[i]],η]];
Clear[Nfe],
{i,NNodos}
];

RESULTADOS

DNfξ

DNfη

MODULO CALCULO FUNCIONES DE FORMA Y JACOBIANO

MODULO GENERICO A COMPLETAR

MODULO GENERICO A COMPLETAR: XX = No. Nodos

(*QuadXXIsoPShapeFunDer[ncoor_,qcoor_]:=Module[{Nf,dNx,dNy,dNξ,dNη,J11,J12,J21,J22,Jdet,ξ,η,x,y},{ξ,η}=qcoor;
Nf=(*Nf*);

dNξ=(*dNξ*);

dNη=(*dNη*);

x=Table[ncoor[[i,1]],{i,XX}];y=Table[ncoor[[i,2]],{i,XX}];
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}]];*)

MODULO COMPLETADO - ###

X;

6. MATRIZ DE RIGIDEZ - INTEGRACION NUMERICA

MODULO DE CALCULO DE LA MATRIZ DE RIGIDEZ

MODULO GENERICO A COMPLETAR

MODULO GENERICO: XX = No. Nodos, YY = Grados Libertad

(*QuadXXIsoPMembraneStiffness[ncoor_,mprop_,fprop_,options_]:=Module[{i,j,k,p=7,numer=False,Emat,th=1,h,qcoor,c,w,Nf,dNx,dNy,Jdet,B,Ke=Table[0,{YY},{YY}]},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}=QuadXXIsoPShapeFunDer[ncoor,qcoor];
If[Length[th]==0,h=th,h=th.Nf];c=w*Jdet*h;
B={Flatten[Table[{dNx[[i]],0},{i,XX}]],Flatten[Table[{0,dNy[[i]]},{i,XX}]],Flatten[Table[{dNy[[i]],dNx[[i]]},{i,XX}]]};
Ke+=Simplify[c*Transpose[B].(Emat.B)];];
Return[Simplify[Ke]]];*)

MODULO COMPLETADO - ###

X;

MODULO REGLAS DE CUADRATURA DE GAUSS

OPCION 1: DEFINICION EN MATHEMATICA

<<NumericalDifferentialEquationAnalysis`;

?GaussianQuadratureWeights

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}=GaussianQuadratureWeights[p1,-1,1][[i1]];
{eta,w2}=GaussianQuadratureWeights[p2,-1,1][[i2]];
info={{xi,eta},w1*w2};
If[numer,Return[N[info]],Return[Simplify[info]]];
];

QuadGaussRuleInfo[{25,False},1]

OPCION 2: DEFINICION DE CARLOS FELIPPA

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]},g5={-Sqrt[5+2*Sqrt[10/7]],-Sqrt[5-2*Sqrt[10/7]],0,Sqrt[5-2*Sqrt[10/7]],Sqrt[5+2*Sqrt[10/7]]}/3,w5={322-13*Sqrt[70],322+13*Sqrt[70],512,322+13*Sqrt[70],322-13*Sqrt[70]}/900,i=point,p=rule,info={Null,0}},
If[p==1,info={0,2}];
If[p==2,info={g2[[i]],1}];
If[p==3,info={g3[[i]],w3[[i]]}];
If[p==4,info={g4[[i]],w4[[i]]}];
If[p==5,info={g5[[i]],w5[[i]]}];
If[numer,Return[N[info]],Return[Simplify[info]]];
];

QuadGaussRuleInfo[{2,False},1]

7. TEST DEL RECTANGULO

DEFINICION DE LA GEOMETRIA

RectanguloT=;

RectanguloTr=Show[RectanguloT, ImageSize->250]

DEFINCION COORDENADAS NODOS ELEMENTO - ANSYS CLASSIC

IMAGEN DEL MODELO EN ANSYS

RectanguloA=;

RectanguloAr=Show[RectanguloA, ImageSize->450]

DATOS COORDENADAS NODALES EN ANSYS

Nodes={{0.`,0.`},{1.`,0.`},{0.5`,0.`},{1.`,0.5`},{1.`,0.25`},{0.`,0.5`},{0.5`,0.5`},{0.`,0.25`},{0.5`,0.25`}};

Dimensions[Nodes][[1]]

ORDENACION DATOS DE ANSYS

CuaT5r

RectanguloAr

ncoor={Nodes[[1]],Nodes[[2]],Nodes[[4]],Nodes[[6]],Nodes[[3]],Nodes[[5]],Nodes[[7]],Nodes[[8]],Nodes[[9]]}

DATOS COORDENADAS NODALES - SEGUN TEMA 23

a=1/2;

ncoor={{0,0},{2*a,0},{2*a,a},{0,a},{a,0},{2*a,a/2},{a,a},{0,a/2},{a,a/2}};

DEFINICION DEL MATERIAL

ClearAll[Em,nu,a,b,e,h,p,num];h=1;

DATOS SEGUN ENUNCIADO

Em=96*39*11*55*7;nu=1/3; (*isotropic material*)

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

Emat//MatrixForm

VERIFICACION DE LA MATRIZ DE RIGIDEZ

NUMERO DE PUNTOS DE GAUSS MINIMO PARA CONSEGUIR SUFICIENCIA DE RANGO

NF=NNodos*2.;

NG=((NF-3)/3)

Se necesitan como mínimo 5 Puntos -- Regla 3 x 3 minima

BUCLE GENERICO A COMPLETAR

BUCLE GENERICO: XX = No. Nodos, ZZ = GRADOS DE LIBERTAD - 3

(*For [p=1, p≤5, p++,
Ke=QuadXXIsoPMembraneStiffness[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p," x ",p];
Print["Ke=",Chop[Ke]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]],10^-7];
If[Valores[[ZZ]]≠0,Break[],Print["Valores proprios matriz Ke=",Valores]]
];
Print["Valores proprios matriz Ke=",Valores];
Print["tenemos la suficiencia de rango para p=",p]*)

DESARROLLO DE LA MATRIZ DE RIGIDEZ - ###

?Chop

X

For [p=4, p≤5, p++,
Ke=X[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p," x ",p];
Print["Ke=",Chop[Ke]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]],10^-7];
Print["Valores proprios matriz Ke=",Valores];
Print["tenemos la suficiencia de rango para p=",p]
];

8. TEST DEL TRAPECIO

DEFINICION DE LA GEOMETRIA

TrapecioT=;

TrapecioTr=Show[TrapecioT, ImageSize->250]

DEFINCION COORDENADAS NODOS ELEMENTO - ANSYS CLASSIC

IMAGEN DEL MODELO EN ANSYS

TrapecioA=;

TrapecioAr=Show[TrapecioA, ImageSize->450]

DATOS COORDENADAS NODALES EN ANSYS

Nodes={{0.`,0.`},{1.`,0.`},{0.5`,0.`},{0.5`,0.5`},{0.75`,0.25`},{0.`,0.5`},{0.25`,0.5`},{0.`,0.25`},{(0.25`+0.5`)/2,0.25`}};

Dimensions[Nodes][[1]]

ORDENACION DATOS DE ANSYS

CuaR4r

TrapecioAr

ncoor={Nodes[[1]],Nodes[[2]],Nodes[[4]],Nodes[[6]],Nodes[[3]],Nodes[[5]],Nodes[[7]],Nodes[[8]],Nodes[[9]]}

DEFINICION DEL MATERIAL

ClearAll[Em,nu,a,b,e,h,p,num];h=1;

DATOS SEGUN ENUNCIADO

Em=96*39*11*55*7;nu=1/3; (*isotropic material*)

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

Emat//MatrixForm

VERIFICACION DE LA MATRIZ DE RIGIDEZ

NUMERO DE PUNTOS DE GAUSS MINIMO PARA CONSEGUIR SUFICIENCIA DE RANGO

NF=NNodos*2.;

NG=((NF-3)/3)

Se necesitan como mínimo 5 Puntos -- Regla 3 x 32 minima

BUCLE GENERICO A COMPLETAR

BUCLE GENERICO: XX = No. Nodos, ZZ = GRADOS DE LIBERTAD - 3

(*For [p=1, p≤5, p++,
Ke=QuadXXIsoPMembraneStiffness[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p," x ",p];
Print["Ke=",Chop[Ke]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]]];
If[Valores[[ZZ]]≠0,Break[],Print["Valores proprios matriz Ke=",Valores]]
];
Print["Valores proprios matriz Ke=",Valores];
Print["tenemos la suficiencia de rango para p=",p]*)

DESARROLLO DE LA MATRIZ DE RIGIDEZ - ###

X

For [p=4, p≤5, p++,
Ke=X[ncoor,{Emat,0,0},{h},{True,p}];
Print["Gauss integration rule: ",p," x ",p];
Print["Ke=",Chop[Ke]//MatrixForm];
Valores=Chop[Eigenvalues[N[Ke]],10^-7];
Print["Valores proprios matriz Ke=",Valores];
Print["tenemos la suficiencia de rango para p=",p]
];