June 21, 2013

1 Introduction

2 Direct method

2.1 Examples using the direct beam stiffness matrix

2.1.1 Example 1

2.1.2 Example 2

2.1.3 Example 3

2.2 Adding more elements

2.2.1 Example 3 redone with 2 elements

3 Generating shear and bending moments diagrams

4 Finding the stiffness matrix using methods other than direct method

4.1 Virtual work method for derivation of the stiffness matrix

4.2 Potential energy (minimize a functional) method to derive the stiffness matrix

5 References

2 Direct method

2.1 Examples using the direct beam stiffness matrix

2.1.1 Example 1

2.1.2 Example 2

2.1.3 Example 3

2.2 Adding more elements

2.2.1 Example 3 redone with 2 elements

3 Generating shear and bending moments diagrams

4 Finding the stiffness matrix using methods other than direct method

4.1 Virtual work method for derivation of the stiffness matrix

4.2 Potential energy (minimize a functional) method to derive the stiffness matrix

5 References

A short review for solving the beam problem in 2D is given. The deflection curve, bending moment and shear force diagrams are calculated for a beam subject to bending moment and shear force using direct stiffness method and finite elements method. The problem is solved first by finding the stiffness matrix using the direct method and then using the virtual work method.

Starting with only one element beam subject to bending and shear forces. There are 4 nodal degrees of freedom. Rotation at the left and right nodes of the beam and transverse displacements at the left and right nodes. The following diagram shows the sign convention used for external forces. Moments are always positive when anti-clockwise direction and vertical forces are positive when in the positive direction.

The two nodes are numbered and from left to right. is the moment at the left node (node 1), is the moment at the right node (node 2). is the vertical force at the left node and is the vertical force at the right node.

The above diagram shows the signs for the applied forces directions when acting in the positive sense. Since this is a one dimensional problem the displacement field (the unknown being solved for) will be function of one independent variable which is the coordinate. The displacement field is in the vertical direction called . This is the vertical displacement of a point on the beam from the . The following diagram shows the notation used for the coordinates

Angular displacement at of the beam is then found using At the left node the degrees of freedom, or the displacements, are called and at the right node called . At an arbitrary location in the beam, the vertical displacement is and the rotation is .

The following diagram shows the displacement field

In the direct method of finding the stiffness matrix, the forces at the ends of the beam are found directly by the use of beam theory. In beam theory the signs are different from is shown in the first diagram above. Therefore, some of the moment and shear forces obtained using beam theory ( and in the diagram below) will have different signs when compared to the external forces. The signs are then adjusted to reflect the convention as shown in the diagram above using and .

For example, the external moment is opposite in sign to and reaction is opposite to . To illustrate this more, a diagram with both sign conventions is shown below.

The goal now is to obtain expressions for external loads and in the above diagram as functions of the displacements at the nodes .

In other words, the goal is to obtain an expression of the form where is the stiffness matrix, is the nodal forces or load vector, is the nodal displacement vector. In this case will be a matrix and is a vector and is a vector.

Starting with , it is in the same direction as shear force . But hence

But , (from beam theory), hence

But and where is radius of curvature, therefore

But and for small angle of deflection hence , and the above becomes

Before continuing, the following diagram illustrates the above derivation. This comes from beam theory.

Now will be obtained. is in the opposite sense of bending moment hence a negative sign is added giving . But hence

Now is calculated. It is in the opposite sense of shear force , hence a negative sign is added giving

But , hence

But and where is radius of curvature therefore

But and for small angle of deflection hence and the above becomes

Finally is in the same direction as so no sign change is needed. hence

The following is a summary of what was found so far. Notice that the above expressions are evaluates at and accordingly to obtain the nodal end forces vector

| (1) |

Now the RHS of the above is expressed as function of nodal displacements . To do that, the field displacement which is the transverse displacement of the beam is assumed to be a polynomial in of degree 3. Hence

Polynomial of degree was used since there are degrees of freedom, and a minimum of free parameters are needed. Hence

| (2) |

and

| (3) |

Assuming the length of the beam is , then

| (4) |

and

| (5) |

eqs (2-5) gives

Solving for gives

Hence the field displacement function from eq. (A) can now be written as function of the nodal displacements

or in matrix form

The above can be written as

Where the are called the shape functions. The shape functions are

Notice that

and

as expected. Also

and

and

and

and

and

The shape functions have thus been verified.

