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
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
Assuming the length of the beam is , then
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
as expected. Also
The shape functions have thus been verified.
The stiffness matrix is now found by substituting eq (5A) into eq. (1), repeated below
and now do the second derivatives
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
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 removed1 . 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
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
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. , , ,
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
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
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
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
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
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
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
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
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
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
Assuming that displacement can be written as a function of the nodal displacements of the element results in
Since then where is the strain displacement matrix , hence
Now from the stress-strain relation , hence
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,
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
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
Substituting Eqs. (2,3,4) into (1) results in
Which is on the form which means that
As was found by the virtual work method.