王亮 张妍 蔡毅鹏 蔡文杰

摘 要：针对全大气层滑翔的高超声速飞行器动特性预示需求， 本文提出了考虑防热层以及烧蚀的高超声速飞行器截面刚度计算模型。 首先， 提出了三层组合的高超声速飞行器截面刚度计算模型， 分别为最内层承力金属壳体、 中间防热层原始层和最外防热层烧蚀层。 其次， 根据刚度线性叠加原理， 给出了各层以及总的截面刚度参数的计算方法。 最后， 针对实际使用， 给出了高超声速飞行器结构动特性具体的计算分析程序和步骤。

关键词：     高超声速;  动特性;  热防护系统;  烧蚀效应;  抗弯刚度

中图分类号：     TJ760; V212.11 文献标识码：    A  文章编号：     1673-5048（2021）06-0072-04

0 引  言

高超声速飞行器是结构组成十分复杂的弹性体， 其动力学稳定性设计是飞行器稳定飞行的前提， 而结构动力学特性是飞行器姿态稳定控制系统设计的关键参数之一， 参数的准确性将直接影响飞行成败。 结构动力学特性预示是理论性和实践性较强的一项工作， 因此其重要性不言而喻。

在分析细长体的高超声速飞行器时， 常采用梁单元进行模态分析。 学者们在梁单元精确分析方面进行了大量的研究[1-9]， 提出了较多梁单元刚度精细化建模方法。 如多种截面剪切系数计算方法， 对梁单元中的剪切效应进行有效修正;  包含自重项的细长梁单元刚度矩阵， 其中考虑了轴力对弯曲刚度的影响及弯矩引起的轴向刚度修正系数;  针对管状复合材料梁提出一种复合梁单元的有限元建模新方法。

典型的高超声速飞行器， 如传统的战术弹道导弹， 其为再入式飞行， 虽然再入速度较高， 遭受严酷的气动加热影响， 但由于再入飞行时间较短， 且再入体外部防热层较厚， 因此气动加热效应不明显， 壳体温升和防热层烧蚀严酷度不高。 再入体的结构动特性预示难度较低， 采用常温材料参数进行动特性预示精度较高， 偏差带可以包络由于短时防热层烧蚀和壳体温升带来的动特性变化。 在此前提下， 学者们在研究高超声速飞行器动特性预示时， 一般采用考虑常温下飞行器的模态特性预示方法， 或考虑高温对承力壳体的影响下的模态预示方法。 何成[10]基于结构动力学理论、 热传导理论、 模型修正与参数辨识方法， 以理论分析与仿真试验验证相结合的方式开展热环境下的高超声速飞行器结构动力学的建模方法研究。 王亮等[11-13]以高超声速战术导弹为应用背景， 研究了各类全弹的梁-质量块建模方法， 给出了常温下多种导弹全弹有限元质量模型、 扭转模态模型和横向模态模型， 对比了各类模型对计算结果的影响， 并与试验结果进行了对比， 给出了全弹有限元建模时的推荐模型。 张琪等[14]针对以往方法0维点单元和3维实体单元等效非结构质量中存在的问题， 提出等效密度法， 即根据密度和刚度等效的原则， 将非结构质量等效均布到壳体单元， 修改壳体材料密度以保证其质量不变。 应用3种方法分别建立了导弹有限元模型， 并进行了导弹模态试验， 将3种方法得到的模态频率与试验值比较， 发现等效密度法建立的模型更为合理， 应用简便。

综上所述， 学者们在研究梁单元精细建模和高超声速飞行器常规建模方面做出了很多工作。 但是考虑高超声速飞行工况下， 氣动加热下全飞行器动特性模型建模的研究偏少， 尤其是考虑防热层烧蚀情况下的动特性模型尚无公开研究。

随着全大气层高超滑翔技术的广泛应用， 与再入飞行模式相比， 气动加热同样恶劣， 而且一方面气动加热时间较长， 另一方面为了追求方案最优， 材料的承载能力也进一步挖潜， 使用温度的上限也逐步提高。 由此出现的热模态问题显著， 对飞行器结构动特性预示带来了难题， 尤其是在对飞行过程中的结构动特性变化趋势的精确预示， 涉及到壳体热模态和考虑部分烧蚀的外部防热层的热模态预示， 其中关键是缺少飞行过程中的热模态预示的模型， 这一模型的研究显得意义重大。 因此， 本文针对全大气层滑翔的高超声速飞行器动特性预示需求， 提出了考虑防热层以及烧蚀的高超声速飞行器截面刚度计算模型， 并进行了详细的推导和分析。

1 模态计算理论