The stiffness matrix is now found by substituting eq (5A) into eq. (1), repeated below

Hence

| (6) |

But

and

and

and

and now do the second derivatives

and

and

and

Hence eq (6) becomes

or in matrix form, after evaluating the expressions above for and

The above now is in the form

Hence the stiffness matrix is

Knowing the stiffness matrix means knowing the nodal displacements given the forces at the nodes. The power of the finite element method now comes after all the nodal displacements are calculated by solving because the polynomial is now completely determined and hence and can now be evaluated for any along the beam and not just at its end nodes. Eq. 5A above can now be used to find the displacement and everywhere

Therefore, these are the steps to obtain

- Find an expression for
- Solve for
- Calculate by assuming is a polynomial. This gives the displacement to use to evaluate the transverse displacement anywhere on the beam and not just at the end nodes.
- Obtain to evaluate the rotation of the beam any where and not just at the end nodes.
- Obtain strain where is the gradient matrix
- Obtain stress from
- Obtain the bending moment diagram from
- Obtain shear force diagram from

Now the beam stiffness matrix is used to solve few beam problems. Starting with simple one span beam

A one span beam, a cantilever beam of length , with point load at the free end

The first step is to make the free body diagram and show all moments and forces at the nodes

is the given force. since there is no external moment at the right end. Hence for this system is

Now is an important step. The known end displacements from boundary conditions is substituted
into , and the corresponding row and columns from the above system of equations are
removed^{1} .
Boundary conditions indicates that there is no rotation on the left end (since fixed). Hence and
. Hence the only unknowns are and . Therefore the first and the second rows and columns are
removed, giving

Now the above is solved for . Let (steel), , , and , hence

P=400;

L=144;

E=30*10^6;

I=57.1;

A=(E*I/L^3)*[12 6*L -12 6*L;

6*L 4*L^2 -6*L 2*L^2;

-12 -6*L 12 -6*L;

6*L 2*L^2 -6*L 4*L^2

];

load=[P;P*L;-P;0];

x=A(3:end,3:end)\load(3:end)

x =

-0.2324

-0.0024

L=144;

E=30*10^6;

I=57.1;

A=(E*I/L^3)*[12 6*L -12 6*L;

6*L 4*L^2 -6*L 2*L^2;

-12 -6*L 12 -6*L;

6*L 2*L^2 -6*L 4*L^2

];

load=[P;P*L;-P;0];

x=A(3:end,3:end)\load(3:end)

x =

-0.2324

-0.0024

Therefore the vertical displacement at the right end is inches (downwards) and radians. Now that all nodal displacements are found, the field displacement function is completely determined.

Hence from eq. 5A

But inches, hence

To verify, let in the above

Here is a plot of the deflection curve for the beam

v=@(x) 3.992*10^-8*x.^3-1.6956*10^-5*x.^2

x=0:0.1:144;

plot(x,v(x),’r-’,’LineWidth’,2);

ylim([-0.8 0.3]);

title(’beam deflection curve’);

xlabel(’x inch’); ylabel(’deflection inch’);

grid

x=0:0.1:144;

plot(x,v(x),’r-’,’LineWidth’,2);

ylim([-0.8 0.3]);

title(’beam deflection curve’);

xlabel(’x inch’); ylabel(’deflection inch’);

grid

Now instead of removing rows/columns for the known boundary conditions, a is put on the diagonal. Starting again with

since and , then

And now the above system is solved as before. , , ,

P=400;

L=144;

E=30*10^6;

I=57.1;

A=(E*I/L^3)*[12 6*L -12 6*L;

6*L 4*L^2 -6*L 2*L^2;

-12 -6*L 12 -6*L;

6*L 2*L^2 -6*L 4*L^2

];

load=[0;0;-P;0]; %put zeros for known B.C.

A(:,1)=0; A(1,:)=0; A(1,1)=1; %put 1 on diagonal

A(:,2)=0; A(2,:)=0; A(2,2)=1; %put 1 on diagonal

A

A =

1.0e+07 *

0.0000 0 0 0

0 0.0000 0 0

0 0 0.0007 -0.0496

0 0 -0.0496 4.7583

sol=A\load %SOLVE

sol =

0

0

-0.2324

-0.0024

L=144;

E=30*10^6;

I=57.1;

A=(E*I/L^3)*[12 6*L -12 6*L;

6*L 4*L^2 -6*L 2*L^2;

-12 -6*L 12 -6*L;

6*L 2*L^2 -6*L 4*L^2

];

