up

Derivation of stiffness matrix for a beam

Nasser M. Abbasi

June 21, 2013

Contents

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

1 Introduction

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.

2 Direct 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 y  direction.

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

PIC

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 x  coordinate. The displacement field is in the vertical direction called v (x )  . This is the vertical displacement of a point x  on the beam from the x - axis  . The following diagram shows the notation used for the coordinates

PIC

Angular displacement at x  of the beam is then found using        dv(x)
θ(x) = -dx-.  At the left node the degrees of freedom, or the displacements, are called v1,θ1   and at the right node called v2,θ2   . At an arbitrary location x  in the beam, the vertical displacement is v (x )  and the rotation is θ (x )  .

The following diagram shows the displacement field v(x)

PIC

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 (MB  and VB  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 M  and V  .

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

PIC

The goal now is to obtain expressions for external loads Mi  and Ri  in the above diagram as functions of the displacements at the nodes {d} = {v1,θ1,v2,θ2}T  .

In other words, the goal is to obtain an expression of the form {p} = [K ]{d} where [K ]  is the stiffness matrix, {p} = {V ,M  ,R  ,M }T
        1  1   2   2  is the nodal forces or load vector, {d} is the nodal displacement vector. In this case [K]  will be a 4 × 4  matrix and {p} is a 4×  1  vector and {d} is a 4× 1  vector.

Starting with V1   , it is in the same direction as shear force VB1   . But VB1 = dMdBx1-  hence

     dM
V1 = ---B1-
      dx

But               I
MB1  = - σ(x) y  , (from beam theory), hence

       Idσ-(x)
V1 = - y  dx

But σ (x) = Eε (x)  and ε(x) = -y-
        ρ  where ρ  is radius of curvature, therefore

          (  )
V1 = EI-d-  1-
       dx   ρ

But        (d2u)
1=  (---dx2)---
ρ    1+(du)2 3∕2
        dx   and for small angle of deflection du-  1
dx  hence 1= ( d2u)
ρ    dx2 , and the above becomes

        d3u(x)
V1 = EI -dx3---

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

PIC

Now M1   will be obtained. M1   is in the opposite sense of bending moment MB1   hence a negative sign is added giving M  =  - M
  1      B1   . But M    = - σ (x) I
  B1          y  hence

pict

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

pict

But              I
MB2  = - σ (x)y  , hence

     I-dσ(x)-
V2 = y  dx

But σ (x) = Eε (x)  and ε(x) = -y-
        ρ  where ρ  is radius of curvature therefore

            (  )
V2 = - EI-d-  1-
         dx   ρ

But        (d2w)
1=  (---dx2)----
ρ    1+(dw)2 3∕2
        dx   and for small angle of deflection dw  1
dx  hence 1 = (d2u)
ρ    dx2 and the above becomes

          d3u(x)
V2 = - EI -dx3---

Finally M2   is in the same direction as MB2   so no sign change is needed. MB2 = - σ (x ) Iy  hence

pict

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

               (          |    )
      (   )    ||  EI d3u(x)||    ||
      ||| V1|||    ||||      dx32 ||x=0 ||||
      { M1}    {  - EI ddxu2|x=0  }
{p} = || V2||  = || - EI d3u(x)||   ||
      |(   |)    ||||       d2x3| |x=L ||||
        M2     (   EI ddxu2||     )
                         x=L
(1)

Now the RHS of the above is expressed as function of nodal displacements v1,θ1,v2,θ2   . To do that, the field displacement v(x)  which is the transverse displacement of the beam is assumed to be a polynomial in   x  of degree 3. Hence

pict

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

v1 = v (x)|x=0 = a0
(2)

and

θ1 = θ (x)|   = a1
         x=0
(3)

Assuming the length of the beam is L  , then

v2 = v(x)|x=L = a0 + a1L + a2L2 + a3L3
(4)

and

                               2
θ2 = θ(x)|x=L =  a1 + 2a2L + 3a3L
(5)

eqs (2-5) gives

       (  )    (                      )    (              )  (  )
       |||v1|||    |||          a0          |||     1  0   0    0    |||a0|||
       {θ1}    {          a1          }    ||0  1   0    0 ||  {a1}
{d} =  ||v2|| =  || a0 + a1L + a2L2 + a3L3|| = |(1  L   L2  L3 |)  ||a2||
       |(θ |)    |(   a + 2a L + 3a L2   |)     0  1   2L  3L2   |(a |)
         2          1     2     3                              3

Solving for ai  gives

pict

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

pict

or in matrix form

pict

The above can be written as

pict

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

(  N  (x) )    ( 1-(L3 - 3Lx2 + 2x3))
|   1    |    |-L13  2       2    3|
||  N2 (x) || =  ||L2 (1L x - 2Lx  + x )||
(  N3 (x))    (   L3 (3Lx2 - 2x3) )
   N4 (x)         L12 (- Lx2 + x3)

Notice that

N1 (0) = 1

and

N1 (L) = 0

as expected. Also

dN2-(x)-|||     -1-( 2           2)|||
  dx   |x=0 = L2  L  - 4Lx + 3x  |x=0 = 1

and

       |                         |
dN2-(x)||    = -1-(L2 - 4Lx + 3x2)||   = 0
   dx  |x=L   L2                 |x=L

and

N3 (0) = 0

and

N3 (L) = 1

and

       |                     |
dN4-(x)||      1--(         2)||
  dx   |   =  L2  - 2Lx + 3x |   = 0
        x=0                   x=0

and

       |                     |
dN4-(x)||   =  1--(- 2Lx + 3x2)||   = 1
  dx   |x=L    L2             |x=L

The shape functions have thus been verified.

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

               (      d3v(x)||    )
      ( V  )      EI  -dx3-||x=0
      |||{  1 |||}   ||  - EI d2v||    ||
{p} =   M1   = ||       d3x2 x|=0 ||
      ||| V2 |||   || - EI ddvx(x3)||   ||
      ( M2 )   (      d2v||  x=L )
                   EI dx2|x=L

Hence

       (   )    (      3                            )
       || V1||       EIddx3 (N1v1 + N2 θ1 + N3v2 + N4θ2)
       |{M1 |}    || - EI-d22 (N1v1 + N2θ1 + N3v2 + N4θ2)||
{p } = |   | =  ||     dxd3-                           ||
       ||( V2||)    ( - EI d2x3 (N1v1 + N2θ1 + N3v2 + N4θ2))
        M2        EI ddx2 (N1v1 + N2θ1 + N3v2 + N4θ2)
(6)

But

pict

and

pict

and

pict

and

pict

and now do the second derivatives

pict

and

pict

and

pict

and

pict

Hence eq (6) becomes

pict

or in matrix form, after evaluating the expressions above for x = L  and x = 0

pict

The above now is in the form

{p} = [K ]{d}

Hence the stiffness matrix is

         (                       )
            12   6L    - 12  6L
      EI-||  6L   4L2   - 6L  2L2 ||
[K] = L3 |( - 12  - 6L   12   - 6L |)
                   2           2
            6L   2L    - 6L  4L

Knowing the stiffness matrix means knowing the nodal displacements {d} given the forces at the nodes. The power of the finite element method now comes after all the nodal displacements v ,θ ,
 1  1  v2,θ2   are calculated by solving {p} = [K ]{d} because the polynomial v (x)  is now completely determined and hence v(x)  and θ (x)  can now be evaluated for any x  along the beam and not just at its end nodes. Eq. 5A above can now be used to find the displacement v(x)  and θ(x)  everywhere

pict

Therefore, these are the steps to obtain v(x)

  1. Find an expression for [K ]
  2. Solve {p} = [K]{d} for {d }
  3. Calculate v(x) = [N ]{d} by assuming v (x)  is a polynomial. This gives the displacement v(x)  to use to evaluate the transverse displacement anywhere on the beam and not just at the end nodes.
  4. Obtain θ (x) = dv(x)= -d [N ]{d}
         dx    dx to evaluate the rotation of the beam any where and not just at the end nodes.
  5. Obtain strain ϵ(x) = - y [B ]{d} where [B ]  is the gradient matrix        d2
[B ] = dx2 [N ]
  6. Obtain stress from σ = Eϵ = - Ey [B ]{d}
  7. Obtain the bending moment diagram from M  (x) = EI [B ]{d}
  8. Obtain shear force diagram from        -d
V (x ) = dxM (x)

2.1 Examples using the direct beam stiffness matrix

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

2.1.1 Example 1

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

PIC

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

PIC

P  is the given force. M2  = 0  since there is no external moment at the right end. Hence {p} = [K ]{d} for this system is

(    )       (                       ) (  )
|||  R |||          12   6L    - 12  6L    ||| v1|||
{ M1 }    EI-||  6L   4L2   - 6L  2L2 || { θ1}
|| - P||  = L3 |( - 12  - 6L   12  - 6L |) || v2||
|(  0 |)          6L   2L2   - 6L  4L2   |( θ|)
                                          2

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

{   }       (           ) {  }
 - P   = EI-   12   - 6L2   v2
  0      L3   - 6L  4L      θ2

Now the above is solved for {v2 }
  θ
   2 . Let            6
E = 30 × 10 psi  (steel), I = 57    4
in   , L = 144 in  , and P = 400lb  , 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

Therefore the vertical displacement at the right end is v  = 0.2324
 2  inches (downwards) and θ = - 0.0024
 2  radians. Now that all nodal displacements are found, the field displacement function is completely determined.

       (|v )|    (|    0   )|
       ||{ 1||}    ||{        ||}
{d } =  θ1   =      0
       |||(v2|||)    |||( - 0.2324|||)
        θ2       - 0.0024

Hence from eq. 5A

pict

But L = 144  inches, hence

v(x) = 3.9920× 10-8x3 - 1.6956 × 10-5x2

To verify, let x = 144  in the above

pict

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

PIC

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

(    )       (                       ) (  )
||  R ||          12   6L    - 12  6L    || v1||
|{ M1 |}    EI ||  6L   4L2   - 6L  2L2 || |{ θ1|}
|    |  = --3|(                       |) |  |
||( - P||)    L    - 12  - 6L2   12  - 6L2   ||( v2||)
   0            6L   2L    - 6L  4L      θ2

since v1 = 0  and θ1 = 0  , then

(    )       (                 ) (   )
|||  0 |||         1  0   0     0    ||| 0 |||
{  0 }    EI-|| 0  1   0     0  || { 0 }
|| - P|| =  L3 |( 0  0   12   - 6L |) || v2||
|(  0 |)         0  0  - 6L  4L2   |( θ2|)

And now the above system is solved as before.            6
E  = 30× 10 psi  , I = 57    4
in   , L = 144 in  , P =  400lb

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  

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.

2.1.2 Example 2

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

PIC

In using stiffness method, all loads must be on the nodes. The vector {p} 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

PIC

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

PIC

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

PIC

Now the stiffness equation is written as

pict

There is no need to determine R  and M1   at this point since these rows will be removed due to boundary conditions v  = 0
 1  and θ  = 0
 1  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 {p} = [K ]{d} . Hence, after removing the first two rows and columns, the system becomes

{     }       (           ) {  }
 - P∕2   = EI-   12   - 6L2  v2
 P L∕8     L3   - 6L  4L      θ2

Now the above is solved for {  }
 v2
 θ2 using the same numerical values for P,E, I,L  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

Therefore

      (   )   (                    )
      ||| v1|||   |||          0         |||
      { θ1}   {          0         }
{d} = || v2|| = || - 0.072630472854641 ||
      |(   |)   |(                    |)
        θ2      - 0.000605253940455

This is enough to obtain v (x )  as before. Now the reactions R  and M1   can be determined if needed. Going back to the full {p} = [K ]{d} , results in

pict

Hence the first equation becomes

              (                       )
R - P∕2 = EI-  0.87157 - 3.6315 × 10-3L
           L3

and since             6
E  = 30× 10 psi  (steel) and I = 57    4
in   and L = 144 in  and P = 400lb  , then

R = (30×1063)57(0.87157 - 3.6315× 10-3 (144 )) + 400∕2
       144  hence

R =  400 lb

and M1 - P L∕8 = EI3 (0.43578L - 1.2105 × 10-3L2)+ P L∕8,
             L  hence

M1 = 28762 lbft

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

pict

Since L =  144  inches , then

pict

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

PIC

2.1.3 Example 3

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

PIC

Using the following values: a = 0.625L, b = 0.375L,E = 30× 106psi  (steel), I = 57  in4   , L = 144 in, P = 1000lb  ,m = 200  lb/in.

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

PIC

And equivalent load for the uniform distributed loading is

PIC

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

PIC

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

pict

Now boundary conditions are applied. v1 = 0,θ1 = 0  ,v2 = 0  hence the first, second and third rows/columns are removed giving

   2      2
Pa-b-+ mL---=  EI4L2 θ
 L2     12     L3     2

Hence

    (    2       2)(     )
θ2 =  P-a-b+ mL---   -L--
       L2     12     4EI

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

pict

Hence the field displacement u (x)  is now found

pict

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

PIC

2.2 Adding more elements

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 v(x)  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)

