Plane Stress Analysis of Cantilever Beam

You will need Aladdin 2.0 to run this problem.

PROBLEM STATEMENT

Figure 1 shows an elevation and cross section view of a cantilever beam subject to an end moment couple of 10 kN forces spaced at 30 cm centers.

Figure 1 : Elevation and Cross Section Views of Cantilever Beam

The cantilever is 200 cm long, 30 cm high, and 5 cm thick, and is attached to a wall with full fixity (i.e., no rotations allowed). It is constructed from a material having Young's modulus E = 20000000 kN/m^2 and Poisson's ratio = 1/3.

This example presents a finite element analysis of the cantilever beam assuming plane-stress behavior. We will:

• Compute and plot the cantilever displacements due to the end-moment.
• Compute and plot components of stress throughout the cantilever beam.

FINITE ELEMENT MODELING

Finite Element Mesh

Figure 2 is a MATLAB plot of the finite element mesh.

Figure 2 : Finite Element Mesh for Cantilever Beam

```    beam_length  =  200 cm;
beam_width   =    5 cm;
beam_height  =   30 cm;

node = 0;
x = 0 cm ; y = 0 cm;
while(x <= beam_length) {
node = node + 1;

y    = y + beam_height/2;
node = node + 1;

y    = y + beam_height/2;
node = node + 1;

x = x + beam_length/10;
y = y - beam_height;
}
```

creates a rectangular block of 33 nodes beginning at (x,y) = (0 cm, 0 cm) and finishing at (x,y) = (100 cm, 30 cm). The block of code:

```    i = 1; elmtno = 1;
while(elmtno <= 19) {
AddElmt( elmtno, [   i, i+3, i+4, i+1 ], "cantilever_elmt_attr");
elmtno = elmtno + 1;

AddElmt( elmtno, [ i+1, i+4, i+5, i+2 ], "cantilever_elmt_attr");
elmtno = elmtno + 1;
i = i + 3;
}
```

systematically attaches twenty plane-stress/plain-strain finite elements to the nodes.

Section and Material Properties

The section and material properties are:

```    ElementAttr("cantilever_elmt_attr") { type     = "PLANE_STRESS";
section  = "cantilever_section";
material = "cantilever_material"; }

SectionAttr("cantilever_section") { depth = beam_height;
width = beam_width; }

MaterialAttr("cantilever_material") { poisson = 1.0/3.0;
E       = 20000000 kN/m^2; }
```

Boundary Conditions

We assume that at its base the cantilever is fully-fixed in the horizontal direction. Node 1 is also fixed in the y-direction. Nodes 2 and 3, however, are allowed to displace in the vertical direction i.e.,

```    FixNode( 1, [1,1] );
FixNode( 2, [1,0] );
FixNode( 3, [1,0] );
```

The bending moment couple is applied to the tip of the cantilever with:

```    Fx = 10.0 kN; Fy = 0.0 kN;

```

DISPLACEMENTS

Displacement Analysis

The block of statements:

```    stiff  = Stiff();
```

form the stiffness and external load matrices, and compute the displacements.

The statement:

```    PrintDispl(displ);
```

generates the output:

```==========================================
Node                  Displacement
No            displ-x           displ-y
==========================================
units                m                 m
1         0.00000e+00       0.00000e+00
2         0.00000e+00      -2.08696e-07
3         0.00000e+00      -4.24311e-19
4         1.66957e-06       1.11304e-06
5         7.93035e-20       9.04348e-07

.... details removed .....

31         1.66957e-05       1.11304e-04
32         1.89735e-19       1.11096e-04
33        -1.66957e-05       1.11304e-04
```

Displaced Shape of Cantilever Beam

Figure 3 shows the cantilever beam after the displacements have been amplified by a factor of 1000.

Figure 3 : Displaced shape of Cantilever Beam

The block of code:

```    scale = 1000.0;
for (ii = 1; ii <= 33; ii = ii + 1 ) {
coord     = GetCoord([ii]);
nodal_dof = GetDof([ii]);

if( nodal_dof[1][1] > 0 ) {
coord[1][1] = coord[1][1] + scale*displ[ nodal_dof[1][1] ][1];
}
if( nodal_dof[1][2] > 0 ) {
coord[1][2] = coord[1][2] + scale*displ[ nodal_dof[1][2] ][1];
}

print ii, coord[1][1], coord[1][2], "\n";
}
```

systematically retrieves the nodal coordinates and d.o.f. associated with each node from the Aladdin database, and then computes the absolute coordinate of the node after the displacements have been scaled by 1000. The abbreviated output from this block is:

```     1          0 cm          0 cm
2          0 cm      14.98 cm
3          0 cm         30 cm
4      20.17 cm     0.1113 cm
5         20 cm      15.09 cm

.... details removed .....

30      178.5 cm      39.02 cm
31      201.7 cm      11.13 cm
32        200 cm      26.11 cm
33      198.3 cm      41.13 cm
```

This array contains the nodal coordinates used to plot Figure 3.

STRESS CONTOURS

Figure 4 shows the distribution of stresses in the horizontal direction throughout the cantilever beam.

Figure 4 : Distribution of sigma_xx stresses

We can use elementary beam theory to verify the accuracy of the results. Assuming that the applied moment is

```    M = 10 kN x 0.3 m.
```

and that the beam inertia

```    I = 1/12 * (0.3 m)^3
```

it follows that the maximum/minimum stresses are:

```    stress = M.y_max/I = 2.0 x 10^5 Pa.
```

which is in the ball-park of stresses for dark blue/red shown along along top/bottom fibers of Figure 4. (Actually, the extrapolated stresses are 1.77 x 10^5 Pa).

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 <= 20; 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";
}
```

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 simply reorders the extrapolated stresses so that they match the nodal points.

The abbreviated output from this block of code is:

```      1  1.774e+05  1.774e+05 -1.043e+04 -1.043e+04
2  1.043e+04  1.043e+04 -1.774e+05 -1.774e+05
3  1.774e+05  1.774e+05 -1.043e+04 -1.043e+04
4  1.043e+04  1.043e+04 -1.774e+05 -1.774e+05

.... details removed .....

18  1.043e+04  1.043e+04 -1.774e+05 -1.774e+05
19  1.774e+05  1.774e+05 -1.043e+04 -1.043e+04
20  1.043e+04  1.043e+04 -1.774e+05 -1.774e+05
```

These are the stresses plotted in Figure 4.