1968年Peter Glaser教授提出了在外空间向太阳“借”能量的有趣设想：通过建立空间太阳能电站向地面传输能量，这一构想有望在一定程度上解决日益严重的能源匮乏难题。受此启发，此后的四十多年间，各国科学家已经提出了多种空间太阳能电站设计的方案^[1]。2012年由多国科学家共同合作，提出了一种全新的空间太阳能电站的构想——相控阵空间太阳能电站，该方案在无须控制的聚光系统、模块化的结构设计等方面受到科学界、工程界的广泛关注。在SPS-ALPHA概念中，由于结构具有面质比高、所处环境恶劣等技术难点，为相关领域提供了极大的研究空间。就结构动力学分析与控制而言，本课题组已经开展了相关的探索性研究工作^[2]，初步验证了辛算法在模拟SPS-ALPHA局部结构自旋展开过程中的适用性。作为一个典型的近似保守系统，仅仅研究辛算法的适用性是远远不够的，还需要进一步研究辛算法在SPS-ALPHA接收装置——太阳帆骨架自旋展开过程中的保结构特性，用于佐证算法的有效性，这正是本文的出发点。

大型柔性结构构件的自旋展开过程是近年来动力学与控制领域的研究热点之一。为了提高自旋展开过程模拟的数值精度，胡海岩院士等采用广义α法模拟黏弹性薄膜太阳帆的自旋展开过程^[3]，并通过试验验证模拟结果及其广义α法的有效性^[4-5]，这一成果标志着太阳帆这类超大柔性结构动力学问题的研究进入了一个崭新时代。然而，广义α法虽然在提高计算精度和模拟速度方面具有一定的优势，但是从本质上讲，广义α法在处理约束问题过程中，采用了Index思想提高计算速度，因此它属于一类耗散数值算法，即：算法存在人工耗散，这一本质缺陷导致广义α法的分析结果并不具有保结构特性，也不具有良好的数值特性。

中国科学家冯康院士于1984年在“Differential Geometry and Differential Equations”会议上提出了“动力学微分方程是保守系统，其差分积分格式应该保辛”的观点^[6]。哈密尔顿力学为提高数值方法的稳定性和保结构性能提供了新的途径^[7-8]，正是在这样的研究背景下，本文将借助保结构思想研究太阳帆骨架结构自旋展开过程中动力学问题分析的数值稳定性和保结构特性。

本文针对SPS-ALPHA中太阳光接收装置——太阳帆骨架，在哈密尔顿体系下建立太阳帆骨架结构的动力学简化模型，利用辛(几何)算法研究其在自旋展开过程中的动力学特性及结构的保持规律，为太阳帆骨架展开方案的设计与论证提供参考依据。需要说明的是：本文的研究没有考虑太阳帆骨架展开过程中的空天复杂环境条件对展开过程的影响，原因是太阳帆骨架结构阻尼较小，由于没有分析太阳帆本身的展开过程，太阳光压影响亦可以忽略不计；同时，由于本文关注的是模拟骨架展开过程中的保结构特性和算法的数值稳定性，因此，没有考虑具体工程实际中关注的指向控制、振动控制等方面的量化指标。

1 太阳帆骨架结构动力学模型

根据SPS-ALPHA类似“鸡尾酒杯”造型的结构(如图 1所示)，并参照网架式可展开空间结构的展开过程^[5]，建立SPS-ALPHA中太阳帆骨架结构自旋展开过程的简化动力学模型。根据骨架结构的回转对称性，同时为了与文献[2]进行对比研究，本文采用不计质量的刚性杆上作用集中载荷的运动模型，并考虑与回转轴垂直的任一截面的运动情况，结构简化平面模型如图 2所示。质点A, B, C, D, E, F的质量均为m，质点1, 2, 3, 4, 5, 6的质量均为M。由骨架结构的周期对称性，可以得到

(1)

(2)

式中，l_1A, l_1B, …, l_6A, θ_1A, θ_1B, …, θ_6A分别为连接半径为R的外圆和半径为r的内圆上相邻质点之间的距离及夹角。

引入广义坐标

当系统以恒定的角速度绕O转动时, 建立与展开系统一同转动的极坐标系, 并忽略科氏力对系统模型的影响, 因此系统模型的动能、万有引力势能可以分别表示为:

(3)

(4)

式中，ω为系统转动角速度。

系统的Lagrange函数可以表示为

(5)

引入广义动量:

根据Lagrange函数^[8], 由展开, 推导得到系统模型广义动量的具体表达式:

(6)

对系统进行勒让德变换, 得到哈密尔顿函数表达式为

(7)

在自旋展开过程中, 结构模型满足如下的约束条件:

(8)

根据哈密尔顿变分原理并完成变分推导, 得到

(9)

式中，λ为Lagrange乘子。

将(9)式展开后, 即得到哈密尔顿正则方程, 如下所示:

(10)

2 辛Runge-Kutta的格式

在哈密尔顿系统下, 使用保结构的辛(几何)算法能够保持系统的固有特性, 例如随着积分时间的增加, 它的相流形能够保持相空间的面积、体积不变以及可以保持系统的能量、动量不变等固有性质^[9-10]。下面以经典Runge-Kutta方法为基础, 简要介绍辛Runge-Kutta方法所需满足的条件, s级经典Runge-Kutta方法的一般格式为:

(11)

式中，i, j=1, 2, …, s, c_i≥0, 。

冯康等人^[6]已经证明格式(11)是辛的, 当且仅当其系数满足如下的条件:

(12)

式中，i, j=1, 2, …, s。

由于辛Runge-Kutta方法程序模块化程度高, 稳定性好, 备受计算数学家的青睐。在(12)式中, 当系数a_ij, b_i取不同的值时, 可以得到不同的辛Runge-Kutta格式, 其中一种常用的2级4阶辛Runge-Kutta格式如下所示:

(13)

3 数值仿真及结果讨论

为了说明前一部分构造的哈密尔顿体系下辛Runge-Kutta方法的优良性能, 本部分将分别采用辛Runge-Kutta方法和经典Runge-Kutta方法对比分析太阳帆骨架结构简化模型的展开过程, 从位移约束、保结构等方面说明辛(几何)算法的优越性。

3.1 基于辛(几何)算法的结构分析

本节采用(13)式对非线性哈密尔顿方程(10)进行数值仿真。其中结构质量及旋转角速度分别为:M=m=1, ω=1, 初始半径R₀=1, r₀=0.8, 且均为无量纲量。当步长分别为h=0.001、h=0.01、h=0.1时, 得到质点广义位移分量如图 3所示。

从图 3可以看出:当步长分别为h=0.001、h=0.01、h=0.1时, 随着模拟过程的进行, 所得到的广义位移分量变化情况基本一致:2个广义位移分量各自的振幅、振动周期均呈现周期性的变化, 并且这种周期性分布不随时间的推移而变化, 这说明使用辛算法模拟太阳帆骨架展开过程可以保持系统的振动特性。图 4所示为利用辛Runge-Kutta方法计算得到的内外层质点的对偶相图。由图 4中可以看出使用辛Runge-Kutta方法得到质点的相空间面积均保持不变, 这也验证了在哈密尔顿系统下使用辛(几何)算法可以保持系统的相空间面积不变这一结论。

3.2 系统位移约束违约分析

为了验证辛(几何)算法对于微分代数方程约束违约问题的适用性, 本文分别采用辛Runge-Kutta方法、经典Runge-Kutta方法对位移约束违约的情况进行模拟分析。令

式中，i=0, 1, 2, …, n, Δf_i表示约束g(R, r)在第i步的数值计算结果与理论值的差值(g(R, r)的理论值为零)。2种方法得到的位移违约情况分别如图 5所示。

由图 5a)可以看出采用辛Runge-Kutta方法得到的位移约束违约值控制在[-5×10^-16, 7×10^-16]范围内, 数值计算结果与理论值的误差基本可以忽略; 由图 5b)可以看出, 采用经典Runge-Kutta方法得到的位移约束违约值在[-5×10^-2, 0]的范围之内, 约束违约情况严重, 方程的解已经趋向于发散。对比图 5, 充分说明了辛算法能够很好地克服系统约束违约的情况。

3.3 系统能量分析

忽略展开过程中外界复杂空天环境的影响，太阳帆骨架自旋展开过程中的动力学系统是一个近似的保守系统，即系统的总能量应该是守恒的，而辛(几何)方法的最大优点就是能够长时间保持系统的整体能量。图 6a)给出了利用辛Runge-Kutta方法计算得到的系统模型的相对能量误差的情况, 从图中可以看出利用辛Runge-Kutta方法计算得到的系统相对能量误差控制在10^-15量级，并且在自旋展开过程中一直保持该数量级。图 6b)给出的是利用经典Runge-Kutta方法计算得到的相对能量误差的情况，从图中可以看出利用经典Runge-Kutta方法计算得到的系统相对能量误差在10^-5量级，其相对能量误差比辛Runge-Kutta方法高10个数量级，因此即便是目前广泛采用的经典Runge-Kutta方法，仍然不能保持系统的能量守恒性质，在长时间的积分运算后将会失去系统总能量不变这一基本性质，从而使计算结果产生较大的偏差。对比图 6，充分说明了辛算法在保持系统能量守恒性质方面具有突出的优势。

4 结论

本文以任意相控阵空间太阳能电站SPS-ALPHA中太阳帆骨架结构的展开过程为背景，建立了简化的动力学模型，并基于变分原理将Lagrange方程导入到哈密尔顿体系下，得到系统的正则方程组，并利用辛Runge-Kutta离散格式模拟分析了太阳帆骨架结构的展开过程，得到的主要结论如下：

1) 在数值模拟过程中，本文构造的2级4阶辛Runge-Kutta格式能够保持太阳帆骨架的振动特性及系统的相空间面积不变；

2) 使用辛Runge-Kutta方法得到的能量相对误差明显低于经典Runge-Kutta方法得到的能量相对误差，本文构造的辛Runge-Kutta格式在保持系统能量守恒方面具有突出优势，并且在数值计算中辛Runge-Kutta方法能够很好地克服系统约束违约的情况。