Point Mechanics in rigid body rotation. Vectors relationship.

[Sylwester L]

Hi

Sorry for my language but I dont speak english very well and probably I have trouble understanding comments.

Over a year ago Professor Jadczyk was interested me Dzanibekov's effect.

https://www.youtube.com/watch?v=BGRWg4aV2mw

Currently I can calculate and simulate a lot. I will use the equations because they are a universal language and they are understandable to all science on over the world.

Vector product.

a x b = c

Perpendicular axis

i ┴ j ┴ k

Proportions of vectors
c_k=a_i*b_j
a_i=c_k/b_j
b_j=c_k/a_i

The inverse of the vector.

a(a_x,a_y,a_z)=√(a_x²+a_y²+a_z²)
1/a=a/a²=(a_x/a², a_y/a², a_z/a²)

Example

d (1,1,1)    1/d (1/3, 1/3, 1/3)
e (1,1,0)    1/e (1/2, 1/2, 0)
f (1,0,0)    1/f (1,0,0)

Easy vector product equations.

c = a x b
a = 1/b x c
b = c x 1/a

Vector product equations for velocity, angular velocity and position. These equations are correct only for a free point.

v = ω x r = 1/s * m
ω = 1/r x v = 1/m * m/s
r = v x 1/ω = m/s * s

In rigid body rotation usually Angular velocity is not stable.

https://youtu.be/CAVXGDMbquk

How to calculate temporary angular velocity for rigid body rotation?

Rigid body elements

m1x`(√2, √2, 0);       v1(0,0,-2)
m2x`(-√2, -√2, 0);     v2(0,0,2)
m3y`(-√2, √2, 0);      v3(0,0,-2)
m4y`(√2, -√2, 0);      v4(0,0,2)

m1,m2 is x` main axis. m3,m4 is y` main axis. The center of mas is the center of the coordinate system.

We calculate angular velocity for main axis

ωx`=(-√2/2, √2/2, 0)
ωy`=(-√2/2, -√2/2, 0)

These are components temporary angular velocity for rigid body rotation

Ω = ωx` + ωy` =(√2, 0, 0)

https://m.salon24.pl/051c2bf7373b20e081ac0f25ababa170,750,0,0,0.jpg

How to calculate velocity?

Property vector product:

a x b = c ---> c=absinα

If

Ω (x,0,0) and r (x,y,0)

this give

ry┴Ω   i    rx║Ω
v = Ω x ry

Angular momentum for the free point

L = r x p

Angular momentum for rigid body show that equations using tensor moment of inertia.

L_x= ω_xΣm_n(r_n²-x_n²) + ω_yΣm_nx_ny_n + ω_zΣm_nx_nz_n
L_y= ω_xΣm_ny_nx_n + ω_yΣm_n(r_n²-y_n²) + ω_zΣm_ny_nz_n
L_z= ω_xΣm_nz_nx_n + ω_yΣm_nz_ny_n + ω_zΣm_n(r_n²-z_n²)

r is the position point to the center of mass.

Easier equations using moment of inertia for main axis

I_x=m_xr² + m_-xr²;     I_y=m_yr²=m_-yr² 

L_x = (r₁ x p₁)_x + (r₂ x p₂)_-x
L_y = (r₃ x p₃)_y + (r₄ x p₄)_-y
L_z = (r₅ x p₅)_z + (r₆ x p₆)_-z
L = L_x + L_y + L_z

Centripetal acceleration a_d for point.

a_d=rω²  --> rω=v
a_d=v²/r  --> v/r=ω
a_d=ωv

a_d║-r

Vector product equations for centripetal acceleration

a_d = ω x v
ω = 1/v x a_d
v =a_d x 1/ω

Equations temporary angular velocity for rigid body rotation

Ω=ωx`+ωy`+ωz`=(1/v x ad)x` + (1/v x ad)y` + (1/v x ad)z`

Centripetal forces  for point in rigid body rotation, three possibilities.

First version

F=am
ω = 1/v x a_d
ω = 1/mv x ma_d
ω = 1/p x F_d
F_d = ω x p
p = F_d x 1/ω

Second version

Ω = (x, 0, 0)
r = r_x +r_y = r┴ + r║
F_d=mω² r_y

Third version

Ω = (x, 0, 0)
r = rx +ry = r┴ + r║
Fd=mω2 r
Fd=m(1/ry x v)2 r

Everything shows my animation

https://youtu.be/oz1uw9x13kA

Another animation shows the angular acceleration vector for the effect

