光对被照射物体的压力称为光压。依靠光压推进的太阳帆航天器无需消耗推进剂，可以实现传统航天器无法实现的各种非开普勒轨道，被认为是最可行的星际探测技术。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 模型与假设

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

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

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

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

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

4) 太阳帆结构对称。

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

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

2.1 瞬态热传导方程

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

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

$ c\rho \frac{{\partial T}}{{\partial t}} - {k_x}\frac{{{\partial ^2}T}}{{\partial {x^2}}} - \frac{{{k_\varphi }}}{{{R^2}}} \cdot \frac{{{\partial ^2}T}}{{\partial {\varphi ^2}}} + \frac{{\varepsilon \sigma }}{h}{T^4} = \frac{q}{h} $	（1）

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

$ q\left( {x,\varphi ,t} \right) = \left\{ {\begin{array}{*{20}{l}} {{\alpha _{\rm{s}}}{S_0}\cos {\theta _0}\cos \varphi }&{\left| \varphi \right| \le {\rm{ \mathsf{ π} }}/2}\\ 0&{\left| \varphi \right| > {\rm{ \mathsf{ π} }}/2} \end{array}} \right. $	（2）

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

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

$ T\left( {x,\varphi ,t} \right) \approx {T_{\rm{a}}}\left( {x,t} \right) + {T_{\rm{p}}}\left( {x,t} \right)\cos \varphi $	（3）

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

$ \frac{{\partial {T_{\rm{a}}}}}{{\partial t}} + \frac{{\varepsilon \sigma }}{{\rho ch}}T_{\rm{a}}^4 = \frac{1}{{\rm{ \mathsf{ π} }}} \cdot \frac{{{\alpha _{\rm{s}}}{S_0}}}{{\rho ch}}\cos {\theta _0} $	（4）

$ \frac{{\partial {T_{\rm{p}}}}}{{\partial t}} + \left( {\frac{k}{{\rho c{R^2}}} + \frac{{4\varepsilon \sigma }}{{\rho ch}}T_{\rm{a}}^3} \right){T_{\rm{p}}} = \frac{1}{2} \cdot \frac{{{\alpha _{\rm{s}}}{S_0}}}{{\rho ch}}\cos {\theta _0} $	（5）

式中：k为热传导系数。其初始条件为：T_a(x, t₀)=T₀和T_p(x, t₀)=0，其中T₀为太阳帆结构初始温度，t₀表示初始时刻。

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

$ {c_{\rm{m}}}{\rho _{\rm{m}}}\frac{{\partial {T_{\rm{m}}}}}{{\partial t}} + \frac{{{\varepsilon _{\rm{m}}}\sigma }}{{{t_{\rm{m}}}}}T_{\rm{m}}^4 = \frac{{{\alpha _{\rm{m}}}{S_0}\cos {\theta _0}}}{{{t_{\rm{m}}}}} $	（6）

其初始条件为：T_m(t₀)=T₀。式中：c_m为帆膜比热；ρ_m为帆膜密度；ε_m为帆膜表面的热辐射系数；α_m为帆膜表面的吸收系数；t_m为帆膜厚度。

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

2.2 结构动力学方程

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

$ {F_{\rm{T}}} = {\alpha _{\rm{T}}}EA\left( {{T_{\rm{a}}} - {T_0}} \right) $	（7）

$ \begin{array}{l} {M_{\rm{T}}}\left( {x,t} \right) = \int_A {{\alpha _{\rm{T}}}E\left( {{T_{\rm{p}}} - {T_0}} \right)y{\rm{d}}A} = \\ \;\;\;\;\;\;\frac{{E{I_{{\alpha _{\rm{T}}}}}}}{R}{T_{\rm{p}}}\left( {x,t} \right) \end{array} $	（8）

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

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

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

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

$ EI\frac{{{\partial ^4}w}}{{\partial {x^4}}} + P\frac{{{\partial ^2}w}}{{\partial {x^2}}} + \rho A\frac{{{\partial ^2}w}}{{\partial {t^2}}} + \mu \frac{{\partial w}}{{\partial t}} = 0 $	（9）

需满足的边界条件为

$ w\left( {0,t} \right) = 0 $	（10a）

$ \frac{{\partial w\left( {0,t} \right)}}{{\partial x}} = 0 $	（10b）

$ - EI\frac{{{\partial ^2}w\left( {l,t} \right)}}{{\partial {x^2}}} = {M_{\rm{T}}}\left( {l,t} \right) $	（10c）

$ - EI\frac{{{\partial ^3}w\left( {l,t} \right)}}{{\partial {x^3}}} = {F_{\rm{e}}}\left( {l,t} \right) $	（10d）

式中：w(x, t)为桅杆沿y向的位移；l为对角线桅杆长度的一半；μ为阻尼；$F_{\mathrm{e}}=\frac{2 P_{0} A_{\mathrm{s}}}{3}$为帆膜光压通过张拉索施加在桅杆端部的y向载荷，其中P₀为光压，A_s为三角形帆膜的面积；P=F_T+ $\sqrt {2} $F_m为由桅杆热轴力F_T和张拉索共同作用引起的桅杆轴向载荷。

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

$ w\left( {x,t} \right) = {W_1}\left( t \right){N_1}\left( x \right) + {W_2}\left( t \right){N_2}\left( x \right) $	（11）

式中：W₁(t)和W₂(t)为广义位移；N₁(x)和N₂(x)为形函数。考虑强制边界条件式(10a)和式(10b)，可取N₁(x)=x²，N₂(x)=x³。

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

$ \begin{array}{l} \int_0^l {{N_1}\left( x \right)\left( {EI\frac{{{\partial ^4}w}}{{\partial {x^4}}} + P\frac{{{\partial ^2}w}}{{\partial {x^2}}} + \rho A\frac{{{\partial ^2}w}}{{\partial {t^2}}} + \mu \frac{{\partial w}}{{\partial t}}} \right){\rm{d}}x} \\ \;\;\;\;\;\;\; + {N_1}\left( l \right)\left[ {EIw''\left( {l,t} \right) + {M_{\rm{T}}}\left( {l,t} \right)} \right] + \\ \;\;\;\;\;\;\;{N_1}\left( l \right)\left[ {EI{w^{\prime \prime \prime }}\left( {l,t} \right) + {F_{\rm{e}}}} \right] = 0 \end{array} $

（12）

$ \begin{array}{l} \int_0^l {{N_2}\left( x \right)\left( {EI\frac{{{\partial ^4}w}}{{\partial {x^4}}} + P\frac{{{\partial ^2}w}}{{\partial {x^2}}} + \rho A\frac{{{\partial ^2}w}}{{\partial {t^2}}} + \mu \frac{{\partial w}}{{\partial t}}} \right){\rm{d}}x} \\ \;\;\;\;\;\;\; + {N_2}\left( l \right)\left[ {EIw''\left( {l,t} \right) + {M_{\rm{T}}}\left( {l,t} \right)} \right] + \\ \;\;\;\;\;\;\;{N_2}\left( l \right)\left[ {EI{w^{\prime \prime \prime }}\left( {l,t} \right) + {F_{\rm{e}}}} \right] = 0 \end{array} $

（13）

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

$ \mathit{\boldsymbol{M\ddot W}} + \mathit{\boldsymbol{C\dot W}} + \mathit{\boldsymbol{KW}} = \mathit{\boldsymbol{F}}\left( t \right) $	（14）

式中：

$ \mathit{\boldsymbol{W}} = {\left[ {\begin{array}{*{20}{c}} {{W_1}}&{{W_2}} \end{array}} \right]^{\rm{T}}},\;\;\;\;\;\;\mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\rho A{l^5}}}{5}}&{\frac{{\rho A{l^6}}}{6}}\\ {\frac{{\rho A{l^6}}}{6}}&{\frac{{\rho A{l^7}}}{7}} \end{array}} \right] $

