fbpx

Building a Mohr’s Circle Calculator for Stress Analysis in Python

by Dr Seán Carroll
|
Updated: 9th July 2021
|
Python Project

1. How we’re going build a Mohr’s Circle Calculator using Python

In this Python project we’re going to build a Mohr’s Circle calculator. By the end of this project, you will have built your own stress analysis Python code. You’ll also have enough knowledge to expand and personalise your code to your own needs. Along the way we’ll cover all of the fundamental topics that lead up to Mohr’s circle of stress. You will learn about:

  • how we use the 2D stress element to represent the state of stress at a point
  • the purpose of stress transformation equations
  • principal stresses and principal planes
  • planes of maximum shear stress
  • and of course Mohr’s circle!

You don’t need to be a Python programmer to work through this project! We’ll be using the Jupyter Notebook development environment to write our Python code. If you’re new to Jupyter Notebooks (or Python) and need a little help getting up and running – check out this getting started video. If you’re new to Python, don’t let it intimidate you! Once you get comfortable with the basics – it will unlock so many doors in your day-to-day work as a student or engineer. 

You can download the complete Jupyter Notebook file for free by enrolling in the associated DegreeTutors course (free)  – it will be handy to have as a reference as you write your own code. 

Mohr’s circle is an elegant graphical method for representing the state of stress at a single point. By constructing a Mohr’s circle – sketching it out and using basic trigonometry – or in our case capturing the trigonometry in code, you can graphically identify

  • the principal stresses,
  • orientation of the principal planes,
  • maximum shear stresses,
  • the orientation of the maximum shear stress planes
  • and the stresses on any other pair of orthogonal planes. 

We’re going to work our way incrementally towards our goal of building a Mohr’s circle calculator. We’ll be building our stress analysis/Mohr’s circle Python code over the course of this project and incrementally adding to our code. This post is accompanied by a full series of video lectures that you can watch in the playlist below. I suggest reading through this post to get a good overview of the theory, then diving into the video series where I’ll walk you through the theory and build process, step by step. So without any further delay, let’s dive in! 

2. Stress at a point and the 2D Stress Element

At any point in a structure under load the state of stress can be represented by a combination of normal stresses (a.k.a. bending stresses) and shear stresses. If we assume a state of plane stress (we discuss plane stress in more detail in video 2 of the playlist), this combination of stresses can be conveniently represented on a 2D stress element. A stress element is simply an infinitesimally small element of area on which the stresses at that point can be shown acting upon. Don’t be confused – all stresses shown acting on the 2D element are acting at a single point in space – just in different directions. 

To ground this discussion in something a little more approachable, let’s consider a simply supported beam subject to a uniformly distributed load. The external load will give rise to a shear and normal stress field within the structure. One simple representation of this varying stress field is the shear force and bending moment diagram shown below.

Stress Analysis-Mohr's Circle 1 | DegreeTutors.com
Fig 1. Simply supported beam subject to a uniformly distributed load and the resulting bending moment and shear force diagram.

If we imagine making a vertical cut through our beam at some distance x_A along the length of the beam we can visualise the (linear) normal and (parabolic) shear stress distribution through the depth of the beam at that point, Fig 2. The normal stress distribution, \sigma_x is obtained using the engineer’s bending equation,

(1)   \begin{equation*}\sigma_x = \frac{-My}{I}\end{equation*}

where M is the bending moment at the location of the cut, y is the height from the neutral axis of the beam to the point at which the normal stress is being determined and I is the second moment of area of the beam cross-section. The negative sign ensures a negative normal stress is obtained for compression stresses above the neutral axis. 

The shear stress distribution \tau is given by the following equation,

(2)   \begin{equation*} \tau = \frac{VQ}{Ib} \end{equation*}

where V is the shear force at the location of the cut, Q is the first moment of area of the area above the level at which the shear stress is being determined, I is the second moment of area for the cross-section and b is the width of the cross-section. 

Stress Analysis-Mohr's Circle 2 | DegreeTutors.com
Fig 2. Vertical section through the beam at position x=x_A and the normal and shear stress distribution down through the depth of the section.

