太阳帆航天器在绕地轨道中的热诱发振动

2019-12-09张军徽佟安武娜刘应华

航空学报 2019年11期

张军徽，佟安，武娜，刘应华

1.北方工业大学土木工程学院，北京 100144 2.清华大学航天航空学院，北京 100084

光对被照射物体的压力称为光压。依靠光压推进的太阳帆航天器无需消耗推进剂，可以实现传统航天器无法实现的各种非开普勒轨道，被认为是最可行的星际探测技术。2010年日本IKAROS号太阳帆[1]和2011年美国Nanosail-D号太阳帆[2]的相继成功发射把人类星际旅行的梦想照进了现实。绕地球飞行轨道是太阳帆进入星际轨道的起点，更是各种近地球任务太阳帆运行的主要场所。太阳帆在绕地轨道中将会周期性地进出地球阴影，经历热环境的剧烈变化，在热辐射冲击作用下，太阳帆这种大柔性空间结构的热致响应值得研究。

热致结构响应的研究最早可以追溯到1956年，Boley在热致结构响应分析中引入惯性项，提出了热诱发振动的概念[3]。热诱发振动影响航天器正常工作的事件屡有发生，其中，最为著名的是太阳翼热诱发振动造成的哈勃太空望远镜成像模糊畸变事故。Thornton和Kim考虑结构变形对热载荷的影响，建立了太阳翼的热-结构耦合动力学模型，数值再现了哈勃望远镜太阳翼在突加热流载荷作用下的弯扭耦合热诱发振动[4]。Gulick和Thornton研究了自旋稳定卫星轴向天线的热诱发振动，研究表明：当卫星的自旋速度等于天线的弯曲振动频率，卫星天线将发生不稳定的热诱发振动，即热颤振[5]。Song等研究了复合材料薄壁梁的热诱发振动稳定性问题[6]。Zhao等推导了考虑热辐射的复合材料板的有限元方程，为分析航天器柔性复合材料薄板附件的热诱发振动提供了工具[7]。Shen等基于绝对节点坐标方法推导了梁的热-结构耦合方程，得益于绝对节点坐标方法的优点，该方程不但能计算昼夜转换热冲击引起的梁热颤振，也能计算由于大转动引起的梁热诱发振动[8]。Zhang等在Graham和Thorton等工作的基础上，建立了悬臂梁在空间热辐射作用下的热颤振准则：若空间热流入射角大于梁自由端的准静态转角，则热颤振不会发生[9]。Li和Yan研究了蜂窝复合材料太阳翼在近地轨道的热诱发振动[10]。Su等在空间环境模拟器中对复杂桁架的热致响应进行了实验研究，成功观测到了空间桁架的热诱发振动[11]。Wang和Li研究了热载荷作用下索-梁结构的参激振动[12]。Shen和Hu研究了热诱发振动对自旋稳定航天器姿态的影响，数值结果表明热诱发振动引起了Ulysses卫星姿态角的跳变[13]。Liu和Pan建立了刚-柔-热耦合的多体动力学模型，考虑太阳热辐射和地球热辐射等的影响，研究了带有柔性太阳翼航天器在近地轨道的热诱发振动[14]。Chen等研究了预应力对复合材料薄板热诱发振动的影响[15]。Azadi等研究了近地轨道卫星的压电智能太阳帆板的热诱发振动[16]。Liu等研究了热结构耦合作用下柔性航天器的姿态控制问题[17]。Xue等提出了傅立叶温度单元方法，大大提高了薄壁结构非线性温度场的求解效率，并将该方法用于大型柔性空间结构的热-结构耦合分析中，促进了热致结构响应分析的工程应用[18-20]。

目前，国内外关于太阳帆的研究集中在总体设计[21-22]、轨道和姿态控制[23-27]、帆膜制备和展开[28-30]等方面。按照展开和维持帆面构型的方式不同，太阳帆总体设计可分为旋转展开无桅杆式和机械展开有桅杆式两大类，图1所示为NASA研制的方形有桅杆式太阳帆[2]。通常认为有桅杆式太阳帆具有展开可靠性高，姿态控制简单和消耗控制能量小等优点[31]。但是，有桅杆式太阳帆在绕地轨道热辐射冲击作用下，桅杆截面温差可能引起太阳帆结构的热致响应，从而影响太阳帆顺利进入星际轨道，或者影响近地球任务太阳帆的工作性能。目前，国内外关于有桅杆式太阳帆在绕地轨道中的热诱发振动问题未见报道。

