计及轴向力作用的柔性梁热-形耦合振动研究

吴磊; 韩飞; 邓子辰; Wu Lei; Han Fei; Deng Zichen

主管单位：中华人民共和国工业和信息化部
主办单位：西北工业大学中国航空学会
地址：西北工业大学友谊校区航空楼

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

计及轴向力作用的柔性梁热—形耦合振动研究 PDF

- ORCID：
吴磊 ¹
- ORCID：
韩飞 ¹
- ORCID：
邓子辰 ²
✉

1. 西北工业大学力学与土木建筑学院，西安 710129； 2. 西北工业大学航空学院，西安 710072

中图分类号： V214； V414

最近更新：2024-12-19

DOI：10.16615/j.cnki.1674-8190.2024.06.20

摘要

飞行器中的梁—索组合式结构在服役期间由于太阳热辐射引起的交变热流，会使结构产生热致变形及振动。杆件的变形又会进一步影响热流的大小和分布，是一类典型的力—形耦合系统。耦合系统的刚度、质量或载荷矩阵通常具有时变特征，通常只能通过多次求解特征方程并借助逐步积分的方法获得系统的模态及响应。如何保证求解过程的计算精度和效率是一大挑战。基于动力刚度方法与Wittrick-Williams算法，进行系统固有频率和模态的精确求解，结合精细积分法计算结构响应。结果表明：前4阶固有频率与有限元解的偏差在1‰以内，响应最大相对偏差为1%；当时变热轴力接近结构失稳临界载荷时，结构振动频率快速减小并趋于0，热致振动将出现发散现象。

关键词

梁—索结构; 热致振动; 时变热轴力; 动力刚度法; Wittrick-Williams算法

0 引言

飞行器运行过程中，其柔性附件由于太阳热辐射引起的交变热流，会致使结构产生热变形及振动等^［

1-3］，特别在航天器进出地影时，受到的热荷载会急剧变化，甚至诱发结构发生颤振。最著名的就是20世纪美国哈勃望远镜的柔性太阳翼弯扭耦合的热颤振事故^{［参考文献 4-5}4-5］。热致振动直接影响了望远镜的成像精度，后来借助航天飞机将其太阳翼更换成半刚性结构，才使望远镜能够正常工作。因此，飞行器在轨运行过程中的热致振动是关乎结构设计和服役期安全性的重要问题，必须予以重视。

柔性梁热致振动的研究，最早要追溯到20世纪50年代，Boley^［

5］从理论上预测了热致振动现象，并提出了一个无量纲参数B来衡量热致振动的剧烈程度。在航天器热致振动力学建模时，主要采用欧拉梁模型，同时计及热荷载与结构变形的耦合效应，即“热—形耦合”模型^{［参考文献 6

百度学术}6］。Shen Z等^{［参考文献 7-8}7-8］基于绝对节点坐标方法推导了梁的热—结构耦合动力学方程，该方程不但能计算昼夜转换热冲击引起的梁热颤振，也能计算由于大转动引起的梁热致振动；薛明德等^{［参考文献 9-10}9-10］提出了傅里叶温度单元方法，显著提高了薄壁结构非线性温度场的求解效率，并将该方法用于大型柔性空间结构的热—结构耦合动力学分析中，促进了热致结构响应分析的工程应用；薛碧洁^{［参考文献 11

百度学术}11］针对非线性索—梁结构，首先通过结构特征时间和热特征时间推导了非线性索—梁结构Boley系数的解析表达式，然后推导了非线性索—梁结构热致振动的有限元方程，并采用 Newmark方法进行求解。

上述研究工作均忽略了热辐射引起的时变轴向载荷对系统几何刚度的影响，通常将其看作常数^［12］或直接忽略，其理由是该部分效应相比于弯曲荷载是“慢变量”^［

4］。然而，热轴力会随时间变化进而影响结构的动力学特性，甚至导致振动形式的改变，当轴力增大至一定程度，热致振动将会发散。采用有限元方法进行求解时，大多采用Newmark方法式进行时间积分，在每一时间步内用Newton-Raphson 法求解^{［参考文献 10

百度学术}10］。这种数值迭代方法计算非常繁琐，同时与切线斜率、步长密切相关，若考虑时变轴向载荷对结构刚度的影响，迭代法的计算难度将进一步增大。