王 亮， 等： 考虑防热层以及烧蚀的高超声速飞行器动特性计算模型研究

飞行器结构动特性求解方程可以表示为

KΦ=ΛMΦ（1）

式中： K， M为飞行器结构的总刚度、 总质量矩阵;  Λ， Φ为飞行器结构的特征值与特征向量矩阵。

因此， 在计算飞行器结构动特性时， 首先要确定飞行器结构动力学模型的总刚度与总质量矩阵。 而结构动力学模型的总刚度矩阵由选取的梁模型特性所决定， 因此截面刚度模型对梁模型刚度特性影响是显而易见的。

2 动特性模型

高超声速飞行器结构动力学建模时， 一般将飞行器结构简化为梁-质量块模型， 并将飞行器结构分为若干个分站（这些分站一般为舱段连接点、 仪器放置点， 或者某些大型部组件的重心）， 给出高超声速飞行器各分站处的集中质量以及各分站之间梁的材料和截面尺寸。

2.1 质量模型

每个节点为一分站， 分站质量与转动惯量添加在对应的自由度上， 对应分站的转动惯量计算公式为

J=me2（R2+46h2），对于端点me2（R2+16h2），对于中间点 （2）

式中： me为该站的集中质量;  R为该站的弹体半径;  h处于端点时， 取与相邻点间距之半， 处于中间点时， 取相邻前后两点间距之半。

2.2 刚度模型

高超声速飞行器总刚度矩阵由各梁单元刚度矩阵组装得到。 一方面由于高超声速飞行器为连续气动外形， 各节点处截面半径不相等， 需采用变截面梁单元技术;  另一方面， 由于高超声速飞行器长时间气动加热效应的影响， 内部金属壳体温度随时间逐渐升高， 外部防热层沿周向和厚度方向温度存在差异， 而且迎风面会由于气动加热产生烧蚀和碳化， 厚度较原始状态会发生改变。 因此， 计算梁单元刚度矩阵时， 使用型函数导数与材料参数矩阵乘积在单元长度上积分的方法， 得到各单元的刚度矩阵。 梁单元的刚度计算公式为

[Ke]=∫l0（EIBTB+kGAε2NTN）dx（3）

式中： G为截面剪切模量;  A为截面面积;  k为剪应力修正因子， 这里取4/3[3];  ε=6EIl/12EIl+GAlk;  [B]=[B]-6εl2xl-1[N]， B为梁单元的型函数对坐标的二次导数;  [N]=[2/l1-2/l1]T， l为梁单元长度;  EI（x）为梁单元的抗弯刚度， 有

N=1-3x2/l2+2x3/l3x-2x2/l+x3/l23x2/l2-2x3/l3-x2/l+x3/l2;

B=dN2dx2=-6/l2+12x/l3-4/l+6x/l26/l2-12x/l3-2/l+6x/l2。

针对不同截面的抗弯刚度， 本文提出了考虑带烧蚀的防热层的飞行器截面抗弯刚度计算模型。

首先以旋成体高超声速飞行器为例， 给出了飞行器典型截面的示意图， 其中包含金属壳体、 防热层原始层和防热层烧蚀层， 如图1所示。

从图中可以看出， 截面共包括三层：

金属壳体位于最内层， 厚度沿周向相同， 其抗弯刚度计算公式中弹性模量为E1， 随温度t变化， 而温度随周向子午角θ变化， 截面惯性矩I1（x）是外径R1和内外径比例α的函数， I1（x） = πR414（1-α4）， 面积A1为外径R1和内外径比例α的函数， A1（x） = πR21（1-α2）。

防热层原始层是中间层， 位于金属壳体外部， 是防热层除烧蚀层以外的部分， 厚度沿周向子午角θ不同， 其中由于迎风面烧蚀严重， 厚度较背风面小， 因此其外径R2（θ）抗弯刚度计算公式中弹性模量为E2， 随温度t变化， 与E1类似。 截面惯性矩I2（x）和面积A2均为外径R2（θ）的函数， 即θ的函数。

防热层烧蚀层是最外层， 是防热层烧蚀碳化后的部分， 厚度沿周向子午角θ不同， 其中由于迎风面烧蚀严重， 厚度较背风面大， 因此其外径R3（θ）抗弯刚度计算公式中弹性模量为E3， 随温度t变化， 与E1和E2类似。 截面惯性矩I3 （x）和面积A3均为外径R3（θ）的函数， 即θ的函数。

因此， 截面的抗弯刚度计算公式为

EI（x）= E1I1（x）+ E2I2（x）+ E3I3（x）（4）

式中：

E1I1（x）=∫AE1[t（θ）]dθ×πR414（1-α4）;

