利用手机传感器测定刚体的惯性张量

1、利用手机传感器测定刚体的惯性张量摘要:利用手机陀螺仪传感器和phyphox软件,记录手机在空中自由旋转时的角速度曲线，根据质心系中的总角动量和机 械能守恒,可以确定手机惯性张量的相对值.在惯性张量的非对角元素为零的假设下，得到了对角元素相对值的解析表达式. 如果考虑非对角元素，则可以用梯度下降算法求解,将结果代入欧拉动力学方程，验证了其正确性.关键词:惯性张量;phyphox;陀螺仪;欧拉方程(4(4)(5)在研究转动问题时，刚体的惯性张量是一个重 要的参量.在大学物理实验中，惯性张量的测量一 般只限于刚体绕某固定轴的转动惯量.对于质量分 布情况未知的刚体，至少需要测量六个绕不同轴的 转动惯量

2、才可以得到完整的惯性张量.目前利用手机传感器测定刚体某个方向上的转 动惯量已是比较成熟的物理实验$，事实上，手机 陀螺仪传感器可以同时记录手机(作为一个刚体)3 个方向的角速度变化曲线，可利用这个便利性研究 刚体的三维转动.本文讨论一种测量手机惯性张量 (相对值)的测量方法.这个方法只需要以不同初始 角速度释放手机，利用phyphox软件记录手机在空 中下落并自由转动过程中的角速度变化曲线，然后 利用质心系中的转动动能守恒与总角动量守恒，就 可以计算出手机惯性张量各矩阵元的相对大小.这 个方法的一大优点是：由于转动动能与总角动量都 是标量，在数据处理中可以避开繁琐的坐标变换 计算.1实验原理1

3、.1简化惯性张量刚体的惯性张量通常表示为以下的3x3对阵 矩阵， TOC o 1-5 h z (II I 444E4C+ = IE4 1( 1)KC4Icy Icc y因为手机几何对称性好，质量分布大致均匀，不妨先 仅假设+的非对角元素为零，即(I00、+/ = 0I0(2)、00L这里为方便起见，已将I.4、妇，、八分别简记为I.、 I、八在自由落体运动过程中，如果忽略空气阻力，手 机在质心系中满足转动动能守恒与总角动量守 恒.设phyphox软件在一次自由落体运动过程中 记录到的角速度矢量为( 3)其中1 = 1，2,3, N为数据点编号.质心系中的 转动动能和总角动量的平方可以表示为)

4、= y#1 Mg =y(I4+Ie)X2=#c MM#1 = I$&4 +I$&2+I2 &2需要找到合适的I.、I,、Ic，使得按上式计算的)和X2 尽量保持不变.注意只能得到I.、I,、Ic的相对值，因 此可以加上约束Tr( M/2) = I2+I2 +I2 = 1.先考虑转动动能守恒.可把极小化的目标函数 定为N+Ke+匕&$,1-2) 21= !容易证明，当该目标函数取极小值时,厂！、/ T 2. T 2. T 2 容易证明，当该目标函数取极小值时,)=崩(K4&4,1+Ke&e，1+Kc&C,1) 2N 1 = !这样目标函数可化简为 ( K,i +九&；,i +匕,1) 2 1=

5、! TOC o 1-5 h z 1奇(/1 +/灰,1 +11) 2 N1= 1(K(K.，4，匕)8 0，/E，/C) T 其中& 4,1 & 4,1& 4,1 & E，14 _ V 22222 一 &E，1 & 4,1&y,1 &e，11=12222、& C,1 & 4,1& C,1 & E，122& 4,1 &C,122& y22& 4,1 &C,122& y,1 & c,122& c,1 & C,1 /N 3； 32 32 3； 3 3、3232 32 32 32 32, 1N3；=葛& B1B N 1=1 B，1(6)(7)(6)(7)(8)(9)(10)Y/ =也,K) A 也,1

6、； ,E)T (11) 如果没有测量误差,Y)和Yx的最小值都可以 取到零，这意味着向量以4，匕，匕)T和M，/;，K)T都 是矩阵A的特征向量，而且对应的特征值等于零. 对于大多数像手机这样的刚体,K.、匕、匕不会完全相 同,0，/7，匕)T和(k2，/; J； ) T线性不相关,因此矩阵 A必然存在两重0特征值.对于实测数据，由于空气 阻力和测量误差等因素，矩阵A只存在两个接近于 0的特征值，设它们对应的特征向量分别为(511鼻；, 513 ) T和(ezne； ) T &因此可以要求0，弓，匕)和 (E J；，K ) T 都位于(511 ,512,5U ) T 和(5；1 ,522,52

7、3 ) T 形成的子空间中，即J( K4，/e，匕)=(511 ,512，513)= 0 ( 521，522，523)1( K4，Ky，KC)= ( 511，512，513)= 0 ( 521，522，523)(12) 等价地,0，匕，匕)T和(K； J； J； ) T都与矩阵A的另 一个特征向量(531 , 532,533 )垂直，(511 + 0521)2531 + ( 512 + 0522)2 5 32 +(513 + 0523) 2533 = 0( 13)由此得到关于,0的齐次二次方程，令5 = 0,得到P52 + 2gy = z = 0( 14)其中222P = 521531 + 5

