Intuicja też jest ważna ale ważne by jej nie wierzyć:)
Od dawna czułem w kościach że przyspieszenia kątowe ze wzorów Eulera które działają w nieinercjalnym układzie BS będą też działać w inercjalnym chociaż myślałem że będzie je trzeba trochę przekształcić.
Zanim jednak wziąłem się za przekształcenia które nie mam pojęcia jak zrobić, najpierw zrobiłem sobie proste symulacje. Umiejąc już symulować trajektorię końca wektora prędkości kątowej ω bardzo łatwo jest wyznaczyć wektor przyspieszenia kątowego
ɛ= (dω/dt) = (ω1-ω0)/dt
Czym się różni takie ɛ od tych ze wzorów Eulera?
ɛx = (Iyωyωz – Izωzωy) /Ix
ɛy = (Izωzωx – Ixωxωz) /Iy
ɛz = (Ixωxωy – Iyωyωx) /Iz
Postanowiłem sprawdzić i w mojej symulacji przeniosłem ɛ z układu nieinercjalnego do inercjalnego a następnie odjąłem uzyskane wektory obydwoma metodami, wynik? Różnice mniejsze niż 10^-14.
Zobaczcie sami.
Prosty wzór na epsilon w układzie inercjalnym
dL/dt= Ĩɛx+Ĩɛy+Ĩɛz=0
Ĩ - tensor momentu bezwładności.
Bez tensora już ani rusz dlatego też muszę go w końcu w pełni rozgryźć. Nikt nie chce pomóc nie szkodzi zrobię to sam, tylko że w ten sposób może to potrwać parę lat:)
Kod do symulacji
from visual import *
osx=vector(1,0,0) #uklad bryly
osy=vector(0,1,0)
osz=vector(0,0,1)
print osx, osy, osz
masax=2
masay=1
masaz=0
rx=1
ry=1
rz=1
#bryla=frame()
#masa11=sphere(frame=bryla, pos = vector(1,0,0), radius =0.1, color = color.red)
#masa12=sphere(frame=bryla, pos = vector(-1,0,0), radius =0.1, color = color.red)
#masa21=sphere(frame=bryla, pos = vector(0,1,0), radius =0.1, color = color.blue)
#masa22=sphere(frame=bryla, pos = vector(0,-1,0), radius =0.1, color = color.blue)
#ramiex=cylinder(frame=bryla, pos = masa11.pos, radius =0.01, axis=masa12.pos-masa11.pos)
#ramiey=cylinder(frame=bryla, pos = masa21.pos, radius =0.01, axis=masa22.pos-masa21.pos)
xy=0. #iloczyn sklalarny osi bryly
xz=0.
yz=0.
Ix=(masay*ry*ry)+(masaz*rz*rz) #moment bezwaladnosci
Iy=(masax*rx*rx)+(masaz*rz*rz)
Iz=(masay*ry*ry)+(masax*rx*rx)
TIxx=Ix #Tensor momentu bezwladnosci
TIyy=Iy
TIzz=Iz
TIxy=0
TIxz=0
TIyx=0
TIyz=0
TIzx=0
TIzy=0
print "momenty bezwladnosci", Ix,Iy,Iz
wx=0.0 #omega startow
wy=1.
wz=0.01
W=vector(wx,wy,wz)
wosx=vector(wx,0,0) #omega w ukladzie bryly
wosy=vector(0,wy,0)
wosz=vector(0,0,wz)
WBS=wosx+wosy+wosz
print "omega bryly", WBS
L=vector(Ix*wx,Iy*wy,Iz*wz)
Lpoczatek=L
print "kret", L
X=arrow(axis=vector(1,0,0), shaftwidth=0.01)
Y=arrow(axis=vector(0,1,0), shaftwidth=0.01)
Z=arrow(axis=vector(0,0,1), shaftwidth=0.01)
omega=arrow(axis=vector(wx,wy,wz), color= color.blue, shaftwidth=0.01)
epsilon=arrow(axis=vector(0,0,0), color= color.green, shaftwidth=0.01)
epsilon2=arrow(axis=vector(0,0,0), shaftwidth=0.01)
kret=arrow(axis=vector(L.x,L.y,L.z), color= color.red, shaftwidth=0.01)
BSx=arrow(axis=vector(1,0,0), color= (0,0.4,0), shaftwidth=0.01)
BSy=arrow(axis=vector(0,1,0), color= (0.2,0.2,0), shaftwidth=0.01)
BSz=arrow(axis=vector(0,0,1), color= (0,0.2,0.2), shaftwidth=0.01)
kropkaw=sphere(pos=vector(wx,wy,wz), radius=0.01, color= color.blue, make_trail=True)
kropkae=sphere(pos=vector(0,0,0), radius=0.01, color= color.green, make_trail=True)
kropkae2=sphere(pos=vector(0,0,0), radius=0.01, make_trail=True)
kropkaL=sphere(pos=vector(L.x,L.y,L.z), radius=0.01, color= color.red, make_trail=True)
kropkaBSx=sphere(pos=vector(0,0,0), radius=0.01, color= color.green, make_trail=True)
t=0
dt=0.01
while t<20000:
rate(200)
#print mag(WBS)*dt
osx = rotate(osx, mag(WBS)*dt, WBS) #obrot ukladu bryly
osy = rotate(osy, mag(WBS)*dt, WBS)
osz = rotate(osz, mag(WBS)*dt, WBS)
# bryla.rotate (angle=mag(WBS)*dt,axis=WBS)
#xy=(osx.x*osy.x)+(osx.y*osy.y)+(osx.z*osy.z) #iloczyn skalarny wektorow osi bryly
#xz=(osx.x*osz.x)+(osx.y*osz.y)+(osx.z*osz.z)
#yz=(osz.x*osy.x)+(osz.y*osy.y)+(osz.z*osy.z)
ex=((Iy-Iz)*wy*wz)/Ix #rowniania Eulera
ey=((Iz-Ix)*wz*wx)/Iy
ez=((Ix-Iy)*wx*wy)/Iz
wx=wx+ex*dt #wspolrzedne omegi
wy=wy+ey*dt
wz=wz+ez*dt
WBS0=wosx+wosy+wosz
wosx=vector(osx.x,osx.y,osx.z) #spolzedne omegi w ukladzie bryly
wosx.mag=wx
wosy=vector(osy.x,osy.y,osy.z)
wosy.mag=wy
wosz=vector(osz.x,osz.y,osz.z)
wosz.mag=wz
WBS=wosx+wosy+wosz
omega.axis=vector(WBS.x,WBS.y,WBS.z)
#print t, WBS
eosx=vector(osx.x,osx.y,osx.z) #spolzedne omegi w ukladzie bryly
eosx.mag=ex
eosy=vector(osy.x,osy.y,osy.z)
eosy.mag=ey
eosz=vector(osz.x,osz.y,osz.z)
eosz.mag=ez
EBS=eosx+eosy+eosz
epsilon2.axis=vector(EBS.x,EBS.y,EBS.z)
eBS=(WBS-WBS0)*(1/dt)
print t, eBS-EBS
TIxx=masax*(osx.y*osx.y+osx.z*osx.z)+masay*(osy.y*osy.y+osy.z*osy.z)+masaz*(osz.y*osz.y+osz.z*osz.z)
TIyy=masax*(osx.x*osx.x+osx.z*osx.z)+masay*(osy.x*osy.x+osy.z*osy.z)+masaz*(osz.x*osz.x+osz.z*osz.z)
TIzz=masay*(osy.y*osy.y+osy.x*osy.x)+masax*(osx.y*osx.y+osx.x*osx.x)+masaz*(osz.y*osz.y+osz.x*osz.x)
TIxy=masax*osx.x*osx.y+masay*osy.x*osy.y+masaz*osz.x*osz.y
TIxz=masax*osx.x*osx.z+masay*osy.x*osy.z+masaz*osz.x*osz.z
TIyx=masax*osx.x*osx.y+masay*osy.x*osy.y+masaz*osz.x*osz.y
TIyz=masax*osx.z*osx.y+masay*osy.z*osy.y+masaz*osz.z*osz.y
TIzx=masax*osx.x*osx.z+masay*osy.x*osy.z+masaz*osz.x*osz.z
TIzy=masax*osx.z*osx.y+masay*osy.z*osy.y+masaz*osz.z*osz.y
L=vector((WBS.x*TIxx-WBS.y*TIxy-WBS.z*TIxz),(-WBS.x*TIyx+WBS.y*TIyy-WBS.z*TIyz),(-WBS.x*TIzx-WBS.y*TIzy+WBS.z*TIzz))
# print t, Lpoczatek-L
kret.axis=vector(L.x,L.y,L.z)
epsilon.axis=vector(eBS.x,eBS.y,eBS.z)
kropkaBSx.pos=vector(osy.x,osy.y,osy.z)
kropkaL.pos=vector(L.x,L.y,L.z)
kropkaw.pos=vector(WBS.x,WBS.y,WBS.z)
kropkae.pos=vector(eBS.x,eBS.y,eBS.z)
kropkae2.pos=vector(EBS.x,EBS.y,EBS.z)
#BSy.axis=vector(osy.x,osy.y,osy.z)
#BSz.axis=vector(osz.x,osz.y,osz.z)
#BSx.axis=vector(osx.x,osx.y,osx.z)
# print t, osx, mag(W)
t=t+1
