Below is an engineering mechanics problem that can be solved with Python. Follow along to see how to solve the problem with code.
An I-beam (also called a W-shape for wide-flange shape) below with dimension $d$, $t_f$, $t_w$ and $b_f$. The $y$-axis runs vertically down though the center of the beam. The $x$-axis run horizontally through the center of the beam. The beam is symmetric about the $x$ and $y$ axes.
yield strength $\sigma_y = 300 \ MPa$
elastic modulus $E = 200 \ GPa$
(a) For values of $y_y$ from $\frac{1}{2}d$ to $\frac{1}{6}d$ in increments of $\frac{1}{2}t_f$, calculate the bending moment $M$ and the radius of curvature $\rho$
(b) Using the dimensions below, determine the bending moment $M$ and radius of curvature of the beam $\rho$ if the plastic zones at the top and bottom of the beam are $40 \ mm$ thick.
$d = 140 \ mm$
$t_f = 10 \ mm$
$t_w = 10 \ mm$
$b_f = 120 \ mm$
We will assume that for part (a), we have a standard W360 $\times$ 44 metric beam with the following measurements:
$d = 351 \ mm$
$t_f = 9.78 \ mm$
$t_w = 6.86 \ mm$
$b_f = 171 \ mm$
We are going to use Python to code the solution to this problem. If you don't already have Python installed on your computer, I recommend installing the Anaconda distribution of Python. See this post to learn how to install Anaconda on your computer.
I am coding this solution in a Jupyter Notebook. Once you install Anaconda, you can open a Jupyter notebook from the Windows Start Menu or the Anaconda Prompt. See this post to learn about 3 ways to open a Jupyter notebook.
Alternatively, instead of using a Jupyter notebook, you could code your solution in a .py file.
Alright. Let's get coding....
Based on our assumption of a standard W360 $\times$ 44 metric beam, and the given elastic modulus and yield strength given in the problem, we can define a number of constants.
$d = 351 \ mm$
$t_f = 9.78 \ mm$
$t_w = 6.86 \ mm$
$b_f = 171 \ mm$
yield strength $\sigma_y = 300 \ MPa$
elastic modulus $E = 200 \ GPa$
In [1]:
d = 351
tf = 9.78
tw = 6.86
bf = 171
ys = 300
E = 200*10**3 #Elastic modulus in MPa
Next we will compute the moment of inertia, $I$ based on the beam cross-section. The beam is made of three rectangular shapes. We will call the upper horizontal plate rectangle 1, the middle vertical plate rectangle 2, and the bottom horizontal plate rectangle 3.
The moment of inertia, $I_x$, of a rectangular shape is equal to:
$I_x = \frac{1}{12}bh^3 + Ad^2$
where $b$ is the width of the rectangle base, $h$ is the rectangle height, $A$ is the rectangle area, and $d$ is the distance between the centroid of the rectangle to the axis we are trying to find the moment of inertia about (in our case the x-axis).
We have three rectangles: The top, middle, and bottom of the beam. We will call the component each of these rectangles contributes to the total moment of inertia I1, I2, and I3.
After the components of the total moment of inertia are calculated, they can be summed together to produce the total moment of inertia $I$.
$I = I_1 + I_2 + I_3$
The Python code below completes these operations. At the end of the code section, a Python f-string is used to print out the value of $I$ after it is calculated. Make sure prepend f-strings with the letter f before the quotation marks.
In [2]:
A1 = tf*bf # area of the top rectangle
d1 = d/2 - tf/2 # distance between the centroid of the top rectangle and the x-axis
I1 = (1/12)*bf*tf**3 + A1*d1**2
A2 = tw*(d-2*tf) # area of the middle rectangle
d2 = 0 # the centroid of the middle rectangle is on the x-axis
I2 = (1/12)*tw*(d-2*tf)**3 + A2*d2**2
A3 = tf*bf # area of the bottom rectangle
d3 = d/2 - tf/2 # distance between the centroid of the bottom rectangle and the x-axis
I3 = (1/12)*bf*tf**3 + A3*d3**2
# sum the individual moments of inertia to caclulate the moment of inertia of the entire beam
I = I1 + I2 + I3
print(f"The moment of inertia I = {I} mm4")
After the moment of inertia, $I$, is calculated, the next step to solve this problem is to compute the maximum elastic moment, $M_y$. The maximum bending moment is limited by the yield strength, $\sigma_y$ of the material.
The equation for maximum elastic moment $M_y$ is below:
$M_y = \sigma_y\frac{I_x}{d/2}$
Where $\sigma_y$ is the yield strength of the material, $I_x$ is the moment of inertia and $d$ is the full-height of the beam.
The code below completes this calculation. A Python f-string is used to print out the calculated $M_y$ value. The units of $M_y$ come out to be $N \cdot mm$
In [3]:
My = ys*(I)/(d/2)
print(f"The maximum elastic moment My = {My} N mm")
Since the entire beam is assumed to be elastoplastic, we need to calculate the maximum value of $y$, which is typically denoted $c$. Since our $x$-axis runs runs horizontally through the center of the beam, the maximum value of $y$ is half of the total height of the beam
$ c = d/2 $
Coding this calculation is pretty simple
In [4]:
c = d/2
print(f"the maximum value of y is c = {c} mm")
Now let's find the radius of curvature $\rho$ if we just reach the yield stress, $\sigma_y$ at the very top of the beam.
The distance up the beam relative to the neutral axis where the yield strength is reached will be called $y_y$. If the yield strength is only reached at the very top of the beam, that means that $y_y = c = d/2$.
For a rectangular elastoplastic material, the distribution of strain across the section remains linear after the onset of yield. Strain is dependent on the bend radius $\rho$ and how far away from the neutral axis the strain occurs $y_y$.
$$ \epsilon = \frac{y_y}{\rho} $$$$ y_y = \epsilon\rho $$Now take the standard definition of elastic modulus $E$, stress $\sigma$ over strain $\epsilon$ and rearrange it to solve for strain $\epsilon$ in terms of $E$ and $\sigma$.
$$ E = \frac{\sigma}{\epsilon} $$$$ \epsilon = \frac{\sigma}{E} $$We can multiply both sides of the equation above by $\rho$.
$$ \epsilon\rho = \frac{\sigma}{E}\rho $$From above $y_y = \epsilon\rho$. Therefore, we can substitute $y_y$ in for $\epsilon\rho$ on the left-hand side of the equation.
$$ y_y = \frac{\sigma}{E}\rho $$Now we can solve the equation above for $\rho$ in terms of $y_y$, $E$, and $\sigma$.
$$ \rho = \frac{y_yE}{\sigma} $$This gives us an equation for radius of curvature $\rho$ in terms of parameters that we know, elastic modulus $E$, yield strength $\sigma_y$ and location in the beam compared to the neutral axis where plastic elongation starts to occur $y_y$.
Remember from before:
yield strength $\sigma_y = 300 \ MPa$
elastic modulus $E = 200 \ GPa$
And we are going to calculate the radius of curvature $\rho$ when the yield strength is only reached at the very top of the beam.
$y_y=c=d/2$
This can be coded in Python is just a few lines
In [5]:
yy=d/2
p = yy*E/ys
print("The radius of curvature if the yield strength is only reached")
print(f"at the very end of the beam (y_y=d/2) is p = {p} mm")
We can iterate through the $y_y$ distances (the distance compared to the neutral axis where plastic deformation begins and the stress in the beam has reached the yield strength) and calculate the radius of curvature $\rho$ for each $y_y$ distance. In Python, this can be accomplished with a for loop.
The code below iterates through $y_y$ values from $\frac{1}{2}d$ to $\frac{1}{6}d$ in increments of $\frac{1}{2}t_f$ and prints the resulting radius of curvature $\rho$ for each $y_y$ value. Remember that Python counting starts at zero and ends at n-1, so we need to tack on an extra $\frac{1}{2}t_f$ to the stop value of our np.arange() function. Since the $y_y$ values are not integers, we need to use NumPy and NumPy's np.arange() function instead of Python's build-in range() function. Make sure import numpy as np before calling np.arange().
In [6]:
import numpy as np
for yy in np.arange((1/2)*d, (1/6)*d + (1/2)*tf, -(1/2)*tf):
p = yy*E/ys
print(f"The radius of curvature at y_y = {round(yy,2)} mm is p = {round(p,1)} mm")
In [ ]: