ME323 - Mechanics of Materials

Quick audit of the Mechanics of Materials with Python to refresh my memory on MoM so I can design a bookcase.

Examples/Problems out of Mechanics of Materials, 3$^{rd}$ edition. by John T.; Johns Ferdinand Pierre; Dewolf

In [1]:
%pylab inline
import pint  # Units.
import sympy # Symbolic math
from IPython.display import Latex # Display latex from Python windows.
Populating the interactive namespace from numpy and matplotlib

Introduction - Concept of Stress

Normal stress in a member under axial loading: $\sigma=\frac{P}{A}$

To define the stress at a given point $Q$, consider a small area $\Delta A$.

$$\sigma=\lim_{\Delta A \to 0} \frac{\Delta F}{\Delta A}$$

Then more hand waving.

$$P = \int dF = \int_{A} \sigma dA$$

Then there's shearing.

$$\tau_{ave}=\frac{P}{A}$$
1.8 Application to the analysis and design of simple structures.

1.8a: Determination of the Normal Stress:

$F_{BC}=$50kN, $A=314\times10^{-6}$m$^2$ and the corresponding average normal stress is $\sigma_{BC}=159$ MPa. Lets verify the math with Python

In [2]:
# Setup units.
u=pint.UnitRegistry()
F_BC=50*u.kN
A=(314*10**-6)*u.m**2 # ^ is ** in Python. Why? Because. Deal with it.
sigma=F_BC/A
display(Latex("$\\sigma=$ {:.2f}".format(sigma)))
display(Latex("$\\sigma=$ {:.2f}".format(sigma.to_base_units())))
display(Latex("$\\sigma=$ {:.2f}".format(sigma.to(u.MPa))))
$\sigma=$ 159235.67 kilonewton / meter ** 2
$\sigma=$ 159235668789.81 gram / meter / second ** 2
$\sigma=$ 159.24 megapascal

1.8b: Determination of Shearing Stress in Various Connections

Shear in the plane is P=50 kN. The pin diameter is 25 mm.

In [3]:
# Some unit aliases.
kN=u.kN
mm=u.mm
MPa=u.MPa
In [4]:
P=50*kN
d=25*mm
r=d/2
A=pi*r**2

tau=P/A

display(Latex(r"$\tau=$ {:.2f}".format(tau)))
display(Latex(r"$\tau=$ {:.2f}".format(tau.to_base_units())))
display(Latex(r"$\tau=$ {:.2f}".format(tau.to(MPa))))
$\tau=$ 0.10 kilonewton / millimeter ** 2
$\tau=$ 101859163578.81 gram / meter / second ** 2
$\tau=$ 101.86 megapascal

Problem 1.16

Determine the diameter of the largest circular hole that can be punched into a sheet of polystyrene 6 mm thick, knowing that the force exerted by the punch is 45kN and that a 55 MPa average shearing stress is required to cause the material to fail.

Lets do some object-oriented programming. It may not look like it now, but later on these will save us some time.

Circles are good shapes to have in toolbox.

In [5]:
class Circle(object):
    def __init__(self,diameter=None):
        self.diameter=diameter
        
    @property
    def radius(self):
        return self.diameter/2
    
    @property
    def circumference(self):
        return 2*pi*self.radius  
    
    @property
    def area(self):
        return pi*self.radius**2
## And test
# Numeric Math.
circle=Circle(diameter=10*u.m)
print("r={:.2f}".format(circle.radius))
print("A={:.2f}".format(circle.area))
print("C={:.2f}".format(circle.circumference))
# Symbolic math
d=sympy.symbols('d')
circle2=Circle(diameter=d)
print("r={}".format(circle2.radius))
print("A={}".format(circle2.area))
print("C={}".format(circle2.circumference))
r=5.00 meter
A=78.54 meter ** 2
C=31.42 meter
r=d/2
A=0.785398163397448*d**2
C=3.14159265358979*d

And back to the question. Shearing stress needed to cut 55MPa. Force exerted by punch 45 kN.

$$\tau_{ave} = \frac{P}{A} $$

A in this case is the area of the 'cylinder' that the punch will cut out, not the area of the circle. It's the circumference of the circle by the height/thickness of the sheet. (6mm)

First, lets create some shortcuts. Eq as shorthand for sympy.Eq, solve as shorthand for sympy.solve and units.

In [6]:
Eq=sympy.Eq
solve=sympy.solve
mm=u.mm
kN=u.kN
MPa=u.MPa

Next set up symbolic variables and create a circle (the punch) with unknown diameter d.

In [7]:
A, d, h, tau, P=sympy.symbols("A d h tau P")
punch=Circle(diameter=d)