2.2.1 Example 3 redone with 2 elements

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

PIC

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.

PIC

The stiffness matrix for each element is found then the global stiffness matrix is constructed, then {pglobal} = [Kglobal]{dglobal} 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 L  becomes L∕2  . The following diagram show the equivalent loading for P

PIC

The equivalent loading for distributed load m  is

PIC

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

PIC

Now {p} = [K]{d} for each element is constructed. Starting with the first element

(            m(L)      )
|      R1 -  -22--     |        (         (L-)             (L) ) (   )
||           m-(L2)2     ||        |  1(2)   6( 2)2    - 1(2 ) 6( 2)2 | ||| v1|||
||      M21L-   12   L   ||   -EI--|| 6 L2-   4 L2(- ) - 6 L2   2 L2(- )|| { θ1}
||  - Pb-(2+23a)-  m(22)- || = (L-)3 |( - 12  - 6 L2-    12    - 6  L2-|) || v2||
|( m (L)(2L∕2)m (L)2     2 |)    2     6(L)   2(L)2  - 6(L-)  4(L)2   |( θ2|)
  --122--- --122---  PaLb2             2      2         2     2
                   (2)

And for the second element

(    Pb2(L2+2a)   m(L2)  )
|  - --(L∕2)3---  -2--- |        (          ( )             ( ) ) (   )
|| m-(L2)2  m-(L2)2   Pab2-||           12    6  L2-    - 12    6 L2-   || v2||
||   12  -   12  -  (L2)2||    EI  || 6(L-)  4(L-)2   - 6 (L) 2 (L)2|| |{ θ2|}
||      Pa2(L2+2b)-  m-(L2)|| =  (L)3|| - 122  - 62(L-)    122   - 6 2(L)|| | v |
|| R3 -   (L2)3   -   2  ||     2  (  (L )   (L 2)2      (L)   (L 2)2) ||(  3||)
|(     m (L2)2   Pa2b     |)          6 -2   2 2-    - 6 2-  4 -2      θ3
      --12--+ (L-)2
                2

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

