CN106096293B - 一种大转动复合材料伸展臂热致振动预测方法 - Google Patents

Abstract

本发明公开了一种大转动复合材料伸展臂热致振动预测方法，该方法采用周向均匀刚度配置的铺层方式，建立复合材料伸展臂耦合热‑结构动力学控制方程；基于Euler‑Bernoulli梁理论和绝对节点坐标方法(ANCF)，建立考虑热效应的ANCF梁单元；使用Fourier温度单元和ANCF梁单元建立复合材料伸展臂耦合热‑结构动力学求解模型；结合Crank‑Nicolson方法、Generalized‑α方法和Newton‑Raphson方法迭代求解耦合热‑结构动力学模型，最终获得复合材料伸展臂热致振动响应。本发明方法可以对大范围运动的伸展臂进行热致振动分析，不仅避免了坐标转换，而且使惯性力的计算变得十分容易。为未来复杂伸展臂结构设计提供有效的指导，具有实际工程应用价值。

Description

一种大转动复合材料伸展臂热致振动预测方法

技术领域

本发明主要是适用于具有大范围运动的复合材料伸展臂耦合/非耦合热-结构动力学响应问题，具体涉及一种大转动复合材料伸展臂热致振动预测方法。

背景技术

复合材料伸展臂因其质量轻、热容小的优点广泛应用于航天器结构设计。然而，伸展臂结构往往尺寸大、刚度低，这类柔性结构在轨运行时容易受到极端热载荷的作用诱发结构振动。这种热致振动会降低系统性能，甚至造成航天器结构损坏，因此成为现代航天器的典型故障。为避免在轨期间出现热致振动，在设计阶段进行航天器结构热致振动研究至关重要。

由于在地面上模拟空间环境十分复杂，相比于试验方法，热致振动分析成为问题研究的重要手段。1956年，Boley将瞬态热弯矩引入受热载荷作用的悬臂梁响应分析，在研究成果的基础上首次提出了热致振动的概念。随后，Thornton等人针对简化梁模型和太阳翼的热致振动问题进行了一系列有意义的研究，并给出了耦合与非耦合情况下求解薄壁管梁横截面温度场分布的解析方法。考虑各向异性本构关系，Song和Yoon研究了复合材料薄壁管梁的热-结构动力学响应，使得分析模型更加接近实际结构。除了上述解析方法，有限元数值方法也被引入热致振动分析。1968年，通过对比理论解，Mason验证了有限元方法进行热致振动分析的可靠性。Chen等人使用有限元方法分析了瞬态非线性热-结构耦合问题，由于采用块体单元求解温度分布使得分析过程非常耗时。结合Fourier级数和有限单元法，清华大学的薛明德等人提出了一种Fourier温度单元，该单元将温度场分解为平均温度和多谐摄动温度，使得平均温度和摄动温度有限元方程相互解耦。该团队的研究工作促进了热致振动分析的工程应用。Zhao等人借助有限元法研究了复合材料层合板的热致振动，利用包含剪切变形的高阶位移场对板结构进行热-结构动力学分析。

从目前的研究可以发现，相比于迅速发展的热分析有限元方法，结构方面仍然主要采用传统有限元方法。然而，对于大转动、大柔性的空间结构，传统有限元法不能准确描述结构的刚体运动，此外，由于节点坐标含有微小转动自由度，从而限制了一些重要单元的应用。

本发明针对复合材料伸展臂结构，提出了一种大转动复合材料伸展臂热致振动预测方法。本发明方法可以对大范围运动的伸展臂进行热致振动分析，不仅避免了坐标转换，而且使惯性力的计算变得十分容易。同时，本发明方法可以进一步延伸应用到多体柔性复合材料结构的热致振动分析中，为未来复杂伸展臂结构设计提供有效的指导。

发明内容

本发明的目的在于克服现有技术的不足，提供一种大转动复合材料伸展臂热致振动预测方法，该发明方法可以对大范围运动的伸展臂进行热致振动分析，不仅避免了坐标转换，而且使惯性力的计算变得十分容易。为解决空间大型复杂伸展臂耦合热-结构动力学响应问题提供了可行有效的预测方法。

