Drawing the Mohr’s Circle in MATLAB

Aim – To draw the Mohr’s Circle in MATLAB.

Objective –

1. To draw the Mohr’s circle for 3-D stress tensor.
2. To draw the Mohr’s circle for 2-D stress state.
3. To calculate the principal stresses in both 3-D and 2-D cases.

Theory –

Mohr's circle, invented by Christian Otto Mohr, is a two-dimensional graphical representation of the transformation law for the Cauchy stress tensor.

In a part subjected to some forces, stress is generally distributed as a continuously varying function within the continuum of material. Every infinitesimal element of the material can conceivably experience different stresses at the same time. Thus, we must look at stresses as acting on vanishingly small elements within the part. These infinitesimal elements are typically modeled as cubes. The stress components are considered to be acting on the faces of these cubes in two different manners. Normal stresses act perpendicular (i.e., normal) to the face of the cube and tend to either pull it out (tensile normal stress) or push it in (compressive normal stress). Shear stresses act parallel to the faces of the cubes, in pairs (couples) on opposite faces, which tends to distort the cube into a rhomboidal shape. These normal and shear components of stress acting on an infinitesimal element make up the terms of a tensor.

For any particular combination of applied stresses, there will be a continuous distribution of the stress field around any point analyzed. There will always be planes on which the shear-stress components are zero. The normal stresses acting on these planes are called the principal stresses. The planes on which these principal stresses act are called the principal planes. The directions of the surface normal to the principal planes are called the principal axes, and the normal stresses acting in those directions are the principal normal stresses. There will also be another set of mutually perpendicular axes along which the shear stresses will be maximal. The principal shear stresses act on a set of planes that are at 45° angles to the planes of the principal normal stresses.

Mohr’s circles have long provided a means to do a graphical solution and find the principal stresses for the plane stress case. Before the advent of programmable calculators and computers, Mohr’s graphical method was a reasonable and practical way to solve equation.

The abscissa and ordinate of each point on the circle are the magnitudes of the normal stress and shear stress components, respectively, acting on the rotated coordinate system. In other words, the circle is the locus of points that represent the state of stress on individual planes at all their orientations, where the axes represent the principal axes of the stress element.

According to Cauchy, the stress at any point in an object, assumed as a continuum, is completely defined by the nine stress components of a second order tensor of type known as the Cauchy stress tensor,

Mohr’s Circle for 2-D state of stress – In two dimensions, the stress tensor at a given material point P with respect to any two perpendicular directions is completely defined by only three stress components. For the particular coordinate system these stress components are: the normal stresses and the shear stress.

Mohr’s Circle for 3-D state of stress – Mohr’s circles can also be drawn for the 3-D stress case as well, but a graphical construction method is not available to create them directly from the applied-stress data except for the special case where one of the principal stresses is coincident with an axis of the xyz coordinate system chosen, i.e., where one plane is a plane of principal stress.

where, σ is the principal stress magnitude and n_x, n_y, and n_z are the direction cosines of the unit vector n, which is normal to the principal plane:

n.n = 1 ;

n = n_xi + n_yj + n_zk ;

However, once the principal stresses are calculated by a suitable root-finding technique, the 3-D Mohr’s circles can be drawn by using the calculated principal stresses.

Governing Equation –

This equation is a cubic polynomial in σ. The coefficients C₀, C₁, and C₂ are called the tensor invariants, because they have the same values regardless of the initial choice of xyz axes in which the applied stresses were measured or calculated. The units of C₂ are MPa, of C₁ MPa², and of C₀ MPa³. The three principal normal stresses σ1, σ2, σ3 are the three roots of this cubic polynomial. The roots of this polynomial are always real and are usually ordered such that σ1 > σ2 > σ3.

The principal shear stresses can be found from the values of the principal normal stresses using

For 2-D case –

For the special case of a two-dimensional stress state, the equations for principal stress reduce to

Code -

For 2-D case – (Inputs are assumed for plotting as 160MPa, -80MPa, 50MPa) –

clear all
close all
clc
 
% Heading
disp('Drawing the Mohr Circle');
fprintf('\n');
 
% Input stresses
sigma_x = input('Enter the value of normal stress on x face (MPa) : ');
sigma_y = input('Enter the value of normal stress on y face (MPa) : ');
tau_xy = input('Enter the value of shear stress on xy plane (MPa) : ');
tau_yx = -tau_xy;
 
% Calculating the parameters of Mohr's circle
% Radius
R = sqrt(((sigma_x-sigma_y)/2)^2 + (tau_xy)^2);
 
