EJEMPLO No 4  == PLACA CON AGUJERO TP = ELEMENTOS 9 NODOS - v2018

1.- DATOS MALLA DE ELEMENTOS FINITOS

NODOS - ###

s={1,0.70,0.48,0.30,0.16,0.07,0.0};
xy1={0,6};xy7={0,1};xy8={2.5,6};xy14={Cos[3*Pi/8],Sin[3*Pi/8]};
xy8={2.5,6};xy21={Cos[Pi/4],Sin[Pi/4]};xy15={5,6};
xy22={5,2};xy28={Cos[Pi/8],Sin[Pi/8]};xy29={5,0};xy35={1,0};
NodeCoordinates=Table[{0,0},X];
Do[NodeCoordinates[[n]]=N[s[[n]]        *xy1  +(1-s[[n]])           *xy7],  {n,1,7}];
Do[NodeCoordinates[[n]]=N[s[[n-7]]  *xy8  +(1-s[[n-7]])     *xy14],{n,8,14}];
Do[NodeCoordinates[[n]]=N[s[[n-14]]*xy15+(1-s[[n-14]])   *xy21],{n,15,21}];
Do[NodeCoordinates[[n]]=N[s[[n-21]]*xy22+(1-s[[n-21]])   *xy28],{n,22,28}];
Do[NodeCoordinates[[n]]=N[s[[n-28]]*xy29+(1-s[[n-28]])   *xy35],{n,29,35}];
PrintPlaneStressNodeCoordinates[NodeCoordinates,"",{6,4}];

numnod=Length[NodeCoordinates];

ELEMENTOS - ###

ElemNodes=Table[{0,0,0,0,0,0,0,0,0},{6}]; 
ElemNodes[[1]]={X};
Do [ElemNodes[[e]]=ElemNodes[[e-1]]+{2,2,2,2,2,2,2,2,2},{e,2,3}];
ElemNodes[[4]]=ElemNodes[[3]]+{10,10,10,10,10,10,10,10,10};
Do [ElemNodes[[e]]=ElemNodes[[e-1]]+{2,2,2,2,2,2,2,2,2},{e,5,6}];

numele=Length[ElemNodes]

ElemTypes=       Table["Quad9",{numele}]; 
PrintPlaneStressElementTypeNodes[ElemTypes,ElemNodes,"",{}];

VISUALIZACION DE LA MALLA DE ELEMENTOS FINITOS - ###

aspect=6/5;
ProcessOptions={True};
Plot2DElementsAndNodes[NodeCoordinates,ElemNodes,aspect,
    "One element mesh - 4-node quad",True,X];

2.- DATOS DEL PROBLEMA TENSION PLANA

MATERIAL

ClearAll[Em,ν,th,Nsub];  
Em=10000; ν=.25;   Nsub=4;
Emat=Em/(1-ν^2)*{{1,ν,0},{ν,1,0},{0,0,(1-ν)/2}};

ASIGNACION DE MATERIAL Y ESPESOR A ELEMENTOS - ###

th=3;

ElemMaterials=   Table[Emat,   {numele}]; 
ElemFabrications=Table[th,     X];
PrintPlaneStressElementMatFab[ElemMaterials,ElemFabrications,"",{}];

ASIGNACION DE CONDICIONES DE CONTORNO EN DESPLAZAMIENTOS

INICIALIZACION

NodeDOFValues=NodeDOFTags=Table[{0,0},{numnod}];

DEFINICION CONDICIONES DE CONTORNO EN DESPLAZAMIENTOS - ###

NodeDOFValues[[1]]=NodeDOFValues[[15]]={0,37.5}; 
NodeDOFValues[[8]]={0,75}; (* nodal loads *)
Do[NodeDOFTags[[n]]=X,{n,1,7}];   (* vroller @ nodes 1-7 *)
Do[NodeDOFTags[[n]]=Y,{n,29,35}]; (* hroller @ node 4 *)

LISTADO DE CONDICIONES DE CONTORNO - ###

PrintPlaneStressFreedomActivity[X,NodeDOFValues,"",{}];

