# Cantilever Beam with Material Nonlinearity

You will need Aladdin 2.0 to run this problem.

### PROBLEM STATEMENT

Figure 1 shows a two dimensional cantilever beam that carries a monotonically increasing point load at its free end.

#### Figure 1 : Cantilever Beam subject to Tip Load

The cantilever is constructed from a material having bi-linear behavior. The section dimensions and material properties are as shown in the lower half of Figure 1. In this example we will compute:

1. The vertical deflection of cantilever tip for loads varying from 0 lbf to 60 lbf.

2. Contours of cantilever deflection, section rotation, and curvature at loading increments 0 lbf, 10 lbf ... 60 lbf.

### FINITE ELEMENT MODELING

Nodes and Elements

The block of code:

```   /* Generate grid of nodes for f.e. model */

x = 0 in; y = 0 in;
for( node=1 ; x<=50in ; node=node+1 ) {
AddNode( node, [x, y] );
x = x + 5 in;
}

/* Attach elements to grid of nodes */

for( elmt_no = 1 ; elmt_no <= total_elmt ; elmt_no = elmt_no + 1) {
AddElmt( elmt_no, [ elmt_no, elmt_no+1 ], "name_of_elmt_attr" );
}
```

generates the line of beam finite elements shown in the upper half of Figure 2.

#### Figure 2 : Finite Element Mesh for Cantilever Beam

The cantilever is modeled with 10 "FIBER_2D" elements, each element contains 40 fibers that stretch along the length of the finite element. Behavior of the element along its length is computed via integration at 5 sections. The overall deformations in each element are described by only six degrees of freedom, however. After the wall support has been fully-fixed, i.e.,

```    FixNode( 1, [ 1, 1, 1] );
```

the structural model has 30 d.o.f.

Finite Element Attributes

The finite element attribute declaration:

```   ElementAttr("name_of_elmt_attr") { type     = "FIBER_2D";
section  = "mysection";
material = "mymaterial";
fiber    = "fibername";
}
```

establishes the finite element type, FIBER_2D, and the section, material, and fiber layout for the finite element.

Section and Material Properties

The section and material properties are given by:

```   h = 4.0 in;  b = 1.0 in;
ks = 20000 psi;  kt = 0.1*ks;  fy = 500 psi;

SectionAttr("mysection") { area  = b*h;
width = b;
depth = h;
shear_factor  = 1.2;
}

MaterialAttr("mymaterial") { density = 2 lbf*sec^2/in/in^3;
poisson = 0.25;
E     = ks;
Et    = kt;
yield = fy;
shear_yield = 10000*fy;
}
```

Fiber Geometry

The nonlinear beam behavior is modeled with a stack of 40 fiber elements, as shown in the cross section view of the FIBER_2D finite element. This fiber layout is generated with:

```   no_fiber = 40; no_fiber_type = 1; dh = h/no_fiber;

fcoord = Matrix([ 1, no_fiber ]);
farea  = Matrix([ 1, no_fiber ]);
fmap   = Matrix([ 1, no_fiber ]);

for( i=1 ; i <= no_fiber ; i=i+1 ) {
fcoord[1][i] = h/2 - dh/2 - (i-1)*dh;
farea[1][i]  = b*h/no_fiber;
fmap[1][i]   = 1;
}

fattr  = Matrix([ 3, no_fiber_type ]);
for( j=1 ; j <= no_fiber_type ; j=j+1 ) {
fattr[1][j] = ks;
fattr[2][j] = kt;
fattr[3][j] = fy;
}

FiberAttr(no_fiber, "fibername" ){ FiberMaterialAttr = fattr;
FiberCoordinate   = fcoord;
FiberArea         = farea;
FiberMaterialMap  = fmap;
}
```

In this particular example, "fcoord [i][1]" stores the y-coordinate of the i-th fiber element. All of the fibers have a cross section area:

```    b*h/no_fiber = 1 in * 4 in / 40 fibers = 0.1 in^2.
```

and the same values of "ks", "fy" and "kt".

### DISPLACEMENTS

Figure 3 shows the vertical displacement of the cantilever tip (in) versus applied load (lbf).

#### Figure 3 : Vertical Tip Displacement (in) versus Applied Load (lbf)

Displacement Analysis

The abbreviated block of code:

```   total_step = 60;
for( step=1 ; step <= total_step ; step=step+1 ) {

dPk = [ 0 lbf, 1 lbf, 0 lbf*in ];
NodeLoad( total_node, dPk );

dP    = P_new - P_old;
P_old = P_new;

/* Newton-Raphson Iteration for State Determination */

while( L2Norm(dP) > 0.001 ) {

dp    = Solve( Ks, dP );
displ = displ + dp;

ElmtStateDet( dp );

Ks = Stiff();               /* Compute new stiffness    */
PR = InternalLoad( displ ); /* Compute internal load    */
dP = P_new - PR;            /* Compute unbalanced force */

}

/* Update element history data */

UpdateResponse();

/* Save external load, and tip displacement and rotation */

tip_response[step+1][1] = P_new[max_dof][1];
tip_response[step+1][2] = displ[max_dof][1];
tip_response[step+1][3] = displ[max_dof+1][1];

/* store all element results every 10 steps */

..... details removed .....

}
```