图1 NASA研制的Nanosail-D太阳帆[2]Fig.1 Solar sail Nanosail-D designed by NASA[2]

本文以方形有桅杆式太阳帆为研究对象，考虑热辐射等非线性因素，建立太阳帆的热-结构耦合动力学模型，分析太阳帆在绕地轨道中遭受热辐射冲击作用时的结构动力学响应，为绕地轨道太阳帆的结构设计、轨道设计和姿态控制提供理论依据。

1 模型与假设

本文以正方形五点连接式太阳帆为研究对象，总体结构如图2所示，三角形薄膜帆面通过张拉索与对角线桅杆连接，帆面薄膜(简称：帆膜)在预应力作用下保持平整，有效载荷安装在太阳帆中心。

图2 五点连接式方形太阳帆Fig.2 Five-point suspension square solar sail

本文太阳帆热-结构耦合分析采用如下假设：

1) 光压很小，太阳帆结构发生小变形。

2) 帆膜褶皱幅值较小，忽略褶皱对光压压心位置和帆膜吸热的影响。

3) 帆膜通过张拉索与桅杆相连，忽略桅杆与帆膜之间的热传导。

4) 太阳帆结构对称。

2 太阳帆的热-结构耦合动力学方程

太阳帆由地球阴影区进入光照区，桅杆-帆膜结构的温度场在突加太阳热流的作用下将发生瞬态变化。由于热弹性效应，薄壁桅杆和帆膜的瞬态温度变化将在桅杆和帆膜内部产生冲击热应力，进而引发太阳帆结构的动力学响应。本节建立太阳帆桅杆-帆膜结构的热-结构动力学方程。

2.1 瞬态热传导方程

太阳帆进入光照区后，薄壁桅杆受到太阳辐射热流S0的作用如图3所示，(O,x,y,z)为桅杆局部坐标系;θ0为太阳辐射热流入射方向与桅杆外表面法线的夹角;h为桅杆截面壁厚。

图3 太阳辐射热流作用下的太阳帆薄壁桅杆Fig.3 Thin-walled boom subjected to solar heat flux

忽略沿桅杆壁厚方向的温度变化，薄壁桅杆的温度场为二维分布：T(x,φ,t)，满足如下非线性偏微分方程：

(1)

式中：c为桅杆比热；ρ为桅杆密度；t为时间；R为桅杆截面的平均半径；φ为桅杆截面壁厚中线的周向角坐标；kx，kφ分别为桅杆沿x向，φ向的热传导系数；ε为桅杆外表面的热辐射系数；σ为Stefan-Boltzman常数；q为桅杆外表面吸收的太阳辐射热量：

(2)

其中：αs为桅杆外表面的热辐射吸收系数；S0为太阳辐射热流S0的幅值。

将T(x,φ,t)在桅杆截面内的分布近似为截面平均温度Ta(x,t)和截面摄动温度Tp(x,t)的叠加[4]：

T(x,φ,t)≈Ta(x,t)+Tp(x,t)cosφ

(3)

将式(3)代入式(1)，并在桅杆截面内进行积分，可将式(1)解耦为如下2个偏微分方程：

(4)

(5)

式中：k为热传导系数。其初始条件为：Ta(x,t0)=T0和Tp(x,t0)=0，其中T0为太阳帆结构初始温度，t0表示初始时刻。

太阳辐射热流S0作用下，帆膜的温度场也将发生变化。考虑到帆膜厚度很小，忽略膜厚方向的温度变化，且不考虑褶皱对吸热的影响，则帆膜温度场在整个帆膜结构内均匀分布：Tm(t)，满足如下微分方程：

(6)

其初始条件为：Tm(t0)=T0。式中：cm为帆膜比热；ρm为帆膜密度；εm为帆膜表面的热辐射系数；αm为帆膜表面的吸收系数；tm为帆膜厚度。

式(4)～式(6)即为描述太阳帆桅杆-帆膜结构瞬态温度场的热传导方程。

2.2 结构动力学方程