(            m(L2)-     )
|       R1 -   2L 2     |        (          ( )                        ( )                     )
||      M1 -  m(2)--    ||           12    6  L2-        - 12           6 L2-         0       0     (| v1)|
||    P b2(L2+2a)12  m(L2) ||        || 6 (L)  4(L-)2     - 6(L-)         2(L-)2        0       0   || ||||   ||||
||  - 2--(L)3-- - 2--2--||        ||    2      2(L)          2         (L-)2  (L-)            (L) || |||{ θ1|||}
||        2  Pab2       || =  E(I)-||  - 1(2) - 6( 2)2     1(2)+ 12( )  -(6 2)2 + 6(  2)2   - 12( ) 6( 2)2||   v2
||        - 2(L2)2       ||     L2-3|| 6  L2-  2  L2-   - 6 L2- + 6 L2-  4  L2- +(4) L2-   - 6 L2-  2  L2(-)|| ||| θ2|||
||      Pa2(L2+2b)   m(L2)||        |(   0      0          - 12          - 6 L2-        12    - 6 L2-|) |||| v3||||
|| R3 -   (L2)3   -   2  ||            0      0         6 (L)          2(L-)2      - 6 (L) 4(L-)2  |( θ3|)
|(      m(L)2     2     |)                                2              2            2     2
       -122--+  P(aL)b2
                2

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

