热-动力学耦合多体系统建模与降阶1)

2022-10-05田青龙於祖庆陆念力

力学学报 2022年9期

关键词：薄板动力学模态

田青龙於祖庆兰朋,,,2) 陆念力

* (哈尔滨工业大学机电工程学院,哈尔滨 150001)

† (河海大学机电工程学院,江苏常州 213002)

** (西安建筑科技大学机电工程学院,西安 710055)

†† (西安建筑科技大学西部绿色建筑国家重点实验室,西安 710055)

引言

太阳能电池板和抛物面天线等典型的航天器板结构具有重量轻、尺寸大的特点.此类航天器的动力学分析与控制问题越来越多地需要考虑到空间环境温度的影响.由太阳辐射导致的极端热载荷激励引起位移场与温度场的强耦合响应,进而导致热致振动或热跳变[1-4].因此,刚柔热耦合系统的精确建模对航天器的安全至关重要.

早期的热动力学耦合研究较多使用浮动坐标方法(FFRF)[5-6],因而其仅限于柔性体小变形问题.绝对节点坐标方法(ANCF)则采用了全局坐标系下定义的位置与梯度向量作为节点坐标,被广泛用于大范围运动和大变形的多柔体系统动力学的研究[7].Wu等[8]在假设梁截面的温度分布为线性的情况下,研究了瞬态热载荷对细长梁结构的影响.Abbas 等[9]使用变厚度ANCF 板单元研究了热载荷下的气动热弹性动力学.Čepon 等[10]建立了基于二维ANCF 梁单元的双金属片模型.在上述研究中,热载荷均以解析解的形式给出.Li 等[11]建立了热柔耦合的大柔性复合太阳能电池板平面航天器模型,并对其进行了动力学和热耦合分析.于香杰等[12]结合负应变率控制法建立了变截面梁构件在空间热载荷作用下的压电抑振抑变结构.邢晓峰等[13]研究了刚-柔耦合、耦合热弹性、耦合热-结构三重耦合效应对航天器结构热致振动响应及稳定性的影响.Shen 等[14-15]研究了薄壁管和多层材料太阳能电池板的热致振动.温度场采用线性插值网格进行离散且传热方程与动力学方程同步求解.这一阶段的研究均采用两种不同的离散格式分别对柔性体位移场和温度场进行插值,会导致数值求解中的网格错配和映射误差问题[16].为此,Cui 等[17-18]将ANCF 单元形状函数应用于位移场和温度场,自然地将传热和连续介质力学统一为一个单元网格.这一思想也将在本研究中得到应用.

随着ANCF 计算规模增大,数值效率低的问题日渐突出,模型降阶成为了重要的研究方向.在文献[19-20]的研究中,本征正交分解(POD)方法被用于ANCF 单元的模型降阶.然而,降阶模型必须通过实验数据库或计算全阶模型仿真结果来获得.模态综合法(CMS) 可以避免这一缺点.Kobayashi 等[21]在轴向变形较小的前提下,使用固定界面模态综合法(也称Craig-Bampton 方法)来缩减ANCF 平面梁单元建模的多体系统.Sun 等[22-23]使用自由界面模态综合法来提高带滑动铰的柔性梁的计算效率.该方法同样采用了相同的轴向小变形的假设.Tang 等[24]将整个柔性体运动过程划分为若干子区间并在每个子区间内对系统方程进行线性化处理,从而使得切线刚度阵定常,基于模态截断的降阶方法得以应用.此外作者也采用了Craig-Bampton 方法对ANCF 梁单元描述的多柔体系统进行子结构划分,称之为CB-ANCF 方法.在此基础上,作者将此方法推广到了具有阻尼的多柔体系统,以减少黏弹性多柔体系统的自由度[25].基于文献[24-25] 的研究,Tian等[26]使用自由界面模态综合法来减少ANCF 梁单元描述的多柔体系统的自由度,即F-ANCF 方法.