load=[0;0;-P;0]; %put zeros for known B.C.

A(:,1)=0; A(1,:)=0; A(1,1)=1; %put 1 on diagonal

A(:,2)=0; A(2,:)=0; A(2,2)=1; %put 1 on diagonal

A

A =

1.0e+07 *

0.0000 0 0 0

0 0.0000 0 0

0 0 0.0007 -0.0496

0 0 -0.0496 4.7583

sol=A\load %SOLVE

sol =

0

0

-0.2324

-0.0024

The same solution is obtained as before, but without the need to remove rows/column from the stiffness matrix. This method might be easier for programming than the first method of removing rows/columns.

The rest now is the same as was done earlier and will not be repeated.

The same example as above, but the vertical load now is placed in the middle of the beam

In using stiffness method, all loads must be on the nodes. The vector is nodal forces vector. Hence equivalent nodal loads are found for the load in the middle of the beam. The equivalent loading is the following

Therefore, the problem is as if it was the following problem

Now that equivalent loading is in place, we continue as before. Make a free body diagram showing all loads (including reaction forces)

Now the stiffness equation is written as

There is no need to determine and at this point since these rows will be removed due to boundary conditions and and hence those quantities are not needed to solve the equations. Remember that rows and columns are removed for the known boundary displacements before solving . Hence, after removing the first two rows and columns, the system becomes

Now the above is solved for using the same numerical values for as in the first example

P=400;

L=144;

E=30*10^6;

I=57.1;

A=(E*I/L^3)*[12 6*L -12 6*L;

6*L 4*L^2 -6*L 2*L^2;

-12 -6*L 12 -6*L;

6*L 2*L^2 -6*L 4*L^2

];

load=[-P/2;P*L/8];

x=A(3:end,3:end)\load

x =

-0.072630472854641

-0.000605253940455

L=144;

E=30*10^6;

I=57.1;

A=(E*I/L^3)*[12 6*L -12 6*L;

6*L 4*L^2 -6*L 2*L^2;

-12 -6*L 12 -6*L;

6*L 2*L^2 -6*L 4*L^2

];

load=[-P/2;P*L/8];

x=A(3:end,3:end)\load

x =

-0.072630472854641

-0.000605253940455

Therefore

This is enough to obtain as before. Now the reactions and can be determined if needed. Going back to the full , results in

Hence the first equation becomes

and since (steel) and and and , then

hence

and hence

Now all nodal reactions are found and the displacement field is found. The deflection curve can be plotted.

Since inches , then

The following is the plot

clear all; close all;

v=@(x) 1.9459*10^-8*x.^3-6.3047*10^-6*x.^2

x=0:0.1:144;

plot(x,v(x),’r-’,’LineWidth’,2);

ylim([-0.8 0.3]);

title(’beam deflection curve’);

xlabel(’x inch’); ylabel(’deflection inch’);

grid

v=@(x) 1.9459*10^-8*x.^3-6.3047*10^-6*x.^2

x=0:0.1:144;

plot(x,v(x),’r-’,’LineWidth’,2);

ylim([-0.8 0.3]);

title(’beam deflection curve’);

xlabel(’x inch’); ylabel(’deflection inch’);

grid

Now assume the beam is fixed on the left end as above, but simply supported on the right end , and the vertical load now at distance from the left end and distance from the right end, and a uniform distributed load of density lb/in on the beam.

Using the following values: (steel), , , lb/in.

In the above, the left end reaction forces are shown as and moment reaction as and the reaction at the right end as . Starting by finding equivalent loads for the point load and equivalent loads for for the uniform distributed load . All external loads must be transferred to the nodes for the stiffness method to work. Equivalent load for the above point load is

And equivalent load for the uniform distributed loading is

Hence, using free body diagram, with all the loads on it gives the following diagram (In this diagram is the reaction moment and are the reaction forces)

Now that all loads are on the nodes, the stiffness equation is applied

Now boundary conditions are applied. , hence the first, second and third rows/columns are removed giving

Hence

Substituting numerical values for the above as given at the top of the problem results in

Hence the field displacement is now found

And a plot of the deflection curve is

clear all; close all;

v=@(x) 3.7229*10^-7*x.^3-5.361*10^-5*x.^2

x=0:0.1:144;

plot(x,v(x),’r-’,’LineWidth’,2);

ylim([-0.5 0.3]);