(          0          )        (                             )
|                     |          1  0    0      0    0    0     ( 0)
||     2 L  0          ||        || 0  1    0      0    0    0  ||  | 0|
||- 2Pb-(2+23a)- 2 m(L2∕2)||        ||                          (L)||  ||  ||
||     (L∕2) Pab2-      || =  E(I)-|| 0  0   24     (0)2  0  6( 2)2||  ||v2||
||       - 2 (L∕2)2      ||     L2-3|| 0  0    0   8  L2-   0  2 L2- ||  || θ2||
||          0          ||        |( 0  0    0      0    1    0  |)  |( 0|)
|(     m(L2)2-  Pa2b    |)          0  0  6(L-) 2 (L)2  0  4(L)2     θ3
       12  +  (L2)2                       2      2         2

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

(                     )
  - 2Pb2(L2+2a)- 2m-(L2)         (                 ( ) ) (   )
||      (L∕2)3  2    2  ||           24      0    6  L2-     v2
||        - 2-PLab2      || =  E(I)-|(   0   8 (L)2  2(L-)2|) |( θ2|)
|(       L 2( 2) 2     |)     L2-3   (L-)   (2L)2   (L2)2
      m(122)-+ PLab2               6  2  2  2    4 2       θ3
              (2)

giving the same solution.