If we now imagine evaluating the normal and shear stress at a height y_A above the neutral axis, we would obtain normal and shear stresses \sigma_A and \tau_A respectively as shown below. This is where our 2D stress element comes in useful. We can represent the state of stress at this point in the structure on a 2D stress element. To be clear, we are showing the state of stress at location or coordinates (x_A, y_A) on the 2D stress element below.

Stress Analysis-Mohr's Circle 3 | DegreeTutors.com
Fig 3. Stresses at horizontal position x_A and at vertical distance y_A from the neutral axis, represented on a 2D stress element.

We talk more about sign conventions for normal and shear stress in video 2, for now it’s sufficient to recognise that the 2D stress element is acting as a ‘stage‘ on which we can ‘display‘ the stresses that exist at this location. The normal stress is shown as a compression stress (squeezing the element) which we typically define as a negative quantity. The shear stress as shown is negative when considered against the sign convention defined in the video – you’ll definitely want to review this video, particularly for when we discuss maximum shear stresses below.

🚨 The big question now is what would happen to the stress magnitudes if we rotated our 2D element?

Stress Analysis-Mohr's Circle 4 | DegreeTutors.com
Fig 4. Normal and shear stresses on a 2D stress element (left), unknown stress magnitudes after 2D element is rotated through an angle \theta (right).

The stresses we have shown on our 2D element above are merely the set of stresses that exist on the planes depicted by the 2D stress element i.e. a set of orthogonal planes that happen to be aligned conveniently with the beam. This is not to say that these are necessarily the maximum values of normal and shear stress at that point. There could exist an alternative pair of mutually perpendicular planes, at the same location but at some other orientation, on which the normal or shear stresses will have a larger or smaller magnitude. As engineers, this is clearly something we’d like to know more about! This is where stress transformation equations enter the picture.

3. Stress Transformation Equations

If we know the state of stress at a point on one pair of mutually perpendicular planes, the question we posed in the previous section was, what are the stress magnitudes on other planes of different orientations at the same point?  It would make our lives a lot easier if we could determine the normal and shear stresses as a function of the angle of the planes. Fortunately this is exactly what the stress transformation equations allow us to do. 

The stress transformation equations can be derived by using simple statics and we cover the process in detail in videos 2 and 3. For now, we’ll simply state the equations and discuss what they’re telling us. We can start by imagining some 2D stress element subject to a normal stress \sigma_x in the horizontal direction, a normal stress \sigma_y in the vertical direction and and shear stress \tau_{xy} acting on the vertical faces and \tau_{yx} acting on the horizontal faces.

You’ll recall from the principle of complimentary shear stresses that the magnitudes of the shear stress acting on all faces are equal and that the direction of any one shear stress will define the directions of the three others (such that equilibrium of the element is maintained). We can visualise our stress element below. 

Stress Analysis-Mohr's Circle 5 | DegreeTutors.com
Fig 5. Stress element subject to orthogonal normal stresses and shear stresses.

Now imagine the element rotated counter-clockwise by some angle \theta. When we talk about rotating the element, what we’re really talking about is considering stresses on a set of rotated planes.  The x_1-y_1 axis now denotes the local axis of the rotated element/planes and this axis system is rotationally offset from the original x-y axis by \theta degrees. The normal and shear stresses on this rotated set of planes, \sigma_{x1}, \sigma_{y1} and \tau_{x1,y1} are shown below.

Stress Analysis-Mohr's Circle 6 | DegreeTutors.com
Fig 6. Stress element orientated at angle \theta with respect to the x and y axis, subject to orthogonal normal stresses and shear stresses.

We can determine the stresses on the rotated set of planes using stress transformation equations (refer to videos 2 and 3 for the derivation),

(3)   \begin{equation*}\boxed{\sigma _{x1} = \frac{\sigma_x + \sigma_y}{2} + \frac{\sigma_x - \sigma_y}{2} \cos 2\theta + \tau_{xy}\sin 2\theta}\end{equation*}

(4)   \begin{equation*}\boxed{\sigma _{y1} = \frac{\sigma_x + \sigma_y}{2} - \frac{\sigma_x - \sigma_y}{2} \cos 2\theta - \tau_{xy}\sin 2\theta}\end{equation*}