为实现上述目的，本发明提供的一种大转动复合材料伸展臂热致振动预测方法，该方法实现步骤如下：

步骤一、基于空间伸展臂的结构特点，采用周向均匀刚度配置的铺层方式建立复合材料伸展臂耦合热-结构动力学控制方程；

步骤二、基于欧拉-伯努利(Euler-Bernoulli)梁理论和绝对节点坐标方法(ANCF)，建立考虑热效应的ANCF梁单元；

步骤三、采用傅里叶(Fourier)温度单元和ANCF梁单元离散步骤一建立的耦合热-结构动力学控制方程，建立复合材料伸展臂耦合热-结构动力学求解模型；

步骤四、在时间t_n上，采用克兰克-尼克尔森(Crank-Nicolson)方法求解步骤三建立的伸展臂传热模型，得到复合材料伸展臂温度结果

步骤五、根据步骤四得到的温度结果，在时间t_n上，采用广义-a(Generalized-a)方法求解步骤三建立的伸展臂结构动力学模型，得到复合材料伸展臂位移结果q_n+1；

步骤六、判断时间条件是否满足t_n+1≥t_tol，如果不满足条件，令t_n＝t_n+Δt，跳转到步骤四，其中，Δt是时间步长；如果满足条件，输出复合材料伸展臂位移响应时间历程q(t)，迭代终止。

其中，步骤一中建立的复合材料伸展臂耦合热-结构动力学控制方程，包括：伸展臂传热控制方程和结构动力学控制方程，具体表示为：

考虑热传导和热辐射的伸展臂传热控制方程：

式中，是伸展臂上任意点的绝对温度；k,ρ和c分别是材料的导热系数，密度和比热；ε和α分别是伸展臂外表面的发射率和吸收系数；σ是Boltzmann常数；δ是与角度相关的量，用于判断伸展臂外表面是否受到太阳辐射；S₀是太阳辐射热流；β是太阳辐射热流入射角；θ_z是伸展臂的转角位移；h是伸展臂薄壁厚度。

采用周向均匀刚度配置的铺层方式建立复合材料伸展臂结构动力学控制方程：

式中，a₃₃，a₅₅是刚度系数，b₁，b₄和b₁₄是惯性系数，m_z是热应力偶，m_z′是热应力偶关于位置坐标x的一阶导数，v₀是伸展臂横向位移，v₀′，v₀″分别是伸展臂横向位移关于位置坐标x的一阶导数和二阶导数，是伸展臂横向加速度，θ_z是伸展臂转角位移，θ_z′，θ_z″分别是伸展臂转角位移关于位置坐标x的一阶导数和二阶导数，是伸展臂角加速度。

其中，步骤二中采用位置矢量和斜率矢量表示考虑热效应的ANCF梁单元的节点坐标，具体表示为：

式中，r_k＝[r_k1 r_k2]^T,(k＝i,j)是单元节点位置矢量，是单元节点斜率矢量；

梁单元上任意一点的位置矢量由单元节点坐标表示为：

r(t)＝[r₁ r₂]^T＝S(x)q(t)

式中，r₁和r₂分别是位置矢量r在x和y方向的分量，q(t)是梁单元节点坐标向量，S(x)是梁单元的形函数，它是关于单元局部坐标的矩阵函数。

其中，步骤二中建立的考虑热效应的ANCF梁单元的刚度矩阵包括纵向刚度矩阵和横向刚度矩阵梁单元的质量矩阵M^e为常数矩阵，具体表示为：

纵向刚度矩阵

式中，E是材料的弹性模量，A是梁单元截面积，l是梁单元长度，是梁单元纵向热应变，q是梁单元节点坐标，S_l是与单元形函数有关的量，表示为：

式中，S_,x是单元形函数关于坐标x的一阶偏导数，是S_,x的转置矩阵，S_,ξ是单元形函数关于局部坐标ξ＝x/l的一阶偏导数，S_,ξ ^T是S_,ξ的转置矩阵，l是梁单元长度；

