邵嘉健，薛鵬飛

基于Gauss偽譜法的高空飛行器再入段軌跡優(yōu)化對(duì)傳熱效應(yīng)的影響分析

邵嘉健，薛鵬飛

（空間物理重點(diǎn)實(shí)驗(yàn)室，北京，100076）

基于Gauss偽譜法和二階有限差分（Gauss Pseudospectral Method，GPM），研究了再入過程中高速飛行器的傳熱問題，并依據(jù)最內(nèi)溫層升溫最小的目的進(jìn)行飛行軌跡的數(shù)值優(yōu)化。主要思路：構(gòu)造傳熱分析模型，并依據(jù)傳熱方程構(gòu)造每一溫層的傳熱微分方程；將各溫層微分方程以及動(dòng)力學(xué)微分方程作為偽譜法中的微分方程約束條件代入，進(jìn)行軌跡數(shù)值優(yōu)化設(shè)計(jì)；利用少量LG點(diǎn)構(gòu)造拉格朗日多項(xiàng)式，再通過一維插值獲得大量LG點(diǎn)的值，獲得更高精度的擬合結(jié)果曲線。以某高超聲速飛行器為對(duì)象用本方法進(jìn)行數(shù)值計(jì)算，結(jié)果驗(yàn)證了方法具有一定的可行性。

Gauss偽譜法；傳熱分析模型；軌跡優(yōu)化設(shè)計(jì)

0 引 言

在飛行器氣動(dòng)加熱的過程中，會(huì)影響乘員艙殼體溫度的因素有兩個(gè)：防隔熱層材料及厚度以及飛行器由于氣動(dòng)加熱形成的表面熱流。而表面熱流的形成與飛行器飛行模式密切相關(guān)，本文的研究?jī)?nèi)容便圍繞彈道優(yōu)化設(shè)計(jì)對(duì)飛行器殼體加熱的影響展開[1,2]。

為實(shí)現(xiàn)本文對(duì)高速飛行器再入軌跡快速優(yōu)化分析的目的，擬利用Gauss偽譜方法展開研究。Gauss偽譜方法是一種利用全局插值多項(xiàng)式構(gòu)建方程的直接配點(diǎn)法，它相對(duì)于一般直接配點(diǎn)法的優(yōu)勢(shì)在于用較少的節(jié)點(diǎn)代入就可以獲得較高的精度[3]。

以往利用Gauss偽譜法開展的高速飛行器軌跡數(shù)值優(yōu)化的研究中，國防科技大學(xué)的雍恩米、唐國金等人在研究高超聲速飛行器滑翔式再入的快速軌跡優(yōu)化問題中，使用駐點(diǎn)熱流密度積分，得到再入過程的熱載作為優(yōu)化目標(biāo)[4]。但這與本文研究?jī)?nèi)壁面升溫最小的問題仍然不同，為了說明問題，做計(jì)算實(shí)例如下。

保持初始條件一致，分別通過積分計(jì)算內(nèi)壁面溫度的上升，保證二者總熱流關(guān)相同，得到結(jié)果如圖1所示。由圖1可以看出，在相同的總熱載下，不同的加熱時(shí)間，內(nèi)壁面溫度的上升差異很大，因此，總熱載并不能有效地反映飛行器結(jié)構(gòu)的溫度變化。

圖1 相同總熱載下溫度上升對(duì)比

結(jié)合傳熱學(xué)和軌跡數(shù)值優(yōu)化，通過二階有限差分的方法，將飛行器殼體上的溫度按照物面法向進(jìn)行離散之后，構(gòu)建有關(guān)于各層溫度關(guān)于時(shí)間的狀態(tài)方程。并將此方程代入動(dòng)力學(xué)微分方程，即將溫度作為狀態(tài)變量考慮，以此開展Guass偽譜法進(jìn)行飛行軌跡數(shù)值優(yōu)化，并且以末點(diǎn)乘員艙壁溫度最低為優(yōu)化設(shè)計(jì)目標(biāo)。

1 優(yōu)化問題的數(shù)學(xué)模型

1.1 飛行器動(dòng)力學(xué)模型

忽略地球自轉(zhuǎn)的影響，臨近空間飛行器無動(dòng)力再入段的動(dòng)力學(xué)方程如下

1.2 傳熱層的模型

1.2.1 傳熱方程