(5)   \begin{equation*}\boxed{\tau _{x1y1} = -\frac{\sigma_x-\sigma_y}{2}\sin 2\theta + \tau_{xy}\cos 2\theta}\end{equation*}

These three transformation equations allow us to determine the stresses (\sigma_{x1}, \sigma_{y1} and \tau_{x1,y1}) on any set of x_1-y_1 planes provided we know the stresses on a set of x-y planes (\sigma_x, \sigma_y and \tau_{xy}) and the angle \theta between the two sets of planes. 

At this point we need to make an important observation. If we add the first two transformation equations together, i.e. equations (3) and (4) above, we get the following,

(6)   \begin{equation*}\boxed{\sigma_{x1} + \sigma_{y1} = \sigma_x + \sigma_y}\end{equation*}

This result tells us that the sum of the orthogonal normal stresses at a point are constant and independent of the angle or orientation of the planes. This is a powerful observation and can be understood by a simple analogy.

Let’s imagine that the state of normal stress at a point could be represented by a single vector and the normal stresses we visualise on our 2D element are just two orthogonal components of that vector. By rotating the planes, we’re merely considering different orthogonal components of the same underlying vector – so the underlying state of stress at that point (the vector in our analogy) is not changing – only the orthogonal components as the reference axis changes. We can visualise this below.

Stress Analysis-Mohr's Circle 7 | DegreeTutors.com
Fig 7. Normal stress components \sigma_x and \sigma_y aligned along axis system x-y (left), normal stress components \sigma_{x1} and \sigma_{y1} aligned along axis system x_1-y_1 (right). Note the resultant vector \sigma is unchanged between the left and right diagrams.

To get a better insight into what the transformation equations are saying, we can plot them. For the purpose of the exercise let’s assume stresses with the following magnitudes,

Stress Analysis-Mohr's Circle 9 | DegreeTutors.com
Fig 8. Assumed stresses on the x-y faces corresponding to an orientation of zero degrees.
  • \sigma_x = 1\: N/mm^2
  • \sigma_y = 0.25\: N/mm^2
  • \tau_xy = 0.75\: N/mm^2

We walk through the Python code to evaluate the stresses and produce the following plot in video 3. We can clearly see that the stress magnitudes vary sinusoidally. The dashed vertical line at 0^\circ identifies the magnitudes of the three stresses indicated in Fig 8. above at zero degrees. The second dashed vertical line indicates the stresses that would be observed if the the element were rotated by 90^\circ .

Stress Analysis-Mohr's Circle 8 | DegreeTutors.com
Fig 9. Variation of normal and shear stress with angle

By observing the blue and green lines representing normal stresses \sigma_x and \sigma_y respectively, we can see that their magnitudes are perfectly out of phase with each other – in other words, when one is a maximum, the other is a minimum. On close inspection we can also see that at the angle at which the maximum and minimum normal stresses occur, the shear stress, shown in red, is zero. These are all important results that are discussed further in lectures 4 and 5. 

Perhaps the biggest takeaway so far from our plot of the transformation equations is that the stresses we started with, indicated in Fig. 8 are clearly not the maximum values of normal or shear stress that occur at this position in the structure. There exists another pair of planes on which the maximum and minimum normal stresses occur and there are yet another pair of planes on which the maximum shear stresses occur. This observation motivates our discussion of Principal Stresses in the next section.

4. Principal Stresses and Principal Planes

As the name suggests the principal stresses are the maximum and minimum values of the normal stress. These stresses occur on mutually perpendicular planes called principal planes. We can confirm this by observing the peak in the blue line at approximately 30^\circ in Fig 9. This peak denotes the maximum principal stress occurring on a principal plane orientated at 30^\circ to the horizontal x axis. If we continue to rotate our element/plane by a further 90 degrees to approximately 120^\circ, we see that the blue line is now at a minimum. This is the minimum principal stress occurring on a principal plane orientated at 90^\circ to the first principal plane. 

We can derive a convenient expression for the magnitudes of the principal stresses and the orientation of the planes on which they act, see video 4 for the derivation and further explanation. We can state the key equations here; the principal stresses are given by,

