Plotting stress-strain curves is a useful skill in mechanical engineering because it allows us to derive mechanical properties of materials such as tensile strength, elastic modulus, yield strength and ductility. Plotting stress strain curves can be accomplished with Excel and MATLAB, but what about plotting stress strain curves with Python? Follow along to find out.
We will use Python, matplotlib, pandas, holoviews and scipy to build the bar plot. I recommend undergraduate engineers use the Anaconda distribution of Python, which comes with matplotlib, pandas, and scipy already installed. For help installing Anaconda, see a previous blog post: Installing Anaconda on Windows 10. If matplotlib, pandas, holoviews, and scipy are not available in your version of Python, open a terminal or the Anaconda Prompt and type:
$ pip install matplotlib
$ pip install pandas
$ pip install holoviews
$ pip install scipy
or
> conda install matplotlib
> conda install pandas
> conda install holoviews
> conda install scipy
The data we are going to plot is stored in two Microsoft Excel Files. You can download the sample data here (clicking the link will start the download):
We'll use pandas to load the data into the notebook.
Note that when I first tried to run the pd.read_excel() function, I was returned an error:
ImportError: Install xlrd >= 0.9.0 for Excel support
To solve this, I went to the Anaconda Prompt and typed:
> conda install xlrdOnce the xlrd module was installed, the pd.read_excel() function worked just fine.
To start the Jupyter notebook, we need to import the required packages:
pandas numpy matplotlib scipy
The %matplotlib inline magic command is add so that we can see our plots right in the jupyter notebook.
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import linregress
%matplotlib inline
In [2]:
dfa=pd.read_excel('Al_6061.xls')
pd.options.display.max_rows = 999
dfa
Out[2]:
We have a couple columns in our pandas dataframe. The columns of interest are FORCE and EXT. The FORCE column conntians the force values in units of pounds (lb). The EXT column contains the extension or strain in units of %. So line 215 of our dataframe shows a force = 8206.19283 lb and an extension or strain of 3.003093% which is the same as a strain of 0.03003093 in/in (unitless strain rather than % strain).
To convert the FORCE column to stress $\sigma$, we need to divide the force $F$ (in lb) by the cross-sectional area $A_0$ of the sample.
$$ \sigma = \frac{F}{A_0} $$The cross-sectional area of the sample $A_0$ is the $\pi r^2$. Since we measured the diameter of the sample, and diamter $d$ is equal to two times the radius $r$
$$ A_0 = \pi (\frac{d}{2})^2 $$In our case diamter $d = 0.506 \ in$
In [3]:
d = 0.506 #units: in
A0 = np.pi*(d/2)**2 #units: in^2
stress_Al = dfa['FORCE'][1:172]*(1/A0)*(0.001)
strain_Al = dfa['EXT'][1:172]*0.01
l0 = 2.0 #units: in
lf = 2.24 #units: in
strain_Al_f = (lf-l0)/l0 #units: in/in
In [4]:
plt.plot(strain_Al, stress_Al)
plt.show()
To zoom in on the elastic region, we'll do something a little fancy, we'll spin up a bokeh plot using holoviews.
In [5]:
import holoviews as hv
hv.extension('bokeh')
import bokeh
In [6]:
hv.Curve((strain_Al, stress_Al),'strain (in/in)','stress (ksi)')
Out[6]:
From this plot of force vs. extension, we can zoom into the linear region. It looks like the linear elastic region is from 1.5 ksi to 39.5 ksi. Pandas has this nice little method .between(). We can set lower and upper limits from the stress series and then use this to index out of both the stress and the strain series in the linear elastic region.
In [7]:
E_stress = stress_Al[stress_Al.between(1.5,39.5)]
E_strain = strain_Al[stress_Al.between(1.5,39.5)]
E_strain
Out[7]:
Now we'll use scipy's linear regression function called linregress() to calculate the slope in the linear elastic region.
In [8]:
res = linregress(E_strain, E_stress)
res
Out[8]:
We see the slope = 9554.61 ksi. The rvalue = 0.999685 is very close to 1, so the linear elastic region is in fact pretty linear, not perfectly linear, but pretty close. Let's make a new variable E_Al for the elastic modulus of Aluminum
In [9]:
E_Al = res[0]
E_Al
Out[9]:
Now we will add a 0.002 or 0.2 % offset line to our plot. We will create a new linear series using the general form:
$$ y = mx + b $$$y$ will be the $\sigma$ value, $x$ will be the offset strain values ($\epsilon - 0.002$), $m$ (the slope) is the elastic modulus $E$, and $b$ is the y-intercept. Adapting $y = mx + b$ to our 0.002 off-set line looks like the equation below.
$$ \sigma_{offset} = E (\epsilon - 0.002) + 0 $$We can code $\sigma_{offset}$ line into a new series based on the equation above.
In [11]:
stress_offset = E_Al*(strain_Al - 0.002)
Now that the stress_offset series is defined, we can add our stress_offset line to the plot. When we call the plt.plot() command, we will pass in two x-y pairs. The first x-y pair is strain_Al, stress_Al, the second x-y pair is strain_Al, stress_offset.
In [12]:
plt.plot(strain_Al, stress_Al, strain_Al, stress_offset)
plt.show()
We can use holoviews again to zoom into the intersection point of the stress-strain curve (the blue line) and the 0.002 offset line (the orange line)
In [13]:
curve = hv.Curve((strain_Al, stress_Al),'strain (in/in)','stress (ksi)')
offset_curve = hv.Curve((strain_Al, stress_offset),'strain (in/in)','stress (ksi)')
curve * offset_curve
Out[13]:
If we use the zoom box (the zoom box button looks like a magnifying glass), we can clearly see the two lines cross at $\sigma_y = 46.28 \ ksi$. This compares well to a published value of around 40 ksi.
Let's define a variable to store yield strength.
In [57]:
sigma_y = 46.28
But what if we want to find this programatically, rather than by zooming into a plotted curve? Let's look at the points (without a line betweeen points) with holoviews. Instead of using the hv.Curve() method that makes line plots, we will use the hv.Scatter() method which builds a scatter plot of points.
In [15]:
scatter = hv.Scatter((strain_Al, stress_Al),'strain (in/in)','stress (ksi)')
offset_scatter = hv.Scatter((strain_Al, stress_offset),'strain (in/in)','stress (ksi)')
scatter * offset_scatter
Out[15]:
We need to find the two points on the blue stress strain curve where the crossover happens, and we need to find the the two red offset curve points where the crossover happens. Then we can do a little algebra and find the crossover point.
In [49]:
for i in range(2, len(strain_Al)):
#when the offset line is first goes above the stress strain curve...
if stress_offset[i] > stress_Al[i]:
# pull out the points from the stress strain curve
sx1 = strain_Al[i-1]
sy1 = stress_Al[i-1]
sy2 = stress_Al[i]
sx2 = strain_Al[i]
# pull out the points from the 0.002 offset line
ox1 = strain_Al[i-1]
oy1 = stress_offset[i-1]
ox2 = strain_Al[i]
oy2 = stress_offset[i]
break
print(ox1)
print(oy1)
print(ox2)
print(oy2)
print(sx1)
print(sy1)
print(sx2)
print(sy2)
In [84]:
# first set of two points, from the 0.002 offset line
x1 = ox1
y1 = oy1
x2 = ox2
y2 = oy2
# second set of two points, from the stress strain curve
x3 = sx1
y3 = sy1
x4 = sx2
y4 = sy2
# formula for the y corrdinate of the intersection two lines defined by two points
YS = ( (x1*y2-y1*x2)*(y3-y4) - (y1-y2)*(x3*y4-y3*x4) ) / ( ( (x1-x2) * (y3-y4) ) - ( (y1-y2)*(x3-x4) ))
In [85]:
YS
Out[85]:
So our yield strength is 46.279. Rounded to four sig figs, this is the same value we found for $\sigma_y$ when we zoomed into the intersection point on the plot.
We can round the value of Py to three sig figs and print it out using an f-string.
In [86]:
print(f'the yield strength of the aluminum is {round(Py,1)} ksi')
In [87]:
TS = max(stress_Al)
We can print out the tensile strength to three sig figs using an f-string
In [88]:
print(f'the tensile strength of the aluminum is {round(TS,1)} ksi')
In [89]:
xf = strain_Al[len(strain_Al)]
yf = stress_Al[len(strain_Al)]
Now using the point-slope formula for a line, we can determine where this line intecects the strain axis (the x-axis). The point-slope formula of a line is below.
$$ y - y_1 = m(x-x_1) $$Since we want to find the x-intercept, we set y = 0.
$$ 0 - y_1 = m(x-x_1) $$If we solve the equation above for x, we are left with
$$ x = x_1 - \frac{y_1}{m} $$Our known points $x_1$ and $y_1$ are the last strain and stress values, and the slope $m$ is the elastic modulus $E$. The $x$ value that we solve for is equal to the ductility ($\%EL$).
$$ \%EL = \epsilon_{last} - \frac{ \sigma_{last}}{E} $$We can code this with the line below
In [90]:
EL = xf - (yf/E_Al)
EL
Out[90]:
We can print out the ductility using an f-string.
In [91]:
print(f'the ductility is {round(EL*100,3)}%')
In [104]:
print(f'The tensile strength is {round(TS,1)} ksi')
print(f'The yield strength is {round(YS,1)} ksi')
print(f'The elastic modulus is {round(E_Al,0)} ksi')
print(f'the ductility is {round(EL*100,3)}%')
In [ ]: