基于传热增广模型的轨迹优化与防热结构分析

2022-03-11张腾飞龚春林薛鹏飞

系统工程与电子技术 2022年3期

张腾飞, 龚春林, 粟华,*, 薛鹏飞

(1. 西北工业大学陕西省空天飞行器设计重点实验室, 陕西西安 710072;2. 空间物理重点实验室, 北京 100076)

0 引言

高超声速飞行器在未来国际局势中有非常重要的战略作用[1],这类飞行器具有突防能力强、运输效率和打击精度高、杀伤威力大等优势,但高速飞行造成的热防护问题是其设计过程中必须着重考虑的因素之一[2-4]。

轨迹优化在飞行器初步方案设计阶段,对方案可行性论证、设计参数分析、指导制导控制系统设计等具有重要意义。飞行器轨迹优化问题是一类包含非线性约束的最优控制问题,数值方法已成为解决复杂轨迹优化问题的主要方法[5]。根据是否直接对性能指标进行寻优将其分为间接法和直接法[6-8],随着伪谱法的兴起[9-12]以及数值优化技术的不断改进,使得在轨迹优化问题中使用更加精细的模型成为可能[13-17]。

在轨迹优化中考虑热约束是为了保证飞行过程中结构不遭受破坏,而结构破坏与否一般是依据是否超出对应材料的耐温极限进行判断的。飞行过程中的热流密度仅是计算结构温度的边界条件,而现有的轨迹优化大多只通过约束热流密度来体现热防护要求[18-21]。需要明确的是,过高的热流密度不是结构破坏的原因,是否超出结构材料的耐温极限才是关键,热流密度本质上属于不必要的“软约束”。随着高超声速飞行器设计中对热防护的重视程度越来越高,在初步轨迹优化设计中就要考虑结构温度约束。如此,约束热流密度的方式便存在明显的不足。

首先,最大热流密度无法与结构最高温度建立等价的简单映射关系,事实上大多数轨迹优化是通过经验或表面热辐射平衡条件给出最大热流密度约束,这显然是不够准确的。在一些着重考虑热约束的研究中[22-23],将热流密度积分值最小作为优化目标,但这样得到的轨迹是十分保守的。因此,为了解决考虑结构温度约束的轨迹设计问题,采用约束热流密度峰值的轨迹优化模型需基于热流密度上限越大,结构温度峰值越高的设想,并且经过多次轨迹优化与传热分析的迭代,才能给出保证结构安全又不十分保守的热流密度上限与飞行轨迹。其次,将传热分析与轨迹优化分步进行,则不能考虑到飞行轨迹与结构传热之间的相互影响,得到的轨迹优化结果必然不是最优的。综上,基于一般模型求解会导致轨迹与热防护结构设计之间反复迭代、互相妥协的问题。此外,为了对后续详细设计提供参考,方案设计阶段需要对不同防热结构尺寸方案进行可行性分析,而由于热流密度约束调节范围有限,导致结构尺寸的可行域较小,不利于详细设计阶段总体层面的协调。

总之,一般约束热流密度的轨迹优化模型不便于解决考虑结构温度约束的轨迹优化设计问题,需建立耦合传热过程的轨迹优化模型,以便对结构温度直接进行约束。文献[24-25]尽管考虑了运动和热传导的耦合关系,但研究主要针对最优轨迹周期性巡航的规律研究,并没有对防热结构作详细分析,且所建立的传热模型未考虑材料属性随温度的变化导致误差较大。

本文通过空间离散的方法,将一维传热方程转化为一组一阶微分方程,考虑到了材料属性随温度的变化,使得计算精度更高,有利于得出更为准确的分析结论。使用该模型对运动方程进行增广,增广之后的系统状态方程引入了结构温度这一状态量,在轨迹优化中可加入更为直接的结构温度约束,通过一次轨迹优化便可得到满足结构温度约束的最优轨迹,而不用反复验证。

以高超声速飞行器CAV-H的再入轨迹优化问题为例,首先,对于基准热防护结构方案,计算结果表明:基于传热增广轨迹优化模型的求解过程更为便捷和高效,且得到的最大航程较一般模型的结果有所提升,最优轨迹在防热结构性能的利用以及姿态变化等方面也更为合理。其次,基于不同轨迹优化模型确定了防热结构尺寸的可行域,结果表明,增广模型对应的可行域远大于一般模型,这为设计提供了更大的选择空间。最后,在结构尺寸可行域内计算最优轨迹,分析总结各层结构尺寸对最大航程的影响规律以及以不同形式再入对应的结构传热机制的影响,对防热结构设计提供一定规律性的参考意见。

1 运动模型

为简化问题和突出重点,仅考虑地球为非旋转均质圆球的二维情况。根据文献[26],给出简化运动方程组:

(1)

式中:r和θ分别为飞行器质心到地心的距离(简称地心距)和航程角;v和γ为速度和当地速度倾角;L和D为升力和阻力;m为飞行器质量;μ为引力常数。

气动力计算如下:

(2)

式中:qd为动压;ρa为大气密度(使用USSA1976大气模型);Sref为参考面积;CL和CD为升阻力系数,是攻角α和速度v的函数。

2 热传导模型

2.1 一维传热模型

2.1.1 一维传热问题

本文不考虑主动热防护。典型的被动热防护结构方案一般为高温防热层+隔热层+承力层的多层方案[27-28],其中防热层一般采用陶瓷基复合材料或C/C复合材料;隔热层一般采用气凝胶、柔性隔热毡等材料;内层承力结构为金属材料,在计算传热时可以忽略,即将三层结构简化为两层,并将内表面视为绝热壁面。

若非详细研究结构的热传导问题,一维模型已经足够精确[29-30],即只考虑温度沿厚度方向的分布。一维传热问题如图1所示,图中qs和qr分别表示表面热流密度和热辐射密度,Tω为外表面温度,δf和δs表示两层结构的厚度。

一维热传导的控制方程和边界条件分别为

(3)

(4)

式中:T为温度;ρ、c、k分别为材料的密度、比热和热导率,后二者为温度T的函数;x为结构厚度方向的位置。qs和qr的计算如下:

qs=Csρ0.5v3

(5)

(6)

式中:Cs为与飞行器外形有关的常数;σ和ε分别为表面辐射率和斯特藩-玻尔兹曼常数;T∞为环境温度。

2.1.2 单元划分

式(3)是一抛物线型二阶非线性偏微分方程,需要对其进行离散化处理。一般热传导问题的求解一般是将热传导方程在时间和空间上都进行离散,而这里则需要保留温度对时间的偏导数项,从而与运动方程保持一致。

首先,沿x方向布置节点:要求在内外表面以及交界面上必须布置节点,而在每层结构中,则沿厚度方向均匀布置若干节点。假设防热层和隔热层分别布置n和m个节点(交界面处的节点是共用的)。接着,作每两个相邻节点连线的中垂线,连同内外表面共同将整个结构划分为N=n+m-1个单元,每个单元对应一个节点,称为该单元的控制点。单元划分结果如图2所示,Δf和Δs分别表示两层结构内部的单元厚度。

2.1.3 有限差分模型

对于划分好单元的防热结构,通过在单元内部对控制方程进行积分,并利用差分代替导数以消除其中关于x的偏导数项。

首先,对于而非交界面处的单元i,对式(3)在单元内部沿x方向进行积分,得

(7)

式中:xi-和xi+表示单元i左右两侧边界的位置。上面处理得到的内外边界处的单元厚度是内部单元厚度的一半,在积分时需要注意积分区间长度的差异。对于式(7)左端项,认为单元内部温度相同,则式(7)可进一步简化为

(8)

其次,对于交界面处的单元n,由于材料不同,需要分段积分,即:

(9)

将式(9)进行简化如下:

(10)

认为交界面处没有热阻,则等号右端中间两项可以消去。进而将式(8)和式(10)表示如下:

(11)

式中:

(12)

(13)

(14)

2.2 传热模型验证

对于式(11),给定热防护结构参数和热流边界条件后,便可以通过数值积分的方法进行求解。下面通过算例对本文所建一维传热模型进行分析验证,具体参数如表1所示。结构整体初始温度为280 K。

表1 传热问题算例

表1中,两种材料的热导率和比热与材料温度的关系如下所示:

(15)

(16)

热流密度和环境温度变化规律(总时间为2 500 s)如下:

(17)

(18)

由于内外表面温度具有代表性,后面分析主要针对这两处的温度。计算发现,随着隔热层节点个数的增加,各节点温度迅速收敛,且收敛速度逐渐变慢;由于防热层厚度较小,热导率较大,导致温度梯度很小,因此较少的节点(但至少有一个内部节点,即防热层节点数目取3即可)就可得到十分准确的外表面温度变化规律。

(19)

图3给出了内外表面节点温度相对误差随隔热层节点数的变化情况,如前面所述,外表面温度误差几乎为0,而随着热量向内表面传递,温度误差也逐渐积累增大(内表面温度相对误差最大,其余各处的相对误差均在图中两条曲线之间)。可以看出,当隔热层节点个数增大到8之前,外表面温度相对误差持续降低,但在此之后基本不再提升,而且节点数为8对应的相对误差已经足够小(不超过2%)。因此,最终选取离散模型的节点个数为:n=3,m=8。

3 基于增广模型构建轨迹优化问题

3.1 系统状态方程

由式(1)和式(11)即可构成传热增广方程,简化表示如下:

(20)

式中:[XT]=[rθvγT1T2…T10],共同构成状态量;F,G分别为4维和10维向量值函数,由式(1)和式(11)的右端项组成;U为控制量,在这里指攻角α。两式通过qs,T∞关联在一起。增广后的系统状态方程组维数由4增加到14,而且考虑到各变量之间数值上差异较大,都使得数值求解变得困难,因此对方程进行无量纲化处理是必要的。

(21)

需要说明的是,由于关于温度的一组等式中,右端包含很多量纲复杂的材料属性,不对每个变量单独进行无量纲化处理,而是对其整体除以Tref/tref。

3.2 路径约束与优化目标

一般轨迹优化会对过载、动压以及热流密度添加路径约束以体现结构强度和防热的要求。除了过载和动压,本文直接对结构温度进行路径约束以代替热流密度约束,即:

(22)

(23)

(24)

式中:qdmax和nmax分别为动压和过载上限;最大节点温度Timax根据对应材料的耐温极限确定,如对于外表面处的T1max即材料PM1000的耐温极限:1 570 K,而对于内表面处的T10max根据内部设备要求,取为620 K。

优化目标为最大航程角(在末端高度一定的情况下与最大航程是等价的),即:

minJ=-θ(tf)

(25)

4 算例分析

4.1 问题描述

本文以CAV-H高超声速再入飞行器的最大航程再入问题为例,构建基于传热增广模型的轨迹优化问题进行求解与分析。首先对飞行器相关参数以及再入轨迹要求作以说明。CAV-H飞行器的相关参数如表2所示。

表2 CAV-H相关参数

气动模型根据文献[31]给出的数据拟合如下所示:

(26)

CAV-H飞行器被运载器运送到指定高度和速度状态之后与分离,开始再入[28]。再入过程一般分为返回段、能量管理段和攻击段。其中,返回段占据大部分时间,最大航程主要取决于这一阶段,本文也主要针对这一阶段的最大航程轨迹优化问题进行分析。动压约束和过载约束分别为70 000 Pa和2.5。具体轨迹要求如表3所示(表中高度H即地心距减地球半径Re=6 378.145 km)。

表3 再入轨迹要求

4.2 基于不同模型的计算结果对比

针对上述轨迹优化问题,使用伪谱法进行求解,作为对比,同样求解了相同条件下的约束热流密度的轨迹优化问题。如前面所述,采用一般模型进行轨迹优化之前无法给定合适的热流密度约束值,只能使用大致确定的约束值进行初步轨迹优化,然后对结构温度进行安全性验证,再根据具体情况对热流密度约束值进行更新,经过迭代计算,最终确定保证结构温度满足要求且热防护性能无过多剩余的热流密度约束值。

由于数值计算噪声的影响,通过梯度算法进行迭代较为困难,本文通过以下爬山法迭代确定热流密度约束:根据表面热辐射平衡条件大致确定初始热流密度约束为10×105W·m-2、经过多次尝试,确定收敛速度最高的初始迭代步长约为-2×105W·m-2。迭代过程中,当节点温度与温度上限之差的符号发生改变时改变搜索方向并将步长减半,由于无法同时使得内外表面温度均恰好达到温度上限,迭代以内表面温度作为收敛依据(低于温度极限不超过0.1 K)。经过17次迭代计算,得到收敛结果。上述计算与基于增广模型的计算对比如表4所示。尽管增广模型变量维度大大增加,但由于飞行状态与结构温度之间的耦合关系并不强,一次轨迹优化耗时并没有过久,且不需要迭代计算,整体计算效率并没有劣势。另外,采用一般模型计算时对迭代算法有一定要求,一般很难顺利进行迭代和快速收敛,整体来看,采用传热增广模型求解考虑结构温度约束的轨迹优化问题是十分高效的。

表4 不同模型计算结果对比

一般模型热流约束为4.047×105W·m-2,内表面最高温度达到温度上限,而外表面的最高温度则低于对应的温度极限,说明外层结构的防热性能存在剩余;基于增广模型的结果,内外表面温度均达到耐温极限,末端航程角为106.15°,相比于前者的103.56°提升了2.5%。

考虑到再入过程中速度随时间是递减的,而且为便于表示动压、热流边界以及最优升阻比攻角,图5和图6中分别给出最优轨迹的高度和攻角随飞行速度的变化规律。可以看出,基于一般模型的最优轨迹在高度第一次降低到极小值点时达到热流上限,此后高度基本持续下降;而由增广模型求解的最优轨迹,取消了对热流的限制,高度第一次可降低到动压约束边界处,由于“再入走廊”的放宽,最终航程有所增加。尽管基于增广模型的轨迹存在更大程度的跳跃,但对于无动力再入来讲严格的平衡滑翔条件并非必要的,跳跃式再入是合理的[2],事实上,尽管高度波动较大,但倾角始终保持在±2.5°以内。

从姿态(攻角)变化的角度来看,基于增广模型结果不存在攻角的大幅变化,是更为稳定的。另外,在无约束的情况下,CAV-H为实现最大航程再入,在整个中期滑翔阶段攻角是在最大升阻比攻角附近作幅度逐渐减小的波动变化,这与基于增广模型的求解结果是类似的,不同的是后者波动略有滞后,且振幅较小,基本保持在最大升阻比攻角曲线之上,主要是为了控制高度下降,降低热流密度峰值从而满足外表面温度约束。反观限制热流密度的轨迹,由于再入过程飞行器的势能没有充分转化为动能,向上跳跃的能力有限,后续则维持最大升阻比攻角飞行。

热流密度、热辐射密度以及结构温度随时间变化曲线如图7所示。图7(a)中,增广模型对应的热流密度最高达到7.27×105W·m-2,远超一般模型对热流密度所限制的4.051×105W·m-2。尽管峰值较高,但与高度类似,热流密度变化也有较大波动。对热流密度曲线积分得到增广模型轨迹总的气动加热密度为5.888×108J·m-2,而一般模型则为5.981×108J·m-2。由此可见,“热流密度上限越大,结构温度峰值越高”的设想并不严谨。另外,一般模型的结果中,热流密度与热辐射密度之间存在唯一的分界点,在该点之前外界向结构内部传热,之后结构才向外散发热量;而增广模型的结果中,热流密度在再入过程前期就存在小于热辐射密度的情况,这使得内表面温度的上升趋势略有减缓,从图7(b)可以看出,最终达到温度限制的时间也有所推迟,前者飞行时间为2 808.2 s,后者为2 855.8 s。

由前面的分析可知,约束热流密度的模型不仅计算过程繁琐,而且会造成飞行性能的损失;而采用增广模型较为便捷且结果更加合理。

4.3 防热结构分析

在初步方案设计阶段,对于一些关键设计参数,一方面需要大致确定其可行域,为后续的详细设计阶段提供参考;另一方面则需要分析与总结其对总体性能的影响规律。下面就利用传热增广轨迹优化模型对防热结构尺寸(这里指厚度)这一关键设计参数进行分析。

4.3.1 防热结构尺寸可行域分析