横向刚度矩阵

式中，是横向结构刚度矩阵，是热弯矩引起横向刚度矩阵，I是梁单元横截面惯性矩，M^T是梁单元的热弯矩，S_t是与单元形函数有关的量，表示为：

式中，S_,ξξ是单元形函数关于局部坐标ξ＝x/l的二阶偏导数，是自定义矩阵，它是与单元形函数有关的量，无具体含义，是的转置矩阵；

质量矩阵M^e：

式中，ρ是材料密度，A是梁单元截面积，l是梁单元长度，S是单元形函数。

其中，步骤三中建立的复合材料伸展臂耦合热-结构动力学求解模型，包括：伸展臂传热模型和结构动力学模型，具体表示为：

采用Fourier温度单元建立伸展臂传热求解模型：

式中，T⁰(t)，T¹(t)分别是节点平均温度向量和摄动温度向量，C是热容矩阵，K⁰，K¹分别是0谐和1谐热传导矩阵，R(T⁰(t))是与节点平均温度相关的热辐射系数矩阵，Q⁰(t,q)，Q¹(t,q)分别是0谐和1谐热流载荷向量，q是t时刻节点坐标向量。

采用考虑热效应的ANCF梁单元建立大转动伸展臂结构动力学求解模型：

式中，M是质量矩阵，D是阻尼矩阵，K_l是纵向刚度矩阵，K_t1是横向结构刚度矩阵，P(T,q,t)是载荷向量，Φ是代数约束方程，Φ_q是约束方程的Jacobian矩阵，λ是Lagrange因子向量，q，和分别是t时刻节点坐标向量，节点速度向量和节点加速度向量。

本发明与现有技术相比的优点在于：本发明提供了一种大转动复合材料伸展臂热致振动预测方法，主要优点如下：

(1)本发明采用绝对节点坐标法描述结构变形，可以准确描述具有大转动、大变形结构的动力学响应，不仅避免了坐标变换，而且由于质量矩阵是常数矩阵，因而使惯性力的计算变得十分容易。

(2)本发明建立的ANCF梁单元考虑了热效应引起的结构应变，可以有效进行伸展臂耦合/非耦合热-结构动力学问题分析，为空间伸展臂热致振动研究提供了重要的分析方法。从某种意义上而言，弥补了复合材料伸展臂热致振动地面试验研究的不足。

(3)由于绝对节点坐标方法对结构大转动、大变形描述的突出优点，本发明方法可以进一步延伸应用到多体柔性复合材料结构的热致振动分析中，为未来复杂伸展臂结构设计提供有效的指导。

附图说明

图1是本发明大转动复合材料伸展臂热致振动预测方法实现流程图；

图2是本发明复合材料伸展臂耦合热-结构分析模型；

图3是本发明复合材料伸展臂圆形截面示意图；

图4是本发明两节点ANCF欧拉-伯努利(Euler-Bernoulli)梁单元；

图5是本发明实例中大转动复合材料伸展臂结构；

图6是本发明实例中在不同太阳入射角情况下伸展臂顶端相对位移的时间历程；

图7是本发明实例中在不同转动角速度情况下伸展臂顶端相对位移的时间历程。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

如图1所示，本发明提出了一种大转动复合材料伸展臂热致振动预测方法，其具体实现步骤是：

(1)基于空间伸展臂的结构特点，采用周向均匀刚度配置的铺层方式建立复合材料伸展臂耦合热-结构动力学控制方程，包括：伸展臂传热控制方程和结构动力学控制方程。

航天器伸展臂结构通常由单根或多根复合材料薄壁管组成，在受到空间辐射热流作用下，伸展臂结构的耦合热-结构分析模型如图2所示。根据复合材料薄壁梁的结构特点，可以得到考虑热传导和热辐射的伸展臂传热控制方程，具体表示为：

式中，是伸展臂上任意点的绝对温度；k,ρ和c分别是材料的导热系数，密度和比热；ε和α分别是伸展臂外表面的发射率和吸收系数；σ是Boltzmann常数；δ是与角度相关的量，用于判断伸展臂外表面是否受到太阳辐射；S₀是太阳辐射热流；β是太阳辐射热流入射角；θ_z是伸展臂的转角位移；h是伸展臂薄壁厚度。

图3给出了具有圆形截面的复合材料薄壁梁模型。其中，sxn是梁上定义的局部坐标系，局部坐标系的x轴沿着梁的轴向，sn定义的平面垂直于梁的轴线，s轴沿着薄壁中面的切向，n指向薄壁中面的法向。

采用周向均匀刚度配置的铺层方式，要求铺层角满足ψ(y)＝ψ(-y)，ψ(z)＝ψ(-z)，此时复合材料薄壁梁的位移场表示为：

式中，θ_y(x,t)，θ_z(x,t)和φ(x,t)分别为绕y轴和z轴的弯曲以及绕x轴的扭转。

复合材料薄壁管梁的本构方程表示为：

式中，N_xx，N_sx和N_xn是横截面上单位宽度上的内力，L_xx和L_sx是横截面上单位宽度上的内力矩，是热引起的内力。K_ij(i≠3,j≠5)是局部修正刚度系数。γ_sx为应变分量，这些应变分量可通过如下应变-位移关系得到：

根据式(2)、(3)和(4)，可以得到复合材料薄壁梁的弹性应变能，具体表示为：

根据动能的表达式，复合材料薄壁梁的动能具体表示为：

式中，r＝[u,v,w]^T表示梁上任意一点的位置矢量，则为相应的速度矢量。

根据哈密顿(Hamilton)变分原理，可以得到复合材料伸展臂结构动力学控制方程，具体表示为：

式中，a₃₃，a₅₅是刚度系数，b₁，b₄和b₁₄是惯性系数，m_z是热应力偶，m_z′是热应力偶关于位置坐标x的一阶导数，v₀是伸展臂横向位移，v₀′，v₀″分别是伸展臂横向位移关于位置坐标x的一阶导数和二阶导数，是伸展臂横向加速度，θ_z是伸展臂转角位移，θ_z′，θ_z″分别是伸展臂转角位移关于位置坐标x的一阶导数和二阶导数，是伸展臂角加速度。

(2)基于欧拉-伯努利(Euler-Bernoulli)梁理论和绝对节点坐标方法，建立考虑热效应的ANCF梁单元。

对于一个两节点欧拉-伯努利梁单元，采用绝对节点坐标方法表示梁单元的节点，如图4所示。梁单元的节点坐标含有位置矢量和斜率矢量，具体表示为：

式中，r_k＝[r_k1 r_k2]^T,(k＝i,j)是单元节点位置矢量，是单元节点斜率矢量；

梁单元上任意一点的位置矢量由单元节点坐标表示为：

r(t)＝[r₁ r₂]^T＝S(x)q(t) (9)

式中，r₁和r₂分别是位置矢量r在x和y方向的分量，q(t)是梁单元节点坐标向量，S(x)是梁单元的形函数，它是关于单元局部坐标的矩阵函数，具体表示为：

S(x)＝[s₁I s₂I s₃I s₄I] (10)

式中，I是2×2的单位矩阵，s_i(i＝1,2,3,4)是形函数分量，l是梁单元的长度，ξ＝x/l，ξ的取值范围为[0,1]。

基于连续介质力学理论，欧拉-伯努利梁单元的弹性应变能含有纵向应变能和横向应变能。引入热效应的梁单元的弹性应变能表示为：

式中，U_l，U_t分别为梁单元纵向应变能和横向应变能，E是材料的弹性模量，A是梁单元截面积，l是梁单元长度，V是梁单元体积，ε_l和分别为纵向结构应变和热应变，ε_t和分别为横向结构应变和热应变。纵向热应变表示为：

式中，α_T是材料的热膨胀系数，是梁横截面上的瞬时平均温度，T_ref是参考温度。

航天器伸展臂通常选用热膨胀系数较小的复合材料制成，并且伸展臂横截面上的温度梯度不会很大，由此引发的横向热应变可视为小量，忽略横向热应变的二阶项，式(12)中的横向应变能进一步表示为：