与动力学方程不同,传热方程为一阶偏微分方程.很多研究表明,传统的动力学降阶方法如模态叠加法[27]、模型识别法[28]、POD 法[29]等依然适用于传热方程的降阶.然而,对于热柔耦合系统的整体降阶研究却较少见于文献.Nachtergaele 等[30]提出了一种弱耦合投影基用于热柔耦合系统的降阶,这是对Craig-Bampton 方法的扩展,温度场隐含在投影基中.Yamashita 等[31]也使用该方法对基于FFRF 建模的热柔耦合系统进行了降阶研究.

本文提出了一种刚柔热耦合系统模型降阶的方法.引入温度梯度,采用高阶插值对温度场进行离散,使位移场和温度场采用统一的ANCF 单元形函数进行插值离散,避免数值求解中网格错配和误差映射的问题.ANCF 参考节点用于描述中心刚体.采用泰勒展开对动力学方程和传热方程进行线性化处理,得到常值切线刚度阵.引入改进的模态综合法分别对动力学方程和传热方程进行自由度缩减.最后,以刚柔热耦合大口径抛物面天线为例验证了该方法的有效性.

1 热-动力学耦合建模

本文拟分别采用热耦合的ANCF 薄板[18]与三维梁单元[17]分别构造星载抛物面天线及其肋梁模型.而卫星本体视为刚体,由ANCF 参考节点[32]描述.由此可将卫星刚柔热耦合系统纳入ANCF 框架内描述和建模.由于采用了全局坐标系下定义的位置与梯度向量,单元内部任意点的空间位置可以定义为r=Se,S为单元形函数矩阵,e为单元节点向量.ANCF 薄板单元基于Kirchhoff-Love 板理论,其节点向量包含中面四个角点的位置和平面内梯度向量,即其中上标P 代指板单元,i=1,2,3,4 为节点号.ANCF 三维梁单元节点为中轴线两端点,每个节点使用全局位置向量和沿三个方向的梯度向量作为节点坐标,可以表示为:eB=上标B 代指梁单元,j=1,2.

ANCF 参考节点是一种具有完备梯度向量但不属于任何单元的孤立节点.它通过连续性条件与柔性体网格相连,并通过非线性约束的形式保证其梯度向量模长始终为1 并且正交.因此它可以在ANCF的表达范式下高效地描述刚体的动力学行为.系统动力学方程的形式不因其的加入而发生改变.现有的研究表明,无论柔性体离散使用全参数单元还是梯度不完备单元,均不影响ANCF 参考单元描述刚体运动的正确性.因此,考虑热应力的系统动力学方程为

其中,M为系统质量矩阵,Qe是弹性力.薄板单元和梁单元的弹性力分别在文献[33-34]中给出.QH是对应于热应变的广义力,表达式在附录中给出.Ce为约束方程,∂Ce/∂e则为约束方程对广义坐标e的偏导数,λ 是拉格朗日乘子,Qext代表广义外力.

由于连续体内温度场采用了与运动学描述相同的离散格式,单元节点处的温度及其梯度被用作传热方程的广义坐标,如图1 所示.薄板单元的温度广义坐标为其中Ti(i=1,2,3,4)为节点i的温度,Txi和Tyi分别表示节点i在x和y方向上的温度梯度.同理,三维梁单元的温度广义坐标可以写成

图1 ANCF 热耦合单元Fig.1 ANCF thermal coupling element

与位置插值同理,单元内任意点的温度可以写为

其中,ST是温度场的单元形状函数,其元素与描述位移场形状函数相同,eT表示单元的节点广义温度坐标向量.同样地,上角标P 和B 分别指代薄板单元和梁单元.

在太空中,影响航天器温度变化的主要原因是太阳热流输入和表面自热辐射,因此在考虑热的边界条件的情况下,传热方程的一般形式如下