(7)   \begin{equation*}\boxed{\sigma_{1,2} = \frac{\sigma_x + \sigma_y}{2}\pm\sqrt{\left(\frac{\sigma_x-\sigma_y}{2}\right)^2 + \tau_{xy}^2}}\end{equation*}

where \sigma_1  is the algebraically larger principal stress and \sigma_2 is the smaller one. The principal angles \theta_p are determined from,

(8)   \begin{equation*}\boxed{\tan 2\theta_p = \frac{2\tau_{xy}}{\sigma_x-\sigma_y}}\end{equation*}

The angle \theta_p has two values, 90^\circ apart. One value identifies the plane on which \sigma_1 acts while the other identifies the plane on which \sigma_2 acts. Both angles can be substituted into our original transformation equation (Eqn 3.) to determine the principal stress for each angle. The values obtained will agree with the principal stresses obtained using equation 7.

Using equations 7 and 8 we can determine the principal stresses for our example to be \sigma_1 = 1.46\:N/mm^2 and \sigma_2 = -0.21\:N/mm^2. We also identify principal angles of 31.7^\circ and 121.7^\circ. We can plot vertical lines indicating the principal angles and horizontal lines indicating the principal stresses. If we overlay these on our original plot of the transformation equations, we can see that equations 7 and 8 do indeed identify the local maxima and minima, i.e. our principal stresses. From the graph below we can see that the maximum principal stress occurs on the plane oriented at 31.7^\circ while the minimum principal stress occurs on the plane orientated at 121.7^\circ

Equations 7 and 8 are valuable tools as they allow us to identify principal stresses and their orientations provided we know the stresses on any pair of planes. So a typical stress analysis task would be to identify the principal stresses and principal planes at a point by first identifying a more conveniently obtainable set of stresses – like for example the normal and shear stress in the beam at the top of the page. 

Stress Analysis-Mohr's Circle 10 | DegreeTutors.com
Fig 10. Variation of normal and shear stress with angle. Red vertical lines indicate principal angles identified from equation 8 while horizontal lines indicate principal stresses identified from equation 7.

5. Maximum shear stresses and planes of maximum shear stress

We can now turn our attention to the maximum shear stresses and the planes on which they occur. The reference video for this discussion is video 5. In that video, through a similar derivation process, we show that the maximum shear stress is given by,

(9)   \begin{equation*}\boxed{\tau_{max} = \sqrt{\left(\frac{\sigma_x-\sigma_y}{2}\right)^2 + \tau_{xy}}}\end{equation*}

and that these stresses occur on planes orientated at 45^\circ to the principal planes,

(10)   \begin{equation*}\boxed{\theta_s = \theta_p\pm 45^\circ}\end{equation*}

If we now use equations 9 and 10 we can identify the maximum shear stress magnitude as \tau_{max}=0.84 \: N/mm^2. The mutually perpendicular planes on which these shear stresses act are orientated at -13.3^\circ and 76.7^\circ. We can confirm this graphically by plotting these values on top of a plot of the shear stress transformation equation, Fig 11. 

Stress Analysis-Mohr's Circle 11 | DegreeTutors.com
Fig 11. Variation of shear stress with angle with max shear stress magnitude and orientation of planes of max positive and negative shear stress identified.

Before reading the next paragraph, make sure you’ve reviewed the discussion of shear stress sign convention in video 2. With reference to Fig 11. we can see that when the element is orientated at -13.3^\circ, i.e. when the outward pointing normal from the positive x-face makes an angle of -13.3^\circ with the global positive x-axis, the shear stress is \tau_{max}=0.84\:N/mm^2. However, when the element (and its reference axis system) is rotated by 90^\circ, the shear stress acting on the face orientated at 76.7^\circ is \tau_{max}=-0.84\:N/mm^2. The change in sign here is because the reference axis for the element has rotated through 90^\circ but all of the shear stress arrow directions remain unchanged. If this last paragraph has totally confused you, make sure to go back and review videos 2 and 5! 

6. Mohr’s Circle and Building our Mohr’s Circle Calculator