式中，κ是梁轴线变形后的曲率，M^T是梁单元的热弯矩，它们分别定义为：

式中，r_,x，r_,xx分别为位置矢量关于局部坐标x的一阶和二阶偏导数，T是梁横截面上的瞬时温度，它是局部坐标x，y和时间t的函数。

式(12)中的弹性应变能关于节点坐标向量求偏导，可以得到考虑热效应的ANCF梁单元的弹性力，具体表示为：

梁单元的纵向应变ε_l和变形后的曲率κ可由纵向变形梯度f分别表示为：

式中，是变形梯度，并且有：

式中，S_,x是单元形函数关于坐标x的一阶偏导数，S_,ξ和S_,ξξ是单元形函数关于局部坐标ξ＝x/l的一阶和二阶偏导数，是自定义矩阵，它是与单元形函数有关的量，无具体含义。

将式(18)、(19)代入式(17)中，考虑热效应的ANCF梁单元的弹性力进一步表示为：

式中，和分别是梁单元纵向弹性力和横向弹性力。

根据式(23)可知，考虑热效应的ANCF梁单元的刚度矩阵包括纵向刚度矩阵和横向刚度矩阵具体表示为：

纵向刚度矩阵

式中，E是材料的弹性模量，A是梁单元截面积，l是梁单元长度，是梁单元纵向热应变，q是梁单元节点坐标，S_l是与单元形函数有关的量。

横向刚度矩阵

式中，是横向结构刚度矩阵，是热弯矩引起横向刚度矩阵，I是梁单元横截面惯性矩，M^T是梁单元的热弯矩，S_t是与单元形函数有关的量。

根据欧拉-伯努利梁单元的动能表达式，可以得到：

根据式(27)可知，考虑热效应的ANCF梁单元的质量矩阵M^e为常数矩阵，具体表示为：

式中，ρ是材料密度，A是梁单元截面积，l是梁单元长度，S是单元形函数。

(3)采用傅里叶(Fourier)温度单元和ANCF梁单元建立复合材料伸展臂耦合热-结构动力学求解模型，包括：伸展臂传热模型和结构动力学模型。

采用傅里叶温度单元对步骤(1)建立的伸展臂传热控制方程式(1)进行离散，可以得到复合材料伸展臂传热求解模型，具体表示为：

式中，T⁰(t)，T¹(t)分别是节点平均温度向量和摄动温度向量，C是热容矩阵，K⁰，K¹分别是0谐和1谐热传导矩阵，R(T⁰(t))是与节点平均温度相关的热辐射系数矩阵，Q⁰(t,q)，Q¹(t,q)分别是0谐和1谐热流载荷向量，q是t时刻节点坐标向量。

采用步骤(2)建立的考虑热效应的ANCF梁单元对步骤(1)建立的伸展臂结构动力学控制方程式(7)进行离散，可以得到复合材料伸展臂结构动力学求解模型，具体表示为：

式中，M是质量矩阵，D是阻尼矩阵，K_l是纵向刚度矩阵，K_t1是横向结构刚度矩阵，P(T,q,t)是载荷向量，T＝T⁰+T¹是节点温度。

对于大范围运动的伸展臂结构，需要施加运动约束方程。使用增广计算将大转动伸展臂运动方程表示为指标-3的微分-代数方程形式：

式中，Φ是代数约束方程，Φ_q是约束方程的Jacobian矩阵，λ是Lagrange因子向量，q，和分别是是t时刻节点坐标向量，节点速度向量和节点加速度向量。

(4)在时间t_n上，采用克兰克-尼克尔森(Crank-Nicolson)方法求解步骤(3)建立的复合材料伸展臂传热模型式(29)，得到复合材料伸展臂温度结果

(5)根据步骤(4)得到的温度结果，在时间t_n上，采用广义-a(Generalized-a)方法求解步骤(3)建立的复合材料伸展臂结构动力学模型式(31)，得到复合材料伸展臂位移结果q_n+1。

