pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
inertiaFormulas.py
Go to the documentation of this file.
00001 ## GNU Lesser General Public License
00002 ## 
00003 ## Program pyNastran - a python interface to NASTRAN files
00004 ## Copyright (C) 2011-2012  Steven Doyle, Al Danial
00005 ## 
00006 ## Authors and copyright holders of pyNastran
00007 ## Steven Doyle <mesheb82@gmail.com>
00008 ## Al Danial    <al.danial@gmail.com>
00009 ## 
00010 ## This file is part of pyNastran.
00011 ## 
00012 ## pyNastran is free software: you can redistribute it and/or modify
00013 ## it under the terms of the GNU Lesser General Public License as published by
00014 ## the Free Software Foundation, either version 3 of the License, or
00015 ## (at your option) any later version.
00016 ## 
00017 ## pyNastran is distributed in the hope that it will be useful,
00018 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
00019 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00020 ## GNU General Public License for more details.
00021 ## 
00022 ## You should have received a copy of the GNU Lesser General Public License
00023 ## along with pyNastran.  If not, see <http://www.gnu.org/licenses/>.
00024 ## 
00025 from sympy import Symbol,simplify,Rational,Matrix
00026 
00027 ## The purpose of this code is to automatically generate
00028 ## the moment of inertia equations.  only the locations of
00029 ## The nodes of the elements are required, as well as the
00030 ## triangular sections that are generated.
00031 ## For example, [1,a,b] means a positive triangular section (+1)
00032 ## with the 3 nodes at a, b, and <0,0>.  
00033 ## This is easily understood if the triangle creates an area.
00034 ## Some are positive, some are negative, but if you sum them all
00035 ## up you get the total area of the section.  This allows you
00036 ## to generate the moments of inertia for complex sections.
00037 
00038 def bar():
00039     """
00040          y
00041          ^
00042          |
00043     1----|----2
00044     |    |    |
00045     |    |    | h
00046     |    o----------> x
00047     |         |
00048     |         |
00049     4---------3
00050          b
00051 
00052     b = w1
00053     h = h1
00054     """
00055     print "---Bar---"
00056     #A = b*h
00057     h = Symbol('h')
00058     b = Symbol('b')
00059     half = Rational(1,2)
00060     A = b*h*half
00061 
00062     p1 = [-b*half, h*half]
00063    #p2 = [ b*half, h*half]
00064    #p3 = [ b*half,-h*half]
00065    #p4 = [-b*half,-h*half]
00066     
00067     p2 = addx(p1,b)
00068     p3 = suby(p2,h)
00069     p4 = subx(p3,b)
00070     print "p1 = ",p1
00071     print "p2 = ",p2
00072     print "p3 = ",p3
00073     print "p4 = ",p4
00074     Atri = half*p1[0]*p1[0]
00075     print "Atri = ",Atri
00076     
00077     sections = [ [1,p1,p2],[1,p2,p3],[1,p3,p4],[1,p4,p1] ]
00078     Ixx(sections)
00079 
00080 def cross():
00081     """
00082            y
00083     b1     ^      b3
00084            |
00085           b2
00086        a---p1--b              
00087        |   |   |               h1
00088     k--l   |   c-----d
00089     |      |         |
00090     |      o---------p2-> x    h2
00091     |                |
00092     j--i       f-----e
00093        |       |               h3
00094        |       |
00095        h-------g
00096     """
00097     #given
00098     #k-l = 0.5 D1
00099     #k-j = D4
00100     #h-g = D2
00101     #b-g = D3
00102     print "---Cross---"
00103 
00104     D1 = Symbol('D1')
00105     D2 = Symbol('D2')
00106     D3 = Symbol('D3')
00107     D4 = Symbol('D4')
00108     o2 = Rational(1,2)
00109     # bases
00110     b2 = D2
00111     b1 = b3 = D1*o2
00112 
00113     # heights
00114     h2 = D4
00115     h1 = h3 = (D3-D4)*o2
00116 
00117     p1 = [0,         D3*o2]
00118     p2 = [(D2+D1)*o2,    0]
00119 
00120     a = subx(p1,D2*o2)
00121     b = addx(p1,D2*o2)
00122     c = suby(b,h1)
00123     d = addx(c,b3)
00124     e = suby(d,h2)
00125     f = subx(e,b3)
00126     g = suby(f,h3)
00127     h = subx(g,b2)
00128     i = addy(h,h3)
00129     j = subx(i,b1)
00130     k = addy(j,h2)
00131     l = addx(k,b1)
00132     A = addy(l,h1)
00133     #print "a = ",a
00134     #print "b = ",b
00135 
00136     #print "d = ",d
00137     #print "e = ",e
00138 
00139     #print "g = ",g
00140     #print "h = ",h
00141 
00142     #print "k = ",k
00143     #print "j = ",j
00144 
00145     #print "A = ",A
00146     sections = [ [1,a,b],[1,b,c],[1,c,d],[1,d,e],
00147                  [1,e,f],[1,f,g],[1,g,h],[1,h,i],
00148                  [1,i,j],[1,j,k],[1,k,l],[1,l,a] ]
00149     Ixx(sections)
00150 
00151 
00152 def box():
00153     """
00154           y
00155           ^
00156           |
00157       w1  | w2  w3
00158           |
00159     a-----p1------b
00160     |     |       |  h1
00161     |  e--|----f  |
00162     |  |  |    |  |
00163     |  |  |    |  |  h2
00164     |  |  o-------|--> x
00165     |  |       |  |
00166     |  h-------g  |
00167     |             |  h3
00168     d-------------c
00169 
00170     # given
00171     #a-d     = D2
00172     #(f-g)_y = D3
00173     #d-c     = D1
00174     #(g-c)_x = D4
00175     """
00176     print "---Box---"
00177 
00178     D1 = Symbol('D1')
00179     D2 = Symbol('D2')
00180     D3 = Symbol('D3')
00181     D4 = Symbol('D4')
00182 
00183     # height
00184     h1 = h3 = D3
00185     h2 = D2-h1-h3
00186 
00187     # bases
00188     w1 = w3 = D4
00189     w2 = D1-w1-w3
00190 
00191     o2 = Rational(1,2)
00192 
00193     # outer
00194     p1 = [0,D2*o2]
00195     a = subx(p1,D1*o2)
00196     b = addx(p1,D1*o2)
00197     c = suby(b,D2)
00198     d = subx(c,D1)
00199 
00200     # inner
00201     e = suby(p1,D3)
00202     f = addx(e,w2)
00203     g = suby(f,h2)
00204     h = subx(g,w2)
00205 
00206     sections = [ [ 1,a,b],[ 1,b,c],[ 1,c,d],[ 1,d,a],
00207                  [-1,e,f],[-1,f,g],[-1,g,h],[-1,h,e] ]
00208     Ixx(sections)
00209     
00210 def Box1():
00211     """
00212           y
00213           ^
00214           |
00215       w1  | w2  w3
00216           |
00217     a-----p1------b
00218     |     |       |  h1
00219     |  e--|----f  |
00220     |  |  |    |  |
00221     |  |  |    |  |  h2
00222     |  |  o-------|--> x
00223     |  |       |  |
00224     |  h-------g  |
00225     |             |  h3
00226     d-------------c
00227 
00228     # given
00229     #a-b = D1
00230     #a-d = D2
00231 
00232     #(b-f)_y = D3
00233     #(g-c)_y = D4
00234     #(g-c)_x = D5
00235     #(d-h)_x = D6
00236     
00237     @warning inertia is wrong b/c 0 is assumed to be <0,0> when it should
00238     be the output of the first inertia calculation
00239     """
00240     print "---Box1---"
00241 
00242     D1 = Symbol('D1')
00243     D2 = Symbol('D2')
00244     D3 = Symbol('D3')
00245     D4 = Symbol('D4')
00246     D5 = Symbol('D5')
00247     D6 = Symbol('D6')
00248     
00249     # height
00250     h1 = D3
00251     h2 = D2-D3-D4
00252     h3 = D4
00253     
00254     # bases
00255     w1 = D6
00256     w2 = D1-D5-D6
00257     w3 = D5
00258 
00259     o2 = Rational(1,2)
00260 
00261     # outer
00262     p1 = [0,D2*o2]
00263     a = subx(p1,D1*o2)
00264     b = addx(p1,D1*o2)
00265     c = suby(b,D2)
00266     d = subx(c,D1)
00267 
00268     # inner
00269     e = suby(p1,D3)
00270     f = addx(e,w2)
00271     g = suby(f,h2)
00272     h = subx(g,w2)
00273 
00274     sections = [ [ 1,a,b],[ 1,b,c],[ 1,c,d],[ 1,d,a],
00275                  [-1,e,f],[-1,f,g],[-1,g,h],[-1,h,e] ]
00276     Ixx(sections)
00277     
00278 
00279 def T():
00280     """
00281            y
00282            ^
00283            |
00284     a---k-p1-l---b
00285     |      |     |  h1
00286     |      |     |
00287     h---g  | d---c
00288         |  | |
00289         |  o----------> x
00290         |    |
00291         |    |  h2
00292     j   f----e   i
00293       w1  w2   w3
00294 
00295     # given
00296     a-b     = D1
00297     (b-e)_y = D2
00298     a-h     = D3
00299     f-e     = D4
00300     @todo Area seems off
00301     @todo Inertias should be relative to the CG
00302     @warning wrong...
00303     """
00304     print "---T---"
00305 
00306     D1 = Symbol('D1')
00307     D2 = Symbol('D2')
00308     D3 = Symbol('D3')
00309     D4 = Symbol('D4')
00310 
00311     o2 = Rational(1,2)
00312     w2 = D4
00313     h1 = D3
00314     h2 = D2-h1
00315     w1 = w3 = (D1-D4)*o2
00316 
00317     p1 = [0,0]
00318     a = subx(p1,D1*o2)
00319     b = addx(p1,D1*o2)
00320     c = suby(b,h1)
00321     d = subx(c,w3)
00322     e = suby(d,h2)
00323     f = subx(e,w2)
00324     g = addy(f,h2)
00325     h = subx(g,w1)
00326     A = addy(h,h1)
00327 
00328     i = suby(c,h2)
00329     j = subx(i,D1)
00330     k = subx(p1,D4*o2)
00331     l = addx(p1,D4*o2)
00332     #print "a = ",a
00333     #print "A = ",A
00334 
00335     sections = [ [ 1,a,b],[1,b,c],[1,d,c],[1,d,e],[1,f,e],
00336                  [ 1,f,g],[1,g,h],[1,h,a],
00337                  [-1,l,d],[-1,d,g],[-1,g,k] ]
00338     Ixx(sections)
00339 
00340 
00341 def Ixx(sections):
00342     """
00343     from http://en.wikipedia.org/wiki/Second_moment_of_area
00344     J_{xx} = \frac{1}{12} \sum_{i = 1}^{n-1} ( y_i^2 + y_i y_{i+1} + y_{i+1}^2 ) a_i
00345     J_{yy} = \frac{1}{12} \sum_{i = 1}^{n-1} ( x_i^2 + x_i x_{i+1} + x_{i+1}^2 ) a_i
00346     J_{xy} = \frac{1}{24} \sum_{i = 1}^{n-1} ( x_i y_{i+1} + 2 x_i y_i + 2 x_{i+1} y_{i+1} + x_{i+1} y_i ) a_i
00347     """
00348     h = Symbol('h')
00349     b = Symbol('b')
00350     Ixx = 0
00351     Iyy = 0
00352     Ixy = 0
00353     half = Rational(1,2)
00354     CG = [0,0]
00355     Area = 0
00356     o3 = Rational(1,3)
00357     for i,section in enumerate(sections):
00358         (sign,p1,p2) = section
00359         #print "p%s = %s" %(i,p1)
00360         #print "p%s = %s" %(i+1,p2)
00361         #Atri = half * abs(p1[0]*p2[1])
00362         AA = Matrix([p1,p2])
00363         #print "AA = \n",AA
00364         detAA = AA.det()
00365         msgAA = '%s'%(detAA)
00366         #print "msgAA = ",msgAA
00367         if '-' in msgAA:
00368             detAA *= Rational(-1,1)
00369         Atri = detAA
00370         #Atri = half*b*h
00371         #Atri = (p1[0]*p2[1] - p2[0]*p1[1])*half
00372         ixxi = p1[1]*p1[1] +   p1[1]*p2[1] +   p2[1]*p2[1]
00373         iyyi = p1[0]*p1[0] +   p1[0]*p2[0] +   p2[0]*p2[0]
00374         ixyi = p1[0]*p1[1] + 2*p1[0]*p1[1] + 2*p2[0]*p2[1] + p2[0]*p1[1]
00375         Ixxi = ixxi*Atri
00376         Iyyi = iyyi*Atri
00377         Ixyi = ixyi*Atri
00378         #print "Atri = ",Atri
00379         #print "Ixxi = ",Ixxi
00380         Ixx += Ixxi
00381         Iyy += Iyyi
00382         Ixy += Ixyi
00383         #print "Ixx = ",simplify(Ixx),'\n'
00384         CG = [CG[0]+p1[0]*o3+p2[0]*o3,CG[1]+p1[1]*o3+p2[1]*o3]
00385         Area += Atri
00386     CG = [CG[0]/Area,CG[1]/Area]
00387     #print "Ixx = ",simplify(Ixx)
00388     o12 = Rational(1,12)
00389     o24 = Rational(1,24)
00390     Ixx = Ixx*o12
00391     Iyy = Iyy*o12
00392     Ixy = Ixy*o24
00393     J = Ixx+Iyy
00394     print "Area = ",simplify(Area*half)
00395     print "CG   = ",simplify(CG)
00396     print "Ixx  = ",simplify(Ixx)
00397     print "Iyy  = ",simplify(Iyy)
00398     print "Ixy  = ",simplify(Ixy)
00399     print "J    = ",simplify(J),'\n'
00400     
00401 def add(p1,p2,s1=1,s2=1):
00402     x = p1[0]*s1+p2[0]*s2
00403     y = p1[1]*s1+p2[1]*s2
00404     return [simplify(x),simplify(y)]
00405 
00406 def sub(p1,p2):
00407     return add(p1,p2,s2=-1)
00408 
00409 def subx(p1,p2x):
00410     return sub(p1,[p2x,0])
00411 
00412 def suby(p1,p2y):
00413     return sub(p1,[0,p2y])
00414 
00415 def addx(p1,p2x):
00416     return add(p1,[p2x,0])
00417 
00418 def addy(p1,p2y):
00419     return add(p1,[0,p2y])
00420 
00421 #bar()
00422 #cross()
00423 #box()
00424 T()
00425 
 All Classes Namespaces Files Functions Variables