高超聲速滑翔飛行器再入氣動系數改進擬合模型

徐 慧，蔡光斌，張勝修

(火箭軍工程大學導彈工程學院，西安 710025)

0 引 言

高超聲速滑翔飛行器憑借著全球快速到達、高空高速突防、大范圍機動等優良特性，成為世界各軍事強國爭相研究的熱門領域[1-4]。近年來，高超聲速滑翔飛行器再入軌跡優化問題成為飛行控制領域內研究熱點[5-7]。基于空氣熱力學外形的通用飛行器[8-9](Common aero vehicle，CAV)成為再入軌跡優化問題研究中常用的研究對象。基于此外形，美國波音公司和洛馬公司分別設計了CAV-L和CAV-H飛行器，最大升阻比分別約為2.5和3.5。

再入機動飛行過程中，高超聲速滑翔飛行器主要利用氣動力在大氣層進行遠距離滑翔，并且受到非線性氣動力和過載、動壓、熱流率等過程約束影響，使得再入軌跡優化成為一個極具挑戰性的問題[10-14]。由于再入過程復雜，氣動升力、阻力會影響到再入過程中飛行器速度、高度、航向角、航跡角的變化，進一步影響飛行器再入過程中的熱流率、動壓和過載變化。因此，建立精確的高超聲速滑翔飛行器再入氣動系數模型，對研究高超聲速滑翔飛行器再入軌跡優化，具有較大的理論意義與應用價值[15]。

在較早的高超聲速滑翔飛行器再入軌跡優化問題研究中，一般采用的是只考慮攻角影響得到的較為簡化的氣動系數模型。實際上，復雜的再入環境導致高超聲速滑翔飛行器氣動系數受到攻角、馬赫數、大氣密度等多種因素的影響。只考慮攻角這一單變量影響的氣動系數模型會帶來較大擬合誤差。李惠峰[16]等對CAV氣動外形進行了參數化建模，并分析了得到模型的氣動特性，但得到的模型縱向靜不穩定，俯仰通道氣動控制效率較低。Gao等[17]考慮攻角和馬赫數的共同影響，采用二元二次多項式的形式，進行高超聲速滑翔飛行器氣動系數的計算。孫勇[18]等通過對升力系數和阻力系數的分析，在擬合模型中使用馬赫數的負指數冪形式，提高了模型精度，減小了擬合誤差，但該模型考慮的攻角和馬赫數影響仍不夠充分，模型擬合精度有進一步提高的空間。

通過對上述多種擬合模型的對比分析，本文提出了一種改進的高超聲速滑翔飛行器再入氣動系數擬合模型。本文采用升阻比較大的CAV-H飛行器作為氣動系數模型的研究對象，利用多元非線性最小二乘的方法對本文提出的改進模型進行了參數辨識。引入一種擬合優度指標，將本文模型與文獻[17-18]中的模型進行對比研究。從升阻比的角度分析了本文氣動系數模型對高超聲速滑翔飛行器再入軌跡優化的影響。擬合仿真結果表明，本文提出的改進模型擬合精度較高。通過進一步的仿真實驗討論了改進氣動系數模型的升阻比特性，并對三類模型進行再入軌跡優化的對比實驗，仿真結果表明，基于本文提出的改進氣動系數模型，得到的優化軌跡更為平穩。

1 CAV-H氣動系數及對比模型分析

1.1 CAV-H氣動數據

準確有效的真實氣動數據和高精度的氣動擬合模型，對于研究高超聲速滑翔飛行器軌跡優化與制導具有十分重要的意義。由于世界范圍內高超聲速飛行器風洞實驗數據保密性高，目前僅能查閱到極少部分的公開數據。因此，本文采用文獻[17-19]中給出的CAV-H的升力系數和阻力系數數據開展此類研究。

CAV-H升力系數、阻力系數和升阻比分別在表1、表2和表3中給出[19-21]。

表1 CAV-H升力系數Table 1 CAV-H lift coefficient

