山东大学学士学位论文
????? 用f1、f2、g1、g2分别表示作用于?1和?2上的外力和外力矩，注意g1是???相对于O点的力矩，g2是相对于M点的力矩（见图10）。此外fM和gM是在P点处?1对?2的作用力和力矩。应用式（52）、（53）有
p1?f1?fM p2?f2?fM
h1?v?p1?g1?gM?b?fM h2?(v?ω?b)?p2?g2?gM
(99) (100) (101) (102)
式（99）和（100）相加得
作用在系统上的总外力定义为
p?f (103)
f?f1?f2
(104)
此两刚体系统中总外力为系统自身重力。 式（101）、（102）相加得
其中
h?v?p?g (105)
g?g1?g2?b?f2
(106)
此式为作用在系统上的总外力矩。 3.2.3 标量运动方程
在航天器力学中一般采用三个坐标系Fi 、F1 、F2 。F1 、F2 是固连在刚体?1和?2
上的。而惯量矩cj和Jj在坐标系Fj（j=1,2）中是常量。下面将矢量和张量按其在不同坐标系中的分量表达为矩阵：
?p?c
f?v??F1?pc1ωhbg??F1????c??fv
??????c1ωhbg
?35
山东大学学士学位论文
?J1?J??F1?J1??TJ?F1
??c2ωMh2gg??F?c?2M2?ω?2Mh?g?22g?M? J?FJ?T22?2?F2
由式（56）得
p?F1?p?F1?(mv?c?ω?c2?ωM)?mv?c?ω?F1?cTT2F2?F2ωM?mv?c?ω?FT1?F2c?2ωM?mv?c?ω?A?12c2ωM?vx???czcy??0?c2zc2y???Mx??m??v?0y????cz0?c???x??x???A?c0?c???2x????vz?????cycx0???y?12?2z????z?????c2yc2x0??My?????Mz??由式（57）得
h?F1?h?F1?(c?v?J?ω?J12?ωM)?c?v?Jω?F1?J12?FT2ωM ?c?v?Jω?J12ωM? 0?czcy?????c???vx?c???x????Mx?z0x??J?????cycx0??vy??y?J??????vz????12??My??z?????Mz??J12部分在F1 ，部分在F2中，
J?F?T121?J12?F2 J21?F?T2?J21?F1
容易看出
JT12?J21
由式（96）得
h?2?c2A21v?J21ω?J2ωM
??0?c2zc2y???c2z0?c??A?vx?????x????????Mx?? 2x?21?vy??J21y??J2??My??c2yc2x0???v??z?????z?????Mz??动能为
36
(107)
(108)
(109)
(110)
山东大学学士学位论文
111T=mvTv?ωTJω?ωMTJ2ωM222?vTc?ω?vTA12c2?ωM?ωTJ12ωM1?m?vxvy?21???Mx?My?2?vx???x??v??1?????J???vz???y?2?xyz??y????vz????z???0?czcy???x???Mx?????????vvv??c?Mz?J0?cx???y??2?My??xyz??z??cycx??0???Mz????z????
????x?y?0?c2zc2y???Mx???Mx?????????vvv?A?c?z?J0?c12Myxyz122z2x???My??????? ??c2yc2x????0?Mz??Mz??????(111)
系统的惯量矩阵
?mE?M=?c???c2?A21?c?JJ21???A12c2?J12? J2??(112)
其中
?p??p??h????h2??;p=Mv
1T?vTMv
2(113) (114)
?v??v??ω????ωM??;?f?? f??g????g2?gM??(115)
系统的惯量矩阵M是对称、正定但是时变的。旋转矩阵A12,J,J12,J21都是时变的。如果把参考点取为系统的质心，那么c?0，但质心在?1中并不是固定的，从式（88）可知
c?c1?A12c2?m2b
(116)
因A12是时变的。从式（114）有
?T?p?v;?T?h?ω;?T?h2 ?ωM(117)
方程（103）、（105）和（102）的矩阵形式为
37
山东大学学士学位论文
p??ω?p?f
?0?????z???y?h??ω?h?v?p?g??z0?x?y??px??fx?????? ??x??py???fy?0??pz????fz????(118)
?0?=-??z???y???z0?x?y??hx??0????v??x??h?y??z?vz0vx0??hz???????vyvy??px??gx?????? ?vx??py???gy?0??pz????gz????(119)
????h2??h2?A21?v?ωb?A12c2???Aω?ωM??g2?gM???21
???v??0?????b????0?czcy?xzyx?????????????A21??v??0??bAc0?cx??y??12?zx? ?y??z???v????????cycx?0?bz?0?z??yx???????????0?h2zh2y?????x???Mx???g2x??gMx?????????????g???g???h2z0?h2x???A12???y??My???2y??My????g??g???h2yh2x0???z?????Mz?????2z??Mz??????(120)
由式（113）和（118）、（119）、（120）组成的方程可解出（p,h,h2）和（v,ω,ωM）,同时还要解旋转矩阵A1i、A21，它是?1和?2在空间的指向。注意这时M是时变的，所以要解这些方程只有在计算机上采用数值解法。
38
山东大学学士学位论文
第4章 MATLAB数值分析
图10 数值分析模型
建立如图所示模型（其中刚体1为长2m宽0.6m高0.8m的长方体，刚体2为长0.2m宽0.1m高0.4m的长方体，它们在点M处连接，密度为1000kg/m3）进行运用MATLAB进行数值分析。将参考坐标系F1的坐标原点O设在刚体1的质心C1上，将参考坐标系F2的坐标原点O设在刚体1与刚体2的连接点M上。在惯性坐标系中设置初始位移参量rx?1m,ry?1m,rz?1m，bx?0,by?1m,bz?0.4m。刚体1的初始转角为
??0.1rad,??0.1rad,??0.1rad，刚体2相对于刚体1的初始转角为
?M?0.1rad,?M?0.1rad,?M?0.1rad。刚体1上受到作用于其质心C1上的外力为
f1x?0,f1y?0,f1z?8N，且不受外力矩，刚体2上受到作用于点M上的外力为f2x?0,f2y?0,f2z?0.07N，且不受外力矩，点
M
处的力矩为
gMx?0.03N?m,gMy?0.05N?m,gMz?0.03N?m。初始时总动量为零，总动量矩为零，刚体2的动量矩为零。由于拉格朗日方程法编程较为困难所以此处选择牛顿-欧拉法编写MATLAB程序得到一下刚体姿态运动随时间变化图（时间t取为0到5秒，步长为
39