為了研究本文的問題，需要建立合適的傳熱分析模型。已有的分析結(jié)果表明，如果不需要準(zhǔn)確分析連接件附近區(qū)域溫度分布，可以不必要建立三維模型。因此采用一維簡(jiǎn)化熱分析模型完全能夠滿足本文優(yōu)化的精度要求。

由于使用有限插分方法的限制，要求各層的厚度相同。建立一維熱分析模型如圖2所示，由一系列連續(xù)的厚度相同、性能參數(shù)不同、對(duì)傳熱起不同作用的材料相互串連而成，最后一層是內(nèi)部冷結(jié)構(gòu)，各層內(nèi)部及各層不同材料之間存在沿厚度方向的連續(xù)導(dǎo)熱。在直角坐標(biāo)系中的無熱源一維瞬態(tài)導(dǎo)熱控制方程為[6]

式中 ，，及分別為各層所使用材料的密度、熱傳導(dǎo)系數(shù)、比熱容、及溫度。

1.2.2 傳熱微分方程

a）中間層的微分方程。

b）外邊界微分方程。

1）外邊界邊界條件。

對(duì)于外邊界，本文采用第3類邊界條件，熱平衡方程為[7]

外壁面的能量平衡方程為

氣動(dòng)加熱產(chǎn)生的熱流如下：

至此，外壁面離散控制方程可以寫為

2）邊界熱流分析。

外邊界熱流為由氣動(dòng)熱產(chǎn)生的冷壁熱流，取飛行器迎風(fēng)面某特征位置為參考點(diǎn)，其計(jì)算公式為

c）內(nèi)邊界微分方程。

對(duì)于內(nèi)邊界，一階離散導(dǎo)熱控制方程為

綜上，微分方程組如下式所示：

2 利用偽譜法轉(zhuǎn)化問題

2.1 Gauss偽譜法簡(jiǎn)介

Gauss偽譜方法通過將狀態(tài)變量和控制變量在一系列Legendre-Gauss（LG）點(diǎn)上進(jìn)行離散，并將這些離散點(diǎn)作為節(jié)點(diǎn)，構(gòu)造Lagrange插值多項(xiàng)式，來逼近原狀態(tài)變量和控制變量[8]。再通過對(duì)插值多項(xiàng)式求導(dǎo)獲得微分矩陣，以矩陣來逼近狀態(tài)變量對(duì)時(shí)間的導(dǎo)數(shù)，這樣就將微分方程約束轉(zhuǎn)換為代數(shù)方程約束。性能指標(biāo)中的積分項(xiàng)和終端狀態(tài)約束都由Gauss積分計(jì)算即可獲得性能指標(biāo)中的積分項(xiàng)和終端狀態(tài)約束。經(jīng)上述變換，可將最優(yōu)控制問題轉(zhuǎn)化為通過代數(shù)約束的一系列參數(shù)優(yōu)化求解問題，稱為非線性規(guī)劃問題（NLP），之后利用SQP算法求解該問題[9]。

2.2 連續(xù)優(yōu)化問題的轉(zhuǎn)化

a）時(shí)域變化。

使用legendre插值多項(xiàng)式需要保證自變量范圍為[-1,1]，因此做時(shí)域變換：

通過上式，時(shí)間區(qū)間從[0,f]轉(zhuǎn)換到[-1,1]。

b）將時(shí)間離散化。

將變換后的時(shí)域離散，獲得個(gè)LG點(diǎn)。LG點(diǎn)的定義為勒讓德多項(xiàng)式的根。

c）狀態(tài)與控制變量的全局插值多項(xiàng)式近似。

式中 Lagrange插值基函數(shù)由下式計(jì)算：

同理，獲得控制變量的近似表達(dá)式：

d）微分方程約束轉(zhuǎn)化。

由式（14），即有：

從而將問題由一開始的微分方程約束轉(zhuǎn)變?yōu)槿缦率剿镜拇鷶?shù)方程約束:

式中=1,…,；= 0,…,。

e）終端狀態(tài)約束。

上式的狀態(tài)變量未包含終端時(shí)刻節(jié)點(diǎn)處的值，終端狀態(tài)由微分方程對(duì)時(shí)間積分再加上初始狀態(tài)值得到，其中，用離散的Guass積分近似原有連續(xù)函數(shù)積分過程，可得：

