mechanizmos58

vpython_spinning_and_flipping_T-handle_5.py

Aug 13th, 2021 (edited)
258
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # 23_08_21
  2.  
  3. '''
  4. Originally typed in from from the video:
  5.  "Modeling a Spinning and Flipping T-Handle With Spring Forces and Python"
  6.   in Rhett Allain's Physics Explained Channel (c. May 2021)
  7.  
  8. tested on: Linux Mint 20.2, Geany 1.36, python 3.8.10,
  9. Web Epiphany, i5-7260U CPU @ 2.20GHz × 2, Iris Plus Graphics 640
  10. vpython version: ['7.6.1', 'jupyter']
  11.  
  12.  
  13. references to Classical Mechanics by John R. Taylor (2005)
  14.  
  15. p. 83
  16.  
  17.   Principle of Conservation of Momentum
  18.  
  19.   Ṗ = Fext
  20.  
  21.     𝐏 = ∑ mα vα   α = 1,..., N
  22.  
  23.  
  24. p. 344 The Centrifugal Force
  25.  
  26.   Fcf = m(𝛀 × r) × 𝛀
  27.  
  28. p. 348 The Coriolis Force
  29.  
  30.   Fcor = 2mV × 𝛀
  31.  
  32.  
  33. p. 394-400 Euler's Equations
  34.  
  35. p. 681 Continuum Mechanics
  36. p. 715 Relation between Stress and Strain: Hooke's Law
  37. p. 718 The Equation of Motion for an Elastic Solid
  38.  
  39.  
  40. An Introduction to the Coriolis Force
  41. by Henry M. Stommel & Dennis W. Moore (1989)
  42. "Particles that are not at rest in the rotating reference frame...
  43. find themselves acted upon by the Coriolis forces." (p.62)
  44.  
  45. Particles that deviate from the equilibrium spring position find
  46. themselves acted upon by the restoring spring forces.
  47.  
  48. #todo:
  49.  
  50. update omega!!
  51.  
  52. move the object with the Coriolis forces and only use the spring force
  53.      to hold the object together (and allow movement to measure velocity)
  54.  
  55. plot Euler's equation output against the object's movement
  56.  
  57. menu to select g.omega
  58.  
  59. '''
  60.  
  61. from vpython import *
  62. # vpython does this:
  63. #   from math import *
  64. #   from numpy import arange
  65. #   from random import random
  66. #   import colorsys
  67. #     use colorsys directly, no need for list_to_vec(L), use vec(*L)
  68.  
  69. import sys
  70.  
  71. def main():
  72.  
  73.     scene.width, scene.height = 1000, 800
  74.     scene.align = 'left'
  75.     scene.resizable = False
  76.     scene.autoscale = False
  77.     scene.autozoom = False
  78.     scene.range = 0.08
  79.  
  80.     PAUSED_AT_START = False
  81.  
  82.     TBAR = False
  83.  
  84.     TRAILS = False
  85.  
  86.     GRAPH = True
  87.     XMAX = 20
  88.  
  89.     RATE = 500
  90.  
  91.     # ~ SPRING_CONSTANT = 2000
  92.     SPRING_CONSTANT = 166700 # MAX
  93.  
  94.     DELTA_T = 0.0017
  95.  
  96.     #todo
  97.     OMEGA_LIST = ((3, 0, 0), (3, 0.2, 0), (0, 0.1, 3), (0, 3, 0.2))
  98.     # ~ g.omega = vec(3, 0, 0) # stable (unperturbed) rotation
  99.     g.omega = vec(3, 0.2, 0) # unstable (perturbed) rotation
  100.     # ~ g.omega = vec(0, 0.1, 3) # stable circular rotation
  101.     # ~ g.omega = vec(0, 3, 0.2) # stable with a wobble
  102.  
  103.  
  104.     # THE MASSES
  105.     posA = vec(0, 0.035, 0)
  106.     posB = vec(0, -0.035, 0)
  107.     posC = vec(0.05, 0, 0)
  108.     # ~ posC = vec(0.07, 0, 0)
  109.  
  110.     # m1 is the mass bottom of the T
  111.     m1 = 0.3
  112.     # m2 are the masses on the T-handle
  113.     m2 = 0.4
  114.  
  115.     # Center of Mass
  116.     CoM = (m2*posA + m2*posB + m1*posC) / (m2 + m2 + m1)
  117.     label(text='CoM', pos=vec(0,0,0), font='bold',
  118.           height=11, color=color.red, box=False, opacity=0)
  119.  
  120.  
  121.     retain = 17
  122.     radius = 0.0002
  123.     massA = sphere(pos=posA-CoM, radius=m2/60, color=color.green,
  124.                    make_trail=TRAILS, retain=retain, trail_radius=radius)
  125.     massB = sphere(pos=posB-CoM, radius=m2/60, color=color.red,
  126.                    make_trail=TRAILS, retain=retain, trail_radius=radius)
  127.     massC = sphere(pos=posC-CoM, radius=m1/60, color=color.blue,
  128.                    make_trail=TRAILS, retain=retain, trail_radius=radius)
  129.  
  130.     massA.m = m2
  131.     massB.m = m2
  132.     massC.m = m1
  133.  
  134.     # v = 𝝎  â¨¯ r   # p. 338
  135.     massA.v = cross(g.omega, massA.pos)
  136.     massB.v = cross(g.omega, massB.pos)
  137.     massC.v = cross(g.omega, massC.pos)
  138.  
  139.     # p = mv   # p. 14
  140.     massA.p = massA.m * massA.v
  141.     massB.p = massB.m * massB.v
  142.     massC.p = massC.m * massC.v
  143.  
  144.     # THE VECTORS
  145.     rAB = massB.pos - massA.pos
  146.     rBC = massC.pos - massB.pos
  147.     rCA = massA.pos - massC.pos
  148.  
  149.  
  150.     if TBAR:
  151.         posD = massA.pos + rAB/2.0
  152.         cyl1 = cylinder(pos=massA.pos, axis=rAB, radius=m2/60)
  153.         cyl2 = cylinder(pos=massC.pos, axis=posD-massC.pos, radius=m1/60)
  154.     else:
  155.         # MAKE THE SPRINGS
  156.         coils = 25
  157.         thick = 0.0005
  158.         radius = 0.0015
  159.         steel = vec(0.33, 0.33, 0.45)
  160.         springAB = helix(pos=massA.pos, axis=rAB, coils=coils,
  161.                          thickness=thick, radius=radius, color=steel)
  162.         springBC = helix(pos=massB.pos, axis=rBC, coils=coils,
  163.                          thickness=thick, radius=radius, color=steel)
  164.         springCA = helix(pos=massC.pos, axis=rCA, coils=coils,
  165.                          thickness=thick, radius=radius, color=steel)
  166.  
  167.     # "UN-STRETCHED" LENGTH OF THE SPRINGS (L₀)
  168.     rAB_0 = rAB.mag
  169.     rBC_0 = rBC.mag
  170.     rCA_0 = rCA.mag
  171.  
  172.     # BUILTIN COLORS:
  173.     # black, white, red, green, blue, yellow, cyan, magenta, orange, purple
  174.     if GRAPH:
  175.         tgraph = graph(xtitle='Time (s)', width=880, height=800,
  176.                        ytitle='',
  177.                        xmin=0, xmax=XMAX,
  178.                        #ymin=-8e-4, ymax=8e-4,
  179.                        fast=False, align='right')
  180.  
  181.         red_pen = gcurve(color=color.red, width=2)
  182.         blue_pen = gcurve(color=color.blue, width=2)
  183.         black_pen = gcurve(color=color.black, width=2)
  184.         sumofforce = gcurve(color=color.black, label='Conbined force', width=2)
  185.         gray_pen = gcurve(color=vec(0.5,0.5,0.5), width=2)
  186.         centrifugal_pen = gcurve(color=color.black, label='centrifugal force', width=2)
  187.         coriolis_pen = gcurve(color=color.red, label='Coriolis force', width=2)
  188.         spring_pen = gcurve(color=vec(0.6,0.6,1), label='Spring force', width=2)
  189.  
  190.         def reset_graph():
  191.             red_pen.delete()
  192.             blue_pen.delete()
  193.             black_pen.delete()
  194.             sumofforce.delete()
  195.             gray_pen.delete()
  196.             centrifugal_pen.delete()
  197.             coriolis_pen.delete()
  198.             spring_pen.delete()
  199.             g.t = 0
  200.         scene.append_to_caption(' '*25)
  201.         button(text="Reset Graph", bind=reset_graph)
  202.  
  203.  
  204.     def pause_run():
  205.         g.paused = not g.paused
  206.         pause.text='  RUN   ' if g.paused else 'PAUSE'
  207.     scene.append_to_caption(' '*10)
  208.     pause = button(bind=pause_run)
  209.     g.paused = not PAUSED_AT_START
  210.     pause_run()
  211.  
  212.     g.t = 0
  213.     dt = DELTA_T
  214.     k = SPRING_CONSTANT
  215.  
  216.     while True:
  217.  
  218.         if g.paused: continue
  219.  
  220.         rate(RATE)
  221.  
  222.         omega = g.omega
  223.  
  224.         # UPDATE VECTORS
  225.         rAB = massB.pos - massA.pos
  226.         rBC = massC.pos - massB.pos
  227.         rCA = massA.pos - massC.pos
  228.  
  229.         # UPDATE SPRING FORCE OF ONE MASS ON ANOTHER
  230.         # Fs = -k * (|L| - L₀) * L.hat
  231.         FAB = -k * (rAB.mag - rAB_0) * rAB.hat
  232.         FBC = -k * (rBC.mag - rBC_0) * rBC.hat
  233.         FCA = -k * (rCA.mag - rCA_0) * rCA.hat
  234.  
  235.         # UPDATE MOMENTUM OF THE MASSES
  236.         # p2 = p1 + Fnet * 𝛥t
  237.         massA.p += (FCA - FAB) * dt
  238.         massB.p += (FAB - FBC) * dt
  239.         massC.p += (FBC - FCA) * dt
  240.  
  241.         # UPDATE POSITIONS OF MASSES
  242.         # r2 = r1 + p2/m * 𝛥t
  243.         massA.pos += massA.p * dt / massA.m
  244.         massB.pos += massB.p * dt / massB.m
  245.         massC.pos += massC.p * dt / massC.m
  246.  
  247.         # UPDATE VELOCITIES
  248.         # v = 𝝎  â¨¯ r   # p. 338
  249.         massA.v = cross(omega, massA.pos)
  250.         massB.v = cross(omega, massB.pos)
  251.         massC.v = cross(omega, massC.pos)
  252.  
  253.         # UPDATE CORIOLIS FORCES
  254.         corA = Coriolis_force(massA.m, omega, massA.v)
  255.         corB = Coriolis_force(massB.m, omega, massB.v)
  256.         corC = Coriolis_force(massC.m, omega, massC.v)
  257.  
  258.         # UPDATE CENTRIFUGAL FORCES
  259.         cenA = centrifugal_force(massA.m, omega, massA.pos)
  260.         cenB = centrifugal_force(massB.m, omega, massB.pos)
  261.         cenC = centrifugal_force(massC.m, omega, massC.pos)
  262.  
  263.         sum_of_forcesA = cenA + corA
  264.  
  265.         if TBAR:
  266.             posD = massA.pos + rAB/2.0
  267.             cyl1.pos = massA.pos
  268.             cyl1.axis = rAB
  269.             cyl2.pos = massC.pos
  270.             cyl2.axis = posD-massC.pos
  271.         else:
  272.             # UPDATE SPRING POSITIONS AND LENGTHS
  273.             springAB.pos = massA.pos; springAB.axis = rAB
  274.             springBC.pos = massB.pos; springBC.axis = rBC
  275.             springCA.pos = massC.pos; springCA.axis = rCA
  276.  
  277.         # ℓ = r ⨯ p  # p.90
  278.         # LA = cross(massA.pos, massA.p)
  279.         # LB = cross(massB.pos, massB.p)
  280.         # LC = cross(massC.pos, massC.p)
  281.  
  282.         # black_pen, red_pen, blue_pen, gray_pen, magenta_pen,
  283.         # cyan_pen, orange_pen, purple_pen
  284.         # centrifugal_pen, coriolis_pen
  285.         if GRAPH:
  286.             # PLOT SELECTED STUFF
  287.             if g.t <= 20.01:
  288.                 # ~ black_pen.plot(g.t, LA.z)
  289.                 # ~ red_pen.plot(g.t, corA.z)
  290.                 # ~ black_pen.plot(g.t, cenA.z)
  291.                 spring_pen.plot(g.t, FAB.z)
  292.                 sumofforce.plot(g.t, sum_of_forcesA.z)
  293.                 coriolis_pen.plot(g.t, corA.z)
  294.  
  295.         g.t += dt
  296.  
  297.  
  298. def centrifugal_force(m, 𝛀, r):
  299.  
  300.     """Find the vector of the centrifugal force given
  301.    mass (m), angular velocity vector (𝛀), and position vector (r).
  302.  
  303.    Fcf = m(𝛀 ⨯ r) ⨯ 𝛀    # Taylor p.344  """
  304.  
  305.     return m * cross(cross(𝛀, r), 𝛀)
  306.  
  307.  
  308. def Coriolis_force(m, 𝛀, v):
  309.  
  310.     """Find the vector of the Coriolis force given
  311.    mass (m), angular velocity vector (𝛀), and velocity vector (v).
  312.  
  313.    Fcor = 2mṙ ⨯ 𝛀  = 2mv ⨯ 𝛀    # Taylor p.348  """
  314.  
  315.     return cross(2.0*m*v, 𝛀)
  316.  
  317.  
  318. class g: ...
  319.  
  320.  
  321. if __name__ == '__main__':
  322.     sys.exit(main())
  323.  
  324. '''
  325.    L = OMEGA_LIST #todo
  326.  
  327.    def M(m):
  328.        val = m.selected
  329.        if val == str(L[0]):
  330.            g.omega = vec(*L[0])
  331.        elif val == str(L[1]):
  332.            g.omega = vec(*L[1])
  333.        elif val == str(L[2]):
  334.            g.omega = vec(*L[2])
  335.        elif val == str(L[3]):
  336.            g.omega = vec(*L[3])
  337.        print(g.omega)
  338.        #todo RESET POSITIONS,...
  339.        g.paused = False
  340.    scene.append_to_caption(' '*10)
  341.    menu(choices=['Select g.omega',
  342.          str(L[0]), str(L[1]), str(L[2]), str(L[3])],
  343.          index=2, bind=M)
  344.  
  345. '''
  346.  
RAW Paste Data