/* * ===================================================================== * 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;