(6)判断时间条件是否满足t_n+1≥t_tol，如果不满足条件，令t_n＝t_n+Δt，跳转到步骤(4)，其中，Δt是时间步长；如果满足条件，输出复合材料伸展臂位移响应时间历程q(t)，迭代终止。

实施例：

1.结构参数及模型介绍

为了更充分的了解该发明的特点及其对工程实际的适用性，本发明以图5所示的大转动伸展臂结构为例说明本发明方法的有效性。图5中伸展臂是周向均匀刚度配置的复合材料薄壁管结构，其结构及材料参数列于表1。太阳热流S₀以入射角β突然作用于伸展臂。假定空间黑体辐射温度为0K，伸展臂初始温度为290K。伸展臂通过根部的转动约束使其带动顶端设备做大范围运动，转动约束规律表示为:

式中，T_s＝15s，ω_s＝3rad/s。

表1复合材料伸展臂属性参数

2.大转动伸展臂热致振动响应分析

采用本发明方法对大转动伸展臂的热致振动响应进行求解分析。当伸展臂的转动速度为ω_s＝3rad/s时，首先比较了太阳热流以两种入射角照射时伸展臂顶端位移动态响应，如图6所示。可以看到，在转动初始阶段，伸展臂顶端会产生较大的结构变形，随后变形逐渐恢复，约在15s以后，顶端位移响应基本达到稳定状态。对比两种情况的位移响应曲线发现，入射角为β＝0°时的顶端位移大于入射角为β＝60°时的顶端位移，并且在稳定阶段出现明显的振动现象。此外，图6中也给出了不考虑热载荷作用时伸展臂顶端位移响应。与考虑热载荷的情况相比，不考虑热载荷作用时伸展臂顶端位移响应明显变小，这是因为忽略了热载荷引起的热弯矩对伸展臂的影响。因此，对于复杂空间环境下的伸展臂，考虑热载荷作用获得的结构响应更加符合真实情况。

当太阳热流入射角为β＝0°时，图7给出了两种角速度时伸展臂顶端相对位移的时间历程。由图可见，当伸展臂转速提高时，顶端位移发生很大变化，在稳定阶段，振动现象更加剧烈。为了清晰显示两种情况下的位移差别，选取了稳定阶段几个典型时刻的顶端相对位移进行比较，如表2所示。综合图6、图7和表2的比较结果可知，相比于太阳热流入射角，转动速度对稳定后的位移响应影响更大。由此可见，合理选择转动速度有助于抑制伸展臂稳定后的振动现象。

表2稳定时顶端位移比较

针对大转动复合材料伸展臂结构，采用本发明方法有效的预测了伸展臂在空间辐射环境下的热-结构动力学响应。以上实例验证了本发明方法对大转动复合材料伸展臂热致振动响应预测的可行性和有效性。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制。

本发明未详细阐述部分属于本领域技术人员的公知技术。

Claims (1)

1.一种大转动复合材料伸展臂热致振动预测方法，其特征在于实现步骤如下：

步骤一、基于空间伸展臂的结构特点，采用周向均匀刚度配置的铺层方式，建立复合材料伸展臂耦合热-结构动力学控制方程，包括：伸展臂传热控制方程和结构动力学控制方程；

步骤二、基于欧拉-伯努利(Euler-Bernoulli)梁理论和绝对节点坐标方法(ANCF)，建立考虑热效应的ANCF梁单元；

步骤三、采用傅里叶(Fourier)温度单元和ANCF梁单元离散步骤一建立的耦合热-结构动力学控制方程，建立复合材料伸展臂耦合热-结构动力学求解模型，包括：伸展臂传热模型和结构动力学模型；

步骤四、在时间t_n上，采用克兰克-尼克尔森(Crank-Nicolson)方法求解步骤三建立的伸展臂传热模型，采用牛顿-拉夫逊(Newton-Raphson)方法求解非线性方程组，得到复合材料伸展臂温度结果

步骤五、根据步骤四得到的温度结果，在时间t_n上，采用广义-a(Generalized-a)方法求解步骤三建立的伸展臂结构动力学模型，采用牛顿-拉夫逊(Newton-Raphson)方法求解非线性方程组，得到复合材料伸展臂位移结果q_n+1；

