EJEMPLO No 1  == PLACA RECTANGULAR TP = UN ELEMENTO DE 4 NODOS - v2018

ClearAll[Em,ν,th];
Em=10000;ν=.25;th=3;aspect=6/5;Nsub=4;

- DEFINICION DEL MODELO DE EF -

NodeCoordinates=N[{{0,6},{0,0},{5,6},{5,0}}];
ElemNodeLists={{1,2,4,3}};
numnod=Length[NodeCoordinates];
numele=Length[ElemNodeLists];
ElemTypes=Table[{"Quad4"},{numele}];
ElemMaterial=Table[{Em,ν},{numele}];
ElemFabrication=Table[{th},{numele}];
FreedomValues=FreedomTags=Table[{0,0},{numnod}];

FreedomTags//MatrixForm

FreedomValues//MatrixForm

FreedomValues[[1]]=FreedomValues[[3]]={0,75};(*nodal loads*)

FreedomTags[[1]]={1,0};(*vroller@node 1*)

FreedomTags[[2]]={1,1};(*fixed node 2*)
FreedomTags[[4]]={0,1};(*hroller@node 4*)

FreedomTags//MatrixForm

ElemOptions={True};

aspect=6/5;

Plot2DElementsAndNodes[NodeCoordinates,ElemNodeLists,aspect,"MALLA UN ELEMENTO - CUADRILATERO 4 NODOS",True,True];

- SOLUCION DEL PROBLEMA E IMPRESION RESULTADOS -

{u,f,sig}=MembraneSolution[NodeCoordinates,ElemTypes,ElemNodeLists,ElemMaterial,ElemFabrication,ElemOptions,FreedomTags,FreedomValues];

Print["Computed displacements=",u];
Print["Recovered forces=",f];
Print["Avg nodal stresses=",sig];

sig //MatrixForm

- VISUALIZACION GRAFICA COMPONENTES DE LOS DESPLAZAMIENTOS -

ux=uy=Table[0,{numnod}];
Do[ux[[n]]=u[[2*n-1]];uy[[n]]=u[[2*n]],{n,1,numnod}];
uxmax=uymax=0;
Do[uxmax=Max[Abs[ux[[n]]],uxmax];uymax=Max[Abs[uy[[n]]],uymax],{n,1,numnod}];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodeLists,ux,uxmax,Nsub,aspect,"Displacement component ux"];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodeLists,uy,uymax,Nsub,aspect,"Displacement component uy"];

- VISUALIZACION GRAFICA DISTRIBUCION TENSIONES NODALES PROMEDIADAS -

sxx=syy=sxy=Table[0,{numnod}];
Do[{sxx[[n]],syy[[n]],sxy[[n]]}=sig[[n]],{n,1,numnod}];
sxxmax=syymax=sxymax=0;
Do[sxxmax=Max[Abs[sxx[[n]]],sxxmax];
syymax=Max[Abs[syy[[n]]],syymax];
sxymax=Max[Abs[sxy[[n]]],sxymax],{n,1,numnod}];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodeLists,sxx,sxxmax,Nsub,aspect,"Nodal stress sig-xx"];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodeLists,syy,syymax,Nsub,aspect,"Nodal stress sig-yy"];
ContourPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodeLists,sxy,sxymax,Nsub,aspect,"Nodal stress sig-xy"];

Cell 11:  Main program for verifying  4-node quad implementation with one element

1.- DATOS MALLA DE ELEMENTOS FINITOS

NODOS - ###

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

numnod=Length[NodeCoordinates];

ELEMENTOS - ###

ElemNodes=  {{1,2,4,X}};

numele=Length[ElemNodes];

ElemTypes=       Table["Quad4",{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,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,ElemFabrications,"",{}];

ASIGNACION DE CONDICIONES DE CONTORNO EN DESPLAZAMIENTOS

INICIALIZACION

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

DEFINICION CONDICIONES DE CONTORNO EN DESPLAZAMIENTOS - ###

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

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,ElemTypes,ElemNodes,
    ElemMaterials,ElemFabrications,
    NodeDOFTags,NodeDOFValues,ProcessOptions];

IMPRESION DE RESULTADOS

PrintPlaneStressSolution[NodeDisplacements,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=30;
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,ux,{uxmin,uxmax,uxinc},
  {True,True,True,False,True,True},{2,2},aspect,"Displacement component ux"];

Print["uymin,uymax,uyinc=",{uymin,uymax,uyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,ElemNodes,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[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,ElemNodes,sxx,{sxxmin,sxxmax,sxxinc},
  {True,True,True,False,True,True},{2,2},aspect,"Stress sigma-xx"];

Print["syymin,syymax,syyinc=",{syymin,syymax,syyinc}];
ContourBandPlotNodeFuncOver2DMesh[NodeCoordinates,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,True},{2,2},aspect,"Stress sigma-xy"];