Building a Beam Deflection Calculator in Python

1.0 Beam Deflection Calculator – Project Overview

In this project, we’ll build a beam deflection calculator that can generate beam deflections by directly integrating the bending moment diagram. Since our approach will require a complete bending moment diagram, we’ll add this code onto the end of our Shear Force and Bending Moment Diagram Calculator. So at the end of this project, the final result will be a complete beam analysis code that calculates beam reactions, shear forces, bending moments and deflections.

In our previous beam analysis project, we were limited to analysing statically determinate beams since we relied on the equations of statics. The technique we’ll use for calculating deflection in this project is not limited to statically determinate structures, although you will need a complete bending moment diagram to integrate.

This Python mini-project is accompanied by a complete HD video series. You can access the full playlist directly below or by enrolling in the project on DegreeTutors.com. In the video series, we’ll walk through the build process line by line. So you can build along with me or if you prefer, just download the finished Python code. I recommend reading through this post to get a good overview of the build process and then working through the tutorial videos in more detail to build your own calculator (or just download the finished code by enrolling in the project).

If you’re new to Python, let me say loud and clear at the start – you don’t need to be a Python programmer to complete this project! Complete beginners can work through this project and pick up the programming basics along the way. If you need help getting your coding environment setup, check out this video (part of another course but covers what you need to get up and running).

2.0 Our Approach to Building a Beam Deflection Calculator

We’ve seen in previous tutorials how we can calculate beam deflections by analytically integrating the differential equation of the deflection curve. We saw first how we could do this by establishing expressions for the internal bending moment as a function of . In a later tutorial,  we implemented Macauley’s Method to speed up the process by reducing the amount of analytical integration required. 

In this project, we’re going to take a different approach. We will still integrate the differential equation of the deflection curve, but we’ll leave the analytical integration to one side and perform our integration numerically. This numerical approach suits our needs quite nicely as we want a method that requires minimal ‘human intervention’. There will be no need to manually determine an equation for the bending moment as a function of because we already have the bending moment magnitude at each position along the beam. As mentioned above, another strength of this approach is that because all we require as an input is the bending moment diagram, our deflection calculation is not limited to statically determinate structures. 

In addition to integrating the differential equation of the deflection curve numerically, we will need to implement some logic along the way that uses our structures’ boundary conditions to ensure we end up with a valid deflected shape, but more on that later.

3.0 Numerical Integration of the Bending Moment Diagram

In this section, we’ll discuss the basic idea of numerical integration and then use this knowledge to build a Python function that integrates the bending moment diagram. Our Python function will return the rotation and displacement at each position along the beam.

3.1 The basic idea – trapezoidal rule

Before we think about numerical integration applied to our specific problem, let’s revisit the basics. We could use several schemes, but one that most people have come across at some point is the trapezoidal scheme or trapezoidal rule. Remember that in the simplest of terms,  integration is a way of determining the area under a line. So, let’s imagine some function of sampled at points , yielding values of the function, , Fig. 1 below. 

The trapezoidal rule says that the area beneath the line made up by joining the sample values of the function is approximated as,

(1)  

where , the interval between sampling points (assumed to remain constant). So, the trapezoidal rule boils down to,

(2)  

And if we only have two sample points, i.e. a ‘first’ and ‘last’, we have,

(3)  

which we can quite clearly see is the area of a rectangle, Fig 2. This is going to be the basic scheme we use to perform our numerical integration.

Now we can return to our specific challenge; let’s rearrange the differential equation of the deflection curve,

(4)  

Now consider plotting values of , which is just a scaled version of our bending moment diagram. 

Before we can obtain deflections, we must perform an integration to get the rotations,

(5)  

We do this by ‘stepping’ through the diagram and applying the trapezoidal rule to each pair of values. So, to obtain for example,

(6)  

After we generalise this equation for any two values of we obtain,

(7)  

Notice that to start this whole process off, we needed to provide some value for . We essentially have to guess this value! If the value is incorrect, then everything that follows, rotations and deflections will be incorrect too. So we’ll discuss how to address this below. For now, let’s assume we’re happy with the initial value of rotation and proceed to calculate deflections. We need to apply the integration process to our rotation values,

(8)  

Notice here we also need to provide a starting value for deflection, . If we start our integration at a support, we can safely assume that this value is zero. So generalising our equation we have,

(9)  