title(’beam deflection curve’);

xlabel(’x inch’); ylabel(’deflection inch’);

grid

v=@(x) 3.7229*10^-7*x.^3-5.361*10^-5*x.^2

x=0:0.1:144;

plot(x,v(x),’r-’,’LineWidth’,2);

ylim([-0.5 0.3]);

title(’beam deflection curve’);

xlabel(’x inch’); ylabel(’deflection inch’);

grid

Finite elements generated displacements are smaller in value than the actual analytical values. To improve the accuracy, more elements are added. To add more element, the beam is divided into 2,3,4 and more beam elements. To show how this work, example 3 above is solved again using 2 elements. It will be found that displacement field becomes more accurate (It was compared this with an exact solution based on using the beam 4th order differential equation and found to be almost the same with only 2 elements)

The first step is to divide the beam into 2 element and number the degrees of freedom and global nodes as follows

There will be 6 total degrees of freedom. 2 at each node. Hence the stiffness matrix for the whole beam (including both elements) will be 6 by 6. For each element however, the same stiffness matrix will be used as above and that will remain 4 by 4.

The stiffness matrix for each element is found then the global stiffness matrix is constructed, then is solved as before. The first step is to move all loads to the nodes as before. This is done for each element. The formulas for equivalent loads remain the same, but becomes . The following diagram show the equivalent loading for

The equivalent loading for distributed load is

Now the above 2 diagrams are put together to show all equivalent loads with the original reaction forces to obtain the following diagram

Now for each element is constructed. Starting with the first element

And for the second element

Now the 2 systems above are added to obtain the global stiffness matrix equation giving

Now the boundary conditions are applied since the first node is fixed, and . Putting 1 on the diagonal of the stiffness matrix corresponding to these known boundary conditions results in

As was mentioned above, another method is to remove the rows/columns which results in

giving the same solution.

Hence, there are 3 unknowns to solve for. Once there are solved, for the first element and for the second element are determined. The following code displays the deflection curve for the above beam

clear all; close all;

P=400;

L=144;

E=30*10^6;

I=57.1;

m=200;

a=0.125*L;

b=0.375*L;

A=E*I/(L/2)^3*[1 0 0 0 0 0;

0 1 0 0 0 0;

0 0 24 0 0 6*L/2;

0 0 0 8*(L/2)^2 0 2*(L/2)^2;

0 0 0 0 1 0;

0 0 6*L/2 2*(L/2)^2 0 4*(L/2)^2];

A

load = [0;

0;

-(m*L/2) - 2*P*b^2*(L/2+2*a)/(L/2)^3;

-2*P*a*b^2/(L/2)^2;

0;

P*a^2*b/(L/2)^2+(m*(L/2)^2)/12]

sol=A\load

A =

1.0e+08 *

0.0000 0 0 0 0 0

0 0.0000 0 0 0 0

0 0 0.0011 0 0 0.0198

0 0 0 1.9033 0 0.4758

0 0 0 0 0.0000 0

0 0 0.0198 0.4758 0 0.9517

load =

0

0

-15075

-8100

0

87750

sol =

0

0

-0.2735

-0.0019

0

0.0076

P=400;

L=144;

E=30*10^6;

I=57.1;

m=200;

a=0.125*L;

b=0.375*L;

A=E*I/(L/2)^3*[1 0 0 0 0 0;

0 1 0 0 0 0;

0 0 24 0 0 6*L/2;

0 0 0 8*(L/2)^2 0 2*(L/2)^2;

0 0 0 0 1 0;

0 0 6*L/2 2*(L/2)^2 0 4*(L/2)^2];

A

load = [0;

0;

-(m*L/2) - 2*P*b^2*(L/2+2*a)/(L/2)^3;

-2*P*a*b^2/(L/2)^2;

0;

P*a^2*b/(L/2)^2+(m*(L/2)^2)/12]

sol=A\load

A =

1.0e+08 *

0.0000 0 0 0 0 0

0 0.0000 0 0 0 0

0 0 0.0011 0 0 0.0198

0 0 0 1.9033 0 0.4758

0 0 0 0 0.0000 0

0 0 0.0198 0.4758 0 0.9517

load =

0

0

-15075

-8100

0

87750

sol =

0

0

-0.2735

-0.0019

0

0.0076

The above solution gives inch (downwards displacement) and radians and radians. Now the polynomial is found for each element