% Centre
C = (sigma_x+sigma_y)/2;
 
% Principal Normal stresses
sigma1 = C+R;
sigma2 = C-R;
 
% Defining theta
theta = linspace(0,2*pi,1080);
 
% Center of the circle
Center = [(sigma1+sigma2)/2,0];
 
% Principal Shear stresses
tau12 = (sigma1-sigma2)/2;
tau21 = (sigma2-sigma1)/2;
 
 
% Plotting
 
O = [0,0];
hold on
plot([O(1) sigma_x],[0 0],'k')
plot([O(1) sigma_y],[0 0],'k')
plot([sigma_x sigma_x],[0 tau_xy],'k')
plot([sigma_y sigma_y],[0 -tau_xy],'k')
 
plot([sigma_x sigma_y],[tau_xy -tau_xy],'g')
plot(C,0,'*','color','r')
 
circle = plot(Center(1)+tau12*cos(theta),Center(2)+tau12*sin(theta));
axis('equal')
grid on
hold on


% Labels 
if sigma1>0
    text(sigma1*1.01,0,'\sigma1','fontsize',7);
else
    text(sigma1*0.95,0,'\sigma1','fontsize',7);
end
 
if sigma2>0
    text(sigma2*1.0,0,'\sigma2','fontsize',7);
else
    text(sigma2*(0.99),0,'\sigma2','fontsize',7);
end
 
text(sigma_x*1.01,0,'\sigma_x','fontsize',8);
text(sigma_y*1.01,0,'\sigma_y','fontsize',8);
text(sigma_x*1.05,tau_xy,'\tau_x_y','fontsize',8);
text(sigma_y*1.05,-tau_xy,'\tau_x_y','fontsize',8);
 
text(C*1.1,0,'C','fontsize',10)
text(C,tau12*0.9,'\tau12','fontsize',10);
 
xlabel('Normal Stress,\sigma(MPa)')
ylabel('Shear Stress,\tau(MPa)')
title('2D Mohr Circle')

Result –

1. Mohr’s Circle –

2. Workspace Output –

Explanation –

• First we use ‘clear all’, ‘close all’, ‘clc’ to clear of any previous data or plot.
• In this case, we have used ‘input’ command so that we don’t have to edit the program if the value for input stresses change.
• For plotting purposes, we have assumed the value of input stresses as 160 MPa, -80 MPa and 50 MPa.
• Then we calculate the center (C) and radius (R) of our desired Mohr’s circle with the help of formulas stated in 2-D stress state condition.
• The value of σ₁ is equal to the sum of C and R while the value of σ₂ is equal the the difference of C and R. The value of τ_max or τ₁₂ is equal to the value of R or we can say that it is half of the difference between σ₁ and σ₂.
• As all the values required are calculate, we can now plot our Mohr’s circle. First we define the coordinates of origin ‘O’ as (0,0).
• Then we draw the lines equal to the value of σ_x and σ_y from the origin parallel to x-axis. After that, we draw the lines for τ_xy from σ_x and σ_y.
• Then we draw a line connecting the coordinates (σ_x , τ_xy) and (σ_y , -τ_xy) which bisects the line σ_x-σ_y at point C, which is the center of our Mohr’s circle. -
• With ‘C’ as center and and ‘R’ as radius we plot the desired Mohr’s circle.
• After that we label all the points and also x-axis and y-axis. ‘Title’ command is used to give the title to the diagram.

For 3-D case – (Values in Stress tensor (T) are assumed in the beginning) –

clear all
close all
clc
 
% Stress tensor
T = [160 50 -100; -50 -80 90; 100 -90 120];
 
sigma_x = T(1,1);
sigma_y = T(2,2);
sigma_z = T(3,3);
 
tau_xy = T(1,2);
tau_yz = T(2,3);
tau_zx = T(3,1);
 
% Coefficients for Mohr's Ciecle
c3 = 1;
c2 = sigma_x+sigma_y+sigma_z;
c1 = sigma_x*sigma_y+sigma_y*sigma_z+sigma_z*sigma_x -tau_xy^2 -tau_yz^2 -tau_zx^2;
c0 = sigma_x*sigma_y*sigma_z + 2*tau_xy*tau_yz*tau_zx -(sigma_x*tau_yz^2+sigma_y*tau_zx^2+sigma_z*tau_xy^2);
 
% Principal stresses
normal_stresses = roots([c3 c2 c1 c0]);
A = sort(normal_stresses,'descend');
 