After processing all values of , we obtain deflection values for each point along our beam.

3.2 Building a numerical integration function

Now that we understand the basic logic involved in the numerical integration and how to apply it in principle to our problem, we can capture this with some Python code (remember we’re just getting an overview here – refer to the video tutorials above for line by line explanation). For now, we’ll just focus on building a function that takes as its inputs,

• the values of bending moment 
• the flexural rigidity,
• the spacing between sample points,
• the initial rotation and deflection, and

Our function will process this information and return a rotation and deflection for each sample point along the beam. The reason we’re capturing this logic in a function is so that we can call the function multiple times as we search for the correct initial rotation (discussed below). The complete function is shown in the code block below. Note that some setup code is required before this function, so refer to the video tutorial for the full setup details.

In lines 4 and 5, we’re simply renaming some of the input arguments; the ‘im1’ in the name stands for . Next, lines 8 to 11 define some data containers to hold our rotation and deflection values. We also assign our initial and values at the location in the array that corresponds to the location of the beam support (again, I dig into this more in the video tutorial).

In lines 14 to 30, we iterate through each value from the array of bending moment values (calculated previously) and perform the numerical integration discussed above. Finally, we return the calculated arrays of rotation and deflection values from the function on line 32. Now that we have a function to perform the numerical integration, we need to write the logic that will call this function.

4.0 Identifying the correct initial rotation

So far, we’ve talked about our deflection calculation in an abstract way; let’s base our remaining discussion on a specific beam. Consider the beam below with a cantilever overhang on each side. For now, we’ll ignore the overhang to the left of support A. 

As we said previously, the biggest weakness in our solution so far is that we need to guess a value for the initial rotation, i.e. the rotation at support A. If the wrong value is selected, then all that follows, including our beam deflections, will be incorrect. In this section, we’ll discuss how to address this and determine the correct initial rotation. Our approach will be iterative. We will initially guess a random value for and progressively refine our estimate until will arrive at an acceptable value.

4.1 The basic idea – minimising the error

After we guess the initial rotation, we will call our function and calculate the beam deflection. Once we have a deflection, we need to judge its accuracy and, therefore the accuracy of our guess. Fortunately, one thing we know for sure is that the deflection is zero at support B. So this becomes our test; the vertical deflection at the right support becomes our error. All we need to do is keep changing our initial rotation value, calculating the beam deflection and checking if the error is zero or close enough that we’re happy with the approximation.  

Next, we need to think about how to alter our guess based on the value of deflection (error) observed at B. If the deflection at B is above the line, we should increase the rotation at A, which will rotate the whole deflected shape such that the deflection at B is reduced back down towards zero. If the deflection at B is below the beam, we should do the reverse and reduce the initial rotation to rotate the deflected shape back up towards zero. 

In terms of coding, we’re going to tackle this in two steps. First, we’ll write some code that checks whether increasing or reducing our initial guess reduces the error. Once we know the direction to sweep through rotation values, we’ll call another function that determines the actual value of that leads to the smallest error at B. 

The code below makes reference to two important parameters that we can tweak to refine our deflection estimates; initRot and deltaRot. initRot is the initial rotation value; the closer it is to the correct value, the fewer iterations are required to calculate the correct initial rotation. deltaRot is the amount by which we change the initial rotation estimate on each iteration of our test. The smaller deltaRot is, the closer we can get to the correct initial rotation value, but the more iterations will be required to get there. You should play around with these values in your final code.  

In lines 3-6 below, we’re writing a for loop that will calculate and store the deflection at support B for three different initial rotation values; 

• our initial guess, 
• slightly below the initial guess,
• slightly above the initial guess. 

This allows us to test whether reducing or increasing the initial rotation leads to a smaller or larger deflection error at support B, which we’re testing on lines 8 to 16. 

Notice that on lines 11 and 16, we’re calling a function, zeroCrossing() that we haven’t written yet. This function is called with a different value of deltaRot, depending on whether or not we need to continue increasing or reducing our initial rotation. This function returns a single value which we assign to a variable called solvedInitRotation. This is our final estimate of the correct initial rotation at support A. 

This allows us to call the deflection calculation function one last time on line 19 and pass in the correct value of the initial rotation. Therefore the deflections and rotations returned on line 19 are the final ‘correct’ values or, in other words, our best approximation. 

4.2 Building a function to find the initial rotation

