~Kindergarten thru Statics

K thru Statics

I had to take a break from Mechanics of Materials notebook to make a notebook for all the prerequisites. The reason being I was missing a Python representation of material I learned from Kindergarten.

Lets start here, but I'm only going to be make what I need for now. (i.e. The classes are no where near complete)

[For those not familiar with object oriented programming or Python classes: Python Classes and Object Oriented Programming]

In [1]:
%pylab inline
import pint  # Units.
import sympy # Symbolic math
from IPython.display import display,Latex # Display latex from Python windows.
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-c6bda6b4dd27> in <module>()
----> 1 get_ipython().magic('pylab inline')
      2 import pint  # Units.
      3 import sympy # Symbolic math
      4 from IPython.display import display,Latex # Display latex from Python windows.

/usr/local/lib/python3.4/site-packages/IPython/core/interactiveshell.py in magic(self, arg_s)
   2334         magic_name, _, magic_arg_s = arg_s.partition(' ')
   2335         magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2336         return self.run_line_magic(magic_name, magic_arg_s)
   2337 
   2338     #-------------------------------------------------------------------------

/usr/local/lib/python3.4/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line)
   2255                 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
   2256             with self.builtin_trap:
-> 2257                 result = fn(*args,**kwargs)
   2258             return result
   2259 

<decorator-gen-107> in pylab(self, line)

/usr/local/lib/python3.4/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    191     # but it's overkill for just that one bit of state.
    192     def magic_deco(arg):
--> 193         call = lambda f, *a, **k: f(*a, **k)
    194 
    195         if callable(arg):

/usr/local/lib/python3.4/site-packages/IPython/core/magics/pylab.py in pylab(self, line)
    154             import_all = not args.no_import_all
    155 
--> 156         gui, backend, clobbered = self.shell.enable_pylab(args.gui, import_all=import_all)
    157         self._show_matplotlib_backend(args.gui, backend)
    158         print ("Populating the interactive namespace from numpy and matplotlib")

/usr/local/lib/python3.4/site-packages/IPython/core/interactiveshell.py in enable_pylab(self, gui, import_all, welcome_message)
   3169         from IPython.core.pylabtools import import_pylab
   3170 
-> 3171         gui, backend = self.enable_matplotlib(gui)
   3172 
   3173         # We want to prevent the loading of pylab to pollute the user's

/usr/local/lib/python3.4/site-packages/IPython/core/interactiveshell.py in enable_matplotlib(self, gui)
   3118         """
   3119         from IPython.core import pylabtools as pt
-> 3120         gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
   3121 
   3122         if gui != 'inline':

/usr/local/lib/python3.4/site-packages/IPython/core/pylabtools.py in find_gui_and_backend(gui, gui_select)
    237     """
    238 
--> 239     import matplotlib
    240 
    241     if gui and gui != 'auto':

/usr/local/lib/python3.4/site-packages/matplotlib/__init__.py in <module>()
    120 # cbook must import matplotlib only within function
    121 # definitions, so it is safe to import from it here.