由于热弹性效应，截面平均温度Ta的升降将引起桅杆沿轴向的伸缩变形，截面摄动温度Tp的变化将引起桅杆的弯曲变形，这些变形可等效为分别由外载荷热轴力FT和热弯矩MT引起的，其计算公式为

众创空间是集聚创客、高校、企业和科研机构等各类主体的协同创新合作体系，可以实现跨行业跨领域的创新资源弹性配置，并形成聚合效应和叠加效应。它不仅为创客团队提供了实现突破性创新的要素和资源，还为协同创新搭建了深度互联的网络服务平台，有助于创客团队进行全向度创新。而且众创空间在协同创新过程中能够不断演化，形成正反馈链条，从而加速产品创新和迭代速度，打造出融合、开放、多元的协同创新生态圈。

FT=αTEA(Ta-T0)

(7)

(8)

式中：αT为桅杆热膨胀系数；A为桅杆截面面积；E为桅杆弹性模量；I为桅杆截面关于z轴的惯性矩。

对于帆膜结构，由于帆膜厚度很小，帆膜厚度方向的温差变化无法形成帆膜热弯矩，因而不会引起帆膜的热诱发振动，但是，帆膜温度Tm将引起帆膜预应力的变化，从而影响张拉索作用在桅杆上的拉力Fm。

根据对称性，太阳帆桅杆可看作为薄壁悬臂梁结构，其受到热轴力FT、热弯矩MT和张拉索拉力Fm的作用，如图4所示。

采用Euler-Bernoulli梁理论，太阳帆桅杆在热载荷作用下的运动方程为

(9)

图4 热应力作用下的太阳帆薄壁桅杆Fig.4 Thin-walled boom subjected to thermally- induced loads

需满足的边界条件为

w(0,t)=0

(10a)

(10b)

(10c)

(10d)

设w(x,t)可近似表示为

w(x,t)=W1(t)N1(x)+W2(t)N2(x)

(11)

式中：W1(t)和W2(t)为广义位移；N1(x)和N2(x)为形函数。考虑强制边界条件式(10a)和式(10b)，可取N1(x)=x2，N2(x)=x3。

采用伽辽金加权余量法，分别取罚函数为N1(x)和N2(x)，考虑到边界条件式(10c)和式(10d)，可得式(9)的等效积分格式：

+N1(l)[EIw″(l,t)+MT(l,t)]+

N1(l)[EIw′″(l,t)+Fe]=0

(12)

+N2(l)[EIw″(l,t)+MT(l,t)]+

N2(l)[EIw′″(l,t)+Fe]=0

(13)

将式(11)分别代入式(12)和式(13)，并在桅杆域内分部积分，可得

(14)

式中：

式(14)和式(4)～式(6)共同构成了描述太阳帆桅杆-帆膜结构热诱发振动的热-结构耦合动力学方程组。

2.3 热-结构动力学方程的求解

太阳帆桅杆-帆膜结构的热-结构耦合动力学方程组可采用直接积分法求解，由于热辐射项的存在，其中式(4)和式(6)是非线性的，在每个时间步还需要迭代求解，其具体求解流程如图5所示。

图5 热-结构耦合动力学方程的求解流程Fig.5 Solving algorithm of thermal-structure coapling dynamic equations

按照上述求解流程，本文使用Visual-Fortran 90编写了太阳帆桅杆-帆膜的热-结构动力学计算程序。

3 太阳帆的热诱发振动分析

文献[4]考虑桅杆静态变形、屈曲载荷、总质量等约束条件，以太阳帆特征加速度达到1 mm/s2为目标给出了边长100 m的方形太阳帆优化设计方案，本文对该太阳帆方案进行热致结构响应分析，其结构尺寸和材料参数见表1和表2。

表1 太阳帆结构尺寸Table 1 Structural dimensions of solar sail

表2 太阳帆材料参数Table 2 Material parameters of solar sail

太阳帆绕地运行，由地球阴影进入光照区，受到突加太阳热流S0=1 350 W/m2和突加光压P0=9.12×10-6Pa的作用。

3.1 太阳帆的热致结构动力学响应

以太阳帆进入光照区为初始时刻t0=0 s，设太阳帆桅杆截面初始平均温度Ta(x,0)=290 K,初始摄动温度Tp(x,0)=0 K，热流入射角度为0°。

