Plane Strain Analysis of Pipe Cross Section

[ Problem Statement ] [ Finite Element Mesh ] [ Displacements ] [ Stress Contours ]
[ Input/Output Files ]

You will need Aladdin 2.0 to run this problem.


PROBLEM STATEMENT

Figure 1 shows an elevation and cross sectional view of a pipe

Figure 1 : Pipe Elevation and Cross Section Views

loaded by two-line loads P = 12 kN/m. The inside and outside diameters of the pipe are 20 cm and 30 cm, respectively. We assume that the pipe is fully-fixed in its longitudinal direction (i.e., strains in the longitudinal direction are zero). The material properties are Young's modulus = 200 GPa and Poisson's ration = 1/3.

This example presents a finite element analysis of the steel pipe cross-section, assuming plane-strain behavior. We will:


FINITE ELEMENT MESH

Because the pipe cross section has axes of symmetry on both the horizontal and vertical axes, we need only model 1/4 of the pipe cross section.

Figure 2 : Finite Element Mesh for Pipe Cross Section

Figure 2 is a MATLAB plot of the 1/4 pipe cross section. Although the finite elements will be represented by a (x,y) cartesian coordinate system, the mesh is most easily generated with a (radius, angle) polar coordinate system. The block of code:

    radius_min  = 10.0 cm;   angle_min  =    0.0; 
    radius_max  = 15.0 cm;   angle_max  =   PI/2;
    radius_incr =  0.5 cm;   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;
   }

generates the circular grid of 99 nodes

    =================================================
    Node                 x (cm)                y (cm)
    =================================================
       1                  10.0                    0.0
       2                  10.5                    0.0
       3                  11.0                    0.0

       ...........

      98                   0.0                   14.5
      99                   0.0                   15.0
    =================================================

We use the test:

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

to ensure that nodes along the horizontal and vertical axes are "exactly" positioned.

Eighty finite elements are attached to the nodes with the block of commands:

    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" );
        }
    }

Section and Material Properties

The block of code:

    ElementAttr("pipe") { type     = "PLANE_STRAIN";
                          section  = "mysection";
                          material = "mymaterial"; }

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

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

defines the material properties and plain-strain modeling assumption for the "pipe" element attribute.

Boundary Conditions

We assume that nodes along the vertical axis (i.e., x = 0) are restrained in the x-direction but can freely move in the y-direction. The block of commands:

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

applies this boundary condition to nodes 1 through 11.

Similarly, nodes along the horizontal axis (i.e., y = 0) are fixed in the vertical direction, but allowed to freely displace in the horizontal direction.

After the boundary conditions have been applied, the finite element model has 176 degrees of freedom.

External Loads

We assume that each half of the pipe section carries 6 kN/cm, exactly, ahd that the loads are carried by equally by nodes in the two uppermost finite elements. The block of commands:

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

applies nodal loads of -3.0 kN to nodes 88 and 89.


DISPLACEMENTS

Only three Aladdin statements:

    stiff = Stiff();
    eload = ExternalLoad();
    displ = Solve(stiff, eload);

are needed to for the stiffness matrix, external load vector, and compute the displacements. The command:

    PrintDispl(displ);

generates the output:

   ==========================================
    Node                  Displacement       
      No            displ-x           displ-y
   ==========================================
    units                 m                 m 
      1        -2.60880e-08       0.00000e+00
      2        -2.63854e-08       0.00000e+00
      3        -2.65135e-08       0.00000e+00
      4        -2.63598e-08       0.00000e+00

      ........ 

     97         0.00000e+00      -4.45024e-08
     98         0.00000e+00      -5.45846e-08
     99         0.00000e+00      -9.30605e-08
   ==========================================

As expected, the nodes along the vertical axis move downwards under the vertically applied loads. Intuition also dictates that the nodes along the horizontal axis will move outwards -- the numerical results indicate, however, a small displacement inwards.


STRESS CONTOURS

Plotting the Stress Contours

Figures 3 and 4 show the distribution of finite element stresses in the horizontal and vertical directions, respectively.

Figure 3 : Distribution of sigma_xx Stresses

Figure 4 : Distribution of sigma_yy Stresses

Rather that produce a smooth countour map of element stress components over the pipe cross section, we have deliberately chosen to plot the stresses within each finite element, and to highlight discontinuities across element boundaries.

While it is probably not a good idea to rely on the accuracy of stresses within any particular element, collectively the distributions of stresses can be used to validate our calculations. You should notice in Figure 3, for example, that the average color of the elements along the vertical axis is almost the same as those elements along the horizontal axiis. Since the pipe wall thickness is the same in both cases, we can deduce that in an overall sense, forces in the horizontal direction sum to zero.

Similarly, the stress contours in Figure 4 can be used to validate equilibrium of forces in the vertical direction. From the color code bar we see that "dark blue" corresponds to a vertical stress of about -1.2 x 10000 Pa. Given that the thickness of the pipe wall is 5 cm, it follows that the external loads should be approximately:

    External load = 1.2 x 10000 Pa * 5 cm 
                  = 6.0 kN / m.

Hence the vertical stresses along the horizontal axis are balanced by the externally applied loads a nodes 88 and 99.

Extrapolation of Stresses to Nodal Coordinates

By default the isoparametric plane-strain/plane-stress element will compute and print stress components at the gaussian points of integration. This is where evaluation of the stresses is most accurate.

The block of code:

   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;
   T = Inverse((1/12)*M);

sets up a (4x4) transformation matrix T that will extrapolate stresses at the 4 internal gauss points out to the nodal coordinates. With the matrix T in place, we can now systematically use

    actions = GetStress( [ii], displ );
to retrieve the stresses for element "ii" and then compute and print the "approximate" stresses at the nodal coordinates. The complete block of code is:
   print "\n";
   print "Element Stresses (sigma_xx at the nodal points)\n";
   print "===============================================\n\n";

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

        /* retrieve 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";
   }

Points to note:

  1. An unexpected problem we encountered with this method is the mismatch in ordering of the element nodes and gauss integration points. The four-by-four matrix of zeros and ones:
                     [ 0, 0, 1, 0;
                       0, 0, 0, 1;
                       1, 0, 0, 0;
                       0, 1, 0, 0 ]
    

    simply reorders the extrapolated stresses so that they match the nodal points.

  2. The Aladdin function QDimenLess() strips the units from the stress components contained in the array "extrap" so that the numerical values are suitable for plotting with MATLAB.


INPUT AND OUTPUT FILES


Developed in August 1997 by Clara Popescu and Mark Austin
Last Modified March 10, 1998
Copyright © 1998, Mark Austin, Department of Civil Engineering, University of Maryland