Hence, there are 3 unknowns to solve for. Once there are solved, v (x )  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

The above solution gives v2 = - 0.2735  inch (downwards displacement) and θ2 = - 0.0019  radians and θ3 = 0.0076  radians. Now the v (x )  polynomial is found for each element

pict

The above polynomial is the transverse deflection of the beam for the region 0 ≤ x  ≤ L ∕2  .

v(x)  for the second element is found in similar way

pict

Hence

           0.0304(  3  1-  2)   0.0076 (1- 2      2    3)    2.188-( 1-3   3-  2    3)
velem2(x) =   L2    x - 2 Lx   -   L2    4 L x- Lx  + x   -   L3    8L -  2Lx  + 2x

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

PIC

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.

PIC

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 w  lbs per unit length is given by v(x) = - wx2-(3L2 - 5Lx + 2x2)
         48EI  , while the analytical deflection for the same beam but when there is a point load P at distance a  from the left end is given by                   3         2    2                2           3
v(x) = - P(⟨L - a⟩ (3L - x)x + L (3(L - a)(L-  x)x + 2L ⟨x- a⟩ ))  , therefore, the analytical expression for deflection is given by the sum of the above expressions, giving

            2
v (x) = --wx--(3L2 - 5Lx + 2x2) - P(⟨L - a⟩3(3L - x)x2 + L2 (3(L - a)(L - x)x2 + 2L ⟨x- a⟩3))
         48EI

Where ⟨x - a⟩ means it is zero when x - a  is negative. In other words ⟨x - a⟩ = (x - a)U nitStep(x-  a)

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

PIC

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.

3 Generating shear and bending moments diagrams

After solving the problem using finite elements and obtaining the field displacement function v (x)  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              d2v(x)-
M  (x) = - EI  dx2   and shear force is given by                    3
V (x) = dMdx-=  - EI ddv(xx3)   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

                            (          )       (          )                }
             (         )2.1L838 2x3 -( 32Lx2 -  0.0L0726-x)3 - 12Lx2(                )    0 ≤ x ≤ L∕2
v(x) =  0.03L024 x3 - 12Lx2  - 0.0L0276- 14L2x - Lx2 + x3 - 2.L1838  18L3 - 32Lx2 + 2x3     L∕2 ≤ x ≤ L