systematically increases the applied load at the cantilever tip from 1 lbf to 60 lbf in increments of 1 lbf, and computes and saves the beam response due the external loads.

The built-in function "ElmtStateDet()" to performs of the element state determinations. Similarly, after each load step finished, the element stress, strain, and material information is updated by the built-in function "UpdateResponse()".

The matrix "tip_response" stores the applied external load, tip displacement, and tip rotation for each step, including the initial non-loaded status. Matices "displacement" and "rotation" store the nodal displacement and rotation once every 10 steps.

Incipient Yielding

We can use elementary structural mechanics theory to estimate the external load that will cause incipient yielding of the beam fiber.

Assuming that the distribution of stresses in the beam can be approximated with:

```            M . y
sigma = -----  and  M = P.L  and  I = 1/12 b.h^3.
I
```

then first yielding of the fibers occurs in the outer layers of the beam at the fixed end when the tip load reaches 26.7 lbf. The coordinates of first yield are x = 0 in and y = 2 in and y = -2 in.

Of course the beam deflection grows rapidly after yielding.

Accuracy of Deflection Calculations

This fiber element includes effects of shear deformation. When P = 20 lbf, for example, the calculated tip deflection is 7.85489 in. Theoretical considerations [1] indicate that the exact solution for the vertical displacement of the cantilever is:

```         P.L^3                       f_s.P.L
v_b = ----- = 7.8125 in     v_s = ------- = 0.0375 in.
3.EI                        G.A.
```

The combined deflection is:

```    v_t = v_b + v_s = 7.85 in.
```

where "f_s" is the shear factor. Neglecting the shear deformation effects results in a relative deflection error of: If the relative shear deformation effects are not considered, like in the FRAME_2D element, the relative numerical error will be

```          v_s
err = ---  = 0.478 %
v_t
```

The relative numerical error for this example is

```          delta - v_t
err = ----------- = 0.0623 %
v_t
```

which is much less than the traditional plane frame element. In this example the numerical error of FRAME_2D is small because the beam is slender (i.e., h/L < 1/10). Shear deformation would make a much larger contribution to the overall displacement for deep beams with large h/L ratios.

Contours of Beam Displacement

Figure 4 shows how the deflected shape of the beam changes as a function of applied load.

#### Figure 4 : Nodal Deflections along Cantilever Beam

We only display the cantilever shape for applied loads 0 lbf, 10 lbf, 20 lbf and so forth. The block of Aladdin code:

```    flag = flag+1;
if( flag == 0 ) {

for( node=2 ; node<=total_node ; node=node+1 ) {
node_displacement         = GetDispl( [node] , displ);
displacement[node][index] = node_displacement[1][2];
rotation[node][index]     = node_displacement[1][3];
}

for(elmt_no=1 ; elmt_no<=total_elmt ; elmt_no=elmt_no+1){
curvature[elmt_no][index] = (rotation[elmt_no+1][index]
- rotation[elmt_no][index])/(5 in);
}
flag  = -10;
index = index+1;
}
```

uses the "flag" variable to detect increments of 10 lbf external loading. The vertical displacement and rotation of nodes 2 through 10 are retrieved and saved in the arrays "displacement" and "rotation", respectively.

### ROTATIONS

Figure 5 shows the tip rotation as a function of applied tip load.

#### Figure 5 : Tip Rotation (rad) versus Applied Load (lbf)

You should observe that this curve is smooth even though the material behavior is bi-linear. This behavior can be attributed to the gradual spread of fiber yielding and from the exterior fibers towards the beam centroid.

### CURVATURE

We can easily calculate the curvature for fiber element "i" by

```                theta_i+1 - theta_i
kappa_i =  -------------------
L_i
```

because curvature along an element is constant. Here "L_i" is the length of element "i", "theta_i" and "theta_i+1" are the end rotations. The results are shown in Figure 6.

#### Figure 6 : Element Curvature along Cantilever Beam

You should observe that in regions of elastic behavior - this includes all of the contours for loading 10 lbf and 20 lbf -- the element curvature will be linearly proportional to the section bending moment. The element curvature will increase significantly when the extreme element fiber yields -- one can therefore estimate which elements are subject to partial and almost total section plasticity.

References

1. Gere J.M., Timoshenko S.P., Mechanics of Materials, Wadsworth Inc, 1984.

### INPUT/OUTPUT FILES

• Click here to visit the input file.
• Click here to visit the output file.

Developed in July 1997 by Wane-Jang Lin and Mark Austin