f）性能函數(shù)指標(biāo)的近似。

原始含有積分方程的性能函數(shù)泛函為

將其中的積分利用Gauss積分方程代換，得到性能指標(biāo)的新的函數(shù)表達(dá)式

3 求解優(yōu)化問題策略

針對(duì)上述問題及其特殊性，本文基于Gauss偽譜法提出以下求解大量約束條件下最優(yōu)飛行軌跡的策略：

a）先利用龍格-庫塔法求取式（16），建立合適的控制策略，并利用粒子群算法以防隔熱層厚度為優(yōu)化變量，進(jìn)行優(yōu)化以期獲得全程再入飛行過程中各控制層合適的溫度變化，要求溫度變化在材料熱受能力之內(nèi)，且具有一定程度的溫度變化，便于優(yōu)化問題的求解。將獲得的結(jié)果以LG點(diǎn)轉(zhuǎn)化后的時(shí)刻離散，作為Gauss偽譜法優(yōu)化的初值。

b）利用Gauss偽譜法計(jì)算LG點(diǎn)個(gè)數(shù)=5時(shí)刻的結(jié)果，獲得其節(jié)點(diǎn)處離散點(diǎn)的數(shù)值之后，利用拉格朗日插值法進(jìn)行全過程的設(shè)計(jì)變量隨時(shí)間的變化曲線的擬合，求解=10的勒讓德多項(xiàng)式的解，獲得=10時(shí)的LG點(diǎn)，通過式（21）將其還原成飛行過程的時(shí)刻，利用對(duì)=5時(shí)刻的節(jié)點(diǎn)數(shù)值進(jìn)行一維插值獲得這些新的節(jié)點(diǎn)處設(shè)機(jī)變量的值，再利用這些值進(jìn)行拉格朗日插值法獲得最終結(jié)果曲線。

c）將獲得的結(jié)果與代入控制策略的龍格-庫塔法的結(jié)果進(jìn)行對(duì)比，觀察二者的差異。

4 數(shù)值計(jì)算實(shí)例

以遠(yuǎn)程高速滑翔式再入飛行器為仿真對(duì)象，氣動(dòng)參數(shù)采用波音公司設(shè)計(jì)研究的錐形體再入機(jī)動(dòng)飛行器CAV-H的氣動(dòng)數(shù)據(jù)，并通過擬合獲得。飛行器最大升阻比約為2.4，選取氣動(dòng)參考面積為0.35 m2，質(zhì)量為907 kg，最大飛行攻角和最大升阻比攻角為30°，最大升阻比攻角為10°。再入飛行器初始數(shù)值計(jì)算初值條件、終端條件及過程相關(guān)約束見表1。

表1 數(shù)值計(jì)算初值條件、終端條件及過程約束

Tab.1 Numerical Calculation of Initial Conditions,Terminal Conditions and Process Constraints

初始條件x/kmy/kmv/（m?s-1）θ/（°） 0805000-2 T0～T2/KT3/K—— 280280—— 終端條件x/kmy/kmv/（m?s-1）T3max/K 240020≧500350 過程約束Q/Pa/(°)/(°)— ≧500030—

圖3 再入軌跡優(yōu)化結(jié)果

續(xù)圖3

5 結(jié)束語

基于Gauss偽譜法求解高速飛行器再入時(shí)的傳熱問題。創(chuàng)新性地將溫度利用有限差分方法作為一個(gè)新的狀態(tài)變量代入偽譜法微分方程組，并整理出各離散點(diǎn)溫度與控制變量的關(guān)系，與速度、高度等狀態(tài)變量共同優(yōu)化，獲得結(jié)果。優(yōu)化計(jì)算結(jié)果表明，彈道軌跡變化切實(shí)引起了飛行器殼體溫度變化；偽譜法獲得結(jié)果與龍格庫塔法求解的結(jié)果基本一致，表明本文使用的方法可靠性較高。本文研究的成果可進(jìn)一步應(yīng)用到飛行器的防熱問題中。