步骤六、判断时间条件是否满足t_n+1≥t_tol，如果不满足条件，令t_n＝t_n+Δt，跳转到步骤四，其中，Δt是时间步长；如果满足条件，输出复合材料伸展臂位移响应时间历程，迭代终止；

其中，步骤一中采用周向均匀刚度配置的铺层方式建立复合材料伸展臂结构动力学控制方程：

式中，a₃₃，a₅₅是刚度系数，b₁，b₄和b₁₄是惯性系数，m_z是热应力偶，m′_z是热应力偶关于位置坐标x的一阶导数，v₀是伸展臂横向位移，v′₀，v″₀分别是伸展臂横向位移关于位置坐标x的一阶导数和二阶导数，是伸展臂横向加速度，θ_z是伸展臂转角位移，θ′_z，θ″_z分别是伸展臂转角位移关于位置坐标x的一阶导数和二阶导数，是伸展臂角加速度；

步骤二中采用位置矢量和斜率矢量表示考虑热效应的ANCF梁单元的节点坐标，具体表示为：

式中，r_k＝[r_k1 r_k2]^T,k＝i,j是单元节点位置矢量，是单元节点斜率矢量；

梁单元上任意一点的位置矢量由单元节点坐标表示为：

r(t)＝[r₁ r₂]^T＝S(x)q(t)

式中，r₁和r₂分别是位置矢量r在x和y方向的分量，q(t)是梁单元节点坐标向量，S(x)是梁单元的形函数，它是关于单元局部坐标的矩阵函数；

步骤二中建立的考虑热效应的ANCF梁单元含有纵向热应变和横向热应变，梁单元的刚度矩阵包括纵向刚度矩阵和横向刚度矩阵梁单元的质量矩阵M^e为常数矩阵，其中：

纵向刚度矩阵具体表示为：

式中，E是材料的弹性模量，A是梁单元截面积，l是梁单元长度，是梁单元纵向热应变，q是梁单元节点坐标，S_l是与单元形函数有关的量，具体表示为：

式中，S_,x是单元形函数关于坐标x的一阶偏导数，是S_,x的转置矩阵，S_,ξ是单元形函数关于局部坐标ξ＝x/l的一阶偏导数，S_,ξ ^T是S_,ξ的转置矩阵，l是梁单元长度；

横向刚度矩阵具体表示为：

式中，是横向结构刚度矩阵，是热弯矩引起横向刚度矩阵，I是梁单元横截面惯性矩，M^T是梁单元的热弯矩，S_t是与单元形函数有关的量，具体表示为：

式中，S_,ξξ是单元形函数关于局部坐标ξ＝x/l的二阶偏导数，是自定义矩阵，它是与单元形函数有关的量，无具体含义，是的转置矩阵；

质量矩阵M^e，具体表示为：

式中，ρ是材料密度，A是梁单元截面积，l是梁单元长度，S是单元形函数；

步骤三中采用Fourier温度单元建立伸展臂传热求解模型，采用考虑热效应的ANCF梁单元建立伸展臂结构动力学求解模型，其中：

伸展臂传热求解模型，具体表示为：

式中，T⁰(t)，T¹(t)分别是节点平均温度向量和摄动温度向量，C是热容矩阵，K⁰，K¹分别是0谐和1谐热传导矩阵，R(T⁰(t))是与节点平均温度相关的热辐射系数矩阵，Q⁰(t,q)，Q¹(t,q)分别是0谐和1谐热流载荷向量，q是t时刻节点坐标向量；

大转动伸展臂结构动力学求解模型，具体表示为：

式中，M是质量矩阵，D是阻尼矩阵，K_l是纵向刚度矩阵，K_t1是横向结构刚度矩阵，P(T,q,t)是载荷向量，T＝T⁰+T¹是节点温度，Φ是代数约束方程，Φ_q是约束方程的Jacobian矩阵，λ是Lagrange因子向量，q，和分别是t时刻节点坐标向量，节点速度向量和节点加速度向量。