表2 CAV-H阻力系數Table 2 CAV-H drag coefficient

表3 CAV-H升阻比Table 3 CAV-H Lift-to-drag ratio

圖1和2是分別是升力系數和阻力系數隨攻角變化的情況，可以看出升力系數幾乎與攻角成線性關系，而阻力系數與攻角的關系接近二次函數，但由于攻角為自變量的數據較少，并不能簡單地確認升力系數、阻力系數與攻角的關系。

圖1 升力系數隨攻角變化Fig.1 Lift coefficient varies with the angle of attack

圖2 阻力系數隨攻角變化Fig.2 Drag coefficient varies with angle of attack

圖3和4分別是升力系數和阻力系數隨馬赫數的變化情況。顯然，升力系數和阻力系數隨馬赫數的增加而逐漸減小，在馬赫數大于Ma10后，兩者的下降趨勢較為平緩。經分析可知，可以采用二次項或者負指數冪的形式描述升力系數及阻力系數與馬赫數之間的關系。

圖3 升力系數隨馬赫數變化Fig.3 Lift coefficient varies with Mach number

圖4 阻力系數隨馬赫數變化Fig.4 Drag coefficient varies with Mach number

1.2 對比典型模型

常見的再入軌跡優化使用傳統的氣動系數模型，僅考慮攻角帶來的影響，忽視了馬赫數的影響，帶來的擬合精度問題已經在文獻[18]中進行了一定討論。這里我們不再對單變量模型做深入研究，本文著重探討綜合考慮攻角和馬赫數影響的擬合模型。文獻[17]采用二元二次多項式形式的模型來擬合CAV-H氣動系數，表現形式如下：

(1)

式中:CL,CD分別為升力系數、阻力系數；α為攻角,M為馬赫數;ai,bi分別為升力系數、阻力系數擬合模型的待辨識參數，其中i=0,1,2,3,4。文獻[18]提出使用馬赫數負指數冪形式，來進行氣動系數擬合，并且通過討論分析表明，在阻力系數擬合模型中可以忽略攻角一次項。文獻[18]中的氣動系數擬合模型如下：

(2)

式中:CLi,CDi分別為升力系數、阻力系數擬合模型的待辨識參數，其中i=0,1,2,3。

高超聲速滑翔飛行器一般具有高升阻比的氣動外形，且升阻比越大，機動性能越好。在研究高超聲速滑翔飛行器再入軌跡優化問題時，升阻比η計算如下：

(3)

上述式(1)和式(2)兩種模型的系數如下：

[a0,a1,a2,a3,a4]T=[0.11139, -0.019871;4.161×10-4, 2.2991, 1.2292]T; [b0,b1,b2,b3,b4]T=[0.23462, -0.02421, 7.089×10-4, -0.17481, 2.7251]T; [CL0,CL1,CL2,CL3]T=[-0.2335, 2.9451, 0.2949, -3.3943×10-4]T; [CD0,CD1,CD2,CD3]T=[0.0234, 2.3795, 0.3983, -1.0794×10-4]T

表4～6中的數據，展現了兩種模型的擬合結果。對表4～5分析來看，文獻[18]中提出的模型升力系數和阻力系數擬合方面誤差較大，但是對表6擬合升阻比結果分析時，文獻[17]中模型的擬合誤差遠遠超過文獻[18]中的模型，可以看出，使用馬赫數負指數冪形式進行擬合具有一定優勢。

2 改進氣動系數擬合模型及分析

2.1 改進氣動系數擬合模型

通過對CAV公開數據的分析，以及對文獻[17-18]中擬合模型結果的對比，可以發現文獻[17]中通過使用攻角和馬赫數的二次項進行擬合，得到的升力系數、阻力系數的擬合精度較高，但升阻比的擬合精度遠小于文獻[18]中采用馬赫數負指數冪形式的模型。因此，綜合考慮攻角和馬赫數對升力系數、阻力系數的影響，本文提出了一種新的改進擬合模型如下：

(4)

