chaos_pendulum.py 2.2 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455
  1. import sympy.physics.mechanics as _me
  2. import sympy as _sm
  3. import math as m
  4. import numpy as _np
  5. g, lb, w, h = _sm.symbols('g lb w h', real=True)
  6. theta, phi, omega, alpha = _me.dynamicsymbols('theta phi omega alpha')
  7. theta_d, phi_d, omega_d, alpha_d = _me.dynamicsymbols('theta_ phi_ omega_ alpha_', 1)
  8. theta_dd, phi_dd = _me.dynamicsymbols('theta_ phi_', 2)
  9. frame_n = _me.ReferenceFrame('n')
  10. body_a_cm = _me.Point('a_cm')
  11. body_a_cm.set_vel(frame_n, 0)
  12. body_a_f = _me.ReferenceFrame('a_f')
  13. body_a = _me.RigidBody('a', body_a_cm, body_a_f, _sm.symbols('m'), (_me.outer(body_a_f.x,body_a_f.x),body_a_cm))
  14. body_b_cm = _me.Point('b_cm')
  15. body_b_cm.set_vel(frame_n, 0)
  16. body_b_f = _me.ReferenceFrame('b_f')
  17. body_b = _me.RigidBody('b', body_b_cm, body_b_f, _sm.symbols('m'), (_me.outer(body_b_f.x,body_b_f.x),body_b_cm))
  18. body_a_f.orient(frame_n, 'Axis', [theta, frame_n.y])
  19. body_b_f.orient(body_a_f, 'Axis', [phi, body_a_f.z])
  20. point_o = _me.Point('o')
  21. la = (lb-h/2)/2
  22. body_a_cm.set_pos(point_o, la*body_a_f.z)
  23. body_b_cm.set_pos(point_o, lb*body_a_f.z)
  24. body_a_f.set_ang_vel(frame_n, omega*frame_n.y)
  25. body_b_f.set_ang_vel(body_a_f, alpha*body_a_f.z)
  26. point_o.set_vel(frame_n, 0)
  27. body_a_cm.v2pt_theory(point_o,frame_n,body_a_f)
  28. body_b_cm.v2pt_theory(point_o,frame_n,body_a_f)
  29. ma = _sm.symbols('ma')
  30. body_a.mass = ma
  31. mb = _sm.symbols('mb')
  32. body_b.mass = mb
  33. iaxx = 1/12*ma*(2*la)**2
  34. iayy = iaxx
  35. iazz = 0
  36. ibxx = 1/12*mb*h**2
  37. ibyy = 1/12*mb*(w**2+h**2)
  38. ibzz = 1/12*mb*w**2
  39. body_a.inertia = (_me.inertia(body_a_f, iaxx, iayy, iazz, 0, 0, 0), body_a_cm)
  40. body_b.inertia = (_me.inertia(body_b_f, ibxx, ibyy, ibzz, 0, 0, 0), body_b_cm)
  41. force_a = body_a.mass*(g*frame_n.z)
  42. force_b = body_b.mass*(g*frame_n.z)
  43. kd_eqs = [theta_d - omega, phi_d - alpha]
  44. forceList = [(body_a.masscenter,body_a.mass*(g*frame_n.z)), (body_b.masscenter,body_b.mass*(g*frame_n.z))]
  45. kane = _me.KanesMethod(frame_n, q_ind=[theta,phi], u_ind=[omega, alpha], kd_eqs = kd_eqs)
  46. fr, frstar = kane.kanes_equations([body_a, body_b], forceList)
  47. zero = fr+frstar
  48. from pydy.system import System
  49. sys = System(kane, constants = {g:9.81, lb:0.2, w:0.2, h:0.1, ma:0.01, mb:0.1},
  50. specifieds={},
  51. initial_conditions={theta:_np.deg2rad(90), phi:_np.deg2rad(0.5), omega:0, alpha:0},
  52. times = _np.linspace(0.0, 10, 10/0.02))
  53. y=sys.integrate()