动力刚度方法是一种精确算法^［13］，其精确度与高效率在于采用了与结构频率相关的精确形状函数。目前，动力刚度方法已被广泛应用于各式连续结构，包括桁架、连续梁、梁板组合结构等^［14-15］。动力刚度法相比于传统有限元方法的优势在于：只需要少数单元就可以计算各阶模态，计算效率大幅提高。动刚度表征的结构特征方程是超越方程，其精确求解可采用Wittrick-Williams（W-W）算法^［15］。W-W算法本质上是一种计数方法，它给出了结构低于用户指定频率的固有频率的数目，结合二分法等不断缩小区间，最终得到满足要求精度的频率和振型解^［16］，结合精细积分方法^［17］可以实现结构响应的高精度稳定求解。

为了更加精确地求解时变系统的热—形耦合问题，本文采用动力刚度法进行计算，基于欧拉梁模型建立薄壁悬臂梁结构在轨运行时的热—形耦合模型，相比之前的问题，考虑了时变热轴力对结构横向振动的影响；针对系统动力学方程，采用动力刚度法和Wittrick-Williams算法计算结构的时变频率，采用精细积分方法计算结构响应并对其进行参数分析。

1 数值方法

1.1 瞬态温度场分析

薄壁桅杆结构是航天器中的常见结构，如图1所示。

图1 柔性梁结构力学模型图

Fig.1 Structural model of beam

薄壁悬臂梁结构受到空间热流S₀的作用，其中α表示初始入射角， $α + \partial w / \partial x$ 表示结构变形后空间热流的入射角，显然，在考虑热—形耦合作用时，结构受到的实际热流与结构横向位移相关。当结构受到入射空间热流时，根据傅里叶定律以及能量守恒定律可以推导出以下非线性方程：

\frac{\partial T}{\partial t} - \frac{k}{ρ c R^{2}} \frac{\partial^{2} T}{\partial ϕ^{2}} + \frac{σ ε}{ρ c h} T^{4} = \frac{α_{s} S}{ρ c h} δ c o s ϕ

（1）

式中：k 为导热系数；ρ为密度；c为比热容；ε为表面发射率；α_s为表面吸收率；T为绝对温度；δ为热流（只作用于管的一侧）；S为实际作用于表面的有效热流（空间热流S₀在结构外表面的投影）。

S = S_{0} c o s (α + \frac{\partial w}{\partial x})

（2）