3.- SOLUCION DEL PROBLEMA Y VISUALIZACION DE RESULTADOS

SOLUCION DEL PROBLEMA - ###

{NodeDisplacements,NodeForces,NodePlateCounts,NodePlateStresses,
    ElemBarNumbers,ElemBarForces}= PlaneStressSolution[
    X,ElemTypes,ElemNodes,
    ElemMaterials,ElemFabrications,
    NodeDOFTags,NodeDOFValues,X];

IMPRESION DE RESULTADOS - ###

PrintPlaneStressSolution[X,NodeForces,NodePlateCounts,
      NodePlateStresses,"Computed Solution:",{}];

4. - VISUALIZACION DE LOS DESPLAZAMIENTOS NODALES

CALCULO DE LOS VALORES MAXIMOS Y MINIMOS DE LOS DESPLAZAMIENTOS - ###

ueps=10.^(-3); nbands=X;
ux=Table[NodeDisplacements[[n,1]],{n,numnod}]; 
uy=Table[NodeDisplacements[[n,2]],{n,numnod}];
{uxmax,uymax}=Abs[{Max[ux],Max[uy]}]+ueps; 
{uxmin,uymin}=Abs[{Min[ux],Min[uy]}]+ueps;
uxmax=Max[uxmax,uxmin];  uxmin=-uxmax;
uymax=Max[uymax,uymin];  uymin=-uymax;
{uxinc,uyinc}={uxmax-uxmin,uymax-uymin}/nbands;

VISUALIZACION DESPLAZAMIENTOS NODALES - X e Y - ###

Print["uxmin,uxmax,uxinc=",{uxmin,uxmax,uxinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,X,{uxmin,uxmax,uxinc},
  {X,True,True,False,True,True},{2,2},aspect,"Displacement component ux"];

Print["uymin,uymax,uyinc=",{uymin,uymax,uyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,X,{uymin,uymax,uyinc},
  {True,True,True,X,True,True},{2,2},aspect,"Displacement component uy"];

5. - VISUALIZACION DE LAS TENSIONES -  NODALES - NORMALES Y TANGENCIALES

CALCULO DE LOS VALORES MAXIMOS Y MINIMOS DE LAS TENSIONES NORMALES Y TANGENCIALES - ###

sigeps=10.^(-3); nbands=X;
sxx=Table[NodePlateStresses[[n,1]],{n,numnod}]; 
X=Table[NodePlateStresses[[n,2]],{n,numnod}];
sxy=Table[NodePlateStresses[[n,3]],{n,numnod}];
{sxxmax,syymax,sxymax}=Abs[{Max[sxx],Max[syy],Max[sxy]}]+sigeps; 
{sxxmin,syymin,sxymin}=Abs[{Min[sxx],Min[syy],Min[sxy]}]+sigeps;
sxxmax=Max[sxxmax,sxxmin];  sxxmin=-sxxmax;
syymax=Max[syymax,syymin];  syymin=-syymax;
sxymax=Max[sxymax,sxymin];  sxymin=-sxymax;
{sxxinc,syyinc,sxyinc}={sxxmax-sxxmin,syymax-syymin,sxymax-sxymin}/nbands;

VISUALIZACION TENSIONES NODALES - NORMALES Y TANGENCIALES - ###

Print["sxxmin,sxxmax,sxxinc=",{sxxmin,sxxmax,sxxinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,X,sxx,{sxxmin,X,sxxinc},
  {True,True,True,False,True,True},{2,2},aspect,"Stress sigma-xx"];

Print["syymin,syymax,syyinc=",{syymin,syymax,syyinc}];
ContourBandPlotNodeFuncOver2DMesh[X,ElemNodes,syy,{syymin,syymax,X},
  {True,True,True,False,True,True},{2,2},aspect,"Stress sigma-yy"];

Print["sxymin,sxymax,sxyinc=",{X,sxymax,sxyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,sxy,{sxymin,sxymax,sxyinc},
  {True,X,True,False,True,True},{2,2},aspect,"Stress sigma-xy"];