The above polynomial is the transverse deflection of the beam for the region .

for the second element is found in similar way

Hence

Which is valid for . The following is a plot of the deflection curve using the above 2 equations

When the above plot is compared to the case with one element, It can be seen the deflection is larger now. Comparing the above to the analytical solution shows that now the deflection is almost exactly the same as the analytical solution. Hence by using 2 elements instead of one element, the solution now agrees with the analytical solution.

The following diagram shows the deflection curve of problem 3 above when using one element and two elements on the same plot to help illustrate the difference in the result more clearly.

The analytical deflection for the beam in problem 3 above (fixed on the left and simply supported at the right end) when there is uniformly loaded with lbs per unit length is given by , while the analytical deflection for the same beam but when there is a point load P at distance from the left end is given by , therefore, the analytical expression for deflection is given by the sum of the above expressions, giving

Where means it is zero when is negative. In other words

The following diagram shows a plot of the analytical deflection with the 2 elements deflection calculated using Finite elements above.

In the above, the blue dashed curve is the analytical solution, and the red curve is the finite elements solution using 2 elements. It can be seen that the finite element solution for the deflection is now in a very good agreement with the analytical solution.

After solving the problem using finite elements and obtaining the field displacement function as was shown in the above examples, the shear force and bending moments along the beam can be calculated. Since the bending moment is given by and shear force is given by then these diagrams can now be readily plotted as shown below for example 3 above using the result from the finite elements with 2 elements. Recall from above that

Hence

using and inch and the bending moment diagram plot is

The bending moment diagram clearly does not agree with the bending moment diagram that can be generated from the analytical solution show below (generated using my other program which solves this problem analytically)

The reason for this is because the solution obtained using the finite elements method is a 3 degree polynomial and after differentiating twice to obtain the bending moment () the result is a linear function in while in the analytical solution case, when the load is distributed, then solution is a 4 degree polynomial. Hence the bending moment will be quadratic function in in the analytical case.

Therefore, in order to obtain good approximation for the bending moment and shear force diagrams using finite elements, more elements will be needed.

There are 3 main methods to obtain the stiffness matrix

- Variational method (minimize a functional). This functional is the potential energy of the structure and loads.
- Weighted residual. requires the differential equation as starting point. Approximated in weighted average. Galerking weighted residual method is the most common for implementation.
- Virtual work method. Make the virtual work zero for arbitrary allowed displacement.

In virtual work method, a small displacement is assumed to occur. Looking at small volume element, the amount of work done by external loads to cause the small displacement is equal to amount of increased internal strain energy. Assuming the field of displacement is given by and assuming the external loads are given by acting on the nodes, hence these point loads will do work given by on that unit volume where is the nodal displacements. In all these derivation, only loads acting directly on the nodes are considered for now. In other words, body forces and traction forces are not considered in order to simplify the derivations.

The increase of strain energy is in that same unit volume.

Hence, for a unit volume

And for the whole volume

| (1) |

Assuming that displacement can be written as a function of the nodal displacements of the element results in

Therefore

Since then where is the strain displacement matrix , hence

Now from the stress-strain relation , hence

| (4) |

Substituting Eqs. (2,3,4) into (1) results in

Since and do not depend on the integration variables they can be moved outside the integral, giving

Since the above is true for any admissible then the only condition is that

Hence this is in the form , therefore

knowing allows finding by integrating over the volume. For the beam element though, the transverse displacement. This means . Recall from the above that for the beam element,

Hence

Hence

Which is the stiffness matrix found earlier.

This method is very similar to the first method actually. It all comes down to finding a functional, which is the potential energy of the system, and minimizing this with respect to the nodal displacements. The result gives the stiffness matrix.

Let the system total potential energy by called and let the total internal energy in the system be and let the work done by external loads acting on the nodes be , then

Work done by external loads have a negative sign since they are an external agent to the system and work is being done onto the system.

The internal strain energy is given by and the work done by external loads is , hence

| (1) |

Now the rest follows as before. Assuming that displacement can be written as a function of the nodal displacements , hence

Since then where is the strain displacement matrix , hence

Now from the stress-strain relation , hence

| (4) |

Substituting Eqs. (2,3,4) into (1) results in

Hence gives

Which is on the form which means that

As was found by the virtual work method.

- A first course in Finite element method, 3rd edition, by Daryl L. Logan
- Matrix analysis of framed structures, 2nd edition, by William Weaver, James Gere