8、22 5 32 + 523 5 33，g = 511 521 531 + 512 522 532 + 513 523 533 ,222Z = 511531 + 512 5 32 + 513 5 33(5至多有两个实数解.此外，由于三个转动惯量平方 和为1,如果用到的3个特征向量已经归一化，则要 求a2 + 02 - 1.考虑到K.、L、L 0,上面的方程组至 多有两个物理上合理的解，容易根据实际情况选择 正确的解.需要指出的是：由于矩阵8的形式并不依 赖于转动动能与总角动量的大小，为了减小测量的 随机误差，可以记录多次落体运动，先计算每条轨道 测量的8，然后把结果相加，用所得矩阵进行后续的 计

9、算.1.2 完整惯性张量如果考虑惯性张量的非对角元,1时刻手机的 转动动能和总角动量平方可以表示为) =刁# T+# 1 = 5Th(1 # T)22( 16)X2 = #T +# 1 = Tr( MM# 1 #T)假设优化目标函数为F =挫柘1( ；Tr(M#1#T)-) +NS/( Tr( MM#1 #T) -X2) 2)( 17)其中S)、S/为权重系数，下面为简单起见，取S)= 2 ,S/ = 1 /2.类似前面的推导，先对)和X2求极值，目标函数 可以表示为v ( Tr M( # 1 #T - 3) ) 2 +2 1=1v ( TrMM( #1#T-3) ) 2( 18)2 1=1其

10、中1N$=n #1 #TN 1=1记#1#：-$ = G，则目标函数可写以表示为( 19)F=T 1(Tr(M31)2 = T 1( Tr(MM31( 19)( 20) 除了对称性MT= M,可要求M满足Tr( M2) = 1.注 意现在可以把不同轨道的数据合并，得到更全面的 目标函数，即每条轨道的数据独立计算出一个目标 函数，再把所有结果相加.为简单起见，下面的公式略去了对不同轨道的求和.目标函数的梯度可以表 示为 Tr( MMC1) - ( MC1+C1 M)1 = 1为使得Tr( M2) = 1,梯度方向应修正为 $F 4 =-Tr(Ml M$M$M(21)(22(21)(22)*个小的

11、正数*,以Mg = M*-万4的形式不断迭代直 至M收敛，得到最佳的惯性矩阵.2实验方法2.1坐标轴的确定打开phyphox软件,将手机分别绕长宽高转动, 监视陀螺仪角速度数据.可以发现，当手机绕短边 转动时,E轴有明显的数据变化，另外两个方向数据 几乎为0,因此判断手机长边的方向对应e轴.类似 地可确定4轴和c轴，结果如图1所示.图1 phyphox与质心系坐标2.2角速度测量打开phyphox软件中的陀螺仪模块，设置定时 运行,0.5 s后开始测量.运行后将手机从高处释放, 使其自由下落到柔软的平面上.为丰富测量数据, 每次尽可能给手机一个不同的初始角速度.最后截 取出手机在空中自由转动过