Now that we have the ground work covered we can introduce Mohr’s circle. As we stated at the outset, Mohr’s circle allows us to capture the complete state of stress at a point in a single graphic. We’ll break this discussion into 2 parts. First we will outline the basic equations; the equation of the circle itself, the radius and the centre of the circle. We’ll also outline the process of building a Mohr’s circle given the necessary starting information. In the second phase of the discussion we’ll use Mohr’s circle to perform an example stress analysis. It’s in solving this example in our Jupyter Notebook that we effectively build our Mohr’s circle calculator. 

6.1 Building Mohr’s circle from the transformation equations

The equation of the circle (Mohr’s circle) is derived from the transformation equations. We walk through the process in video 8. Here we’ll summarise the key points. We can show that the equation of the circle is,

(11)   \begin{equation*}\boxed{(\sigma_{x1} - \sigma_{avg})^2 + \tau_{x1y1}^2 = R^2}\end{equation*}

where R represents the radius of the circle and is given by, 

(12)   \begin{equation*}\boxed{R = \sqrt{\left(\frac{\sigma_x-\sigma_y}{2}\right)^2 + \tau_{xy}^2}}\end{equation*}

and \sigma_{avg} = 0.5(\sigma_x + \sigma_y). If we compare equation 11 to the standard equation of a circle with radius r and centre at coordinates (h,k),

(13)   \begin{equation*}(x-h)^2 + (y-k)^2 = r^2\end{equation*}

we see that the centre of Mohr’s circle is located at coordinates (\sigma_{avg}, 0). Now that we know the coordinates of the centre of the circle and the radius, we have enough information to actually build or draw the circle. Since every point on the circle represents the stresses on a particular set of planes – the circle captures the complete state of stress at the point. Let’s assume we know the normal and shear stresses (\sigma_x, \sigma_y and \tau_{xy}) on a stress element at a known orientation – say our default orientation where the planes align with the global x-y axis system. From here, the process of actually building and using Mohr’s circle is pretty much the same every time:

  • Establish an axis system for the circle where the horizontal axis measures \sigma_{x1} and is positive to the right and the vertical axis measures \tau_{x1y1} and is positive downwards.
  • Locate the centre of the circle at (\sigma_{avg},0).
  • Locate point ‘A’ on the circle representing the stresses you started with at the default orientation, (\sigma_x, \sigma_y and \tau_{xy}). 
  • Locate a point ‘B’ on the circle which represents the stresses on the adjacent perpendicular face of the element. Point ‘B’ will therefore have coordinates (\sigma_y, -\tau_{xy}) on the circle.

🚨 It’s important to note that a rotation of \theta degrees on our 2D stress element corresponds to a rotation of 2\theta around Mohr’s circle.

6.2 Mohr’s Circle Example

Now we can put Mohr’s circle to use and complete an example stress analysis. We’re going to implement our calculations and plot the Mohr’s circle for this analysis in our Jupyter Notebook. The beauty of coding this is that once complete we will have built a Mohr’s circle calculator! The process is covered in video 9. I’ll briefly discuss the output of our analysis here. We start with a candidate question to answer; consider the following 2D stress element subject to the stresses shown. 

Stress Analysis-Mohr's Circle 12 | DegreeTutors.com
Fig 12. 2D stress element subject to normal and shear stresses.

The first task is to work out \sigma_{avg} to identify the centre of the circle and then to calculate the radius of the circle,

(14)   \begin{equation*}\sigma_{avg} = \frac{\sigma_x + \sigma_y}{2} = \frac{-32.5+8.2}{2} = -12.15\:N/mm^2\end{equation*}

(15)   \begin{equation*}R = \sqrt{\left(\frac{\sigma_x-\sigma_y}{2}\right)^2 + \tau_{xy}^2} = \sqrt{\left(\frac{-32.5-8.2}{2}\right)^2 + (-10)^2} = 22.67\:N/mm^2\end{equation*}

Now we can draw the circle, or plot it using Python in our case and immediately identify the principal stresses, maximum shear stresses and their associated planes. 

Stress Analysis-Mohr's Circle 13 | DegreeTutors.com
Fig 13. Mohr’s circle showing the principal stresses (red crosses), maximum positive and negative shear stresses (red diamonds), stresses on the planes depicted by the 2D stress element in Fig 12. above (green dot for x face and yellow dot for y face)

