双轴太阳帆板卫星在轨惯性特性与柔性模态频率计算方法与流程

1.本发明涉及一种航天器在轨惯性特性与柔性频率计算分析方法，特别是根据卫星中心体、挠性天线与太阳翼的质量惯量和柔性特性，以及附件安装的位姿关系计算系统质心、惯量以及挠性附件频率变化情况的方法。


背景技术：

2.大型空间挠性结构，如大型sar卫星、通讯卫星等，在轨呈现轻量化、变构型的特征。在轨运行过程中，因太阳帆板转动，引起星体质心、转动惯量、非约束模态频率发生变化，为星上姿态控制系统设计带来不确定因素。为解决上述问题，需要根据活动附件的变化情况，计算卫星在轨质量特性与柔性模态频率的变化情况。
3.公开号为cn106295196a，名称为一种带有旋转挠性太阳阵卫星的在轨模态计算方法的中国专利，针对上述问题，研究了单翼太阳帆板绕单轴转动条件下，卫星转动惯量与非约束模态频率的计算方法。
4.然而现有大型挠性卫星为增加太阳翼充电效率，会采用双轴太阳帆板驱动机构控制卫星太阳翼转动，由于卫星太阳帆板驱动机构存在b轴的转角，使卫星构型发生较大变化，同时引起卫星质心、转动惯量以及挠性振动频率等参数发生较大变化，动力学参数变化更为显著，增加了系统设计的不确定性。采用上述专利技术，无法对这种带有双轴太阳帆板卫星的在轨惯性特性与柔性频率进行计算分析。


技术实现要素：