12、程中的数据用于后面的 分析.2.3阻尼修正按照2.1节中的计算原理，在Matlab中编写程 序进行计算简化的惯性矩阵.求出三个对角元后, 画出实验的动能和角动量随时间变化曲线.图2为 一条轨道的结果.相比于3个角速度曲线中存在的 振荡,动能和角动量几乎是不变的，这说明它们的确 近似是守恒的.不过也可以发现动能和角动量存在 逐渐变小的趋势，这说明在自由下落过程中，空气阻 力也对转动产生衰减.为此本文进行简单的修正:= #1 eak.根据10条轨道的结果，取a = 0.05.修正 后的结果参见图2.为避免过拟合，实验中对不同的 轨道使用相同a.不过在一定范围内，最终M对a 的选择并不敏感.此外，我

13、们注意到图2中有动能 和总角动量在一小段时间内增大的现象.这种现象 与衰减修正无关，而且在多条轨道中出现.这也许 是手机陀螺仪的系统误差，但也不排除是某种空气 动力学效应所导致.3实验结果3.1简化的惯性张量经过阻尼修正后，求得10组对角元素的相对 值，结果见表1.可以看出不同轨道得到的结果还是 一致性比较好的.将这十条轨道的8矩阵相加，用 所得矩阵进行计算综合值，结果可得? = 0.221, ?= 1.229.图2修正前后的能量，角动量平方随时间变化图表1惯量矩阵对角元素相对值实验次数10.63120.13890.763120.62840.14460.764430.62120.13820.7

14、71440.62560.14010.767550.62200.14690.769160.62630.13820.7673实验次数70.62990&13380&7650实验次数70.62990&13380&765080.63020&15070&761790&62030&15210&7695100&61700&15710&7711平均值0&62520.14410&7670标准差0&00490&00750&0034综合值0&62490&13860&7682续表计算得到将记录的18条角速度曲线进行阻尼修正，按照 2. 2中的方法，在Matlab中编写程序进行计算.可以将 上一步中计算得到的对角元作为初

15、始值，以加快收敛18组数据得出的惯性张量结果为(0.6265-0.0022 0.0022 )+= -0.00220.1515-0.0020、0.0022-0.0020 0.7647 /为评估各个矩阵元的误差大小，将这18组数据分为 三组，分别利用相同方法求解出相对惯量矩阵，结果 如下.(0.6293-0&00160.0028 )+(11 =-0&00160&1606-0.0011、0.0028-0&00110.7604 /( 0&6268-0&00560.0003 )+(2)=-0&00560&1221-0.0037、0.0003-0&00370.7695 /( 0&6237-0&00300.

16、0024 )+(3)=-0&00300&1480-0.0056、0.0024-0&00560.7676 /基于这3次结果，估计各个矩阵元的不确定度， 结果为(0.0016 0.0014 0.0010)( M) = 0.0014 0.0126 0.0017.0.0010 0.0017 0.0029,从以上结果可以看到，非对角元远小于对角元.说 明简化惯性张量仅考虑对角元素的假设合理.由此 方法得到的更为准确的对角元的比值关系为亍=0. 242, ： =1.221，与简化模型的结果很接近.3.3结果验证为验证结果的准确性，可以将手机近似认为质 量均匀分布的长方体刚体，根据测量得到的长度参 数,计算

17、三个对角元的相对值.手机的长宽高分别 为 155.3 mm，75.0 mm，7.5 mm，根据K，=二(E2)匕=二(42+C2)4i.4i.+c2 . 2E +c= 0.235ic =.2= 0.235ic =.2+E2 寸 +c21.230在上述测量、计算方法过程中，将角动量的平方 作为一个守恒量.事实上，角动量本身也是一个守 恒量,而欧拉动力学则描述了这个过程:/.&.+( Kc-IE) &= 0e &E +(I,-Ic)= 0( 23)& c +( 4-1. ) &. &e = 0注意角加速度只与惯性矩阵的相对大小有关.将测 量的惯性矩阵代入欧拉动力学方程的三个分量式 中，计算出角加速度曲线.同时，对测得的角速度- 时间曲线，取相邻的5个点进行二次多项式拟合，求图3角加速度-时间曲线1得中间点的斜率作为此时的角加速度，得到实测角 加速度曲线.图3和图4给出在两条轨道上得到的 两种