https://youtu.be/exwM5bTuO6Q

https://youtu.be/v2kwwzLA3aM

Another animation shows the vectors in no inertial frame

https://youtu.be/LZ9YwG9cVBE

code for the first simulation Vpython.

from visual import *

mx=0.5 #masy x`,y`

my=1.

x1=1 #pozycja m1,m2 na osi x`

y1=0

z1=0

x2=0 #pozycja m3,m4 na osi y`

y2=1

z2=0

r=vector(x1,y1,z1) #promien

R=mag(r)

#print R

W=vector(0.9,0,0) #omega

v1=W.x*y1 #predkosci

v2=W.x*-y1

v3=W.x*y2

v4=W.x*-y2

p1=mx*v1 #ped

p2=mx*v2

p3=my*v3

p4=my*v4

ax=(v1*v1) #przyspieszenia a=(v^2)/r; r=1

ay=(v3*v3)

#print "v",v1,v2,v3,v4

TIxx=(2*mx*((y1*y1)+(z1*z1)))+(2*my*((y2*y2)+(z2*z2))) #mx1 i mx2 --> 2*mx

TIyy=(2*mx*((x1*x1)+(z1*z1)))+(2*my*((x2*x2)+(z2*z2))) #elementy tensora

TIzz=(2*mx*((x1*x1)+(y1*y1)))+(2*my*((x2*x2)+(y2*y2)))

TIxy=(2*mx*x1*y1)+(2*my*x2*y2)

TIxz=(2*mx*x1*z1)+(2*my*x2*z2)

TIyx=(2*mx*y1*x1)+(2*my*y2*x2)

TIyz=(2*mx*y1*z1)+(2*my*y2*z2)

TIzx=(2*mx*z1*x1)+(2*my*z2*x2)

TIzy=(2*mx*z1*y1)+(2*my*z2*y2)

Fdx=mx*W.x*W.x*y1 #sily dosrodkowe na osiach glownych z Fd=mW^2ry

Fdy=my*W.x*W.x*y2

#print "F", Fd1,Fd2

Fdax=mx*ax #sily dosrodkowe na osiach glownych z F=ma

Fday=my*ay

#L=vector((W.x*TIxx-W.y*TIxy-W.z*TIxz),(-W.x*TIyx+W.y*TIyy-W.z*TIyz),(-W.x*TIzx-W.y*TIzy+W.z*TIzz))

#print "W",W

omega=arrow(axis=vector(W.x,W.y,0), color= color.blue, shaftwidth=0.05) #omega startowa

#kret=arrow(axis=vector(0,0,0), color= color.red, shaftwidth=0.04)

kret2=arrow(axis=vector(0,0,0),color=vector(1,0.4,0.4), shaftwidth=0.04)

#kret2x=arrow(axis=vector(0,0,0),color=vector(1,1,0.3), shaftwidth=0.04)

#kret2y=arrow(axis=vector(0,0,0),color=vector(1,1,0.3), shaftwidth=0.04)

dv1=arrow(axis=vector(0,0,0), color=vector(0.3,0.6,0), shaftwidth=0.05) #przyspieszenie punktow

dv2=arrow(axis=vector(0,0,0), color=vector(0.3,0.6,0), shaftwidth=0.05)

dv3=arrow(axis=vector(0,0,0), color=vector(0.3,0.6,0), shaftwidth=0.05)

dv4=arrow(axis=vector(0,0,0), color=vector(0.3,0.6,0), shaftwidth=0.05)

dvg=arrow(pos=vector(0,1,0),axis=vector(0,0,0), color=vector(0.6,0.6,0), shaftwidth=0.05)

dvd=arrow(pos=vector(0,-1,0),axis=vector(0,0,0), color=vector(0.6,0.6,0), shaftwidth=0.05)

#omegax=arrow(axis=vector(0,0,0), color= vector(0,0,0.01), shaftwidth=0.02)

#omegay=arrow(axis=vector(1,0,0), color= vector(0,0,0.01), shaftwidth=0.02)

masa1x=sphere(pos=vector(x1,y1,0),radius=0.05) #bryla sztywna

masa2x=sphere(pos=vector(-x1,-y1,0),radius=0.05)

masa1y=sphere(pos=vector(x2,y2,0),radius=0.05)

masa2y=sphere(pos=vector(-x2,-y2,0),radius=0.05)

promien1=arrow(pos=masa2x.pos, axis=masa1x.pos-masa2x.pos, color= color.yellow, shaftwidth=0.005)