式中:li,di分別為升力系數、阻力系數擬合模型的待辨識參數，其中i=0,1,2,3,4,5。

本文模型得到升阻比η計算公式為：

(5)

2.2 多元非線性最小二乘參數辨識

本文采用多元非線性最小二乘的方法辨識模型參數。非線性最小二乘法是以各數據點誤差平方和最小為目標函數，來估計非線性模型中的參數。

本文需要辨識參數的模型，自變量為攻角α和馬赫數M，因變量為氣動系數C。則可以將氣動系數C表示如下：

C(α,M)≈f(α,M,p)=p0+p1α+p2α2+p3M+p4ep5M

(6)

式中:p=(p0,p1,p2,p3,p4,p5),f(·)為擬合模型。現已知數據(αi,Mi,C(αi,Mi))，可以將多元非線性最小二乘模型參數辨識問題轉化為尋找參數p，使得下式最小的優化問題：

(7)

式中:p為待優化參數，F(·)為擬合模型誤差。

根據表1～3中的CAV-H的升力系數阻力系數數據，采用多元非線性最小二乘的方法對本文提出的模型進行參數辨識。得到CAV-H的新的氣動系數模型的參數。

CAV-H氣動系數模型參數辨識結果如下：

(8)

2.3 擬合模型的評價指標

常用的模型擬合評價指標，包括和方差(Sum square error，SSE)、均方根誤差、均方差和平均絕對誤差。但這四種常見的性能指標，僅體現模型輸出與真值之間的誤差大小，不能衡量擬合模型輸出與真實數據變化的趨勢是否一致。擬合優度能夠較好地刻畫模型輸出與真實數據的動態變化趨勢[22]。本文在和方差擬合性能指標之外，引入擬合優度指標，進一步評價模型擬合好壞。擬合優度記作τ,τ值越接近1，表明提出的模型擬合效果越好，越能刻畫真實的數據；反之，越接近0，表明提出的模型效果擬合越差，用該模型擬合出來的氣動系數結果越不準確。擬合優度的計算如下：

(9)

注1.擬合優度是回歸問題中常用到的表述模型對數據解釋程度的指標，相對于僅使用SSE類型的指標，能夠更準確評價模型的好壞。本文將其引入到CAV-H氣動系數擬合模型的評價問題中，后續的仿真實驗證明本文改進的氣動系數擬合模型擬合優度更大，即能更準確的擬合和解釋CAV-H氣動系數數據。

2.4 模型升阻比特性分析

高升阻比飛行是高超聲速滑翔飛行器的重要特點。對高超聲速滑翔飛行器機動飛行的升阻比特性的充分討論，有利于后續軌跡優化問題研究。因此，針對本文提出的改進氣動系數擬合模型的升阻比特性進行分析。

1)改進氣動模型的最大升阻比分析

首先討論最大升阻比與攻角、馬赫數的關系，可將式(5)改寫為：

(10)

式中:lM=l0+l3M+l4el5M,dM=d0+d3M+d4ed5M。

令式(10)對攻角α求導且等于0，可得:

(11)

將α*代入式(5)，即可得最大升阻比與馬赫數的關系為：

(12)

通過式(12)，可以得到本文氣動系數模型的最大升阻比對應的攻角隨速度變化的關系，在后續的仿真實驗中會展示相應的仿真結果。

注2.高升阻比氣動特性是高超聲速再入滑翔飛行器一個重要特點，對本文改進氣動系數模型最大升阻比的討論，有利于掌握該模型最大升阻比條件下攻角與速度之間的關系，對后續應用該模型的再入軌跡優化與制導問題研究具有重要意義。

2)不同升阻比對再入飛行軌跡的影響分析

在2.4節的基礎上，討論本文改進氣動模型不同升阻比對高超聲速飛行器再入飛行的影響。

對于不同升阻比，可以得到相應的攻角-速度剖面。定升阻比情況下，式(10)可改寫為：

(13)