$ \mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{l}} {\frac{{\mu {l^5}}}{5}}&{\frac{{\mu {l^6}}}{6}}\\ {\frac{{\mu {l^6}}}{6}}&{\frac{{\mu {l^7}}}{7}} \end{array}} \right] $

$ \mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{l}} {4EIl + \frac{{2P{l^3}}}{3}}&{6EI{l^2} + \frac{{3P{l^4}}}{2}}\\ {6EI{l^2} + \frac{{P{l^4}}}{2}}&{12EI{l^3} + \frac{{6P{l^5}}}{5}} \end{array}} \right] $

$ \mathit{\boldsymbol{F}} = \left[ {\begin{array}{*{20}{l}} {{F_{\rm{e}}}{l^2} - 2l{M_{\rm{T}}}}\\ {{F_{\rm{e}}}{l^3} - 3{l^2}{M_{\rm{T}}}} \end{array}} \right] $

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

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

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

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

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

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

太阳帆绕地运行，由地球阴影进入光照区，受到突加太阳热流S₀=1 350 W/m²和突加光压P₀=9.12×10^－6 Pa的作用。

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

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

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

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

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

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

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

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

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

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

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

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

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

4 结论

本文首次建立了太阳帆桅杆-薄膜结构的热-结构动力学模型，给出了热-结构耦合动力学方程组的求解流程，分析了太阳帆桅杆-薄膜结构热诱发振动的特点和影响因素。得到了以下结论：

1) 绕地运行的方形桅杆支撑式太阳帆，由地球阴影进入光照区将会发生明显的热诱发振动，振动形式为围绕准静态变形的等幅振动。

2) 太阳帆由阴影区进入光照区的姿态角会影响太阳帆结构的热诱发振动，太阳热流垂直入射时，桅杆截面的平均温度和摄动温度最大，热诱发振动最剧烈，随着热流入射角的增大，桅杆截面摄动温度减小，热诱发振动减弱，频率不变。

3) 随着桅杆壁厚增加，刚度增大，桅杆热诱发振动的振幅和准静态变形减小，频率增大，桅杆刚度增大有利于抑制太阳帆结构发生热诱发振动。

4) 帆膜预应力将造成桅杆轴力的变化，影响桅杆的热致振动响应，随着桅杆轴力增大，桅杆热诱发振动的振幅增大，频率减小。当轴力增大至一定程度，热诱发振动将会发散。