图6为太阳帆桅杆截面平均温度随时间变化的历程图，桅杆截面平均温度Ta在太阳辐射热流作用下缓慢上升，在t=450 s后达到稳态温度413 K。

图6 桅杆自由端截面平均温度时间历程Fig.6 Variation of average temperature of boom cross-section with time

图7为太阳帆桅杆截面摄动温度随时间变化的历程图，桅杆截面摄动温度Tp在太阳辐射热流作用下快速上升，在t=250 s后达到稳定，桅杆截面最大温差达到了92.98 K。

图7 桅杆自由端截面摄动温度时间历程Fig.7 Variation of perturbation temperature of boom cross-section with time

图8给出了桅杆自由端截面内的温度变化，可以看出桅杆温度在截面内的正弦分布，光照面的最大温度达到了496 K，背光面的温度为329.9 K，桅杆截面内的温差达到了166 K。

图8 桅杆自由端截面温度分布Fig.8 Temperature distribution on cross-section of boom tip

太阳帆桅杆截面内的温差将引起热弯矩，其作用在桅杆上会引起太阳帆桅杆的动态响应，图9给出了冲击热弯矩引起的桅杆自由端位移时间历程。

图9 桅杆自由端的位移-时间历程Fig.9 Displacement-time history of boom tip

由图9可以看出，在太阳辐射热冲击作用下，太阳帆桅杆发生了明显的热诱发振动，振动形式为围绕准静态变形的等幅振动，其准静态变形为3.16 m，平均振幅为0.75 m，最大动态位移达到了4.10 m。

3.2 太阳帆热诱发振动的参数分析

太阳帆发生热诱发振动的条件和特点取决于太阳帆结构刚度、进入光照区的姿态角和帆膜预应力等因素，本小节对这些影响太阳帆热诱发振动的因素进行参数分析。

太阳帆进入光照区的姿态会造成太阳辐射热流相对于结构的入射角θ0不同，从而影响结构吸收热量的不同。图10为不同热流入射角下，桅杆平均温度随时间变化的历程图，可以看出随着辐射热流由垂直入射变为斜射，桅杆表面吸收热量减少，太阳帆桅杆截面平均温度Ta的稳态值变小：由垂直入射θ0=0°的时的413.1 K，逐渐减小为θ0=75°的时的294.5 K。图11为不同热流入射角下，桅杆截面摄动温度随时间变化的历程图，同样，随着热流由垂直入射变为斜射，桅杆截面的温差Tp也逐渐减小：由垂直入射θ0=0°时的83.2 K，逐渐减小为θ0=75°时的32.0 K。

图10 热流入射角度对桅杆截面平均温度的影响Fig.10 Effect of incident angle of heat flux on average temperature of boom cross-section

图11 热流入射角度对桅杆截面摄动温度的影响Fig.11 Effect of incident angle of heat flux on perturbation temperature of boom cross-section

如前所述，摄动温度Tp会对桅杆的热弯矩造成影响，图12为不同热流入射角下，桅杆自由端位移随时间变化的历程图，可以看出热流入射角的增大造成热弯矩减小，使得桅杆热诱发振动的准静态变形和振幅均有所减小。热流垂直入射θ0=0°时的太阳帆桅杆热诱发振动最为剧烈，当热流入射角θ0=75°时桅杆的热诱发振动可近似看作静态的热致变形。

图13为桅杆壁厚参数变化时，桅杆自由端位移随时间变化的历程图，随着桅杆刚度增大，桅杆热诱发振动的振幅和准静态变形减小，频率增大，即太阳帆桅杆刚度的增大有利于抑制热诱发振动的发生。

提高帆膜预应力可以有效减小帆面的褶皱，提高光压的反射效率，但是，帆膜预应力也会造成桅杆轴力的变化，从而影响桅杆的热致振动响应。图14为帆膜预应力引起桅杆轴力变化时，桅杆自由端位移随时间变化的历程图，其中Pcr为桅杆的一阶屈曲载荷值。可知随着轴力增大，桅杆热诱发振动的振幅增大，频率减小。当轴力增大至一定程度，热诱发振动将会发散，因此，帆膜预应力对于热诱发振动是不利因素，须适当设计。