其中,DT是热容矩阵,Kc是热传导矩阵,Kr是表面辐射矩阵[17-18].eT∞代表热辐射环境温度,在太空中自热辐射环境温度近似为0 K.QT是外热载荷,与热边界条件有关.当受到太阳热流照射时,外热载荷通常与柔性体当前构型和太阳热流入射角度有关:QT=∫H·ndS,H=H0nT,n为受热面法线方向,nT为热流入射方向,H0为太阳热流密度.CT传热方程的约束方程. ∂CT/∂eT是传热约束方程对广义温度坐标eT的偏导数,λT是传热约束的拉格朗日乘子.

2 动力学和传热方程的线性化

基于ANCF 建立的柔性体动力学方程具有常值质量阵的优势,但其弹性力却是节点坐标向量的高度非线性函数,切线刚度阵亦然.这就导致基于模态截断的降阶方法无法直接应用.为此,将系统的整个运动过程划分为多个区间,认为柔性体在每个区间内部的动力学行为是准静态的.利用一阶泰勒展开将基于初始构型定义的全量平衡方程变换为基于当前区间起始构型的增量平衡形式.每个区间内的切线刚度矩阵可视为定常,进而采用改进的模态综合法对系统进行降阶.

根据泰勒展开式,每个准静态区间的初始构型处的动力学方程可展开为

其中,e0是每个准静态区间初始构型处的广义坐标.同理,eT0是初始广义温度坐标.JQe是弹性力雅可比矩阵(表达式见文献[33-34]).JQH为热应变广义力的雅可比矩阵,其表达式在附录中给出.约束方程也可以线性化为

将式(4)和式(5)代入式(1)中,线性化的动力学方程形式可以写成

其中,K0是广义刚度矩阵,F0是线性化的广义外力,表达式为

用同样的方法,传热方程可以展开如下

传热约束方程可展开为

将式(7)和式(8)代入式(3)中,线性化的传热方程形式可以写成

式中,FT0是线性化的广义热载荷.

3 改进的模态综合法

根据子结构界面约束的不同,模态综合方法可分为固定界面模态综合法[24]和自由界面模态综合法[35].如图2 所示,多体系统ABCD被连接节点EF分为两个子结构.对于固定界面模态综合法,子结构之间连接的自由度被全部保留.然后固定每个子结构的界面自由度,并分别对子结构进行模态分析.最后,子结构之间通过共享边界节点的形式连接起来,如图2(b)所示.而自由界面模态综合法则无需保留连接界面自由度,即每个子结构的模态都包含了刚体位移信息.结构的组装通过界面位移和界面力相等的双协调条件实现,如图2(c)所示.而本文在自由界面模态综合法的基础上提出了一种改进的模态综合法.通过施加界面位置和梯度的约束方程保证子结构间的连接精度,如图2(d)所示.

图2 系统子结构划分Fig.2 The mesh of the system and substructures

在动力学方程部分,子结构j的动力学方程齐次形式如下

其中,[]j表示子结构j.通过模态分析获得当前子结构的特征值和特征向量

其中,ωn是第n个固有频率,[Φ]j为子结构振型,是子结构j的自由度数.

子结构j的广义坐标增量可以表示为

以两个子结构组成的系统为例,其广义坐标可以表示为

显然,每个子结构的坐标是独立的,除了连接处的自由度.在本研究中,直接施加界面位置和梯度的约束方程,以确保子结构连接处的精度,即

将式(14)和式(15)代入式(6),系统动力学方程可以表示为

其中N=diag([Φk]1,[Φk]2)

传热方程是一阶微分方程.子结构j的特征值和特征向量可以通过热特征值分析获得

式中,ωm为第m热模态频率,物理意义为热衰减系数[27],[Ψ]j是对应于热特征值的特征向量是子结构j的温度坐标数.传热方程的特征值的求解在文献[36]中给出.

子结构j的广义温度坐标增量可以表示为