5.本发明解决的技术问题是：克服现有技术的不足，综合考虑卫星质量特征信息与柔性特征信息以及星体变构型信息的数据交互规律，关注卫星柔性动力学分析与控制系统设计中重点分析参数，提供了一种带有双轴太阳帆板卫星在轨惯性特性与柔性频率的计算分析方法，根据卫星双轴太阳帆板驱动机构a、b轴转角的变化情况，以及卫星的惯性与柔性特征参数，计算卫星质心、转动惯量、非约束模态频率。同时，结合卫星挠性附件的阻尼，计算非约束模态复频率的变化情况，用于解决变构型卫星惯量特性参数与柔性特性参数计算分析的问题，为卫星柔性动力学分析与控制系统设计提供设计依据。
6.本发明的技术解决方案是：双轴太阳帆板卫星在轨惯性特性与柔性模态频率计算方法，包括卫星质心计算、卫星转动惯量计算、卫星附件非约束模态频率计算、以及卫星附件非约束模态复频率计算，其中：
7.卫星质心在整星机械坐标系下的位置计算方法如下：
8.9.其中：(
·
)
t
为矩阵转置，mc为卫星中心体质量，m
sa
为南太阳帆板质量，m
na
为北太阳帆板质量，m
at
为柔性天线质量，为卫星中心体质心在整星机械坐标系下的位置，为南太阳帆板安装坐标系的原点在整星机械坐标系下的位置，为北太阳帆板安装坐标系的原点在整星机械坐标系下的位置，c
ssj
为南太阳帆板安装坐标系相对于整星机械坐标系的转移矩阵，c
nsj
为北太阳帆板安装坐标系相对于整星机械坐标系的转移矩阵，为南太阳帆板质心在南太阳帆板安装坐标系下的位置，为北太阳帆板质心在北太阳帆板安装坐标系下的位置，c
atj
为柔性天线安装坐标系相对于整星机械坐标系的安装矩阵，为柔性天线质心在柔性天线安装坐标系下的位置，为柔性天线安装坐标系的原点在整星机械坐标系下的位置；
10.卫星相对于其质心的转动惯量计算方法如下：
[0011][0012]
其中为卫星中心体相对于其质心的转动惯量，为南太阳帆板相对于卫星质心的转动惯量，为北太阳帆板相对于卫星质心的转动惯量，为柔性天线相对于卫星质心的转动惯量；
[0013]
卫星附件非约束模态频率ω
sati
计算方法如下：
[0014][0015]
其中，i＝1、2、
…
、n
sat
，n
sat
为对矩阵进行特征值分解后特征值的个数，其中特征值大于0的记为λ
sati
，m
sat
为柔性卫星系统的质量矩阵，k
sat
为柔性卫星系统的刚度矩阵；
[0016]
卫星附件非约束模态复频率ω'
sati
计算方法如下：
[0017][0018]
其中，i等于1、2、
…
、n
sat
，对矩阵a
sat
进行特征值分解，得到按照从小到大排列的特征值依次为-n
sat1
±
jω
dsat1
、-n
sat2
±
jω
dsat2
、
…
、选取不等于0的特征值记为-n
sati
±
jω
dsati
，a
sat
为柔性卫星系统的动力矩阵。
[0019]
进一步的，所述的南太阳帆板安装坐标系相对于整星机械坐标系的转移矩阵c
ssj
为：
[0020][0021]
其中α
ss
、β
ss
为南太阳帆板的a轴和b轴转角。
[0022]
进一步的，所述的北太阳帆板安装坐标系相对于整星机械坐标系的转移矩阵c
nsj
为：
[0023][0024]
其中α
ns
、β
ns
为北太阳帆板的a轴和b轴转角。
[0025]
进一步的，所述的南太阳帆板相对于卫星质心的转动惯量为为为南太阳帆板相对于其质心的转动惯量，(
·
)
×
为由矢量(
·
)的分量构成的斜对称矩阵。
[0026]
进一步的，所述的北太阳帆板相对于卫星质心的转动惯量为为为北太阳帆板相对于其质心的转动惯量，(
·
)
×
为由矢量(
·
)的分量构成的斜对称矩阵。
[0027]
进一步的，所述的柔性天线相对于卫星质心的转动惯量为为为柔性天线相对于其质心的转动惯量，(
·
)
×
为由矢量(
·
)的分量构成的斜对称矩阵。
[0028]
进一步的，所述的柔性卫星系统的质量矩阵m
sat
为：
[0029][0030]
其中，n
ss
为南太阳帆板柔性模态的阶数，n
ns
为北太阳帆板柔性模态的阶数，n
at
为柔性天线柔性模态的阶数，为n
ss
×nss
阶的单位矩阵，为n
ns
×nns
阶的单位矩阵，为n
at
×nat
阶的单位矩阵，为n
ss
×nns
阶的零矩阵，为n
ss
×nat
阶的零矩阵，为n
ns
×nat
阶的零矩阵，
[0031][0032][0033]
和分别为南太阳帆板在南太阳帆板安装坐标系中相对于连接点的平动耦合系数和转动耦合系数，和分别为北太阳帆板在北太阳帆板安装坐标系中相对于连接点的平动耦合系数和转动耦合系数，和分别为柔性天线在柔性天线安装坐标系中相对于连接点的平动耦合系数和转动耦合系数。
[0034]
进一步的，所述的柔性卫星系统的刚度矩阵k
sat
为：
[0035][0036]
其中03×3为3
×
3阶的零矩阵，为3
×nss
阶的零矩阵，为3
×nns
阶的零矩阵，为3
×nat
阶的零矩阵，λ
ss
为n
ss
×nss
阶的南太阳帆板的频率矩阵，其对角线元素依次为2πω
ss1
、2πω
ss2
、
…
、其余元素为0，λ
ns
为n
ns
×nns
阶的北太阳帆板的频率矩阵，其对角线元素依次为2πω
ns1
、2πω
ns2
、
…
、其余元素为0，λ
at
为n
at
×nat
阶的北太阳帆板的频率矩阵，其对角线元素依次为2πω
at1
、2πω
at2
、
…
、其余元素为0，南太阳帆板各阶振动频率为ω
ss1
、ω
ss2
、
…
、北太阳帆板各阶振动频率为ω
ns1
、ω
ns2
、
…
、柔性天线各阶振动频率为ω
at1
、ω
at2
、
…
、
[0037]
进一步的，所述的柔性卫星系统的动力矩阵a
sat
为：
[0038][0039][0040]
其中，s
ss
为n
ss
×nss
阶的南太阳帆板的阻尼矩阵，其对角线元素依次为ξ
ss1
、ξ
ss2
、
…
、其余元素为0，s
ns
为n
ns
×nns
阶的北太阳帆板的阻尼矩阵，其对角线元素依次为ξ
ns1
、ξ
ns2
、
…
、其余元素为0，s
at
为n
at
×nat
阶的北太阳帆板的阻尼矩阵，其对角线元素依次为ξ
at1
、ξ
at2
、
…
、其余元素为0，南太阳帆板各阶阻尼为ξ
ss1
、ξ
ss2
、
…
、北太阳帆板各阶阻尼为ξ
ns1
、ξ
ns2
、
…
、柔性天线各阶振动阻尼为ξ
at1
、ξ
at2
、
…
、
[0041]
本发明与现有技术相比的优点在于：本发明根据卫星中心体、挠性天线与太阳帆板质量、质心位置、转动惯量、附件安装坐标系原点位置，计算卫星在不同转角下的质心和转动惯量的变化情况；再根据挠性附件的振动频率、阻尼系数、平动耦合系数与转动耦合系数计算卫星在不同转角下的非约束模态频率以及考虑附件阻尼特性的复频率。本发明主要用于卫星姿态控制系统设计过程中计算大型挠性卫星在轨惯量特性与柔性特性，输入参数为卫星设计中的常用参数，相关设计方法简便易行，具有较强的实用性。该方法有效解决了
变构型卫星惯量特性参数与柔性特性参数计算分析的难题，可以直接用于计算分析带有双轴太阳帆板卫星在轨惯性特性与柔性频率，是带有大型挠性附件变构型遥感卫星进行高精度高稳定度控制设计的重要一环，可以推广应用于其它有质心、转动惯量、挠性振动频率计算分析需求的大型挠性遥感卫星。
附图说明
[0042]
图1为双轴太阳帆板卫星组成与连接方式示意图；
[0043]
图2为本发明方法的计算原理图；
[0044]
图3为本发明实施例中卫星质心在整星机械坐标系下的相对位移(相对于帆板零位状态)示意图；
[0045]
图4为本发明实施例中卫星转动惯量的相对变化(相对于帆板零位状态)示意图；
[0046]
图5为本发明实施例中太阳帆板非约束模态频率示意图；
[0047]
图6为本发明实施例中太阳帆板非约束模态复频率示意图。
具体实施方式
[0048]
本发明涉及一种带双轴柔性太阳帆板和柔性天线的卫星，该卫星由卫星中心体(不含南太阳帆板、北太阳帆板和柔性天线，下同)、北太阳帆板、南太阳帆板、柔性天线组成。双轴太阳帆板卫星组成与各单体之间的连接方式如图1所示。
[0049]
图1中主要坐标系定义如下：
[0050]
(1)卫星机械坐标系ojxjyjzj：坐标原点oj位于星箭分离面对接环理论圆心；ojzj轴过坐标原点oj，垂直于星箭分离面，沿卫星的纵轴方向，正方向指向对地面方向；ojxj轴过坐标原点oj，位于星箭分离面内，与卫星东板理论法线方向平行，正方向与东板外法线方向一致，指向卫星正常飞行方向；ojyj轴过坐标原点oj，位于星箭分离面内，与ojxj轴、ojzj轴构成右手系。
[0051]
(2)卫星中心体坐标系ocxcyczc：坐标原点oc位于卫星中心体质心，坐标系的ocxc、ocyc、oczc轴依次与卫星机械坐标系的ojxj、ojyj、ojzj轴重合。
[0052]
(3)南太阳帆板安装坐标系osx
syszs
：坐标原点os位于南太阳帆板根部与卫星中心体连接处，南太阳帆板未转动时，坐标系的osxs、o
sys
、o
szs
轴依次与卫星机械坐标系的ojxj、ojyj、ojzj轴重合。双轴帆板转动时，南太阳帆板安装坐标系osx
syszs
在卫星机械坐标系ojxjyjzj的基础上按照21转序进行转动，先绕o
sys
轴转动α
ss
角，再绕osxs轴转动β
ss
角。其中，南帆板的a轴对应于o
sys
轴，b轴对应于osxs轴。
[0053]
(4)北太阳帆板安装坐标系onx
nynzn
，坐标原点on位于北太阳帆板根部与卫星中心体连接处，北太阳帆板未转动时，坐标系的onxn、o
nyn
、o
nzn
轴依次与卫星机械坐标系的ojxj、ojyj、ojzj轴重合。双轴帆板转动时，北太阳帆板安装坐标系onx
nynzn
在卫星机械坐标系ojxjyjzj的基础上按照21转序进行转动，先绕o
nyn
轴转动α
ns
角，再绕onxn轴转动β
ns
角。其中，北帆板的a轴对应于o
nyn
轴，b轴对应于onxn轴。
[0054]
(5)柔性天线安装坐标系o
at
x
atyatzat
，坐标原点o
at
位于柔性天线根部与卫星中心体的连接处，坐标系的o
at
x
at
、o
atyat
、o
atzat
轴依次与卫星机械坐标系的ojxj、ojyj、ojzj轴重合。
[0055]
(6)卫星本体坐标系obxbybzb，坐标原点ob位于卫星整星(含卫星中心体、南太阳帆板、北太阳帆板和柔性天线，下同)质心，坐标系的obxb、obyb、obzb轴依次与卫星机械坐标系的ojxj、ojyj、ojzj轴重合。
[0056]
上述卫星结构及其坐标系定义关系，具体可参考太阳翼不同转角的卫星在轨模态频率计算方法，航天器工程，2017-08；cn106295196a，一种带有旋转挠性太阳阵卫星的在轨模态计算方法；以及挠性卫星在轨非约束模态计算研究，宇航学报，2014-04。
[0057]
如图2所示，为本发明方法的计算原理示意图，需要提前获取的量包括：
[0058]
南太阳帆板的a轴和b轴转角α
ss
、β
ss
，单位rad
[0059]
北太阳帆板的a轴和b轴转角α
ns
、β
ns
，单位rad
[0060]
卫星中心体质量mc，单位kg
[0061]
南太阳帆板质量m
sa
，单位kg
[0062]
北太阳帆板质量m
na
，单位kg
[0063]
柔性天线质量m
at
，单位kg
[0064]
卫星中心体质心在整星机械坐标系下的位置单位m
[0065]
南太阳帆板安装坐标系的原点在整星机械坐标系下的位置单位m
[0066]
北太阳帆板安装坐标系的原点在整星机械坐标系下的位置单位m
[0067]
南太阳帆板质心在南太阳帆板安装坐标系下的位置单位m
[0068]
北太阳帆板质心在北太阳帆板安装坐标系下的位置单位m
[0069]
柔性天线安装坐标系的原点在整星机械坐标系下的位置单位m
[0070]
柔性天线安装坐标系相对于整星机械坐标系的安装矩阵c
atj
[0071]
柔性天线质心在柔性天线安装坐标系下的位置单位m
[0072]
卫星中心体相对于其质心的转动惯量单位kgm2[0073]
南太阳帆板相对于其质心的转动惯量单位kgm2[0074]
北太阳帆板相对于其质心的转动惯量单位kgm2[0075]
柔性天线相对于其质心的转动惯量单位kgm2[0076]
南太阳帆板在南太阳帆板安装坐标系中相对于连接点的平动耦合系数单位kg
1/2
和转动耦合系数单位kg
1/2m[0077]
北太阳帆板在北太阳帆板安装坐标系中相对于连接点的平动耦合系数单位kg
1/2
和转动耦合系数单位kg
1/2m[0078]
柔性天线在柔性天线安装坐标系中相对于连接点的平动耦合系数单位kg
1/2
和转动耦合系数单位kg
1/2m[0079]
南太阳帆板各阶振动频率ω
ss1
、ω
ss2
、
…
、单位hz
[0080]
北太阳帆板各阶振动频率ω
ns1
、ω
ns2
、
…
、单位hz
[0081]
柔性天线各阶振动频率ω
at1
、ω
at2
、
…
、单位hz
[0082]
南太阳帆板各阶阻尼ξ
ss1
、ξ
ss2
、
…
、
[0083]
北太阳帆板各阶阻尼ξ
ns1
、ξ
ns2
、
…
、
[0084]
柔性天线各阶振动阻尼ξ
at1
、ξ
at2
、
…
、
[0085]
一、卫星质心计算
[0086]
此处的卫星质心计算，是指计算卫星整星(含卫星中心体、南太阳帆板、北太阳帆板和柔性天线)质心在整星机械坐标系下的位置单位m，具体计算流程如下：
[0087]
(1)计算南太阳帆板安装坐标系相对于整星机械坐标系的转移矩阵c
ssj
，具体计算公式如下：
[0088][0089]
(2)计算北太阳帆板安装坐标系相对于整星机械坐标系的转移矩阵c
nsj
，具体计算公式如下：
[0090][0091]
(3)计算卫星质心在整星机械坐标系下的位置具体计算公式如下：
[0092][0093]
(
·
)
t
为矩阵转置。
[0094]
二、卫星转动惯量计算
[0095]
在第一部分的基础上，可以计算卫星整星相对于其质心的转动惯量单位kgm2，具体计算流程如下：
[0096]
(1)计算南太阳帆板相对于卫星质心的转动惯量单位kgm2，具体计算公式如下：
[0097][0098]
(
·
)
×
为由矢量(
·
)的分量构成的斜对称矩阵。
[0099]
(2)计算北太阳帆板相对于卫星质心的转动惯量单位kgm2，具体计算公式如下：
[0100][0101]
(3)计算柔性天线相对于卫星质心的转动惯量单位kgm2，具体计算公式如下：
[0102][0103]
(4)计算卫星相对于其质心的转动惯量具体计算公式如下：
[0104][0105]
三、卫星附件非约束模态频率计算
[0106]
计算由中心体、南太阳帆板、北太阳帆板、柔性天线构成的柔性卫星系统中的附件(含南太阳帆板、北太阳帆板、柔性天线)的非约束模态频率ω
sat1
、ω
sat2
、
…
、单位hz，具体计算流程如下：
[0107]
(1)计算本体坐标系中南太阳帆板相对于质心的转动耦合系数单位kg
1/2
m，具体计算公式如下：
[0108][0109]
(2)计算本体坐标系中北太阳帆板相对于质心的转动耦合系数单位kg
1/2
m，具体计算公式如下：
[0110][0111]
(3)计算本体坐标系中柔性天线相对于质心的平动耦合系数单位kg
1/2
和转动耦合系数单位kg
1/2
m，具体计算公式如下：
[0112][0113]
(4)计算柔性卫星系统的质量矩阵m
sat
和刚度矩阵k
sat
，具体计算公式如下：
[0114][0115][0116]
其中，n
ss
为南太阳帆板柔性模态的阶数，n
ns
为北太阳帆板柔性模态的阶数，n
at
为
柔性天线柔性模态的阶数；为n
ss
×nss
阶的单位矩阵、为n
ns
×nns
阶的单位矩阵、为n
at
×nat
阶的单位矩阵；03×3为3
×
3阶的零矩阵，为3
×nss
阶的零矩阵，为3
×nns
阶的零矩阵，为3
×nat
阶的零矩阵，为n
ss
×nns
阶的零矩阵，为n
ss
×nat
阶的零矩阵，为n
ns
×nat
阶的零矩阵；λ
ss
为n
ss
×nss
阶的南太阳帆板的频率矩阵，其对角线元素依次为2πω
ss1
、2πω
ss2
、
…
、其余元素为0；λ
ns
为n
ns
×nns
阶的北太阳帆板的频率矩阵，其对角线元素依次为2πω
ns1
、2πω
ns2
、
…
、其余元素为0；λ
at
为n
at
×nat
阶的北太阳帆板的频率矩阵，其对角线元素依次为2πω
at1
、2πω
at2
、
…
、其余元素为0。
[0117]
(5)对矩阵进行特征值分解，得到按照从小到大排列的特征值依次为λ
sat1
、λ
sat2
、
…
、选取大于0的特征值λ
sati
，按照如下公式计算卫星非约束模态频率：
[0118][0119]
其中，i＝1、2、
…
、n
sat
，n
sat
为柔性卫星系统特征值的个数，等于n
ss
+n
ns
+n
at
+3。
[0120]
四、卫星附件非约束模态复频率计算
[0121]
计算由卫星中心体、南太阳帆板、北太阳帆板、柔性天线构成的柔性卫星系统中的附件(含南太阳帆板、北太阳帆板、柔性天线)非约束模态复频率ωs′
at1
、ωs′
at2
、
…
、单位hz，具体计算流程如下：
[0122]
(1)计算柔性卫星系统的阻尼矩阵c
sat
，具体计算公式如下：
[0123][0124]
其中，s
ss
为n
ss
×nss
阶的南太阳帆板的阻尼矩阵，其对角线元素依次为ξ
ss1
、ξ
ss2
、
…
、其余元素为0；s
ns
为n
ns
×nns
阶的北太阳帆板的阻尼矩阵，其对角线元素依次为ξ
ns1
、ξ
ns2
、
…
、其余元素为0；s
at
为n
at
×nat
阶的北太阳帆板的阻尼矩阵，其对角线元素依次为ξ
at1
、ξ
at2
、
…
、其余元素为0。
[0125]
(2)计算系统的动力矩阵a
sat
，具体计算公式如下：
[0126][0127]
(3)对矩阵a
sat
进行特征值分解，得到按照从小到大排列的特征值依次为-n
sat1
±jω
dsat1
、-n
sat2
±
jω
dsat2
、
…
、选取不等于0的特征值-n
sati
±
jω
dsati
，按照如下公式计算卫星非约束模态复频率：
[0128][0129]
其中，i等于1、2、
…
、n
sat
。
[0130]
实施例
[0131]
为验证本发明的有效性，采用数学仿真对本发明方法进行验证。
[0132]
图3给出了卫星南、北太阳帆板的a轴由0
°
变化至360
°
，b轴依次为0
°
、20
°
、40
°
工况下，卫星整星质心在整星机械坐标系下相对于太阳帆板处于零位时卫星质心的变化情况。结果分析表明，卫星整星质心因帆板b轴的变化引起的质心位移量约为2.721
×
10-6
m(δy
fcj
)。
[0133]
图4给出了卫星南、北太阳帆板的a轴由0
°
变化至360
°
，b轴依次为0
°
、20
°
、40
°
工况下，卫星转动惯量相对于太阳帆板处于零位时卫星转动惯量的变化情况。结果分析表明，卫星整星转动惯量因帆板b轴的变化引起的惯性变化量约为2000kgm2(δi
zzsat
)。
[0134]
图5给出了卫星南、北太阳帆板的a轴由0
°
变化至360
°
，b轴依次为0
°
、20
°
、40
°
工况下，卫星太阳帆板非约束模态频率的变化情况。结果分析表明，卫星太阳帆板非约束模态频率因帆板b轴的变化引起的变化量约为0.0105hz(ω
sat4
)。
[0135]
图6给出了卫星南、北太阳帆板的a轴由0
°
变化至360
°
，b轴依次为0
°
、20
°
、40
°
工况下，卫星太阳帆板非约束模态复频率的变化情况。结果分析表明，卫星太阳帆板非约束模态复频率因帆板b轴的变化引起的变化量约为0.0118hz(ω'
sat4
)。
[0136]
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。