大型航天器结构-控制-光学一体化建模与微振动响应快速计算方法

李建宏宇，庞贺伟

（1.北京卫星环境工程研究所 可靠性与环境工程技术重点实验室，北京 100094；2.中国空间技术研究院，北京 100094）

微振动是影响航天器有效载荷特性的重要因素[1]。 大型航天器微振动涉及到结构系统、控制系统、光学系统等。传统的集中在单系统的设计和评估方法已经不能满足实际工程的需求，建立结构-控制-光学一体化模型是目前比较有效的分析手段[2]。21 世纪初，国外在SIM[3]、JWST[4]、TPF[5]等高分辨率空间望远镜的研制中开发了IMOS[6]、DOCS[7-8]和IME[9]等一体化建模软件。国内一些学者对一体化建模方法进行了介绍和综述性的研究[10-11]，并应用到基于整星级的微振动仿真[12-13]。

目前，结构-控制-光学一体化建模方法主要应用于概念设计阶段，航天器结构系统相对简单，建立的一体化模型的状态空间维数较低，不存在数值求解困难的问题。如SIM 中一体化模型状态空间的维数为316[3]；JWST 中状态空间的维数为320[4]；国内的整星级一体化模型状态空间的维数达到 1×103量级[12-13]。在实际大型航天器详细设计阶段，由于结构复杂，且自由度数高，一体化模型的状态空间维数将达到1×104量级，其计算成本将增长到千倍以上。在建立一体化模型计算微振动响应时，需要考虑模型和数值算法的计算效率。

文中使用状态空间法和附加刚度法建立一体化模型，根据两种模型的特点分别选取合适的数值计算方法，并制定了微振动频率响应计算策略，极大地提高了微振动响应的计算效率。

1 大型航天器微振动的分系统模型

1.1 结构动力学系统模型

结构动力学模型是基于有限元模型建立的。结构动力学方程为：

由于航天器结构复杂，自由度数高，一般利用模态叠加法对结构进行频率响应分析。考虑前r 阶模态，其中r 由需要考察的频率范围决定。一般情况下，考察前n Hz 的频率响应需要至少提取模态频率小于2n Hz 的所有模态。

有限元模型由航天器结构设计单位提供，内容包括特征值信息和特定节点的特征向量信息。一般情况下，不会要求输出全部自由度上的特征向量信息。这是由于大型航天器结构有限元模型的自由度达到百万以上，导出全部自由度的特征向量信息存在着存储成本过大和数据存储和读取时间过长的问题。

以模态坐标及其一阶导数作为线性系统的状态量，则结构动力学的线性控制系统为：

线性控制系统的输入分为两部分：扰动力（力矩）和控制系统的控制力（力矩）。线性控制系统的输出分为两部分，分别是控制系统中传感器自由度上的位移（转动）和光学系统中光学结构自由度上的位移（平动和转动）。

1.2 控制系统模型

控制系统模型由姿态控制设计单位以线性控制系统的状态空间形式提供。该系统的输入为控制系统中传感器自由度上的位移（转动），输出为控制系统的控制力（力矩）。

其中输入 xs与结构动力学系统模型中的状态量的关系为：

控制系统在频域的传递函数为：

1.3 光学系统模型

光学系统模型由光学设计单位提供。在微振动的力学环境下，光学镜面之间的相对运动为微小运动，可以用光学系统的一阶差分计算出光学灵敏度矩阵。

式中： yo是光轴绕三个坐标轴的转角；为光学灵敏度矩阵； xo为光学系统局部模型中光学结构自由度上的位移（平动和转动）。

1.4 选择矩阵

控制模型和光学模型都是建立在局部模型中。与结构动力学模型建立一体化模型时，需要将局部的自由度扩展到整个结构的全局自由度上。

设一个自由度集合 P = { p1, p2, … , ps}。自由度 pj在局部模型中为第个自由度，在全局模型中为第个自由度。局部坐标系的总自由度数为d，全局坐标系的总自由度数为n。则选择矩阵为n ×d 矩阵，其中第行第列的元素为1，其他元素均为0。可以通过选择矩阵 βP将自由度集合P 在局部模型和全局模型之间相互转换。

在一体化模型中，选择矩阵与特征值矩阵是以βTΦ 的形式成对出现。根据选择矩阵的定义可得出，对于自由度集合P：