[1] 閆曉東, 唐碩. 基于偽譜法的亞軌道飛行器返回軌跡優(yōu)化設(shè)計(jì)[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2010, 28(5): 748-751.

Yan Xiaodong, Tang Shuo. Optimal design of return trajectory of suborbital vehicle based on pseudospectral method[J]. Journal of Northwestern Polytechniacal University, 2010, 28(5): 748-751.

[2] 宗群, 田栢苓, 竇立謙. 基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化[J]. 宇航學(xué)報(bào), 2010, 31(7): 1775-1781.

Zong Qun, Tian Bailing, Dou Liqian. Trajectory optimization of near space vehicle ascent stage based on gauss pseudospectral method[J]. Journal of Astronautics, 2010, 31(7): 1775-1781.

[3] 周文雅, 楊滌, 梁新剛. 利用高斯偽譜法求解最優(yōu)氣動(dòng)輔助軌道轉(zhuǎn)移[J].吉林大學(xué)學(xué)報(bào)(工學(xué)版), 2010, 40(5): 454-1459.

Zhou Wenya, Yang Di, Liang Xingang. Use Gauss pseudospectral method to solve the optimal aerodynamic assisted orbit transfer[J]. Journal of Jilin University(Engineering and Technology Edition), 2010, 40(5): 1454-1459.

[4] 雍恩米, 唐國金, 陳磊. 基于Gauss偽譜法的高超聲速飛行器再入軌跡快速優(yōu)化[J]. 宇航學(xué)報(bào), 2008, 29(6): 1766-1772.

Yong Enmi, Tang Guojin, Chen Lei. Rapid optimization of reentry trajectory of hypersonic vehicle based on Gauss pseudospectral method[J] Journal of Astronautics, 2008, 29(6): 1766-1772.

[5] 雍恩米, 陳磊, 唐國金. 飛行器軌跡優(yōu)化數(shù)值方法綜述[J]. 宇航學(xué)報(bào), 2008, 29(2): 397- 406.

Yong Enmi, Chen Lei, Tang Guojin. A review of numerical methods for aircraft trajectory optimization[J]. Journal of Astronautics, 2008, 29(2): 397- 406.

[6] John D, et al. Hypersnoic and high temperature gas dynamics[M]. NewYork: McGraw-Hill Book Company. American Institute of Aeronautics and Astronautics, 2002.

[7] 薛鵬飛, 龔春林, 谷良賢. 基于軸對(duì)稱比擬的氣動(dòng)熱計(jì)算方法研究[J].計(jì)算機(jī)仿真, 2012, 29(6): 97-100.

Xue Pengfei, Gong Chunlin, Gu Liangxian. Aerodynamic thermal calculation method based on Axisymmetric comparison[J]. Computer Simulation, 2012, 29(6): 97-100.

[8] Huntington G T. Advancement and analysis of a gauss pseudospectraltranscription for optimal control problems[D]. MA: Massachusetts Institute of Technology, 2007.

[9] David Benson. Direct trajectory optimization and costate estimation via an orthogonal collocation method[J]. Journal Of Guidance, Control, and Dynamics, 2006, 29(6): 1435-1440.

Influence Analysis of Heat Transfer Effect on Reentry Trajectory Optimization of Adjacent Space Vehicles based on Gauss Pseudospectral Method

Shao Jia-jian, Xue Peng-fei

(Science and Technology on Space Physics Laboratory, Beijing, 100076)

Based on Gauss Pseudospectral and the second order finite difference method, to study the heat transfer problem in the process of hypersonic flight vehicle reentry, and according to the purpose of the lowest temperature rise in the innermost layer to perform numerical optimization of flight path. The main train of thought. constructing the heat transfer analysis model andthe heat transfer differential equation of each temperature layer according to the heat transfer equation. Substituting differential equation of every temperature layer and dynamic differential equation into the constraint conditions of the differential equation in the pseudo-spectral method, for the trajectory numerical optimization design; constructing lagrangian polynomials with a small number of LG points, and a large number of LG points are obtained through one-dimensional interpolation, to obtain the fitting result curve with higher accuracy. This method is used for numerical calculation of a hypersonic vehicle, and the results verify the feasibility of it.

gauss pseudo-spectral method; heat transfer analysis model; trajectory optimization design

2097-1974(2023)01-0011-05

10.7654/j.issn.2097-1974.20230103

V412.4+4

A

2018-12-14；

2022-12-01

邵嘉健（1994-），男，工程師，主要研究方向?yàn)轱w行力學(xué)。

薛鵬飛（1987-），男，博士，高級(jí)工程師，主要研究方向?yàn)轱w行力學(xué)。