由于防热层尺寸一般较小且要考虑诸如结构变形、维持外形等问题,这里给定防热层厚度的设计范围为0.012～0.030 m。另外,考虑内部空间的制约,要求整体结构厚度不超过0.13 m,防热结构尺寸基本设计空间如图8灰色区域所示。首先,针对不同防热层厚度,基于增广模型对不同隔热层厚度方案进行轨迹优化,根据是否可行大致确定隔热层最小厚度构成图8中的尺寸下界,黄色区域表示基于增广模型的防热结构尺寸可行域。其次,对一般模型,计算确定了热流密度约束值最低为2.5×105W·m-2,基于此约束下最优轨迹的热流密度规律,确定不同防热层厚度下的最小隔热层厚度,图8中紫色区域则表示基于一般模型的可行域。由于热流密度约束可调节的范围十分有限,即便对于热流密度峰值最小的轨迹,仍需较厚的隔热层,导致基于一般模型的尺寸可行域较增广模型有很大的收缩。

以黄色区域内(紫色区域外)一点:[0.025 m,0.075 m]为例,其轨迹高度与结构温度变化规律如图9所示。相比于基准方案,为满足结构温度要求,这里高度“跳跃”的次数和幅度都有所增加,外表面降温的次数也更多。尽管航程稍短,但通过一般模型判定为不可行的方案,使用增广模型计算却是可行的。

4.3.2 防热结构尺寸对飞行性能的影响

4.3.3 不同再入形式对结构传热的影响

可以看出,能够满足一定航程的防热结构尺寸方案并不唯一,而在相同的航程下,不同的结构尺寸对应的再入轨迹形式是有所差别的。下面选取尺寸可行域内航程相同的两点:a=[0.012 m,0.099 m],b=[0.030 m,0.089 m],分析不同形式的再入轨迹对防热结构传热的影响。对应的轨迹以及结构内外表面温度变化规律如图12所示(a和b对应的曲线分别用1和2表示)。

可以看出,轨迹1只在前期进行了2次跳跃,后面一直保持滑翔飞行,而轨迹2则持续跳跃,但最终的航程角是相同的。两条不同形式的再入轨迹,均受到结构温度约束的影响:对于轨迹1来讲,由于防热层较薄,为了保证防热层温度不超出限制,需要在飞行前期进行跳跃,减缓外表面温度持续升高,而较厚的隔热层使得内表面温度升高较慢,允许外表面温度长时间保持较高温度。对于轨迹2来讲,防热层较厚,相比之下飞行初期外表面温度升高的趋势十分缓慢,但由于隔热层太薄,即使外表面温度较低,再入飞行中后期需要始终保持跳跃以周期性降低外表面温度,延缓内表面温度的升高。图13给出了两组结果对应的结构温度分布随飞行过程的变化情况对比(左侧为a点对应的结果,右侧为b点对应的结果),δx表示距离内表面的厚度。

可以看出,对于防热层更薄的a,从30%航程开始一直到70%航程的飞行历程中,防热层结构温度都保持在较高值,而对于b,由于防热层很厚,结构温度远低于极限温度,造成飞行前中期的防热结构性能浪费,但尽管外层结构温度较低,而由于长时间的飞行,热量持续地向内传递,也使得内表面温度最终达到极限。另外,在再入结束时刻,a对应的温度分布比较均匀,而b对应的结构温度分布则存在较大梯度,说明b在再入飞行结束时,结构有更多的热冗余,存在一定的安全隐患。因此,对于长时间的再入飞行,在一定的结构尺寸约束下,应尽可能增大隔热层厚度,防热层只需满足在再入初期飞行轨迹的协助下,外层结构温度不超出极限即可。

4.3.4 轨迹与防热结构相互影响规律总结

上面首先对防热结构尺寸可行域进行了计算和对比,基于增广模型的可行域较一般模型有很大扩展,有利于避免由于结构尺寸可行域过小导致设计过程不能继续进行的隐患,体现出飞行轨迹对结构传热过程的良好调节作用。

其次,利用增广模型分析了防热结构尺寸改变对飞行性能的影响规律,防热结构足够厚时,热约束变为消极约束,而在热约束为积极约束的结构尺寸区域内,随着厚度的增大航程角有所提升,其中隔热层厚度影响较大而隔热层厚度影响很小;结合结构总尺寸以及总质量密度随各层尺寸的变化趋势,发现在有限的结构尺寸和质量内,尽可能厚的隔热层是使得航程最大的防热结构尺寸方案。

最后,从更为细致传热过程的角度分析了相同航程下,较厚的隔热层方案相比于较厚的防热层方案的优势,即后者对结构防热性能利用不充分,而且最终时刻的热冗余较大。