2 大型航天器结构-控制-光学一体化模型

2.1 通过状态空间法建立一体化模型

将结构动力学模型和控制系统模型整合，建立新的状态方程：

令一体化模型的状态量为动力学系统中的模态坐标及其一阶导数和控制系统中的状态量：

则有：

将光学系统的输出作为一体化模型的输出：

则一体化模型在频域的传递函数为：

光轴三个方向的转动欧拉角与扰动力（力矩）在频域的关系为：

这是通过状态空间法建立一体化模型的微振动响应计算公式。

尾矿的组成较为复杂，我国尾矿主要成分有石英、长石、石灰石、滑石、白云石、云母、高岭石、石榴石和绿泥石等。按照尾矿的主要组成成分，尾矿可分为以下5类：

2.2 通过附加刚度建立一体化模型

在控制系统模型中，通过传递函数计算控制系统的控制力（力矩）与控制系统中传感器自由度上的位移（转动）的关系为：

通过选择矩阵将局部模型中控制系统的控制力（力矩）与传感器自由度上的位移（转动）扩展到全局模型中。在全局模型中，控制力（力矩）与位移（转动）的关系为：

将控制力（力矩）作为附加刚度加入到动力学方程中，动力学方程变为：

使用模态叠加法，将动力学方程缩聚到Φ 构成的r 阶线性子空间后，转换到频域：

对于光学系统模型，光轴三个方向的转动欧拉角与模态坐标的关系为：

这是通过附加刚度法建立一体化模型的微振动响应计算公式。

2.3 微振动响应求解过程中的数值方法

对于两种建模方法计算微振动响应，主要的计算成本均为求解线性方程组。其基本形式为：