Our basic shear formula displayed in $\LaTeX$ and again with python's print function.

In [8]:
eqn1=Eq(tau,P/(A))
display(Latex(r"$$\tau_{ave} = \frac{P}{A}$$"))
print(eqn1)
$$\tau_{ave} = \frac{P}{A}$$
tau == P/A

Substitute the calculation for A, displayed in $\LaTeX$ and again as represented by Python symbolic math.

In [9]:
eqn2=eqn1.subs(A,punch.circumference*h)
display(Latex(r"$$\tau_{ave} = \frac{P}{circ*h}$$"))
print(eqn2)
$$\tau_{ave} = \frac{P}{circ*h}$$
tau == 0.318309886183791*P/(d*h)

Solve for the diameter:

In [10]:
diameter=sympy.solve(eqn2,d)
print("d={}".format(diameter))
d=[0.318309886183791*P/(h*tau)]

Substitute values to get our answer:

In [11]:
diameter[0].subs(P,45*kN)
---------------------------------------------------------------------------
SyntaxError                               Traceback (most recent call last)
/usr/local/lib/python3.4/site-packages/sympy/core/sympify.py in sympify(a, locals, convert_xor, strict, rational, evaluate)
    316         a = a.replace('\n', '')
--> 317         expr = parse_expr(a, local_dict=locals, transformations=transformations, evaluate=evaluate)
    318     except (TokenError, SyntaxError) as exc:

/usr/local/lib/python3.4/site-packages/sympy/parsing/sympy_parser.py in parse_expr(s, local_dict, transformations, global_dict, evaluate)
    819 
--> 820     return eval_expr(code, local_dict, global_dict)
    821 

/usr/local/lib/python3.4/site-packages/sympy/parsing/sympy_parser.py in eval_expr(code, local_dict, global_dict)
    732     expr = eval(
--> 733         code, global_dict, local_dict)  # take local objects in preference
    734 

SyntaxError: invalid syntax (<string>, line 1)

During handling of the above exception, another exception occurred:

SympifyError                              Traceback (most recent call last)
<ipython-input-11-b60ca2e87b77> in <module>()
----> 1 diameter[0].subs(P,45*kN)

/usr/local/lib/python3.4/site-packages/sympy/core/basic.py in subs(self, *args, **kwargs)
    841         for i in range(len(sequence)):
    842             o, n = sequence[i]
--> 843             so, sn = sympify(o), sympify(n)
    844             if not isinstance(so, Basic):
    845                 if type(o) is str:

/usr/local/lib/python3.4/site-packages/sympy/core/sympify.py in sympify(a, locals, convert_xor, strict, rational, evaluate)
    317         expr = parse_expr(a, local_dict=locals, transformations=transformations, evaluate=evaluate)
    318     except (TokenError, SyntaxError) as exc:
--> 319         raise SympifyError('could not parse %r' % a, exc)
    320 
    321     return expr

SympifyError: Sympify of expression 'could not parse '45 kilonewton'' failed, because of exception being raised:
SyntaxError: invalid syntax (<string>, line 1)

Well that wasn't supposed to happen. It looks like sympy and pint have some issues playing nice. So this won't quite be as pretty as it is on the TI-89, so lets do it the 'hard' way. Copying the output d=[0.318309886183791*P/(h*tau)] of the solver and setting values

In [12]:
P=45*kN
h=6*mm
tau=55*MPa

d=0.318309886183791*P/(h*tau)
print("The diameter of the punch needs to be {:.1f} or less".format(d.to(mm)))
The diameter of the punch needs to be 43.4 millimeter or less

Solution in the back of the book. 43.4 mm.

Now that I've done a proof of concept on MoM with Python lets skip ahead to where I need for a bookshelf, shear and bending-moment diagrams of a simply supported beam. Chapters 1-4 will be seft as an exercize to the reader.

Chapter 5. Analysis and Design of Beams for Bending

Yada, yada, yada:

$$\frac{dV(x)}{dx}=w(x)$$$$\frac{dM(x)}{dx}=V(x)$$$$\Delta V(x) = \int w(x) dx$$$$\Delta M(x) = \int V(x) dx$$

Sample Problem 5.1