注3.高升阻比的氣動特性直接關系到高超聲速再入飛行器的機動能力，對改進模型的不同升阻比再入飛行討論，得到不同升阻比對應的攻角-速度剖面解析式形式，為了保持一定的機動能力，可以由式(13)得到攻角-速度剖面內的飛行走廊，能夠對考慮禁飛區約束的高超聲速滑翔飛行器再入軌跡優化與機動制導問題研究中攻角-速度剖面的設計，帶來一定的指導意義。

2.5 高超聲速滑翔飛行器再入軌跡優化模型

本文采用的高超聲速再入軌跡優化模型與文獻[19-20]中模型一致。高超聲速滑翔飛行器再入軌跡優化問題的動力學方程為：

(14)

式中:r是飛行器的地心距;V是飛行器的地球相對速度;ψ與γ分別是飛行器的航向角與航跡角;飛行器所處的經度θ緯度φ是判斷飛行器路徑約束的主要因素;m和g是飛行器的質量和當前地心距的重力加速度;飛行器的側滑角σ與攻角α分別控制飛行器制導策略中的橫向、縱向制導剖面。D=ρV2SrefCD/2與L=ρV2SrefCL/2是飛行器在飛行過程中的氣動阻力與升力，其中ρ是飛行器當前高度的空氣密度，S是飛行器的參考橫截面積，CL與CD分別是與飛行器攻角有關的空氣動力學參數。該模型中，側滑角σ與攻角α是控制量。

在高超聲速飛行器再入軌跡優化問題中，飛行器飛行過程還需滿足路徑約束、終端約束、控制量約束和狀態約束。

1) 路徑約束

在飛行過程中需要滿足熱流率約束、過載約束、動壓約束等常規路徑約束:

(15)

2) 控制量約束

(16)

式中:σmin、σmax分別為傾側角的最小值和最大值；αmin、αmax分別為攻角的最小值和最大值。

3) 狀態約束

(17)

式中:下標“min”和“max”分別表示各狀態量最小取值和最大取值。

4) 終端約束

(18)

式中:tf指終端時刻；下標“f”指期望的終端狀態。

5) 性能指標

本文進行討論的仿真實驗中，3.3.2節中的選用最大航程作為性能指標，對本文模型的三種不同升阻比下飛行性能進行考察，即本節仿真實驗的性能指標為：

J1=min(-θ)

(19)

3.3節中的選用最小時間作為性能指標，對三種不同模型，在完成任務情況下對飛行性能進行考察，即本節仿真實驗的性能指標為：

J2=mintf

(20)

3 仿真實驗與結果分析

采用多元非線性最小二乘法對本文提出的模型進行參數辨識，得到辨識結果后，將本文模型同文獻[17-18]中模型進行對比。仿真實驗包括：1)本文改進模型對CAV-H氣動系數數據的擬合，并與文獻[17-18]中的模型擬合結果進行對比；2)從升阻比的角度對本文氣動系數模型的飛行特性進行討論；3)應用本文改進氣動系數模型的再入軌跡優化仿真實驗，并與其它兩種模型進行對比。

本文仿真實驗環境CPU為Intel i7 9700處理器，運行速度3GHz。

3.1 模型對CAV-H氣動系數數據的擬合

利用多元非線性最小二乘方法對本文所提模型進行參數辨識。并將本文模型與文獻[17-18]模型的擬合結果進行對比，如圖5(a)、(b)、(c)分別是升力系數、阻力系數和升阻比的擬合結果。

圖5 三種模型對CAV-H數據擬合Fig.5 Three models fit CAV-H data

表4～6中數據分別展示了三種模型的升力系數、阻力系數與升阻比相對于真值的誤差大小，表中粗體的數據為誤差較小的數據。比較仿真實驗結果可得，本文模型升力系數擬合平均誤差相對于文獻[18]降低10%，相對于文獻[17]降低5%；本文模型阻力系數擬合誤差相對于文獻[18]降低11.4%，相對于文獻[17]同樣降低11.4%；特別地，本文模型升阻比擬合誤差相對于文獻[18]降低了26.6%，相對于文獻[17]降低52.6%，擬合精度明顯提升。可以看出本文提出的新的CAV-H氣動系數擬合模型，能夠準確地得到氣動系數，明顯地降低擬合誤差。