promien2=arrow(pos=masa2y.pos, axis=masa1y.pos-masa2y.pos, color= color.yellow, shaftwidth=0.005)

vm1x=arrow(pos=masa1x.pos, axis=vector(0,0,v1), shaftwidth=0.01) #wektory predkosci punktow

vm2x=arrow(pos=masa2x.pos, axis=vector(0,0,v2), color=color.green, shaftwidth=0.01)

vm1y=arrow(pos=masa1y.pos, axis=vector(0,0,v3), color=color.green, shaftwidth=0.01)

vm2y=arrow(pos=masa2y.pos, axis=vector(0,0,v4), color=color.green, shaftwidth=0.01)

#orbita1=ring(pos=vector(x1,0,0), axis=vector(1,0,0), radius=y1, thickness=0.01) #orbita

os=arrow(pos=vector(-2,0,0), axis=vector(4,0,0),color=vector(0.3,0.3,0.3),shaftwidth=0.005) #os obrotu

sila1=arrow(pos=masa1x.pos,axis=-vector(x1,y1,0)*Fdx, color= vector(1,1,0), shaftwidth=0.05) #wektory sil dosrodkowych punktow

sila2=arrow(pos=masa1x.pos,axis=-vector(-x1,-y1,0)*Fdx, color= vector(1,1,0), shaftwidth=0.05)

sila3=arrow(pos=masa1y.pos,axis=-vector(x2,y2,0)*Fdy, color= vector(1,1,0), shaftwidth=0.05)

sila4=arrow(pos=masa1y.pos,axis=-vector(-x2,-y2,0)*Fdy, color= vector(1,1,0), shaftwidth=0.05)

sumasilag=arrow(color= vector(0.8,0.5,0), shaftwidth=0.05) #suma sil dosrodkowych gora z Fd=mW^2ry

sumasilad=arrow(color= vector(0.8,0.5,0), shaftwidth=0.05) #suma sil dosrodkowych dol z Fd=mW^2ry

#sila1a=arrow(pos=masa1x.pos,axis=-vector(x1,y1,0)*Fdax, color= vector(1,1,0.5), shaftwidth=0.05) #wektory sil dosrodkowych punktow

#sila2a=arrow(pos=masa1x.pos,axis=-vector(x1,y1,0)*Fdax, color= vector(1,1,0.5), shaftwidth=0.05) #wektory sil dosrodkowych punktow

#sila3a=arrow(pos=masa1x.pos,axis=-vector(x1,y1,0)*Fday, color= vector(1,1,0.5), shaftwidth=0.05) #wektory sil dosrodkowych punktow

#sila4a=arrow(pos=masa1x.pos,axis=-vector(x1,y1,0)*Fday, color= vector(1,1,0.5), shaftwidth=0.05) #wektory sil dosrodkowych punktow

#sumasilga=arrow(color= vector(0.8,0.5,0), shaftwidth=0.05) #suma sil dosrodkowych dol F=am

#sumasilda=arrow(color= vector(0.8,0.5,0), shaftwidth=0.05) #suma sil dosrodkowych dol F=am

#sumaF=arrow(axis=-vector(0,Fd1+Fd2,0), color= vector(1,0.5,0), shaftwidth=0.03)

t=0

while t<20:

rate(3)

# print t,L

# print TIxx-TIxy-TIxz,-TIyx+TIyy-TIyz,-TIzx-TIzy+TIzz

# print TIxx,TIxy,TIxz," I ",TIyx,TIyy,TIyz," I ",TIzx,TIzy,TIzz

x1=x1-0.1 #nowe pozycje punktow

y1=sqrt(1-(x1*x1))

y2=y2-0.1

x2=-sqrt(1-(y2*y2))

v1=W.x*y1 # v = W x ry

v2=W.x*-y1

v3=W.x*y2

v4=W.x*-y2

ov1=1/v1 #1/v

ov2=1/v2

ov3=1/v3

ov4=1/v4

# print t,v1*ov1

r1=vector(x1,y1,0) #promienie

r2=vector(-x1,-y1,0)

r3=vector(x2,y2,0)

r4=vector(-x1,-y1,0)

# print mag(r1),mag(r2),mag(r3),mag(r4)

p1=mx*v1 #ped

p2=mx*v2

p3=my*v3

p4=my*v4

op1=1/p1 #1/p

op2=1/p2

op3=1/p3

op4=1/p4

# print "v",v1,v2,v3,v4, "p",p1,p2,p3,p4

wx=vector(y1*v1,-(x1*v1),0) # wx` = r x v12; r=1