Hence

                     -0.0076(6x-L)   2.188(12x-3L)-       }
M  (x) = - EI   0.0076(6x- 2LL)2  0.03+04(6x-L)L32.188(12x-3L)    0 ≤ x ≤ L ∕2
              - ----L2-----+ ----L2----- ----L3-----   L ∕2 ≤ x ≤ L

using E = 30×  106psi  and I = 57  inch4   and L = 144  the bending moment diagram plot is

PIC

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)

PIC

The reason for this is because the solution v(x)  obtained using the finite elements method is a 3rd  degree polynomial and after differentiating twice to obtain the bending moment (             2
M  (x ) = - EI ddxv2   ) the result is a linear function in x  while in the analytical solution case, when the load is distributed, then solution v (x)  is a 4th  degree polynomial. Hence the bending moment will be quadratic function in x  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.

4 Finding the stiffness matrix using methods other than direct method

There are 3 main methods to obtain the stiffness matrix

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

4.1 Virtual work method for derivation of the stiffness matrix

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 u = {u,v,w } and assuming the external loads are given by {p} acting on the nodes, hence these point loads will do work given by {δd}T {p} on that unit volume where {d} 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 {δε}T {σ} in that same unit volume.

PIC

Hence, for a unit volume

    T          T
{δd}  {p} = {δε} {σ}

And for the whole volume

            ∫
{δd}T {p} =   {δε}T {σ} dV
             V
(1)

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

u = [N ]{d}

Therefore

pict

Since {ε} = ∂{u } then {ε} = ∂ [N ]{d} =  [B]{d} where B  is the strain displacement matrix [B ] = ∂ [N ]  , hence

pict

Now from the stress-strain relation {σ} = [E]{ε} , hence

{σ} = [E][B ]{d}
(4)

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

    T      ∫      T   T
{δd}  {p}-    {δd}  [B ] [E][B]{d} dV = 0
            V

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

      (          ∫               )
{δd}T  {p} - {d}   [B ]T [E ][B ]dV   = 0
                  V

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

          ∫
{p} = {d}   [B ]T [E ][B]dV
           V

Hence this is in the form P = K Δ  , therefore

      ∫
            T
[K ] = V [B ] [E][B]dV

knowing [B ]  allows finding [k]  by integrating over the volume. For the beam element though, u =v (x)  the transverse displacement. This means       d2-
[B] = dx2 [N]  . Recall from the above that for the beam element,

      (-1   3      2     3  -1   2       2    3   -1     2     3   -1      2   3 )
[N] =  L3 (L - 3Lx  + 2x )  L2 (L x-  2Lx  + x )  L3 (3Lx  - 2x  )  L2 (- Lx  + x )

Hence

        2       (                                                            )
[B] = -d--[N ] =  -13 (- 6L + 12x)  -12 (- 4L + 6x)  13 (6L - 12x)  12 (- 2L + 6x)
      dx2        L               L               L              L

Hence

pict

Which is the stiffness matrix found earlier.

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

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    U  and let the work done by external loads acting on the nodes be Ω  , then

Π = U  - Ω

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 1 ∫ {σ}T {ε}dV
2 V  and the work done by external loads is {d}T {p} , hence

      ∫
Π = 1-   {σ}T {ε}dV - {d}T {p}
    2  V
(1)

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

u = [N ]{d}

Since {ε} = ∂{u } then {ε} = ∂ [N ]{d} =  [B]{d} where B  is the strain displacement matrix [B ] = ∂ [N ]  , hence

{ε} = [B ]{d}

Now from the stress-strain relation {σ} = [E]{ε} , hence

    T      T   T
{σ } =  {d} [B ] [E ]
(4)

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

    1 ∫     T   T                   T
Π = 2-   {d} [B]  [E ][B ]{d}dV  - {d} {p }
       V

Hence -∂Π-
∂{d} = 0  gives

       ∫
0 = {d}   [B ]T [E][B]dV -  {p}
        V

Which is on the form P =  [K ]D  which means that

      ∫
[K ] =   [B ]T [E][B]dV
       V

As was found by the virtual work method.

5 References

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