直接施加界面连接处温度和温度梯度的约束方程

然后,将式(21)和式(22)代入式(9),系统传热方程可以表示为

4 求解策略

为准确模拟热-动力学耦合效应对系统的影响,需要同步求解传热方程和动力学方程.本文采用广义-α 积分来求解传热方程和动力学方程.图3 给出了算法的流程图.算法的数值参数 αf,αm,β ,γ 等的取值见参考文献[37-38].

图3 算法求解流程Fig.3 The overall flow-chat of the computation algorithm

在第2 节中,使用泰勒展开法将动力学方程和传热方程线性化,将非线性方程转化为线性方程.需要指出的是,线性方程等式成立的前提是广义位置坐标的增量应该很小.当第j个子结构的最大位移增量过大时,第j子结构应重新线性化.同样,温度的升高也会引起热应变的变化,同样会影响切线刚度阵,因此,当第j子结构的温度增量过大时,第j子结构也应重新线性化.保证 δk和 δT足够小,将累计误差控制在一定范围内.本文给出阈值的经验计算公式

其中L是梁结构的长度或者薄板结构长度.通常,对于动力学方程,保留的自由度nred建议如下nred=ntal/3～4ntal/5.对于传热方程,保留的自由度建议如下

需要指出的是,t=0 时的初始模态速度和加速度均为0,任意子结构j第一次缩减保留的模态为前k阶低阶固有频率对应的主模态.当 ‖Δr‖2＞δk或者 ‖ΔT‖2＞δT时,动力学方程和传热方程需要重新线性化,更新系统的切线刚度阵.之后,用新的切线刚度阵和质量阵去重新进行子结构模态分析,然后根据模态速度的绝对值大小排序,保留模态速度最大的前k阶对应的主模态.更多细节可参见文献[26].

5 数值算例

本节给出了四个数值算例来验证所提出方法的准确性和有效性.为方便起见,完整自由度模型记为FOM-TANCF,降阶模型记为ROM-TANCF.

5.1 半圆薄板的纯导热

在半圆形薄板上测试纯导热问题,一侧受到热流输入,并考虑了表面自热辐射.如图4(a)所示,将半圆板划分为四个ANCF 热薄板单元.图4(b)是参考文献[14]中给出的一维线性传热单元,将半圆划分十个单元.总仿真时间设置为80 s,初始温度为290 K,热流为1350 W/m2.其他参数见表1.

图4 半圆板传热模型Fig.4 The heat conduction of a cylinder thin plate

表1 半圆板模型参数Table 1 Parameters of the semicircular thin plate

由于半圆形结构简单,因此不进行子结构划分.图5 给出了FOM-TANCF、线性传统有限元(FEM)和解析解的结果.可以看出,本文给出的全自由度模型分析结果与传统有限元结果较为吻合.而解析解[2]由低阶傅里叶展开得到,温度按照三角函数沿圆周分布,故而与数值解之间存在一定误差.图6 给出了FOM-TANCF 和ROM-TANCF 的温度分布.图7 给出了FOM-TANCF 和ROM-TANCF 之间的温度误差.温度均方根误差(TRMS),用来评估所提出的缩减后的模型和全自由度模型之间温度的误差,表达方式如下

图5 不同时刻圆周温度分布对比Fig.5 Comparison of section temperature distribution at different times

其中n是节点个数.

温度坐标的总自由度为30,ROM-TANCF 中保留了10 个温度坐标.时间步长为0.01 s.计算时间从321.3 s 减少到119.8 s,降幅为62.71%.如图5所示,FOM-TANCF 和FEM 的结果一致.因此,FOM-TANCF 的结果作为基准.如图6 和图7 所示,FOM-TANCF 和ROM-TANCF 的计算结果吻合良好,证明该方法在保证足够精度的前提下能有效提高求解传热方程的效率.

