EJEMPLO No 2  == PLACA RECTANGULAR TP = UN ELEMENTO 9 NODOS - v2018

Cell 12:  Main program for verifying  9-node quad implementation with one element

1.- DATOS MALLA DE ELEMENTOS FINITOS

NODOS - ###

NodeCoordinates=N[{{0,6},{0,3},{0,0},X,{5/2,3},
     {5/2,0},{5,6},X,{5,0}}];  
PrintPlaneStressNodeCoordinates[NodeCoordinates,"",{6,4}];

numnod=Length[NodeCoordinates];

ELEMENTOS - ###

ElemNodes={{1,3,9,X,2,6,X,4,5}};

numele=Length[ElemNodes];

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

VISUALIZACION DE LA MALLA DE ELEMENTOS FINITOS - ###

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

2.- DATOS DEL PROBLEMA TENSION PLANA

MATERIAL

ClearAll[Em,ν,th];  
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,     {numele}];
PrintPlaneStressElementMatFab[ElemMaterials,X,"",{}];

ASIGNACION DE CONDICIONES DE CONTORNO EN DESPLAZAMIENTOS

INICIALIZACION

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

DEFINICION CONDICIONES DE CONTORNO EN DESPLAZAMIENTOS - ###

NodeDOFValues[[1]]=NodeDOFValues[[7]]={0,25}; 
NodeDOFValues[[4]]=X;              (* nodal loads *)
NodeDOFTags[[1]]=NodeDOFTags[[2]]={1,0}; (* vroller @ nodes 1,2 *)
NodeDOFTags[[3]]={1,1};                  (* fixed node 3 *)
NodeDOFTags[[6]]=NodeDOFTags[[9]]={0,1}; (* hroller @ nodes 6,9 *)

LISTADO DE CONDICIONES DE CONTORNO

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

3.- SOLUCION DEL PROBLEMA Y VISUALIZACION DE RESULTADOS

SOLUCION DEL PROBLEMA - ###

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

IMPRESION DE RESULTADOS - ###

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

4. - VISUALIZACION DE LOS DESPLAZAMIENTOS NODALES

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

ueps=10.^(-3); nbands=30;
ux=Table[NodeDisplacements[[n,1]],{n,X}]; 
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,ux,{uxmin,uxmax,uxinc},
  {True,True,True,False,True,True},{2,2},X,"Displacement component ux"];

Print["uymin,uymax,uyinc=",{uymin,uymax,uyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,X,uy,{uymin,uymax,uyinc},
  {True,True,True,False,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=30;
sxx=Table[NodePlateStresses[[n,1]],{n,numnod}]; 
syy=Table[NodePlateStresses[[n,2]],{n,numnod}];
sxy=Table[X[[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,ElemNodes,X,{sxxmin,sxxmax,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,syyinc},
  {True,True,True,False,True,True},{2,2},aspect,"Stress sigma-yy"];

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