δ = \{\begin{array}{l} 1 (0 \leq ϕ \leq π) \\ 0 (- π \leq ϕ < 0) \end{array}

（3）

式中：α为入射角。

杆件承受热载荷的表达式为

M (t) = \frac{E I α_{T} T_{m} (t)}{R}

（4）

P_{T} (t) = E α_{T} (Δ \bar{T}) A

（5）

式中：M（t）为热弯矩；P_T（t）为热轴力，详细计算过程详见文献［

11］；T_m（t）为截面扰动温度；ΔT为结构平均温度差；E为弹性模量；α_T为热膨胀系数；A为截面面积。

1.2 结构时变模态分析

计及轴向载荷影响的杆件热—形耦合振动方程为

E I \frac{\partial^{4} w (x, t)}{\partial x^{4}} + ρ A \frac{\partial^{2} w (x, t)}{\partial t^{2}} - P (t) \frac{\partial^{2} w (x, t)}{\partial x^{2}} = M (x, t)

（6）

由于空间热流S与结构横向位移相关，因此方程（6）中的P、M均与结构横向位移相关，这将大幅增加方程的求解难度。由于轴向荷载P相比弯矩M是缓慢变化的，因此在结构响应分析时可对其采用准静态的方式处理，即仅考虑其对系统几何刚度的影响。此时，系统的动刚度矩阵D，刚度矩阵K，动质量矩阵M以及几何矩阵G分别为

[M (ω, P)] = \int_{0}^{l} ρ A {[N (ω, P, x)]}^{T} [N (ω, P, x)] d x

（7）

\begin{array}{l} [K (ω, P)] = \\ \int_{0}^{l} E I {[N " (ω, P, x)]}^{T} [N " (ω, P, x)] d x \end{array}

（8）

\begin{array}{l} [G (ω, P)] = \\ \int_{0}^{l} {[N' (ω, P, x)]}^{T} [N' (ω, P, x)] d x \end{array}

（9）

[D (ω, P)] = [K (ω, P)] + [G (ω, P)] -

ω^{2} [M (ω, P)]

（10）

式中： $[N (ω, P, x)]$ 为与频率相关的形函数。

[N (ω, P, x)] = [χ (ξ)] Q

（11）

χ (ξ) = [\begin{matrix} c o s α ξ & s i n α ξ & c o s h β ξ & s i n h β ξ \end{matrix}]

（12）

Q = \frac{1}{α^{2} + β^{2}} [\begin{matrix} β^{2} - Q_{4} & Q_{2} l & - Q_{3} & Q_{1} l \\ - Q_{6} / α & (α + Q_{4} / α) / l & - Q_{5} / α & - Q_{3} l / α \\ α^{2} + Q_{4} & - Q_{2} l & Q_{3} & - Q_{1} l \\ Q_{6} / β & (β - Q_{4} / β) / l & Q_{5} / β & Q_{3} l / β \end{matrix}]

（13）

Q_{1} = (β s i n α - α s i n h β) (α^{2} + β^{2}) / δ

（14）

Q_{2} = (α c o s α s i n h β - β s i n α c o s h β) (α^{2} + β^{2}) / δ

（15）

Q_{3} = (c o s α - c o s β) α β (α^{2} + β^{2}) / δ

（16）

Q_{4} = [(β^{2} - α^{2}) (c o s α c o s h β - 1) +

2 α β s i n α s i n h β] α β / δ

（17）

Q_{5} = (β s i n h β + α s i n α) (α^{2} + β^{2}) α β / δ

（18）

Q_{6} = - (α c o s h β s i n α + β s i n h β c o s α) \times

(α^{2} + β^{2}) α β / δ

（19）

δ = 2 α β (c o s α c o s h β - 1) + (α^{2} - β^{2}) s i n α s i n h β

（20）

α^{2} = - \frac{σ^{2}}{2} + \sqrt[]{\frac{σ^{4}}{4} + λ^{4}}

（21）

β^{2} = \frac{σ^{2}}{2} + \sqrt[]{\frac{σ^{4}}{4} + λ^{4}}

（22）

λ^{4} = \frac{ρ A l^{4} ω^{2}}{E I}

（23）

σ^{2} = \frac{P l^{2}}{E I}

（24）

ξ = \frac{x}{l}

（25）

动刚度矩阵、刚度矩阵以及动质量矩阵存在以下关系：

M (ω, P) = - \frac{\partial^{2} D (ω, P)}{\partial ω^{2}}

（26）

K (ω, P) = \frac{\partial D (ω, P)}{\partial ω}

（27）

得到动刚度矩阵D后，通过W-W算法求解结构频率。

W-W算法是一种频率计数方法，它没有提供频率方程的直接解，而是通过计算所选频率阶数下的模态频率数J（ω^*）来给出模态频率的上下限^［17］。通过该过程的迭代并对频率近似解ω^*的不断改进，得到了第k阶频率的上下界ω_u、ω_l，使得满足条件：

\{\begin{array}{l} J (ω_{u}) \geq k \\ J (ω_{l}) \leq k - 1 \end{array}

（28）

在得到粗略的范围后，采用牛顿法进行计算，缩小上下界范围，直至：

ω_{u} - ω_{l} \leq T_{o l}

（29）

式中： $T_{o l}$ 为提前设置的误差限。

对于大多数结构来说，低于给定值ω^*的频率个数J很难直接给出。通常由以下两项的和来表示：

J = J_{0} + J_{K} = J_{0} + s \{D (ω^{*})\}

（30）

式中：J₀为低于ω^*的单元固端频率数；J_K为用Gauss消元法将D（ω^*）消成上三角矩阵D^Δ（ω^*）后，主对角线上负元素的个数。

采用牛顿法缩小所求频率区间时，也可以得到结构振型，详细过程详见文献［18］。

1.3 热致振动响应分析

结构的运动方程

M \ddot{u} + C \dot{u} + K u = F

（31）

式中：M、C、K分别为质量矩阵、阻尼矩阵以及刚度矩阵；F为载荷向量；u为位移向量。其增量形式为

M_{t} Δ {\ddot{u}}_{t} + C_{t} Δ {\dot{u}}_{t} + K_{t} Δ u_{t} = Δ F_{t}

（32）

ΔF_t=F_t_+Δ_t-F_t，求出第n步的位移增量Δu_t=Δu_n后，第（n+1）步的位移为

u_{n + 1} = u_{n} + Δ u_{n}

（33）

式中： $Δ u_{n}$ 为第n步的初始位移。

假设每一时间步长内系统具有线性特性，因此在每一时间步长内可以使用振型叠加法，每一步长内的位移增量可以表示为

Δ u_{t} = ϕ_{t} Δ η_{t}

（34）

振型矩阵 $ϕ_{t} = [ϕ_{t}^{(1)}, ϕ_{t}^{(2)}, \dots, ϕ_{t}^{(N_{n})}]$ 。

根据模态的正交性，增量方程可以去耦为

Δ {\ddot{η}}_{t}^{(i)} + 2 ζ_{i} ω_{t}^{(i)} Δ {\dot{η}}_{t}^{(i)} + [ω_{t}^{(i)}]^{2} Δ η_{t}^{(i)} = Δ N_{t}^{(i)}

（35）

式中：Δη为局部模态坐标增量；ζ为模态阻尼比；ω^（ⁱ^）为第i阶频率；ΔN^（ⁱ^）为局部模态广义力。

本文采用精细积分法求解方程（35），具体过程如下：

1）引入对偶变量

\{\begin{array}{l} Z (t) = [\begin{matrix} Δ η (t) \\ q (t) \end{matrix}] \\ q (t) = Δ \dot{η} (t) + \frac{ζ Δ η (t)}{2} \end{array}

（36）

将式（36）转化为Hamilton系统下一阶常微分方程组：

\dot{Z} (t) = H Z (t) + F (t)

（37）

式中： $H = [\begin{matrix} E_{11} & E_{12} \\ E_{21} & E_{22} \end{matrix}], F (t) = \{\begin{matrix} 0 \\ Δ N_{t}^{(i)} \end{matrix}\}$ ， $E_{11} = - \frac{ζ}{2}, E_{12} = I^{- 1}, E_{21} = \frac{ζ I^{- 1} ζ}{4} - W, E_{22} = - \frac{ζ I^{- 1}}{2}$ 。

选择适当的时间步长Δt，引入微小时段Δτ=Δt/2^N。一般取N=20即可满足计算精度。计算：

T_{a, 0} = H Δ τ + \frac{{(H Δ τ)}^{2}}{2} [I + \frac{H Δ τ}{3} + \frac{{(H Δ τ)}^{2}}{12}]

（38）

再根据下式求出指数矩阵T：

T = I + T_{a, N} = {(I + T_{a, N - 1})}^{2} = \dots = {(I + T_{a, 0})}^{2^{N}}

（39）

其中：

T_{a, i} = 2 T_{a, i - 1} + T_{a, i - 1} T_{a, i - 1} (i = 1,2, \dots, N)

（40）

2）当时间步长较短时，可假设激励的变化是线性的

f (t) = f_{1} + f_{2} (t - t_{i})

（41）

式中：f₁和f₂为常向量。

3）计算每一时刻的局部模态坐标增量

Z_{t + 1} = T [Z_{t} + H^{- 1} (f_{1} + H^{- 1} f_{2})] -

H^{- 1} (f_{1} + H^{- 1} f_{2} + Δ t f_{2})

（42）

再计算下一步的响应。具体计算流程如图2所示。

图2 计算流程图

Fig.2 Calculation process

2 算例与分析

2.1 动响应分析

本文数值算例考虑重载飞艇中的梁—索结构（如图3所示），将拉索作用考虑为约束力，其结构尺寸与材料参数如表1所示。

图3 热应力下梁—索结构示意图

Fig.3 Schematic representation of a beam-cable structure subjected to thermal stress.

表1 材料参数

Table 1 Material parameters

参数	数值
桅杆长度L/m	40
桅杆外径R/m	35.55×10^-3
桅杆截面壁厚/m	R/200
桅杆弹性模量E/GPa	689.5
桅杆密度ρ/（kg·m^-3）	1 660
桅杆比热c/［J·（kg·K）^-1］	502
桅杆表面吸收率α_S	0.5
桅杆热膨胀系数α_T/K^-1	1.69×10^-8
表面发射率ε	0.13
材料导热系数k/［W·（kg·K）^-1］	16.61
太阳辐射热流S₀/（W·m^-2）	1 350
波尔兹曼系数σ/［（W·kg）·K^-4］	5.67×10^-8
张拉索力F_c/N	0.141

考虑飞行器在进入光照区为初始t=0时刻，假设梁截面初始平均温度为290 K，截面摄动温度为0 K，初始角度为0°，忽略轴向载荷。

桅杆整体平均温度随时间变化的历程图如图4所示，可以看出：相比于截面扰动温度，平均温度的变化较慢，这符合文献［

3］中描述的“慢变量”；在600 s时间后达到稳定值413 K，最大温差可以达到123 K。

图4 自由端截面平均温度曲线

Fig.4 Average temperature curve of free end cross-section

动刚度法与有限元软件计算的1~4阶频率对比如表2所示，可以看出：精度可以达到小数点后四位，说明本方法针对频率计算的准确性。

表2 1~4阶频率对比

Table 2 1st~4th order frequency comparison

方法	频率/（rad·s^-1）
方法	第一阶	第二阶	第三阶	第四阶
动刚度法	1.119 06	7.033 39	19.702 23	38.612 29
有限元	1.119 06	7.033 40	19.702 23	38.612 30

动刚度法与有限元方法的计算耗时以及相对误差如表3所示，可以看出：动刚度法在单元数为1的情况下，计算精度已经达到10^-6；而有限元法，在单元数较少时，难以满足精度要求；在增多划分单元数后，计算时间是动刚度法的两倍以上。因此，可以说明动刚度法在计算效率上要优于传统有限元法。

表3 两种方法计算耗时及相对误差

Table 3 Calculation time and relative error of the two methods

方法	计算耗时/s	单元数	相对误差
动刚度法	0.26	1	10^-6
有限元法	0.39	1	0.58
	0.44	2	0.12
	0.65	10	10^-4

梁—索结构在空间热流作用下，自由端横向响应如图5所示，可以看出：在空间热流作用下，结构出现了明显的热致振动现象，最大位移达到了0.03 m；在150 s后，位移响应逐渐趋于稳定。当忽略热轴力时，计算结果与文献［

9］中的计算结果一致。

图5 无轴力自由端横向响应

Fig.5 Lateral response of free end without axial force

考虑轴向力的影响，轴向载荷P（t）为

P (t) = P_{0} + P_{T} (t)

（43）

式中：P₀为张拉绳索对梁结构的作用力，取为-0.2 N；P_T（t）为时变热轴力。

结构受到的热轴力曲线如图6所示，可以看出：热轴力在400 s左右趋于稳定，峰值达到56 N。

图6 轴向热载荷曲线

Fig.6 Axial thermal load curve

考虑时变热轴力对横向振动影响下的尖端横向响应如图7所示，可以看出：在0~5 s内，轴力对结构振动的影响不大，但随着时间推移，轴向载荷使得结构位移变大，不再出现稳定现象；在62 s时，考虑时变热轴力的情况下结构自由端横向位移已经发散。

图7 结构自由端横向响应

Fig.7 Lateral response at the free end of the structure

结构在空间热流入射角为0°时的应力曲线图如图8所示，表示从固定端向自由端应力变化。动应力包络线由 $σ = E R u^{″} (x, t)$ 计算。

图8 0°时结构应力包络线

Fig.8 Structure stress envelope at 0°

从图8可以看出：热应力远小于结构动应力。因此可以说明，时变热轴力对结构振动影响不可忽略，若忽略可能会出现安全隐患。

结构振动一阶频率如图9所示，可以看出：随着热轴力的增加，结构频率在62 s时接近0。

图9 结构一阶振动频率

Fig.9 First order vibration frequency of the structure

根据压杆稳定理论，结构承受的临界压力为 $F_{c r} = \frac{π^{2} E I}{{(μ l)}^{2}} = 26.48$ N。

热轴力在62 s时，达到了26.2 N，也就是说，在考虑时变热轴力的情况下，在62 s时，结构已经接近失稳载荷的临界值。

2.2 参数分析

空间热辐射的入射角直接影响结构承受的实际热辐射，因此，考虑不同初始入射角对结构振动的影响，θ分别取0°、30°、45°、70°以及85°，不同入射角条件下结构受到的热轴力如图10所示，可以看出：随着入射角的增大，轴向载荷峰值随之减小；从0°~70°，热轴力峰值从56 N缩小至12 N，缩小了78%，同时，热轴力达到稳定的时间也有所延长。值得注意的是，当入射角达到85°时，结构所承受的热轴力方向发生改变。

图10 不同入射角下结构受到轴向热载荷

Fig.10 Structures subjected to axial thermal loads at different incident angles

不同入射角下考虑时变热轴力的结构自由端横向位移如图11所示，可以看出：随着入射角的增大，在出现失稳之前，结构位移的幅值逐渐减小，同时，结构达到失稳破坏阶段的时间也有所增长。

图11 不同入射角下结构自由端横向响应

Fig.11 Lateral response of the free end of the structure at different incident angles

不同入射角下结构自由端横向响应如图12所示，可以看出：当入射角超过60°时，结构横向位移趋于稳定。

图12 不同入射角下结构自由端横向响应

Fig.12 Lateral response of the free end of the structure at different incident angles

不同入射角下结构一阶振动频率如图13所示，可以看出：随着入射角的改变，结构的频率减小速率也随之变缓。

图13 不同入射角情况下结构一阶振动频率

Fig.13 First order vibration frequency of structure under different incident angles

考虑长度缩短分别为40、30、20与10 m，入射角为0°，分析考虑热轴力时结构横向响应。不同长度下结构自由端横向响应如图14所示，可以看出：随着长度的缩短，结构振动发散时间明显延长，在没有发散阶段，结构横向振动振幅也明显减小；随着长度进一步缩短，时变热轴力已经不足以使结构振幅出现发散现象。

图14 不同长度下结构自由端横向响应

Fig.14 Lateral response of structural free end under different lengths

3 结论

1）时变热轴力对航空结构热致振动影响较大，结构频率会随之发生改变，热轴力可能使飞行器结构出现频率下降接近至0的现象，结构振动发散；结构动应力明显大于热应力，若忽略时变热轴力，可能使结构出现安全隐患。

2）当飞行器结构入射角较小时，结构受到的空间热流较大，此时结构的平均温度与摄动温度较大，结构振动更容易发散。因此，可以通过调整航天器姿态来改变入射角大小，使结构热致振动减弱。

3）当飞行器结构刚度增大（长度缩短、截面积增大）时，结构受热轴力的影响显著降低，因此，如何在增大结构尺度时，增大结构的刚度，是结构避免失稳的重要考虑方向。

参考文献

HU Weipeng， DENG Zichen. A review of dynamic analysis on space solar power station［J］. Astrodynamics， 2023（2）： 115-130. [百度学术]

张新榃，吴志刚，杨超. 大变形状态机翼振动试验与气动弹性分析［J］. 航空工程进展， 2010， 1（1）： 76-79. [百度学术]

ZHANG Xintan， WU Zhigang， YANG Chao. Vibration test and aeroelastie analysis of wing in large static deformation［J］. Advances in Aeronautical Science and Engineering， 2010， 1（1）： 76-79.（in Chinese） [百度学术]

毛一青，杨飞，谷迎松. 15米翼展太阳能飞机机翼颤振分析和刚度设计［J］. 航空工程进展， 2019， 10（4）： 536-541. [百度学术]

MAO Yiqing， YANG Fei， GU Yingsong. Research on flutter characteristics and stiffness design for 15 meters span wing of a solar-powered aircraft［J］. Advances in Aeronautical Science and Engineering， 2019， 10（4）： 536-541.（in Chinese） [百度学术]

薛明德，向志海. 大型空间结构热—动力学耦合分析方法综述［J］. 力学学报， 2022， 54（9）： 2361-2376. [百度学术]

XUE Mingde， XIANG Zhihai. Review of thermal-dynamical analysis methods for large space structures［J］. Chinese Journal of Theoretical and Applied Mechanics， 2022， 54（9）： 2361-2376.（in Chinese） [百度学术]

BOLEY B A. Approximate analyses of thermally induced vibrations of beams and plates［J］. Journal of Applied Mechanics， 1972， 39（1）： 787-796. [百度学术]

THORNTON E A， KIM Y A. Thermally induced bending vibrations of flexible rolled-up solar array［J］. Journal of Spacecraft & Rockets， 1993， 1（4）： 2138-2150. [百度学术]

SHEN Z， HU G. Thermoelastic-structural analysis of space thin-walled beam under solar flux［J］. AIAA Journal， 2019， 57（4）： 1784-1788. [百度学术]

SHEN Z， LI H， LIU X， et al. Thermal-structural dynamic analysis of a satellite antenna with the cable-network and hoop-truss supports［J］. Journal of Thermal Stresses， 2019， 42（6）： 1-18. [百度学术]

程乐锦，薛明德. 大型空间结构热—动力学耦合有限元分析［J］. 清华大学学报（自然科学版）， 2004（4）： 681-684，688. [百度学术]

CHENG Lejin， XUE Mingde. Thermal-dynamic analysis of large scale space structures by FEM［J］. Journal of Tsinghua University （Science and Technology）， 2004（4）： 681-684，688.（in Chinese） [百度学术]

FAN L， XIANG Z， XUE M， et al. Robust optimization of thermal-dynamic coupling systems using a Kriging model［J］. Journal of Spacecraft & Rockets， 2015， 47（6）： 1029-1037. [百度学术]

薛碧洁. 索梁结构的热致振动及动刚度分析［D］. 西安：西安电子科技大学， 2014. （下转第254页） [百度学术]