对于状态空间法建立的一体化模型， A = ( jω E-；对于附加刚度建立的一体化模型，。在两种方法中，矩阵A 不属于对称矩阵、三角矩阵、埃尔米特（Hermitian）矩阵、黑森贝格（Hessenberg）矩阵中的任意一种矩阵，可采用LU 矩阵分解法进行求解。

对于状态空间法建立的一体化模型，A 为稀疏矩阵，可以使用稀疏矩阵的LU 矩阵分解法进行进一步加速；对于附加刚度建立的一体化模型，A 由两部分组成，其中一部分为对角矩阵，可以采用迭代法进行求解。

以一个典型的大型航天器一体化模型为例。状态空间法建立的一体化模型中，方阵A 的维数为16 702，根据估算[14]，使用LU 分解的计算复杂度在1×1012量级；方阵A 中非零元素个数为1×106量级，根据估算[15]，使用稀疏矩阵LU 分解的计算复杂度在1×109量级。

附加刚度建立的一体化模型中，方阵A 的维数为8348，根据估算，使用LU 分解的计算复杂度在1×1011量级。使用迭代法进行求解线性方程组的步骤为：令；选取 x0=D 作为初始近似 值；迭 代 公 式 为；当迭代次数超过设置的最大迭代次数或者时，迭代结束。

迭代过程中，使用矩阵乘法的结合律，可以减少计算复杂度。迭代过程中的计算为矩阵相乘。经过估算，迭代一次的计算复杂度在1×106量级。迭代初始近似值的物理意义是无控制系统下系统的频率响应。由于控制系统仅在低频起到抑制振动的效果，而对高频振动影响不大。从物理角度分析，在低频段需要迭代次数较多，而在高频段需要迭代次数较少。

迭代收敛的必要条件是矩阵G 的任意一种范数小于1。在整个频段的计算中，当频率逼近0 或者频率逼近固有频率时，可能存在某些频率点使用该迭代公式无法收敛的情况。

选取了5 个典型的频率点对两种建模方法使用四种数值计算方法进行微振动响应计算。在典型的大型航天器微振动预示计算中，需要计算的频率点在1×103量级。以2800 频率点为例，预估四种数值计算方法在全频段的时间成本，结果见表1。

表1 典型频率点计算时间成本Tab.1 Time costs on typical frequency point calculation

通过表1 可以验证估算的计算复杂度的正确性。使用迭代法求解附加刚度建立的一体化模型微振动响应的计算效率远高于其他方法，但在某些频率点可能出现迭代发散的情况。使用稀疏矩阵LU 分解求解状态空间建立的一体化模型微振动响应的计算效率虽然不是最高，但仍然在可接受范围，且适用于所有频率点。

结合两种算法的特点，可以制定全频段微振动响应计算策略：首先使用迭代法求解附加刚度建立的一体化模型微振动响应，当在某些频率点，迭代法发散或者收敛困难时，使用稀疏矩阵LU 分解求解状态空间建立的一体化模型微振动响应。这样既可保证整个频段微振动响应的计算效率，又可避免特殊频率点求解失败。

3 数值算例

3.1 某大型航天器一体化建模与微振动预示

某大型航天器有限元模型由400 000+个节点和400 000+个单元组成，自由度数为2 000 000+。利用Nastran 软件对有限元模型进行模态分析，取前8348阶（小于400 Hz 的所有模态）特征值，并输出特征向量在扰动力、控制力、传感器、光学系统节点自由度上的分量。

航天器的主要扰源有CMG、制冷机等。扰动输入为150 个自由度上的扰动力和扰动力矩。扰动力和扰动力矩由谐波叠加组成。控制模型的控制输入为3个转动自由度，控制输出为3 个控制力矩。光学模型的输入为5 个节点共计30 个平动自由度和转动自由度，输出为光轴三个方向的转动欧拉角。

分别通过状态空间法和附加刚度法建立一体化模型。考察频率从0.02～200 Hz 的微振动频率响应。在0.02～20 Hz 的低频段，取1000 个频率点（间隔0.02 Hz）；在20～200 Hz 的中高频段取1800 个频率点（间隔0.1 Hz）。根据2.3 小节中的计算策略求解整个频段共计2800 个频率点的微振动频率响应。

从表1 可知，使用迭代法求解附加刚度法的一体化模型中，迭代法一次的计算成本是稀疏矩阵LU 分解法求解状态空间法计算成本的1/1000 左右。据此设置最大迭代次数设为100，设置迭代收敛的收敛精度 ε = 10-4x0，保证迭代法求解附加刚度法的一体化模型与稀疏矩阵LU 分解法求解状态空间法的相对误差在0.01%左右。

3.2 计算结果

低频段和中高频段的光轴绕X 轴转动角的频率响应如图1 所示，同时绘制了无控制系统下的微振动频率响应。在低频段，相比于无控制系统，有控制系统明显抑制了微振动响应；在中高频段，控制系统的影响减弱，有控制系统的频率响应和无控制的频率响应几乎重合。

图1 光轴绕X 轴转动角的频率响应Fig.1 Frequency response of rotation angle of optical axis around X axis

在整个频段共2800 个频率点中，共有6 个频率点使用迭代法求解附加刚度建立的一体化模型微振动响应存在困难，其中4 个频率点（0.02、0.04、0.06、0.10 Hz）接近0 Hz，2 个频率点（0.58 Hz 接近第8阶模态0.5807 Hz，1.12 Hz 接近第10 阶频率1.1163 Hz）逼近固有频率。总耗时10.35 s，其中使用迭代法计算2794 个频率点耗时7.20 s，使用稀疏LU 分解法计算6 个频率点耗时3.15 s。微振动的计算成本得到了很好的控制，使未来针对微振动响应的优化提供了保障。

迭代法的迭代次数如图2 所示。可以看出，在小于40 Hz 以下的中低频段，迭代次数较高；在大于40 Hz 以上的中高频段，迭代次数较低。这说明了迭代初始值在中高频段与最终迭代值相近，其物理意义为无控制系统的频率响应与有控制系统的频率响应几乎相等。即从数学角度和物理角度相互印证了控制系统在大型航天器微振动抑制中对中高频段的效果较弱的结论。

图2 数值计算迭代次数Fig.2 Number of iterations for numerical computation

4 结语

大型航天器由于结构复杂自由度数高，在建立结构-控制-光学一体化模型求解微振动响应时，需要考虑计算成本。文中使用状态空间法和附加刚度法建立一体化模型并推导出微振动响应计算公式，并根据计算公式的特性分别采用稀疏矩阵LU 分解法和迭代法两种数值计算方法。使用迭代法求解附加刚度法的一体化模型效率最高，但在特殊频率点存在迭代发散的现象，这时可以换用稀疏矩阵LU 分解法求解状态空间法的一体化模型。数值算例证明了方法的高效性。此方法是适用于实际工程中大型航天器微振动预示的高效算法。