Now we need to write the function zeroCrossing(), which will take as its inputs the deflection values, the deltaRot value (positive or negative as determined by the code above), and the initial rotation and deflection values at support A. The complete function is shown below.

This function’s basic approach is to calculate the deflection at support B and test whether it has changed sign since the previous calculation. The while loop will run, successively altering the initial rotation value until the error changes sign. This indicates that between the previous and current loop iterations, the error passed through zero. When this happens, the value of initRot is saved and returned as the ‘solved for’ initial rotation. The smaller the step size, deltaRot, the closer we’re likely to get to the theoretically correct initial rotation. Zero-crossing is a pretty common and simple algorithm we can use to identify when some quantity is zero, or very close to zero.

5.0 Reversing integration direction for deflection of the left side cantilever

At this point, we have produced a reasonable approximation of the beam’s deflection to the right of support A, including the cantilever region to the right of support B. However, we also need our calculator to evaluate cantilever defections to the left of support A. Fortunately, we can use the initial rotation we’ve already calculated for support A and simply integrate in the reverse direction from A out to the left end of the beam.

5.1 Writing the code to identify the left overhang deflection

The code block below shows the reverse integration. This code is almost identical to the code we wrote in our calcDeflection() function above. In particular, note the moment value assignment on lines 13 and 14. Because we are integrating to the left, as we iterate through our moment values, we’re access values with smaller and smaller indices in the array. It’s also worth noting that we only need to access this code if there is a cantilever overhang to the left of the leftmost support. This condition is captured in line 2 by checking if dimension A (specified previously) is not equal to zero.

6.0 Plotting the beam deflection

All that remains is to plot the deflections. We’ve discussed the plotting code in our previous shear and moment diagram project, and I review it in the video tutorial above, so for brevity, I won’t discuss it here. We’ll simply plot the deflection below. To obtain these deflections, I’ve used a value of for Young’s Modulus and assumed a rectangular beam cross-section wide and deep.

7.0 Comparing our numerical deflection with an analytical solution

The beam we introduced above was the subject of our Macauley’s Method analysis in a previous post. This presents us with a convenient validation opportunity. It’s always a good idea to perform some kind of validation exercise whenever you develop any new code. So that’s what we’ll do here by comparing the deflections calculated using Macauley’s Method with those obtained from our numerical deflection estimate. We’ll actually take the opportunity to perform two validation exercises, one on each of the beams analysed in our Macauley’s Method post.

7.1 Example 1 – Simply supported beam without cantilevers

Let’s start with the more straightforward beam analysis.

The code block below shows the code that we wrote in the original Macauley’s Method post to plot the deflected shape using the derived deflection equation, 

(10)  

If we input the beam and loading details for this structure into our beam analysis notebook, we can plot both deflections on the same set of axes, Fig. 7. We can see that the numerical approximation of deflection obtained using our beam deflection calculator is in very close agreement with the analytical solution obtained using Macauley’s Method; so close in fact that when we plot the numerical deflection line (yellow) over the line for the analytical solution (green), it completely obscures it.

Considering the effort required to implement Macauley’s method, the advantage offered by our new calculator is obvious. With minimal effort, simply inputting the problem parameters, we’ve produced an excellent approximation of the deflection. This gives us great confidence in the code we’ve written.

7.2 Example 2 – Simply supported beam with left and right cantilevers

We’ll compare the numerical and analytical deflections for the beam we started with above (Fig. 4) to conclude our discussion. Again, we’ve derived the Macauley’s Method deflection equation previously,

(11)  

For completeness, the code used to generate the deflected shape using this equation is shown below.

Comparing the numerical and analytically obtained deflections, Fig. 8, we can again see close agreement, strengthening our confidence in the validity of our calculator output.

At this stage, you should have a good sense of how we’ve pieced this beam deflection calculator together. Armed with the roadmap in this post, go ahead and work your way through the tutorial videos at the top of the post. These will cover all of the details that we didn’t have space to go into here. 

I hope this post has given you a taste for combining basic coding with structural analysis. Programming is a huge force multiplier for your structural analysis knowledge. If you want to dive further into code-based structural analysis, consider the stiffness method course bundle below. In these courses, I walk you through building a range of Python analysis tools for analysing 2D and 3D truss systems as well as 2D beam and frame structures. What we discussed in this post is only the tip of the iceberg compared to what you can build!

That’s all for now, see you in the next one!