E2I2（x）=∫AE2[t（θ）]R22（θ）sin2（θ）R2（θ）dR2（θ）dθ;

E3I3（x）=∫AE3[t（θ）]R23（θ）sin2（θ）R3（θ）dR3（θ）dθ。

截面的抗剪刚度计算公式为

GA（x）= G1A1（x）+ G2A2（x）+ G3A3（x）  （5）

式中：

G1A1（x）=∫AE1[t（θ）]2（1+μ）dθ×πR21（1-α2）;

G2A2（x）=∫AE2[t（θ）]2（1+μ）R2（θ）dR2（θ）dθ;

G3A3（x）=∫AE3[t（θ）]2（1+μ）R3（θ）dR3（θ）dθ。

從而求得各截面的抗弯刚度和抗剪刚度， 代入式（3）中， 通过积分可求得梁单元的刚度矩阵， 再通过对相应节点的刚度矩阵进行组装可以获得飞行器结构的总刚度矩阵。 实际使用时， 需要已知截面各层的温度随周向子午角的分布， 以及各层材料随温度变化的力学性能。

3 应用程序

根据上述单元刚度和质量矩阵的计算模型和计算公式， 具体使用的计算程序和流程如图2所示。

4 算  例

由于烧蚀对飞行器整体质量的变化影响较小， 因此模态的变化主要取决非于飞行器整体的刚度分布。 本节采用给出的新模型计算截面抗弯刚度， 并与传统方法进行对比。

假设截面的参数如表1所示。 其中， 由于金属壳体的传热系数较大， 算例中假设整个壳体温度空间分布是一致的， 另外传统模型仅为原始层模型， 因此传统模型中防热层刚度为考虑温度和烧蚀对刚度的影响。

新模型采用本文的计算方法， 与传统方法相比， 金属壳体刚度计算方法一致。 考虑温度效应的影响， 防热层将原始层和烧蚀层合并， 计算参数采用原始层材料参数。 两种方法计算得到的截面抗弯刚度和抗剪刚度的结果如表2所示。

从计算结果可以看出， 新模型计算的截面抗弯刚度和抗剪刚度较传统模型偏小， 相对差值分别为3.38%和2.69%。 一方面， 这是由于新模型对防热层原始层和烧蚀层尺寸建模细化， 这两层的物理参数存在一定的差异; 另一方面， 防热层细化模型考虑了温度的影响， 这说明了有必要对防热层的刚度进行细化建模。

5 结  论

针对全大气层滑翔的高超声速飞行器动特性预示需求， 本文提出了考虑防热层以及烧蚀的高超声速飞行器截面刚度计算模型， 通过分析得出以下结论：

（1） 提出了3层组合的高超声速飞行器截面刚度计算模型， 分别为最内层承力金属壳体、 中间防热层原始层和最外层的防热层烧蚀层。

（2） 根据刚度线性叠加原理， 给出了各层以及总的截面刚度参数的计算方法和具体计算公式。

（3） 针对实际使用， 给出了高超声速飞行器结构动特性具体的计算分析程序和步骤， 从材料参数获取直到飞行器结构动特性计算。

（4） 由于新模型对防热层进行细化建模， 而且考虑了尺寸和温度对防热层刚度的影响， 新模型计算的截面抗弯刚度和抗剪刚度较传统模型偏小。

通过本研究， 细化了全大气层滑翔的高超声速飞行器动特性预示模型。 不仅可以获得高超声速飞行器在各飞行时间点下的动特性数值， 而且可以建立随飞行弹道和轨迹的动特性模型， 为姿态控制提供了精细的飞行器弹性模型。

参考文献：

[1] Cleghorn W L， Tabarrok B. Finite Element Formulation of a Tapered Timoshenko Beam for Free Lateral Vibration Analysis[J]. Journal of Sound and Vibration， 1992， 152（3）： 461-470.

[2] Al-Gahtani H J. Exact Stiffnesses for Tapered Members[J]. Journal of Structural Engineering， 1996， 122（10）： 1234-1239.

[3] Berman A. Mass Matrix Correction Using an Incomplete Set of Mea-sured Modes[J]. AIAA Journal， 1979， 17（10）： 1147-1148.

[4] 陆念力， 张宏生. 计及二阶效应的一种变截面梁精确单元刚度阵[J]. 工程力学， 2008， 25（12）： 60-64.

Lu Nianli， Zhang Hongsheng. Exact Elemental Stiffness Matrix of a Tapered Beam Considering Second-Order Effects[J]. Engineering Mechanics， 2008， 25（12）： 60-64.（in Chinese）