表4 升力系數擬合誤差Table 4 Lift coefficient fitting error

表5 阻力系數擬合誤差Table 5 Resistance coefficient fitting error

表6 升阻比擬合誤差Table 6 Lift-to-drag ratio fitting error

表7和圖6展示了三種模型擬合優度的對比，表中粗體的數據為擬合優度較大的數據。可以看出本文提出的模型，與文獻[17-18]中的模型相比，對CAV-H升力系數、阻力系數和升阻比的數據擬合優度更接近1。特別地，本文模型在升阻比擬合方面，擬合優度相對于文獻[17]提高了52.4%。由此可知，本文模型對樣本數據隨自變量的變化趨勢刻畫更為精確，對樣本數據擬合更為準確。

表7 擬合優度對比Table 7 Goodness of fit comparison

綜上所述，本文提出的改進氣動系數模型，能夠更好地描述CAV-H高超聲速飛行器的氣動系數變化，得到的氣動數據比傳統模型結果更為精確。

3.2 本文改進模型的升阻比特性分析

3.1節的仿真結果表明，本文提出的模型對數據的解釋性更好，擬合精度更高。利用式(10)～(12)對此模型的升阻比特性進行分析，有利于此模型更好的應用于高超聲速飛行器再入軌跡優化問題研究。

3.2.1改進氣動系數模型最大升阻比分析

對于本文提出新的氣動系數擬合模型，進行升阻比特性的分析。圖6和圖7是利用式(5)進行的升阻比特性分析得到的結果。圖6是升阻比隨攻角變化的曲線，由圖可知，當馬赫數一定時，隨著攻角的增大，飛行器升阻比先增后減。由圖6可知當馬赫數約為Ma14時，升阻比能取得較大的值。圖7是升阻比隨馬赫數變化的曲線，由圖可知，當攻角一定時，隨著馬赫數的增大，飛行器升阻比先增后減。由圖7可知當攻角約10°時，升阻比能一直保持較大的值。圖8～9是由式(11)、(12)，進行仿真實驗得到的結果。分析圖8和圖9曲線可知，當馬赫數為Ma13.6，攻角為9.583°時，本文提出的氣動系數擬合模型能夠達到最大升阻比約為3.679。這與圖7和圖8的結果分析相一致。需要指出的是，CAV公開的資料中最大升阻比約為3.5，本文提出的模型最大升阻比基本與此數據一致。

3.2.2定升阻比對再入飛行軌跡的影響分析

為分析本文模型在不同升阻比情況下對飛行軌跡產生的影響，選取下面三種定升阻比情況進行仿真：1)η=3.2，升阻比稍大；2)η=3.0，升阻比中等大小；3)η=2.7，升阻比稍小。

本文提出的氣動系數模型理論最大升阻比為3.679，但是從圖6～9中可以發現，當定升阻比大于3.5時，對應的馬赫數、攻角的變化范圍非常小，不利于對本文氣動系數模型的飛行特點進行分析，因此本文并未采用最大升阻比進行仿真實驗。

圖6 本文模型升阻比與攻角的關系Fig.6 Relationship between lift-drag ratio and angle of attack

圖7 本文模型升阻比與馬赫數的關系Fig.7 Relationship between lift-drag ratio andMach number

圖8 最大升阻比與馬赫數的關系Fig.8 Relationship between maximum lift-to-drag ratio and Mach number

圖9 最大升阻比對應攻角與馬赫數的關系Fig.9 Relationship between the angle of attack of maximum lift-to-drag ratio and Mach number

