Input file for Plane Strain Analysis of Pipe Cross Section

/*
 *  =====================================================================
 *  Plane strain analysis of circular pipe cross section.
 *
 *  Written By : Mark Austin                                 Spring, 1998
 *  =====================================================================
 */ 

/* Define problem specific parameters   */

NDimension         = 2;
NDofPerNode        = 2;
MaxNodesPerElement = 4;

StartMesh();

/* Generate grid of finite element nodes */

radius_min  = 10.0 cm;
radius_max  = 15.0 cm;
radius_incr =  0.5 cm;
angle_min   =     0.0; 
angle_max   =    PI/2;
angle_incr  =   PI/16;

node  = 0;
angle = angle_min;
while( angle <= angle_max ) {

   radius = radius_min;
   while( radius <= radius_max ) {
      node = node + 1;
      x = radius*cos(angle);
      y = radius*sin(angle);

      if ( abs(x) <= 0.0000001 cm ) { x = 0.0 cm; }
      if ( abs(y) <= 0.0000001 cm ) { y = 0.0 cm; }

      AddNode(node, [ x, y]);
      radius = radius + radius_incr;
   }

   angle = angle + angle_incr;
}

/* Attach elements to grid of nodes */

nodeno = 0; elmtno = 0;
for ( i = 0; i < 8; i = i + 1 ) {
    nodeno = 11*i;
    for ( j = 1; j <= 10; j = j + 1 ) {
        nodeno = nodeno + 1;
        elmtno = elmtno + 1; 
        AddElmt( elmtno, [ nodeno, nodeno + 1, nodeno + 12, nodeno + 11 ], "pipe" );
    }
}

/* Define element, section and material properties */

ElementAttr("pipe") { type     = "PLANE_STRAIN";
                      section  = "pipesection";
                      material = "pipematerial"; }

SectionAttr("pipesection") { depth = 30 cm;
                             width = 30 cm; }

MaterialAttr("pipematerial") { poisson = 1/3;   
                               E       = 200 GPa; }

/* Setup boundary conditions */

   for (ii = 1; ii <= 11; ii = ii + 1 ) {
        FixNode( ii , [0,1] );
   }
   for (ii = 89; ii <= 99; ii = ii + 1 ) {
        FixNode( ii , [1,0] );
   }

/* Add point nodal loads to end of cantilever */

   Fx =  0.0 kN; Fy = 3.0 kN;
   NodeLoad( 88, [  Fx, -Fy ]);
   NodeLoad( 99, [  Fx, -Fy ]);

/* Compile and print finite element mesh */ 

   EndMesh();
   PrintMesh();

/* Compute stiffness and external load matrices */

   stiff = Stiff();
   eload = ExternalLoad();

/* Solve static analysis problem */

   displ  = Solve(stiff, eload);

/* Print displacements in tidy format */

   PrintDispl(displ);

/* Setup matrix for exptrapolation of stresses to nodal coordinates */

   M = Zero([4,4]);

   M[1][1] = (sqrt(3) + 1)^2;
   M[1][2] = 2;
   M[1][3] = (sqrt(3) - 1)^2;
   M[1][4] = 2;
   M[2][1] = 2;
   M[2][2] = (sqrt(3) + 1)^2;
   M[2][3] = 2;
   M[2][4] = (sqrt(3) - 1)^2;
   M[3][1] = (sqrt(3) - 1)^2;
   M[3][2] = 2;
   M[3][3] = (sqrt(3) + 1)^2;
   M[3][4] = 2;
   M[4][1] = 2;
   M[4][2] = (sqrt(3) - 1)^2;
   M[4][3] = 2;
   M[4][4] = (sqrt(3) + 1)^2;

   PrintMatrix(M);
   T = Inverse((1/12)*M);
  
/* Systematically retrieve stresses from individual elements */

   print "\n";
   print "Element Stresses (sigma_xx at the nodal points)\n";
   print "===============================================\n\n";

   for( ii = 1; ii <= 80; ii = ii + 1 ) {

        /* get stresses and extrapolate out to nodal coordinates  */

        actions = GetStress( [ii], displ );
        extrap  = T*[ actions[1][3];
                      actions[2][3];
                      actions[3][3];
                      actions[4][3] ];

        /* map extrapolated stresses onto nodal coordinate system */

        extrap = [ 0, 0, 1, 0;
                   0, 0, 0, 1;
                   1, 0, 0, 0;
                   0, 1, 0, 0 ] * extrap;
       
        /* print extrapolated stresses in format for MATLAB plotting */

        print ii;
        print QDimenLess(extrap[1][1]);
        print QDimenLess(extrap[2][1]);
        print QDimenLess(extrap[3][1]);
        print QDimenLess(extrap[4][1]);
        print "\n";
   }

quit;