图6 FOM-TANCF 和ROM-TANCF 温度分布对比Fig.6 Comparison of section temperature distribution between FOMTANCF and ROM-TANCF

图7 FOM-TANCF 和ROM-TANCF 温度误差Fig.7 Error of temperature between FOM-TANCF and ROM-TANCF

5.2 薄板热膨胀

本节提出了一个长方形薄板模型来测试热膨胀.如图8(a)所示.将系统在长度方向上划分为6 个热柔耦合薄板单元.系统分为两个子结构,每个子结构划分为3 个单元,如图8(b)所示.初始温度为0 K.总仿真时间为6 s.其中前5 s 是恒定的内热,最后1 s 内热停止,其他参数见表2.均方根误差(RMS),用来评估所提出的缩减后的模型和全自由度模型之间的误差,表达方式如下

图8 薄板热膨胀模型Fig.8 Thermal deformation of thin plate mode

表2 薄板热膨胀板模型参数Table 2 Parameters of the thermal expansion thin plate

给定的阈值 δk和 δT由式(24)获得,分别为2.5 mm和0.01 K.时间步长为0.1 ms.计算时间从5089.7 s减少到1422.5 s,减少了72.05%.动力学自由度从126 个减少到70 个,其中子结构1 保留26 个,子结构2 保留26 个,子结构连接边界约束分别为18 个.为了方便起见,可以将记为70 (26+26+18).同样,对于传热方程,自由度从42 个减少到26 个(10+10 +6).采用ANSYS 中的SOLID90单元(网格划分为200×50×2),采用间接耦合求解.如图9 和10 所示,FOM-TANCF,ROM-TANCF 和ANSYS 的结果一致.图11 中给出了RMS 误差.计算结果表明,该方法能在有效提高求解热膨胀问题的计算效率的同时保证较高的计算精度.

图9 检测点温度Fig.9 The test point temperature

图10 检测点在长度方向位置Fig.10 The test point position in length direction

图11 FOM-TANCF 和ROM-TANCF 误差Fig.11 Error of deformation between FOM-TANCF and ROM-TANCF

5.3 柔性太阳能电池板

如图12(a)所示,最初水平放置的柔性太阳能电池板一端铰接,并在微重力的作用下下落.太阳能电池板由中间的薄板结构和两侧的梁结构组成.薄板单元和梁单元节点通过约束方程连接.太阳热流水平射入,初始时不受热流照射,随着微重力下摆时,受热面积逐渐增大,温度逐渐升高.薄板结构的网格为4×2 (长×宽).每个梁离散为8 个圆形截面ANCF梁单元[39].系统分为两个子结构,如图12(b)所示.总模拟时间设置为10 s.初始温度为0 K,柔性太阳能电池板的其他参数如表3 所示.

图12 柔性太阳能电池板模型Fig.12 The flexible solar panel model

表3 柔性太阳能电池板模型参数Table 3 Parameters of the flexible solar panel

给定的阈值 δk和δT分别为1 mm 和0.01 K.时间步长为0.2 ms.计算时间从30 363.6 s 减少到9293.2 s,减少了69.39%.动力学的自由度从351 个减少到227 个(100+100+27).对于传热方程,自由度从117 个减少到75 个(33+33+9).

X方向上的测试点位置如图13 所示.梁和薄板的热膨胀系数不同,如表3 所示.由于大变形和旋转,柔性太阳能电池板的下摆导致不均匀温升和热应变.与没有温度变化的情况相比,不均匀的热应变会导致测试点轨迹的不同.RMS 误差如图14 所示,FOM-TANCF 和ROM-TANCF 的计算结果一致.图15给出了由ROM-TANCF 获得的柔性太阳能电池板在不同时间的运动过程和温度分布.计算结果表明,该方法能有效提高大变形、大转动的热柔耦合系统的效率.

图13 检测点X-方向位置对比Fig.13 Comparison of test point position in X-direction