wy=vector(y2*v3,-(x2*v3),0) # wy` = r x v34; r=1

Wk=wx+wy # omega koncowa Wk=wx`+wy`

# print t, "Ws=", W, "Wk=", Wk

a1=(v1*v1)/mag(r1) #przyspieszenia a=(v^2)/r; r=1

a2=(v2*v2)/mag(r2)

a3=(v3*v3)/mag(r3)

a4=(v4*v4)/mag(r4)

a1v=vector(x1,y1,0)*-a1 #wektory przyspieszen

a2v=vector(-x1,-y1,0)*-a2

a3v=vector(x2,y2,0)*-a3

a4v=vector(-x2,-y2,0)*-a4

# a1r=-r1*(v1*v1)

# f1am=a1*mx

# print t, f1am

# print t,a1,a2

if t<10: #suma par przyspieszen dosrodkowych gora, dol z Fd=mW^2ry

adg=a1v+a3v

add=a2v+a4v

else:

adg=a1v+a4v

add=a2v+a3v

Wax=vector(-ov1*a1v.y,ov1*a1v.x,0) # w = 1/v x a

Way=vector(-ov3*a3v.y,ov3*a3v.x,0)

Wa=Wax+Way #W=wx`+wy`

# print t, Wa

# r1=sqrt((x1*x1)+(y1*y1))

TIxx=(2*mx*((y1*y1)+(z1*z1)))+(2*my*((y2*y2)+(z2*z2))) #mx1 i mx2 --> 2*mx

TIyy=(2*mx*((x1*x1)+(z1*z1)))+(2*my*((x2*x2)+(z2*z2)))

TIzz=(2*mx*((x1*x1)+(y1*y1)))+(2*my*((x2*x2)+(y2*y2)))

TIxy=(2*mx*x1*y1)+(2*my*x2*y2)

TIxz=(2*mx*x1*z1)+(2*my*x2*z2)

TIyx=(2*mx*y1*x1)+(2*my*y2*x2)

TIyz=(2*mx*y1*z1)+(2*my*y2*z2)

TIzx=(2*mx*z1*x1)+(2*my*z2*x2)

TIzy=(2*mx*z1*y1)+(2*my*z2*y2)

L=vector((W.x*TIxx-W.y*TIxy-W.z*TIxz),(-W.x*TIyx+W.y*TIyy-W.z*TIyz),(-W.x*TIzx-W.y*TIzy+W.z*TIzz))

if y1<0: #wartosc sily dosrodkowe Fd=mW^2ry

Fdx=mx*W.x*W.x*y1

else:

Fdx=-mx*W.x*W.x*y1

if y2<0:

Fdy=my*W.x*W.x*y2

else:

Fdy=-my*W.x*W.x*y2

# print t, Fdy, x2,y2

Lpr1=vector(y1*p1,-x1*p1,0) # L = r x p

Lpr2=vector(-y1*p2,x1*p2,0)

Lpr3=vector(y2*p3,-x2*p3,0)

Lpr4=vector(-y2*p4,x2*p4,0)

Lprx=Lpr1+Lpr2

Lpry=Lpr3+Lpr4

Lpr=Lprx+Lpry #suma kretow

# print t, L-Lpr

Fd1=vector(x1,y1,0)*Fdx #wektory sil dosrodkowch z Fd=mW^2ry

Fd2=vector(-x1,-y1,0)*Fdx

Fd3=vector(x2,y2,0)*Fdy

Fd4=vector(-x2,-y2,0)*Fdy

if t<10: #suma par sil dosrodkowych gora, dol z Fd=mW^2ry

Fdg=Fd1+Fd3

Fdd=Fd2+Fd4

else:

Fdg=Fd1+Fd4

Fdd=Fd2+Fd3

Fd1a=a1v*mx #F=am

Fd2a=a2v*mx

Fd3a=a3v*my

Fd4a=a4v*my

WFpx=vector(-op1*Fd1a.y,op1*Fd1a.x,0) #w = (1/p) x F

WFpy=vector(-op3*Fd3a.y,op3*Fd3a.x,0)

WFp=WFpx+WFpy

# print t, WFp, W

if t<10: #suma par sil dosrodkowych gora, dol z F=ma

Fdag=Fd1a+Fd3a

Fdad=Fd2a+Fd4a

else:

Fdag=Fd1a+Fd4a

Fdad=Fd2a+Fd3a

print t,Fdag, Fdad

# kret.axis=vector(L.x,L.y,L.z+0.01)