圖10～圖15是在本文模型基礎上進行不同升阻比的再入飛行仿真結果。可以看出，采用本文氣動系數模型進行定升阻比飛行時，升阻比越大，高超聲速滑翔飛行器再入時高度下降越慢；在相同飛行仿真時間內，升阻比越大，航程越遠，速度下降越慢，動能越大。從圖13可知，本文改進模型高升阻比飛行時，航跡角波動相對較小。另外，如圖14 所示大升阻比飛行時，對應攻角較小，這與圖9中展現的結果相一致。圖15結果可以發現，升阻比η分別取3.2、3.0、2.7時，路徑約束都能滿足條件，但大升阻比對應的這些路徑約束值相對較大。

圖10 不同升阻比飛行時高度變化Fig.10 Altitude changes with different lift-to-drag ratio

圖11 不同升阻比飛行時經度變化Fig.11 Longitude changes with different lift-to-drag ratios

圖12 不同升阻比飛行時速度變化Fig.12 Speed changes different lift-to-drag ratios

圖13 不同升阻比飛行時航跡角變化Fig.13 Track angle changes with different lift-to-drag ratios

圖14 不同升阻比飛行時攻角變化Fig.14 Angle of attack changes with different lift-to-drag ratios

圖15 不同升阻比飛行時路徑約束變化Fig.15 Path constraints changs with different lift-to-drag ratios

從上述對定升阻比仿真結果分析來看，大升阻比有利于大航程飛行，但對應的攻角較小，會導致動壓、過載和熱流率較大，因此，在再入飛行的開始階段，采用小升阻比，可以在滿足約束的條件下，使得攻角的控制范圍變的更大，有利于高超聲速滑翔飛行器再入軌跡的優化。

3.3 本文氣動系數擬合模型再入軌跡優化

采用本文提出的氣動系數擬合模型，針對CAV-H飛行器，進行高超聲速滑翔飛行器再入軌跡優化仿真實驗，并與文獻[17-18]氣動系數模型進行對比。CAV-H飛行器質量為907.2 kg，參考橫截面積為0.484 m2，再入滑翔飛行實驗的任務參數見表8，路徑約束中最大熱流率為4000 kW/m2，最大負載為2.5g(g為重力加速度)，最大動壓為60 kPa。實驗在Matlab2018b編程環境下進行，采用GPOPS2軟件進行軌跡求解。

表8 CAV-H飛行器再入飛行任務數據Table 8 CAV-H aircraft reentry mission data

圖16～19為三種模型在CAV-H再入滑翔飛行任務中得到的結果。圖16～18可以看出，三種模型都可以按預定飛行狀態到達終點，且到達終點時間幾乎一致。但仔細觀察圖16和圖17，可以發現，在各飛行階段，采用本文氣動系數模型的飛行器跳躍飛行的高度變化幅值，相對其他兩種模型較小，速度的下降也比文獻[17-18]中的模型更為平穩，從而帶來動壓、過載、熱流率的變化整體上幅值波動相對較小，對飛行器結構強度、熱防護措施等設計具有參考價值。

圖16 高度隨時間變化Fig.16 Altitude changes with time

圖17 速度隨時間變化Fig.17 Speed changes with time

圖18 航跡角隨時間變化Fig.18 Track angle changes with time

圖19 路徑約束隨時間變化Fig19 Path constraint changes with time

4 結 論

1)本文融合攻角二次項和馬赫數負指數冪項，提出一種新的氣動系數擬合模型，使得氣動系數擬合誤差大大降低，本文改進模型對數據的解釋程度更好。

2)通過對改進模型升阻比特性的討論，采用本文改進模型的高超聲速滑翔飛行器大升阻時，飛行航程有所提高，但對應的攻角較小，會導致動壓、過載和熱流率增大。

3)給出本文改進模型在不同升阻比飛行任務時，攻角-速度剖面的解析形式。在后續再入軌跡優化和機動制導問題研究時提供一定的參考。

4)與其它兩種模型的再入軌跡優化的對比實驗表明，本文模型得到的再入軌跡優化更為平穩，具有較好的潛力和應用前景。

下一步可以將本文模型應用到再入飛行走廊討論、準平衡滑翔條件制導、阻力加速度剖面設計等牽涉到氣動問題的研究中。