1.   PRE-PROCESSING

The Pre-Processing involve the building up of finite element model by nodal coordinates, user have to prepare nodal and element listing in order to be read in by program.

Table of nodal and elemental listing will be produced for user to counter check.  Nodal plot in colour will be displayed as well, shape functions, element properties could be read in for pre-processing.

1.1   Input Packages

<<readqm6.m;

1.2  Reading in and tabulation of Nodal Coordinates

ReadCo["node"];
OutputForm[TableForm[co,TableHeadings->{Automatic,
{"X - Coord","Y - Coord"}}]] >> resfile.m;

1.3   Reading in and tablution of Element Connectivity

ReadElem["elem"];
OutputForm[TableForm[elem,TableHeadings->{Automatic,
{"Node - 1","Node -2","Node - 3","Node - 4"}}]] >>>resfile.m

1.4  Reading in Shape Functions, Elasticity Matrix

(*  Reading in Shape Functions  -- dkl  *)
<<shqm6.m;

(*  Reading in Elasticity Matrix  --  dpss *)
e=1500;v=0.3; (* e is the young's modulus and v poisson ratio*)
<<dpss.m;

(* Declare thickness for Plane Stress *)
thick=1.  ; (* default thickness  *)

"Elasticity Matrix">>>resfile.m
OutputForm[MatrixForm[dpss]]>>>resfile.m

1.5 Generation and Assembly of Stiffness Matrix

FormGlobalStifMat[]

1.6 Specification of Boundary Conditions

(*  Reading in Displacement and Mechanical  
    Boundary conditions   *)

ReadDispAnsys["disp"];  ReadLoadAnsys["load"];
paload=Table[0,{i,totdof},{j,1}];
dadisp=Table[0,{i,totdof},{j,1}];

ApplyDispAnsys[];  ApplyLoadAnsys[];

1.7  Nodal Plot and Element Plot

(* Routine to do element plot  *)
ModelPlot2D[elem, co, ielem, inode, nnode];

ModelMesh2D[elem, co, ielem, inode, nnode];

2.  SOLUTION PHASE

In this phase, the assembled stiffness matrix with the known boundary conditions are called into user programmed routine to solve for the displacement unknowns.  Gauss Elimination Method is used for the solution of simultaneous equations. Pivoting and row interchange procedures are used to treat ill-conditioned matrix.

SolutionLinear[];

3. POST-PROCESSING

Post-Processing involved the print out of displacement and stress values of the analysis.  Displacement plots are produced to give a better picture of final shape with magnification factor adjustable.

3.1 Stress Computation and Print-out

StressCompu[];x=.;y=.;

StressFunction[];

3.2 New coordinates and Print-out

newco=Table[0,{i,inode},{j,2}];

(*  Set displacement manification *)
mag=1.;

Do[newco[[i,1]]=co[[i,1]]+dadisp[[i*2-1,1]]*mag;
   newco[[i,2]]=co[[i,2]]+dadisp[[i*2,1]]*mag;
,{i,inode}];


OutputForm[TableForm[newco-co,TableHeadings->
{Automatic,{"X - Disp",
 "Y - Disp"}}]]  >>> resfile.m

3.3 Displacement Plots

DispPlot2D[elem,newco,ielem,inode,nnode];

Show[pldisp,wiremesh,AspectRatio->Automatic];

3.4 Stress Plots

Plnssx[]

Show[%,wiremesh,Frame->False]

Plnssy[]

Show[%,wiremesh,Frame->False]

Plnssxy[]

Show[%,wiremesh,Frame->False]

3.5 Stress at Locations within Model

StressAtLocation[1.,2.]

4. THEORETICAL COMPARISION

From the theory of beam bending, the stress at node 8 could be calculated through the following formulae :

Moment at node 7 = M = force x distance
 	                     = 5 x 10 Nm
	                     = 50 Nm

Stress = ( M * y) / I          where   y = 1 m    and    I = (b d^3) / 12  =  0.6666 m^4

Theoretical Stress = ( 50 * 1 ) / 0.6666  N/m^2   =  75  N/m^2

theorystress=75;polystress=xpoly[0.,2.];
errorstress1=(75 - calstress )*100/(75);
errorstress2=(75 - polystress )*100/(75);

TableForm[{{theorystress,calstress,errorstress1,
polystress,errorstress2}},
TableHeadings -> {{"Node 8"}, {"Theoy Stress",
"Cal. Stress","% of Error","Poly. Stress","% of Error"}}]