连续系统和刚体

这篇文章主要从连续体到刚体，介绍一些物理量在姿态的计算和坐标变换的方法。

一、一些简单物理量计算

这里对一个物体分成了四个部分，分别给出质量、位置矢量、速度矢量。计算整个物体的平动能、旋转动能、动量、相对原点的角动量、相对质心的角动量。

对一个物体而言，其运动可以分为平动与旋转。平动只需要计算质心，然后用中学方法代入计算即可。唯一的区别就是这里采用的是矩阵计算。对于旋转，需要用物体的总量，减去平动的量，才能得到相应的旋转的量。具体可以看代码。

值得一提的是，对于角动量的计算，相对于原点需要用总量进行计算，而相对于质心因为只需要考虑自身旋转，故只需要用旋转部分的量进行计算。

%定义质量、位置矢量、速度矢量
m1=1; m2=1; m3=2; m4=2; %质量
M = m1+m2+m3+m4; %总质量
R1=[1 -1 2];R2=[-1 -3 2];R3=[2 -1 -1];R4=[3 -1 -2]; %位置矢量
dR1=[2 1 1]; dR2=[0 -1 1]; dR3=[3 2 -1]; dR4=[0 0 1]; &速度矢量
R = [R1; R2 ;R3 ;R4];
dR =[dR1;dR2;dR3;dR4];

%计算质心位置矢量、质心速度矢量
Rc=(R1*m1+R2*m2+R3*m3+R4*m4)/(m1+m2+m3+m4);
dRc=(dR1*m1+dR2*m2+dR3*m3+dR4*m4)/(m1+m2+m3+m4);

%计算物体旋转分量，定义速度矢量
r=R-Rc;
r1 =r(1,:);r2 =r(2,:);r3 =r(3,:);r4 =r(4,:);
dr = dR-dRc; 
dr1=dr(1,:); dr2=dr(2,:); dr3=dr(3,:); dr4=dr(4,:);

%计算
T_m= 1/2*M*dot(dRc,dRc); %平移动能
T_r= 1/2*(m1*(dot(dr1,dr1)+m2*dot(dr2,dr2)+m3*dot(dr3,dr3)+m4*dot(dr4,dr4))); %旋转动能
p= M*dRc; %线动量
H0= cross(R1, m1*dR1)+ cross(R2, m2*dR2)+ cross(R3, m3*dR3)+ cross(R4, m4*dR4); %关于原点的角动量
Hc= cross(r1, m1*dR1)+ cross(r2, m2*dR2)+ cross(r3, m3*dR3)+ cross(r4, m4*dR4); %关于质心的角动量

%输出
disp(T_m)
disp(T_r)
disp(p)
disp(H0)
disp(Hc)

二、惯性张量

通过引入惯性张量的概念，可以简化角动量的计算。如图，括号内即表示惯量。

但需要注意，所有计算要基于同一个框架。一般有body frame 和inertial frame，分别用B、N表示。如要计算，需要转换为一个框架下的物理量。而转换方法是之前文章提到的DCM。可以参考前面的转换方法。

这里给出一段代码，给出了B的惯量和N的角速度。因此要利用欧拉角转换为DCM，通过转换得到B的角速度，再利用上面的公式得到B的角动量。

format long g;
%给定body frame惯性张量、inertial frame角速度。通过欧拉角转DCM，转换角速度至body frame速率
%EA转DCM
EA=[-10,10,5];
phi = deg2rad(EA(1,1));
theta = deg2rad(EA(1,2));
psi = deg2rad(EA(1,3));
BN=angle2dcm( phi,theta,psi,'ZYX' );

%给定惯性张量
I_B=[10,1,-1;
     1,5,1;
     -1,1,8];

%角速度变换
omega_N=[0.01,-0.01,0.01];
omega_B=BN * omega_N';

%用惯性张量计算角动量
H=I_B*omega_B;
disp(H)
    

三、平行轴定理

如果想将质心的惯量转换为原点的惯量，可以利用下面公式进行转换。即原点的惯量等质心的惯量加上总质量乘以质心位置矢量的反对称矩阵再乘以其转置。

 同样的，这里需要在同一个框架下。用一道例题解释：

format long g;
%已知欧拉角、body frame下质心的惯性张量、inertial frame的p点位置矢量，求body frame下p点的惯性张量
%欧拉角转DCM变换p点位置向量为body frame
EA=[-10,10,5];
phi = deg2rad(EA(1,1));
theta = deg2rad(EA(1,2));
psi = deg2rad(EA(1,3));
BN=angle2dcm( phi,theta,psi,'ZYX' );
Rcp_N=[-0.5, 0.5, 0.25];
Rcp_B=BN*Rcp_N';

%将Rcp_B转换为反对称矩阵
Rcp_B_H=[0, -Rcp_B(3), Rcp_B(2);
         Rcp_B(3), 0, -Rcp_B(1);
         -Rcp_B(2), Rcp_B(1),0];

%计算p点惯性张量
M=12.5;
Ic_B=[10,1,-1;
      1,5,1;
      -1,1,8];
Ip_B=Ic_B+M.*Rcp_B_H*Rcp_B_H';
disp(Ip_B)

四、坐标变换

那有没有对不同框架的任意转换呢？有的，如公式：

如图，要将B的惯量转换为F，最重要的还是找到FB的旋转矩阵。

这里需要介绍一个主惯量的概念。因为利用这个方法可以进行一一对应无限变换，因此会有一个特殊框架下的惯量I：是对角矩阵。会有这个形式：

那如何求解这个问题呢？很简单，看等式左半边，实际上，C的行，实际上就是I_B的特征向量，即C=V' 

但这里有一点要注意，因为涉及左右手系的判断。方法很简单，最V1和V2做叉乘，看他等于V3还是等于-V3。如果等于V3，那C就这样；如果和V3相反，那C使用需要加一个负号

这里给出一道题目，已知一个惯量。

1、将其转换参考系

2、将主惯量算出来并且从大到小输出

3、将C根据主惯量从大到小输出

代码如下：

format long g
Ic_B=[10,1,-1;
      1,5,1;
      -1,1,8];
%1、已知B惯性张量，坐标变换为D
%将MRP转换为DCM
sigma=[0.1;0.2;0.3];
sigma1=[0, -sigma(3), sigma(2);
        sigma(3), 0, -sigma(1);
       -sigma(2), sigma(1), 0];
DB=eye(3)+(8*sigma1^2-4*(1-sigma'*sigma)*sigma1)/(1+sigma'*sigma)^2;
%坐标变换
Ic_D=DB*Ic_B*DB';
disp(Ic_D)

%2、已知惯性张量，求主惯性张量（对角）,再从大到小排列
[V, D]=eig(Ic_B);
if isequal(cross(V(:,1), V(:,2)),V(:,3))
    C=V';
else
    C=-V';
end
I=C*Ic_B*C';
x=[I(1,1),I(2,2),I(3,3)];
[m,n]=sort(x,'descend');
disp(m)

%3、主惯量大小向下排列C
C_sort=[C(n(1,1),:);C(n(1,2),:);C(n(1,3),:)];
if isequal(cross(C_sort(1,:),C_sort(2,:)),C_sort(1,:))
    C_0=C_sort;
else
    C_0=-C_sort;
end
disp(C_0)

 五、结语

这部分大学物理和之前提到的DCM等相关内容。主要通过例题学习思想与方法。希望有所帮助