kret2.axis=vector(Lpr.x,Lpr.y,Lpr.z+0.01)

# kret2x.axis=vector(Lprx.x,Lprx.y,Lprx.z+0.01)

# kret2y.axis=vector(Lpry.x,Lpry.y,Lpry.z+0.01)

masa1x.pos=vector(x1,y1,0)

masa2x.pos=vector(-x1,-y1,0)

masa1y.pos=vector(x2,y2,0)

masa2y.pos=vector(-x2,-y2,0)

promien1.pos=masa2x.pos

promien1.axis=masa1x.pos-masa2x.pos

promien2.pos=masa2y.pos

promien2.axis=masa1y.pos-masa2y.pos

vm1x.pos=masa1x.pos

vm1x.axis=vector(0,0,v1)

vm2x.pos=masa2x.pos

vm2x.axis=vector(0,0,v2)

vm1y.pos=masa1y.pos

vm1y.axis=vector(0,0,v3)

vm2y.pos=masa2y.pos

vm2y.axis=vector(0,0,v4)

# orbita1.pos=vector(x1,0,0)

# orbita1.radius=y1

# omegax.axis=wx

# omegay.axis=wy

# dv1.pos=vector(x1,y1,0)

# dv1.axis=vector(a1v.x,a1v.y,a1v.z)

# dv2.pos=vector(-x1,-y1,0)

# dv2.axis=vector(a2v.x,a2v.y,a2v.z)

# dv3.pos=vector(x2,y2,0)

# dv3.axis=vector(a3v.x,a3v.y,a3v.z)

# dv4.pos=vector(-x2,-y2,0)

# dv4.axis=vector(a4v.x,a4v.y,a4v.z)

# dvg.axis=vector(adg.x,adg.y,adg.z)

# dvd.axis=vector(add.x,add.y,add.z)

sila1.pos=vector(x1,y1,0)

# sila1.axis=vector(0,-mag(Fd1),0)

sila1.axis=vector(Fd1.x,Fd1.y,Fd1.z)

sila2.pos=vector(-x1,-y1,0)

# sila2.axis=vector(0,mag(Fd2),0)

sila2.axis=vector(Fd2.x,Fd2.y,Fd2.z)

sila3.pos=vector(x2,y2,0)

sila3.axis=vector(Fd3.x,Fd3.y,Fd3.z)

sila4.pos=vector(-x2,-y2,0)

sila4.axis=vector(Fd4.x,Fd4.y,Fd4.z)

# if t<10:

# sila3.axis=vector(0,-mag(Fd3),0)

# sila4.axis=vector(0,mag(Fd4),0)

# else:

# sila3.axis=vector(0,mag(Fd3),0)

# sila4.axis=vector(0,-mag(Fd4),0)

sumasilag.axis=vector(Fdg.x,Fdg.y,Fdg.z)

sumasilad.axis=vector(Fdd.x,Fdd.y,Fdd.z)

# sila1a.pos=vector(x1,y1,0)

# sila1a.axis=vector(Fd1a.x,Fd1a.y,Fd1a.z)

# sila2a.pos=vector(-x1,-y1,0)

# sila2a.axis=vector(Fd2a.x,Fd2a.y,Fd2a.z)

# sila3a.pos=vector(x2,y2,0)

# sila3a.axis=vector(Fd3a.x,Fd3a.y,Fd3a.z)

# sila4a.pos=vector(-x2,-y2,0)

# sila4a.axis=vector(Fd4a.x,Fd4a.y,Fd4a.z)

# sumasilga.axis=vector(Fdag.x,Fdag.y,Fdag.z)

# sumasilda.axis=vector(Fdad.x,Fdad.y,Fdad.z)

# sumaF.axis=sila1.axis+sila2.axis

# LxF=sumaF.axis.x*L.x+sumaF.axis.y*L.y+sumaF.axis.z*L.z

# print t,"L x F",LxF

t=t+1

[Bruce Sherwood]

What is your question?

The purpose of this forum is to answer questions about VPython. We cannot help with math or debug a large program. 

Bruce

[Sylwester L]

I have established some dependencies between vectors this is first half a note, the second half shows the code in VPython which I use to visualize these vectors. I know that you have hundreds of similar codes but in my opinion this is another example how can use the VPython. I hope some professional will see this code and maybe show me how some things make better.

Regard

I have problem with vector product. Currently to do full action every time I must add or subtract six pairs elements of vectors. If I have to do it for three axes this gives action on eighteen pairs. I know it should be done to do a lot easier.