图14 FOM-TANCF 和ROM-TANCF 误差Fig.14 Error of deformation between FOM-TANCF and ROM-TANCF

图15 柔性太阳能电池板不同时刻构型Fig.15 Motion process of the flexible solar panel

5.4 刚柔热耦合抛物线天线

本节使用刚柔热耦合天线验证所提出算法的有效性.刚柔热耦合天线模型如图16 所示.采用薄板单元模拟抛物面天线,三维梁单元模拟抛物面天线上的肋梁.两种单元在公共节点处通过约束方程连接.如图16 所示,柔性天线分为两个子结构.ANCF参考节点用于建模中心刚性轮毂.中心轮毂的质量为300 kg,中心轮毂绕主轴的转动惯量分别为Jxx=100 kg·m2,Jyy=100 kg·m2,Jzz=500 kg·m2.太阳热流的入射方向为[5,0,-2].仿真时长为50 s,初始温度为0 K,其他参数见表4.驱动扭矩施加在中心刚性轮毂上,驱动力矩表达式如下

表4 刚柔热耦合天线模型参数Table 4 Parameters of the rigid-flexible-thermal coupling satellite

图16 刚柔热耦合天线模型Fig.16 The rigid-flexible-thermal coupling satellite schematic model

阈值 δk和 δT分别为3 mm 和0.01 K.时间步长为0.1 ms.测试点A和B在Y方向上的位置分量如图17 所示,温度如图18 所示.RMS 误差如图19 所示.不同时间的温度分布如图20 所示.可以看出相较于直径达10 m 的天线,RMS 最大值不超过3 cm,与此同时,计算时间从337 303.2 s 减少到67 466.5 s,减少了79.99%.动力学的自由度从708 个减少到554 个(250+250+54).对于传热方程,自由度从232 个减少到148 个(65+65+18).尽管本文提出的方法需要在位移或温度增量过大时,对子结构进行重新线性化,更新模态基,但是ANCF全自由度模型每个时间步长内弹性力及其雅可比求解都非常耗时,降阶模型从整体上体现出更高的仿真求解效率.

图17 A 和B 点Y 方向位置对比Fig.17 Comparison of test points A and B position in Y-direction

图18 A 和B 点温度Fig.18 Temperature of test points A and B

图19 FOM-TANCF 和 ROM-TANCF 误差Fig.19 Error of deformation between FOM-TANCF and ROM-TANCF

图20 不同时刻温度分布Fig.20 Temperature distribution at different times

6 结论

本文提出了一种基于改进模态综合法的刚柔热耦合多体系统模型降阶方法.ANCF 单元形函数被同时用于描述具有大变形和大范围运动的多柔体系统及柔性体内温度场差值离散.采用ANCF 参考节点用于描述刚体运动,得到了统一描述的刚柔热耦合多体系统模型.采用泰勒展开法对系统动力学和传热方程进行线性化处理,使得线性化区间内,柔性体的切线刚度矩阵可视为常值矩阵.采用了改进的模态综合法对线性化后的动力学方程和传热方程进行降阶.子结构之间通过约束方程连接,以确保连接精度及边界连续性.最后,通过四个数值算例表明,降阶后的模型计算结果与全自由度模型计算结果一致,并可以有效提高计算效率.半圆形薄板的纯热传导算例考虑了热流输入和表面自热辐射.计算结果与以往文献中的结果吻合,计算时间减少了62.71%,保留了33.33%的自由度.薄板的热膨胀算例中,计算结果与ANSYS 计算结果吻合较好,验证了该方法的准确性.对于柔性太阳能电池板和刚柔热耦合天线,计算时间可分别减少69.39%和79.99%.算例结果表明,该方法能有效地提高刚柔热耦合系统的计算效率,并保持较高的精度.此外,该模型降阶方法可以推广到阻尼系统和其他多个物理耦合的领域.

附录

与热应变对应的广义力可以写成