均匀锤子形刚体的贾尼别科夫效应.docx

建立坐标轴与锤子形刚体（如图1所示）Ω=x,y,zx1≤x≤x2,y2+z2≤R2∪{(x,y,z)|x1-3R≤x≤x1,-4.5R≤y≤4.5R,-1.5R≤z≤1.5R}惯量主轴重合的坐标系Ox'y'z'，设ρ是刚体的密度，经过计算得到下面绕三个主轴转动的转动惯量：Ix'=Ωρy2+z2dV=ρR42πx2-x1+1215R,Iy'=Ωρx2+z2dV=26ρR2x133+ρR2x233+πρR4x2-x14-9ρR2x1-3R3+2434ρR5Iz'=Ωρx2+y2dV=26ρR2x133+ρR2x233+πρR4x2-x14-9ρR2x1-3R3+21874ρR5图1刚体瞬时角速度矢在三个主轴上的分量为ωx',ωy',ωz'，由Euler动力学方程可得：dωx'dt=Iy'-Iz'Ix'ωy'ωz'dωy'dt=Iz'-Ix'Iy'ωx'ωz'dωz'dt=Ix'-Iy'Iz'ωx'ωy'我们接下来用z-y-z顺序的欧拉角（章动角θ，进动角ψ，自转角φ）来描述刚体的转动状态图2可得转动矩阵：M=cosψ-sinψ0sinψcosψ0001cosθ0sinθ010-sinθ0cosθcosφ-sinφ0sinφcosφ0001=cosθcosψcosφ-sinψsinφ-cosθcosψsinφ-sinψcosφsinθcosψcosθsinψcosφ+cosψsinφ-cosθsinψsinφ+cosψcosφsinθsinψ-sinθcosφsinθsinφcosθ这个矩阵用于对刚体进行旋转。

如图2所示，我们可将角速度矢ω分解：ω=dθdtnθ+dψdtnψ+dφdtnφ=ωx'i'+ωy'j'+ωz'k'i',j',k' 是坐标系Ox'y'z'的基令nψ=nψx'i'+nψy'j'+nψz'k'，在坐标系Ox'y'z'看来，有：nψx'nψy'nψz'=M-1001=MT001=-sinθcosφsinθsinφcosθnθ=sinφi'+cosφj'nφ=k'得：ωx'ωy'ωz'=sinφ-sinθcosφ0cosφsinθsinφ00cosθ1dθdtdψdtdφdt由Cramer法则得：dθdt=sinφωx'+cosφωy'dψdt=-cscθcosφωx'+cscθsinφωy'dφdt=cotθcosφωx'-cotθsinφωy+ωz'联立所有的微分方程：dωx'dt=Iy'-Iz'Ix'ωy'ωz',dωy'dt=Iz'-Ix'Iy'ωx'ωz',dωz'dt=Ix'-Iy'Iz'ωx'ωy',dθdt=sinφωx'+cosφωy',dψdt=-cscθcosφωx'+cscθsinφωy',dφdt=cotθcosφωx'-cotθsinφωy+ωz'这是个一阶非线性微分方程组，可用Runge-Kutta法求得z-y-z顺序的欧拉角的数值解以及角速度矢在坐标系Ox'y'z'下的坐标。

这里要注意的一点是，θ=0,π时会遇到奇点，实际上θ=0,π时已经无法区分ψ和φ了，我们可以通过调整初值来避免遇上奇点除了z-y-z顺序的欧拉角，我们当然可以使用其他二十三种顺序的欧拉角，笔者已经用习惯了使用z-y-z顺序的欧拉角这是用GeoGebra模拟的刚体出现贾尼别科夫现象的一个结果。