--> 122 from matplotlib.cbook import is_string_like, mplDeprecation, dedent, get_label
    123 from matplotlib.compat import subprocess
    124 from matplotlib.rcsetup import (defaultParams,

/usr/local/lib/python3.4/site-packages/matplotlib/cbook.py in <module>()
     31 from weakref import ref, WeakKeyDictionary
     32 
---> 33 import numpy as np
     34 import numpy.ma as ma
     35 

/usr/local/lib/python3.4/site-packages/numpy/__init__.py in <module>()
    178         return loader(*packages, **options)
    179 
--> 180     from . import add_newdocs
    181     __all__ = ['add_newdocs',
    182                'ModuleDeprecationWarning',

/usr/local/lib/python3.4/site-packages/numpy/add_newdocs.py in <module>()
     11 from __future__ import division, absolute_import, print_function
     12 
---> 13 from numpy.lib import add_newdoc
     14 
     15 ###############################################################################

/usr/local/lib/python3.4/site-packages/numpy/lib/__init__.py in <module>()
      6 from numpy.version import version as __version__
      7 
----> 8 from .type_check import *
      9 from .index_tricks import *
     10 from .function_base import *

/usr/local/lib/python3.4/site-packages/numpy/lib/type_check.py in <module>()
      9            'common_type']
     10 
---> 11 import numpy.core.numeric as _nx
     12 from numpy.core.numeric import asarray, asanyarray, array, isnan, \
     13                 obj2sctype, zeros

/usr/local/lib/python3.4/site-packages/numpy/core/__init__.py in <module>()
     12         os.environ[envkey] = '1'
     13         env_added.append(envkey)
---> 14 from . import multiarray
     15 for envkey in env_added:
     16     del os.environ[envkey]

ImportError: /lib/libgcc_s.so.1: version GCC_4.6.0 required by /usr/local/lib/gcc48/libgfortran.so.3 not found

Somewhere along the way I learned about 2D shapes and their properties.

In [61]:
class Shape2D(object):
    def __init__(self,area=None,perimeter=None):
        self.area=area
        self.perimeter=perimeter

Test, test, test. Lets make sure we're calculating things correctly.

In [62]:
class Rectangle(Shape2D):
    def __init__(self,width=None,height=None):
        self.width=width
        self.height=height
        self.sides=4
        
    
    def parallel_axis(self,area_moment,distance):
        return area_moment+self.area*distance**2
    
    @property
    def area(self):
        return self.height*self.width
    @property
    def perimeter(self):
        return 2*self.height+2*self.width   
    @property
    def Ix_c(self):
        # Ix about the centroid.
        return 1/12*self.width*self.height**3
    @property
    def Ix_0(self):
        # Ix about 0,0
        return 1/3*self.width*self.height**3
    @property
    def Iy_c(self):
        # Ix about the centroid.
        return 1/12*self.width**3*self.height
    @property
    def Iy_0(self):
        # Iy about 0,0
        return 1/3*self.width**3*self.height
    
    def __retr__(self):
        return "{}({}x{})".format(self.__class__.__name__,
                                      self.width,
                                      self.height)
    def __str__(self):
        return self.__retr__

# Square is just an special instance of rectangles.
class Square(Rectangle):
    def __init__(self,width=None,height=None,area=None,perimeter=None):
        if width is not None:
            super().__init__(width=width,height=width)
            return
        if height is not None:
            super().__init__(width=height,height=height)
            return
        if area is not None:
            side=math.sqrt(area)
        if perimeter is not None:
            side=perimeter/4
        super().__init__(width=side,height=side)
        return

width=5
height=10
my_rectangle=Rectangle(width, height)
assert my_rectangle.area==width*height
assert my_rectangle.perimeter==width+height+width+height

side=4
my_square=Square(side)
assert my_square.area==side*side
assert my_square.perimeter==side+side+side+side

And some classes to represent physical quantities.

In [63]:
class Force(object):
    def __init__(self,magnitude,location):
        self.magnitude=magnitude
        self.location=location
        
    def __repr__(self):
        return "{}<{}@{}>".format(self.__class__.__name__,self.magnitude,self.location)

    def __str__(self):
        return self.__repr__()
        
class Moment(object):
    def __init__(self,magnitude,location):
        self.magnitude=magnitude
        self.location=location
        
    def __repr__(self):
        return "{}<{}@{}>".format(self.__class__.__name__,self.magnitude,self.location)

    def __str__(self):
        return self.__repr__()
In [64]:
class Beam(object):
    # Types of beams, as defined on page 308, Mechanics of Materials, 3rd edition, ISBN: 0-07-365935-5
    # Statically determinate beams
    SIMPLY_SUPPORTED=1
    # TODO: The rest
    # Statically determinate beams
    def __init__(self,length,beam_type=SIMPLY_SUPPORTED,cross_section_shape=None):
        self.forces=list()
        self.moments=list()
        pass 

Lets solve the statics for Sample Problem 5.1

In [65]:
u=pint.UnitRegistry()
_kN=u.kN
_mm=u.mm
_m=u.m

P# to denote force magnitude. F# for the Force class instance of that force.

In [97]:
_kN=u.kN
_N=u.kN
P_b, P_d = sympy.symbols("P_b P_d")
P_a=-20*_kN
P_c=-40*_kN

## See edit below
P_b*=u.lbf
P_d*=u.lbf
In [98]:
F_a=Force(P_a,0*_m)
F_b=Force(P_b,2.5*_m)
F_c=Force(P_c,5.5*_m)
F_d=Force(P_d,7.5*_m)
In [99]:
beam=Beam(7.5*_m,cross_section_shape=Rectangle(80*_mm,250*_mm))
In [100]:
beam.forces.append(F_a)
beam.forces.append(F_b)
beam.forces.append(F_c)
beam.forces.append(F_d)
In [101]:
beam.forces
Out[101]:
[Force<-20 kilonewton@0 meter>,
 Force<P_b force_pound@2.5 meter>,
 Force<-40 kilonewton@5.5 meter>,
 Force<P_d force_pound@7.5 meter>]
In [102]:
sum_forces=0*_kN
sum_moments=0*_kN*_m
moment_about=0*_m
for force in beam.forces:
    sum_forces+=force.magnitude
    if force.location != moment_about:
        sum_moments+=force.magnitude*(force.location-moment_about)
In [103]:
for i, force in enumerate(beam.forces):
    if any([symbol in [force.magnitude.magnitude] for symbol in statics]):
        beam.forces[i].magnitude._magnitude=beam.forces[i].magnitude._magnitude.subs(statics)

~Well that's disappointing. The TI-89 can handle this just fine. But I don't have time to fix those toolboxes now.~ Original is at the bottom.

Edit 1:

/u/ecgite on reddit caught the error. I forgot to multiply the

In [ ]:
 
In [ ]:
Lets redo it assuming kN and meters.
In [73]:
kN=u.kN
P_b, P_d = sympy.symbols("P_b P_d")
P_a=-20
P_c=-40

Create each of the forces as a Force

In [74]:
F_a=Force(P_a,0)
F_b=Force(P_b,2.5)
F_c=Force(P_c,5.5)
F_d=Force(P_d,7.5)

Add each of the forces to the beam.

In [75]:
beam=Beam(7.5,cross_section_shape=Rectangle(0.080,0.250))
beam.forces.append(F_a)
beam.forces.append(F_b)
beam.forces.append(F_c)
beam.forces.append(F_d)

Lets do statics.

Statics = static.

F = Mass * acceleration.

Torque = Inertia * rotational acceleration.

Static = 0 acceleration.

$\sum F=0$ $\sum M=0$

In [76]:
sum_forces=0
sum_moments=0
moment_about=0
unknowns=set()
for force in beam.forces:
    sum_forces+=force.magnitude
    if isinstance(force.magnitude, sympy.Symbol):
        unknowns=unknowns | force.magnitude.free_symbols
    # Sum all of the forces
    sum_moments+=force.magnitude*(force.location-moment_about)
# Sum all of the moments.
for moment in beam.moments:
    sum_moments+=moment
In [77]:
def tex(eqn):
    display(Latex("$"+sympy.latex(eqn)+"$"))
In [78]:
tex(sympy.Eq(0,sum_forces))
tex(sympy.Eq(0,sum_moments))
$0 = P_{b} + P_{d} - 60$
$0 = 2.5 P_{b} + 7.5 P_{d} - 220.0$
In [79]:
unknowns
Out[79]:
{P_b, P_d}
In [80]:
statics=sympy.solve((sum_forces,sum_moments))
print(statics)
{P_b: 46.0000000000000, P_d: 14.0000000000000}
In [81]:
for i in range(len(beam.forces)):
    # If the beam force is symbolic, solve for it.
    if isinstance(beam.forces[i].magnitude, sympy.Symbol):
        beam.forces[i].magnitude=beam.forces[i].magnitude.subs(statics)
In [104]:
beam.forces
Out[104]:
[Force<-20 kilonewton@0 meter>,
 Force<46.0000000000000 force_pound@2.5 meter>,
 Force<-40 kilonewton@5.5 meter>,
 Force<14.0000000000000 force_pound@7.5 meter>]

And there we have, Python can almost do what the TI-89.

And just solved a 1D statics problem with python.

Future versions of each of the classes will live in their own modules.

Todo:

  • Add more 2D shapes.
  • Extend Forces for more than point loads.