sigma1 = A(1);
sigma2 = A(2);
sigma3 = A(3);
sigma = [sigma1;sigma2;sigma3];
 
sigma_max = max(sigma)
sigma_min = min(sigma)
 
tau13 = (sigma1-sigma3)/2;
tau12 = (sigma1-sigma2)/2;
tau23 = (sigma2-sigma3)/2;
tau = [tau13;tau12;tau23];
 
tau_max = max(tau)
 
% Plotting
theta = 0:0.01:2*pi;
 
C1 = [(sigma1+sigma3)/2 0];
C2 = [(sigma1+sigma2)/2 0];
C3 = [(sigma2+sigma3)/2 0];
 
circle1 = plot(C1(1)+tau13*cos(theta),C1(2)+tau13*sin(theta));
hold on
circle2 = plot(C2(1)+tau12*cos(theta),C2(2)+tau12*sin(theta));
hold on
circle3 = plot(C3(1)+tau23*cos(theta),C3(2)+tau23*sin(theta));
hold on
axis('equal')
grid on


% Labels
if sigma1>0
    text(sigma1*1.01,0,'\sigma1','fontsize',10);
else
    text(sigma1*0.95,0,'\sigma1','fontsize',10);
end
 
if sigma2>0
    text(sigma2*1.0,0,'\sigma2','fontsize',10);
else
    text(sigma2*0.99,0,'\sigma2','fontsize',10);
end
 
if sigma3>0
    text(sigma3*1.1,0,'\sigma3','fontsize',10);
else
    text(sigma3*0.99,0,'\sigma3','fontsize',10);
end
 
designation = ['sigma1';'sigma2';'sigma3'];
values_MPa = [sigma1;sigma2;sigma3];
Principal_stresses = table(sigma,values_MPa);
    
text(C1(1),tau13*0.9,'\tau13','fontsize',10);
text(C2(1),tau12*0.9,'\tau12','fontsize',10);
text(C3(1),tau23*0.9,'\tau23','fontsize',10);
xlabel('Normal Stress,\sigma (MPa)')
ylabel('Shear Stress,\tau (MPa)')
title('3D Mohr Circle')

Result –

1. Mohr’s Circle –

2. Work-space Output –

Explanation –

• First we use ‘clear all’, ‘close all’, ‘clc’ to clear of any previous data or plot.
• Instead of putting the input values every time as in case of 2-D stress condition, we have a given some random values to stress tensor 'T'.
• With the use of indexing of matrix, we calculate σ_x, σ_y, σ_z, τ_xy, τ_yz and τ_zx.
• Then we write the formula for coefficients (C₀, C₁, C₂ and C₃) of cubic polynomial of σ.
• Since normal stresses are roots of the cubic polynomial equation, hence to find out the roots (i.e., normal stresses), we use ‘roots’ command with the coefficients of the polynomial. Then we sort them out in descending order.
• The biggest value numerically is our principal normal stress and in this case, it is σ₃ (-268.73 MPa).
• Then we find out the values of τ₁₃, τ₁₂ and τ₂₃. The value of τ_max is the maximum value between these three.
• Now, as we have found all the values so we can start plotting the Mohr’s circle. First, we define the value of theta (θ) from 0 to 2π radian with an interval of 0.01.
• Then we define the coordinates of center of the circles named C1, C2 and C3.
• Then we use the plot command to draw the circles as theta (θ) varies with center as specified above. To make the plot look like a circle instead of an ellipse, we use equal axis and ‘grid on’ command to clearly view the values in Mohr’s circle.
• Our Mohr’s circle has been drawn correctly. Now we have to label the important points in the diagram. To label the values of σ₁, σ₂ and σ₃, we are going to use the if-else loop. Then we label the values of τ₁₃, τ₁₂ and τ₂₃.
• Then we label the x-axis as ‘Normal stress (MPa)’ and y-axis as ‘Shear stress (MPa)’. ‘Title’ command is used to give the title to the figure.

Conclusion –

Thus the study of Mohr’s circle is completed and analysis is performed with the help of MATLAB. We drew the Mohr’s circle for both 2-D and 3-D cases with two slightly different approaches which we learnt during the completion of course and previous challenges.

Reference –

1. https://en.wikipedia.org/wiki/Mohr%27s_circle#:~:text=The%20Mohr%20circle%20is%20used,Pole%20or%20origin%20of%20planes.

2. Textbook – Machine Design, An integrated approach by Robert L. Norton