In [55]:
from sympy.functions.special.delta_functions import DiracDelta, Heaviside
In [14]:
x= sympy.symbols('x')
# These are all point loads, use DiracDelta https://en.wikipedia.org/wiki/Dirac_delta_function
w=-20*DiracDelta(x)+46*DiracDelta(x-2.5)+-40*DiracDelta(x-5.5)+14*DiracDelta(x-7.5)
In [15]:
V=sympy.integrate(w,x)
M=sympy.integrate(V,x)
In [16]:
V
Out[16]:
-20*Heaviside(x) + 14*Heaviside(x - 7.5) - 40*Heaviside(x - 5.5) + 46*Heaviside(x - 2.5)
In [49]:
pyplot.xkcd()
sympy.plot(V,(x,0,8), 
             adaptive=False, 
             xlabel="x (m)",
             ylabel="V (kN)",
             title="Shear Diagram for Sample Problem 5.1",
             ylim=(-30,30))
pyplot.xkcd()
sympy.plot(M,(x,0,7.5), 
             adaptive=False, 
             xlabel="x (m)",
             ylabel="V (kN)",
             title="Bending Diagram for Sample Problem 5.1",

           ylim=(-30,30))
/usr/local/lib/python3.4/site-packages/sympy/plotting/experimental_lambdify.py:165: UserWarning: The evaluation of the expression is problematic. We are trying a failback method that may still work. Please report this as a bug.
  warnings.warn('The evaluation of the expression is'
/usr/local/lib/python3.4/site-packages/matplotlib/font_manager.py:1288: UserWarning: findfont: Font family ['ComicNeue', 'Comic Neue', 'Humor Sans', 'Comic Sans MS'] not found. Falling back to Bitstream Vera Sans
  (prop.get_family(), self.defaultFamily[fontext]))
Out[49]:
<matplotlib.rc_context at 0x8159ecf28>

There's a warning and the bending moment failed. Apparently whatever V(x) integrated to divides by zero. If you give up the first time code doesn't work, you'll never make anything. And I'm sure as hell not doing this by hand again. I need a tool in my toolbox to do it for me.

Lets look at M.

In [51]:
M
Out[51]:
-20*x*Heaviside(x) + 14*Piecewise((0, 0.133333333333333*Abs(x) < 1), (1.0*x, 7.5*Abs(1/x) < 1), (7.5*meijerg(((2, 1), ()), ((), (1, 0)), 0.133333333333333*x), True)) - 40*Piecewise((0, 0.181818181818182*Abs(x) < 1), (1.0*x, 5.5*Abs(1/x) < 1), (5.5*meijerg(((2, 1), ()), ((), (1, 0)), 0.181818181818182*x), True)) + 46*Piecewise((0, 0.4*Abs(x) < 1), (1.0*x, 2.5*Abs(1/x) < 1), (2.5*meijerg(((2, 1), ()), ((), (1, 0)), 0.4*x), True))

Thats ugly. Lets see if sympy.pretty_print looks better.

In [52]:
sympy.pretty_print(M)
                        ⎛⎧                      0                         for 
                        ⎜⎪                                                    
                        ⎜⎪                                                    
                        ⎜⎪                    1.0⋅x                           
-20⋅x⋅Heaviside(x) + 14⋅⎜⎨                                                    
                        ⎜⎪                                                    
                        ⎜⎪    ╭─╮0, 2 ⎛2, 1       │                    ⎞      
                        ⎜⎪7.5⋅│╶┐     ⎜           │ 0.133333333333333⋅x⎟      
                        ⎝⎩    ╰─╯2, 2 ⎝      1, 0 │                    ⎠      

0.133333333333333⋅│x│ < 1⎞      ⎛⎧                      0                     
                         ⎟      ⎜⎪                                            
           │1│           ⎟      ⎜⎪                                            
   for 7.5⋅│─│ < 1       ⎟      ⎜⎪                    1.0⋅x                   
           │x│           ⎟ - 40⋅⎜⎨                                            
                         ⎟      ⎜⎪                                            
                         ⎟      ⎜⎪    ╭─╮0, 2 ⎛2, 1       │                   
      otherwise          ⎟      ⎜⎪5.5⋅│╶┐     ⎜           │ 0.181818181818182⋅
                         ⎠      ⎝⎩    ╰─╯2, 2 ⎝      1, 0 │                   

    for 0.181818181818182⋅│x│ < 1⎞      ⎛⎧               0                  fo
                                 ⎟      ⎜⎪                                    
                   │1│           ⎟      ⎜⎪                                    
           for 5.5⋅│─│ < 1       ⎟      ⎜⎪             1.0⋅x                fo
                   │x│           ⎟ + 46⋅⎜⎨                                    
                                 ⎟      ⎜⎪                                    
 ⎞                               ⎟      ⎜⎪    ╭─╮0, 2 ⎛2, 1       │      ⎞    