We can see that the principal stresses can be identified from knowledge of the circle centre and radius,

(16)   \begin{equation*}\boxed{\sigma_{1,2} = \sigma_{avg}\pm R = -12.15 \pm 22.67 \:N/mm^2}\end{equation*}

Similarly the maximum positive and negative shear stress is simply \tau_{max} = \pm R = \pm 22.67\:N/mm^2

The angles of the principal planes and the planes of maximum positive and negative shear stress can be easily identified using basic trigonometry and using the green dashed line as the reference line that represents \theta=0^\circ i.e. the orientation of the element shown in Fig. 12. Just remember that an angle of 2\theta on Mohr’s circle corresponds to a rotation of \theta degrees on our element. The fact that our vertical axis is positive downwards also means that the sense of rotation of the element matches the sense of rotation on our circle, i.e. a counter-clockwise rotation of the element by \theta degrees corresponds to a counter-clockwise rotation of 2\theta degrees from our reference axis on the circle (green dashed line in our case). 

As I’ve said – the beauty of building this solution using Python in our Jupyter Notebook (video 9) means that simply by changing the input parameters, i.e. the stresses shown in Fig. 12, we can automatically generate the corresponding Mohr’s circle from which we can then determine the stresses at any orientation.

So there you have it – we’ve covered everything from the basics of normal and shear stress all they way up to Mohr’s circle. If you’ve been following along and watched the videos series accompanying this post, you’ll also have built yourself a handy Mohr’s circle calculator. Aside from developing knowledge of the underlying theory – the real win with this project is that you should feel confident enough to start tinkering with your code and making additions to build a stress analysis calculator that suits you. 

Now you’ve seen how we can automate a standard engineering calculation – what other repetitive calculations can you automate? Why not start building your own library of analysis notebooks?!

Finite Element Analysis of Continuum Structures in Python

Use the Isoparametric Finite Element Method to build an analysis tool for 2D structures in Python.

Finite-Element-Analysis-of-2D-Solid-Structures-in-Python (M) | DegreeTutors.com

After completing this course…

  • You will have the tools to analyse continuum structures using your own Isoparameteric Finite element Python code, developed from the ground up.
  • You will understand how the plane stress and plane strain approximations allow us to analyse 3D structures accurately with 2D planar models.
  • You will be able to use open source tools to generate structural models and mesh data that can be analysed with your FE code.

If you want to dive deeper into stress analysis and build your own finite element analysis code, take a look at the course above. In that course, you’ll put the knowledge and coding skills you picked up here to good use and build a complete 2D finite element analysis tool, right inside another Jupyter Notebook!

That’s all for this one – see you in the next one. 👍

Tutorial by:

Dr Seán Carroll

Hi, I’m Seán, the founder and lead tutor at DegreeTutors.com. I hope your found whatever I wrote up there helpful! Feel free to get in touch or follow DegreeTutors on any of the social accounts.

Finite Element Analysis of Continuum Structures in Python

Use the Isoparametric Finite Element Method to build an analysis tool for 2D structures in Python.

Play Video

In this course we’re going to build a full Isoparametric Finite Element solver. This unlocks our ability to model and analyse any 2D structural form. 

When we combine our finite element analysis solver with knowledge of the theory of plane stress and plane strain, this course will equip you with the ability to model 3D structures using 2D finite element models. When you complete this course you will have built your own 2D Finite Element solver, but more importantly, you’ll understand exactly how it works and what every single line of code does! 

Once complete, you can use your solver to show:

  • deflected shapes
  • normal stress and strain fields
  • shear stresses and strains fields
  • principal stress magnitudes fields
  • principal plane orientations
  • von Mises stress fields

We’ll build in the ability to simulate the influence of point load forces, distributed forces and body or self-weight forces. Once you complete this course you’ll have the knowledge, experience and confidence to extend your solver and add the new features that matter to you.

You DO NOT need to be a Python programming guru to take this course. If you’ve taken any of the prerequisite courses – or even if you’re just familiar with basic programming ideas like functions, loops and variables that will be plenty to get you started. 👍