[5] 周凌远， 李彤梅， 李乔. 基于网格截面法的非线性梁单元[J]. 西南交通大学学报， 2009， 44（5）： 726-731.

Zhou Lingyuan， Li Tongmei， Li Qiao. Nonlinear Beam Element with Meshed Sections[J]. Journal of Southwest Jiaotong University， 2009， 44（5）： 726-731.（in Chinese）

[6] 周世軍， 朱晞. 一组新的Timoshenko梁单元一致矩阵公式[J]. 兰州铁道学院学报， 1994， 13（2）： 1-7.

Zhou Shijun， Zhu Xi. A Set of New Consistent Matrix for Formulations of Timoshenko Beam Element[J]. Journal of Lanzhou Railway Institute， 1994， 13（2）： 1-7.（in Chinese）

[7] 李华， 武兰河， 黄羚. 一种新的考虑剪切变形影响的梁单元[J]. 石家庄铁道学院学报， 1999， 12（4）： 58-61.

Li Hua， Wu Lanhe， Huang Ling. A New Element in Consideration  of the Effect of the Shear Deformation on Beam Deflection[J]. Journal of Shijiazhuang Railway Institute， 1999， 12（4）： 58-61.（in Chinese）

[8] 刘文炼. 一种新的Timoshenko梁单元[J]. 广东土木与建筑， 2007， 14（12）： 13-14.

Liu Wenlian. A New Timoshenko Beam Element [J]. Guangdong Architecture Civil Engineering， 2007， 14（12）： 13-14.（in Chinese）

[9] 王乐， 王亮. 一种新的计算Timoshenko梁截面剪切系数的方法[J]. 应用数学和力学， 2013， 34（7）： 756-763.

Wang Le， Wang Liang. A New Method of Obtaining Timoshenko’s Shear Coefficients[J]. Applied Mathematics and Mechanics， 2013， 34（7）： 756-763.（in Chinese）

[10] 何成. 高温环境下结构动力学建模关键技术研究[D]. 南京： 南京航空航天大学， 2014.

He Cheng. Key Technology of Dynamic Modeling in High-Tempera-ture Environment[D]. Nanjing： Nanjing University of Aeronautics and Astronautics， 2014. （in Chinese）

[11] 王亮， 商霖， 牛智玲， 等. 集中质量法与一致质量法对导弹模态计算的影响研究[J]. 强度与环境， 2013， 40（2）： 28-32.

Wang Liang， Shang Lin， Niu Zhiling， et al. Study on the Influence of the Consistent Mass Method and Lumped Mass Method in the Mode of the Rocket[J]. Structure & Environment Engineering， 2013， 40（2）： 28-32.（in Chinese）

[12] 王亮， 商霖， 王樂. 战术导弹结构动力学建模研究[J]. 战术导弹技术， 2013（1）： 24-27.

Wang Liang， Shang Lin， Wang Le. Study on the Structure Dynamic Modeling of Tactical Missile[J]. Tactical Missile Technology， 2013（1）： 24-27.（in Chinese）

[13] 王亮， 周晓丽， 廖选平， 等. 飞行器扭转模态建模分析研究[J]. 战术导弹技术， 2015（6）： 53-56.

Wang Liang， Zhou Xiaoli， Liao Xuanping， et al. Study of the Structure Dynamic Modeling of Aircraft’s Torsion Mode[J]. Tactical Missile Technology， 2015（6）： 53-56.（in Chinese）

[14] 张琪， 刘莉. 导弹等效密度有限元模型的模态分析[J]. 兵工学报， 2008， 29（11）： 1395-1399.

Zhang Qi， Liu Li. Modal Analysis of Missile’s Equivalent Density Finite Element Model[J]. Acta Armamentarii， 2008， 29（11）： 1395-1399.（in Chinese）

Study on  Dynamic Model of  Hypersonic

Aircraft Considering TPS with Ablation Effects

Wang Liang*， Zhang Yan， Cai Yipeng， Cai Wenjie

（China Academy of Launch Vehicle Technology， Beijing 100076， China）

Abstract： To indicate the dynamic characteristics of  hypersonic aircraft gliding in the whole atmosphere， a calculation model of section stiffness considering thermal protection system（TPS） and ablation effects is proposed. Firstly， the section stiffness model with three-layer combination is put forward， which include the inside main backbone， the in-between original layer of TPS， and the outmost ablative layer of TPS. Secondly，

based on the linear superposition principle，

the calculation methods of each layer and total section stiffness parameters are given.  At last， the analysis strategy of the hypersonic aircraft’s dynamic characteristics is given for actual use.

Key words： hypersonic; dynamic characteristics; thermal protection system（TPS）; ablation effect; flexural stiffness