x⎟            otherwise          ⎟      ⎜⎪2.5⋅│╶┐     ⎜           │ 0.4⋅x⎟    
 ⎠                               ⎠      ⎝⎩    ╰─╯2, 2 ⎝      1, 0 │      ⎠    

r 0.4⋅│x│ < 1⎞
             ⎟
      │1│    ⎟
r 2.5⋅│─│ < 1⎟
      │x│    ⎟
             ⎟
             ⎟
 otherwise   ⎟
             ⎠

That's not much better. Thankfully sympy can render to $\LaTeX$.

In [54]:
display(Latex("$"+sympy.latex(M)+"$"))
$- 20 x \theta\left(x\right) + 14 \begin{cases} 0 & \text{for}\: 0.133333333333333 \left\lvert{x}\right\rvert

Looks like there's a divide by zero at 0 at 0, lets graph just after.

In [63]:
pyplot.xkcd()
sympy.plot(M,(x,0.001,7.5), 
             adaptive=False, 
             xlabel="x (m)",
             ylabel="V (kN)",
             title="Shear Diagram for Sample Problem 5.1")
/usr/local/lib/python3.4/site-packages/sympy/plotting/experimental_lambdify.py:165: UserWarning: The evaluation of the expression is problematic. We are trying a failback method that may still work. Please report this as a bug.
  warnings.warn('The evaluation of the expression is'
/usr/local/lib/python3.4/site-packages/matplotlib/font_manager.py:1288: UserWarning: findfont: Font family ['ComicNeue', 'Comic Neue', 'Humor Sans', 'Comic Sans MS'] not found. Falling back to Bitstream Vera Sans
  (prop.get_family(), self.defaultFamily[fontext]))
Out[63]:
<sympy.plotting.plot.Plot at 0x814a25c18>

That doesn't look like what's in the book. So lets try another way. This is M by inspection of the answer. Lets work backwards.

In [64]:
M2=-20*x*Heaviside(x)+46*(x-2.5)*Heaviside(x-2.5)-40*(x-5.5)*Heaviside(x-5.5)+14*(x-7.5)*Heaviside(x-7.5)
sympy.plot(M2,(x,0,8))
/usr/local/lib/python3.4/site-packages/matplotlib/font_manager.py:1288: UserWarning: findfont: Font family ['ComicNeue', 'Comic Neue', 'Humor Sans', 'Comic Sans MS'] not found. Falling back to Bitstream Vera Sans
  (prop.get_family(), self.defaultFamily[fontext]))
Out[64]:
<sympy.plotting.plot.Plot at 0x815afcf28>
In [65]:
# Cheap and dirty function to print an equation in LaTeX.
def L(M):
    display(Latex("$"+sympy.latex(M)+"$"))

Now calculate V & w(x)

In [77]:
L2=sympy.simplify(sympy.diff(M2, x))
w2=sympy.simplify(sympy.diff(M2, x, x))
L(L2)
L(w2)
$- 20 x \delta\left(x\right) + \left(14 x - 105.0\right) \delta\left(x - 7.5\right) - \left(40 x - 220.0\right) \delta\left(x - 5.5\right) + \left(46 x - 115.0\right) \delta\left(x - 2.5\right) - 20 \theta\left(x\right) + 14 \theta\left(x - 7.5\right) - 40 \theta\left(x - 5.5\right) + 46 \theta\left(x - 2.5\right)$
$- 20 x \delta^{\left( 1 \right)}\left( x \right) + \left(14 x - 105.0\right) \delta^{\left( 1 \right)}\left( x - 7.5 \right) - \left(40 x - 220.0\right) \delta^{\left( 1 \right)}\left( x - 5.5 \right) + \left(46 x - 115.0\right) \delta^{\left( 1 \right)}\left( x - 2.5 \right) - 40 \delta\left(x\right) + 28 \delta\left(x - 7.5\right) - 80 \delta\left(x - 5.5\right) + 92 \delta\left(x - 2.5\right)$
In [79]:
sympy.plot(L2,(x,0,8))
/usr/local/lib/python3.4/site-packages/matplotlib/font_manager.py:1288: UserWarning: findfont: Font family ['ComicNeue', 'Comic Neue', 'Humor Sans', 'Comic Sans MS'] not found. Falling back to Bitstream Vera Sans
  (prop.get_family(), self.defaultFamily[fontext]))
Out[79]:
<sympy.plotting.plot.Plot at 0x814bc2128>

Eh, close enough. Surely there has to be a better way. Lets create pythonic class objects to represent everything needed to solve this.

However to do that, I have to go all the way back to higshchool and start making a toolbox with what I learned along the way. (And why you'll need this later on)... So let us detour to Kindergarten to College Statics, just so we can continue here.