超聲速/高超聲速飛行器氣動力快速估算平臺設計及應用

程鋒, 唐碩, 張棟

(1.西北工業大學 航天學院, 陜西 西安 710072; 2.陜西省空天飛行器設計重點實驗室, 陜西 西安 710072)

高超聲速飛行器氣動性能評估影響著其設計的成敗,是基礎設計要素[1]。一方面,高超聲速飛行器因其復雜的系統特性,不能完全在飛行條件下或者風洞環境里測得飛行器性能評估所需要的所有氣動數據[2];另一方面,CFD計算所需要的硬件設備和耗時都很高,難以滿足初步設計階段大量設計方案快速性能評估的需求[3]。基于氣動理論的快速氣動性能估算平臺的開發為此提供了一個可行的解決方案,其不僅可以降低計算對硬件和時間的要求,還能計算飛行實驗和風洞實驗所不能涵蓋的狀態點氣動數據,是高超聲速飛行器初步設計階段設計構型氣動性能快速評估和設計優化的有力工具。

美國NASA Langley研究中心在20世紀80年代開發了APAS(aerodynamic preliminary analysis system)[4],用以快速評估飛行器氣動設計性能,其中的算法包含牛頓方法[5]、激波法以及Pandtl-Meyer膨脹波方法[5]等。Cruz等[6]將其應用于計算HL-20的氣動性能,并和風洞實驗結果進行對比,發現其有較好的吻合程度。另外,基于構型的氣動工具CBAERO(configuration based aerodynamics tool)是一個基于非結構網格求解Euler方程的的求解器,使用三角形面元定義飛行器的外模線構型,而不需要體網格,是NASA第二代運載器氣動性能快速評估工具[7]。其顯著的特點是:不需要體網格,計算速度快;使用Euler方程,計算精度較高;構型適應性好等。最近,The Aerospace Corporation的Lobbia等[8-11]開發了適用于多學科設計優化(MDO)的飛行器氣動性能快速評估工具FAAT(fast aerodynamics analysis tool),其使用牛頓碰撞理論等方法基于C/C++、OpenMPI開發,具有較好的構型適應能力。

國內關于氣動力快速計算研究的公開文獻方面,黃志澄[12]基于美國開發的超聲速/高超聲速任意物體程序(S/HABP),介紹了一種超聲速和高超聲速一致適用的壓力計算方法。李治宇等[13]應用快速氣動分析軟件,通過求解Euler方程得到高超聲速飛行器氣動力特性,并以此為基礎對飛行器使用遺傳算法等優化方法進行了氣動外形的優化。郝佳傲等[14]改進和發展了一套適于有翼再入飛行器氣動布局的部件劃分策略和壓強計算選取準則,并對航天飛機和類X-43飛行器進行了計算,獲得了與實驗結果較吻合的結果。龔安龍等[15]以經典參考溫度法為基礎,基于無黏Euler方程CFD方法獲得的壁面流場參數發展了一種較精確的高超聲速飛行器壁面黏性力計算方法。

高超聲速氣動力快速估算程序對于高超聲速飛行器設計以及性能評估具有重要的影響,而由于絕大部分程序的閉源和非開放特性,限制了高超聲速飛行器初步氣動估算的發展;另一方面,隨著高超聲速飛行器的發展需求,前文所述的工具難以滿足高超聲速飛行器氣動力快速評估的需求,例如摩阻估算、三維效應的考慮等。本文開發的高超聲速氣動力快速估算平臺(aerodynamic forces preliminary evaluation,AFPE)采用了基于流線的壓力和摩阻估算策略,在獨立于研究對象的要求下,保證計算速度的同時,較傳統的基于自由流的策略有較高的精度和可靠性,為高超聲速飛行器初步設計階段構型遴選和性能評估提供了基礎。

1 算法理論

AFPE最主要的核心在于:沿流線參數計算和摩阻的估算,為了實現這2點要求,需要提供任意機體表面上流線的追蹤方法和基于傳統方法的組合計算模型。由于傳統壁面氣動力計算方法存在各自的應用范圍限制,例如牛頓法和切楔/切錐法不能用于背風面氣動力計算,激波方法在鈍頭體頭部計算能力不足等等,需要一種策略來合理選擇算法。

1.1 流線追蹤技術

相比于基于自由來流的高超聲速飛行器全機氣動力估算方法,基于流線的氣動力估算方法在計算當前點的氣動力的數值時,不僅會考慮自由來流條件的影響,而且還會考慮經過當前點的流線受機體壁面的累積氣動影響,如圖1所示,因而后者計算結果會更接近實際飛行過程中的狀況。

圖1 沿流線流動示意圖

基于流線計算策略的重點是流線的確定,流線是計算壁面氣動壓力分布和摩擦力分布的基礎。當前面元上單位速度矢量為

Vi=ni×V∞×ni

(1)

式中,Vi為當地面元上的單位速度矢量,ni為當地面元的單位外法向矢量,V∞自由流單位速度矢量。在不考慮邊界層的影響下,無黏流線是計算壁面氣動參數的合理選擇,當考慮黏性邊界層影響時,此處計算的流線是真實流線的一種近似解,但為了簡化計算過程,在此假設(1)式獲得的流線即為壁面上的流線。

1.2 修正牛頓-激波/膨脹波模型

激波膨脹波理論存在著最大壁面傾角限制,當壁面與來流夾角大于這個限制角度時,激波將會離體,二維激波模型和Taylor-Maccoll方程都不能處理激波離體的情況,而牛頓理論在大傾角時的使用不受激波離體的限制,因而將牛頓理論與激波/膨脹波理論結合起來使用可以覆蓋幾乎所有的高超聲速飛行器外流場壁面。修正牛頓理論采用如下系數和指數組合修正的方式

Cp=Cp,maxsinNθ

(2)

式中,Cp,max為滯止點最大壓力系數,N為修正指數。實際上滯止點的壓力系數是由Ma∞=0時的1增長到Ma∞=1時的1.28,當γ=1.4,Ma∞→∞時,Cp=1.86(γ=1,M∞→∞時,Cp=2)[16]。另外,牛頓碰撞理論在較小的碰撞角使用時,容易引起較大的誤差,通過系數的修正很難匹配滯止點附近和小碰撞角附近的壓力分布。牛頓理論中,如果指數小于2將會在保持滯止點壓力分布的同時,獲得較好的小碰撞角區域壓力分布,此時指數逼近1.86。

激波/膨脹波模型使用典型的二維/三維激波理論和Pandtl-Meyer膨脹波理論,分別計算迎風面和背風面的壓力系數分布,其中三維激波理論由Taylor-Maccoll方程組給出：

(3)

式中,θ為當前點的傾角,vr和vθ分別為徑向和法向速度。

1.3 修正牛頓-切楔/切錐模型

將牛頓理論運用于切楔/切錐理論不能適用的鈍頭體頭部區域,在壁面傾角降低到一定值之后再使用切楔/切錐方法求解高超聲速飛行器壁面壓力分布,能夠擴大牛頓和切楔/切錐方法的應用范圍。圖2為修正牛頓-激波/膨脹波模型和修正牛頓-切楔/切錐模型對一鈍頭體壁面壓力計算結果與高精度CFD計算結果的對比,由圖 2中可以看出,不論是修正牛頓-激波/膨脹波模型還是修正牛頓-切楔/切錐模型都比傳統的單一模型有著更寬的應用范圍和更好的整體逼近程度。

圖2 修正牛頓-激波/膨脹波模型、修正牛頓-切楔/切錐模型和CFD計算結果對比

1.4 底壓模型

對于暴露在高超聲速流中的物體,我們期望物體上任何底部面積都承受著總真空壓,但是真實氣體效應的黏性會使有些壓力作用在底部,而實驗數據也顯示這時的壓力系數大約為70%的大氣真空壓[4],可以使用以下理論修正公式計算底部壓力

(4)

另外,Gaubeaud底壓計算模型可用來計算飛行器背風面角度大于45°時的壓力系數[17]

(5)

圖3為使用2種底壓計算方法對HL-20的底部壓力進行計算并與實驗測量值[18]進行對比的結果。由圖3中可以發現,在較低馬赫數理論修正方法和實驗測量值很吻合,而在較高馬赫數理論修正方法和Gaubeaud模型并沒有太大的差別。

圖3 底壓計算結果與實驗測量值對比

1.5 不可壓層流-van Driest II摩阻計算模型

van Driest II摩阻計算模型[19-20]是基于使用von Karman混合長度來積分動量方程的,因為是對使用Prandtl混合長度積分的摩阻系數方法van Driest方法的修正,所以本方法被稱為van Driest II方法。van Driest表面摩擦因數公式為

(6)

動量厚度雷諾數是由動量方程沿無黏表面流線積分所得,進而,不可壓變形動量厚度雷諾數由下式計算

(7)

式中，Fe=μe/μw,使用Karman-Schoenherr公式,這個變形動量厚度雷諾數用來計算不可壓變形表面摩擦因數,變形無黏摩擦因數通過下式轉變為可壓表面摩擦因數

(8)

式中

(9)

Fc表達式的求解由理想氣體方程封閉,而Crocco邊界層溫度分布和溫度恢復因子0.9可用來估算Fc的值。

通過推導,將可壓流和不可壓流之間的轉換公式概括如下:

(10)

(11)

(12)

(13)

(14)

不可壓層流模型在平板前部對摩阻因數有很高的重現能力,但在湍流區,van Driest II模型能夠獲得更好的計算結果,將2種方法以一定的方式結合起來,在轉捩之前使用不可壓層流模型,轉捩之后使用van Driest II模型。在求得轉捩雷諾數和轉捩區長度以后,就可以對轉捩前的層流區域和轉捩后的湍流區域分別使用相應的壁面摩擦力計算方法進行計算,并在轉捩區使用插值算法。其中轉捩雷諾數由下式來確定[21]

(15)

式中,Rext為當地轉捩雷諾數;Mae為邊界馬赫數。轉捩區長度由以下關系式計算[22]

(16)

式中,xt,s和xt,e分別為沿流線的轉捩區開始和結束位置;Rext是轉捩區雷諾數。

2 平臺開發

本文研究的高超聲速飛行器氣動力快速計算平臺AFPE的架構如圖4所示,總體上由平臺部分、基礎和核心計算部分以及可自定義的算法模型數據庫部分構成。在假設飛行器為剛體結構的情況下,數據計算流程如圖5所示。首先為快速計算平臺準備好算法數據庫、非結構三角形網格文件和計算條件及設置文件。程序在讀入這些文件后判斷文件的完整性,然后進行網格轉化,將原始網格文件輸出為可讀性更強的格點文件和面元文件。在這之后,程序開始計算面元的幾何參數,如面心、外法向矢量、面積、面元上的速度矢量等,將計算結果存入面元文件。根據計算所得的面元幾何參數給定的計算條件計算飛行器外表面流線。計算流線的要求是:外表面所有面元上必須至少有一條流線經過,將計算所得的流線數據存入流線文件。接著由計算所得的流線和流線上的流動參數,利用算法庫中提供的算法和算法選擇策略計算沿流線的氣動參數以及壓力和摩擦力,將計算結果存入流線文件。最后將計算所得流線氣動數據重新分布到格點和面元上,并以此為基礎自動生成MATLAB網格和流線數據顯示腳本和數據、Tecplot數據文件,判斷是否完成所有指定點的計算。程序中所有的輸入輸出文件以ANSI編碼,具有很好的可讀性和修改性,便于后續程序功能升級和與優化設計平臺、彈道仿真平臺的協同數據交互。

圖4 高超聲速飛行器氣動力快速估算平臺架構

圖5 高超聲速飛行器氣動力快速計算平臺數據流程圖

3 驗證與分析

以公開文獻中實驗數據較多的HL-20升力體高超聲速跨大氣層飛行器為參考,分析所建立的高超聲速飛行器氣動力快速估算平臺計算結果。HL-20屬于面對稱飛行器,取飛行器一半作為離散化對象構型。機體頭部和座艙前部氣動參數變化梯度較大,在這些地方將網格加密有助于快速估算平臺更準確地評估飛行器所受的氣動力。在較為平坦的區域,如機體下部和背部靠后位置,由于氣動力沿縱向變化梯度較小,而沿橫向梯度相對較大,故為了節約計算量,這些區域的網格橫向密度大于縱向密度。所得的簡化非結構網格如圖6所示。

圖6 HL-20簡化網格和流線

圖6中還提供了使用AFPE計算所得的HL-20在Ma4.5時2種不同迎角狀態下的沿壁面流線圖。由圖6中可以看出,流線在小迎角或者較小負迎角時流線都集中于機體頭部和翼前緣,而在較大迎角,流線明顯地發源于頭部和機體下表面。在V翼下

表面,大迎角時流線在翼面下部非常集中,也就是說此時V翼翼面下部壓力較大,為飛行器提供足夠的低頭力矩,這和對HL-20的氣動分析結果一致,說明快速計算平臺所采用的流線計算方法能夠真實的反映物理現象本質。

根據HL-20的風洞實驗數據,分別計算Ma4.5和Ma10.07情況下的氣動受力情況,計算過程中機體雷諾數保持與實驗值相同,圖7和圖8為本文開發的高超聲速飛行器氣動力快速估算平臺AFPE對HL-20在2種不同的馬赫數和與實驗對應的條件下計算的各氣動力系數。

對于Ma4.5(見圖7),AFPE估算所得的升力系數稍低于風洞實驗值[6,18],這主要是因為HL-20構型的下部較為平坦,可以使用二維壁面壓力計算方法,而圖7的計算中全部使用了三維方法,三維釋壓效應使得HL-20下部壓力平均水平低于風洞測量值,造成升力系數偏低。阻力系數在很大范圍內和風洞實驗測量值很吻合,甚至于優于APAS估算的結果。由于預測了較低的升力系數和很吻合的阻力系數,AFPE得到的升阻比低于風洞測量值。俯仰力矩方面,相比于實驗測量值,APAS預測較吻合,CBAERO預測較小,而AFPE估算得到了較大的力矩系數。另外,需要注意的是AFPE預測的俯仰力矩系數在零迎角附近呈略增大的趨勢,造成這一現象的原因可以歸結于HL-20壁面壓力計算都使用了三維模型、幾何建模誤差和網格密度。

圖7 Ma4.5 AFPE氣動力系數估算結果

Ma10.07的氣動系數曲線與Ma4.5類似,如圖8所示。由于機體幾何建模誤差和網格密度的影響,而致使升力因數和升阻比略小,俯仰力矩與實驗測量值略有差別,但阻力系數與實驗測量值[23-24]很接近,說明摩阻系數模型至少能夠滿足馬赫數4～10范圍內的氣動力估算。Micol[23]提供了Langley開發的LAURA(langley aerothermodynamic upwind relaxation algorithm,基于有限體積法,使用隱式點松弛算法求解包含化學反應和熱非平衡的高超聲速三維NS黏性流方程)工具,在文獻[23]中提供了基于Euler方程的LAURA計算結果,如圖8所示。相比于CFD計算結果,AFPE在升力系數和升阻比計算時與實驗有一些誤差,而阻力系數計算結果與LAURA CFD計算結果都很接近實驗測量值。

圖8 Ma10.07 AFPE氣動力系數估算結果

通過對比圖7和圖8,我們可以看出,隨著馬赫數的增加,升力系數與實驗測量值的誤差逐漸減小,這是因為隨著馬赫數的增加,氣動系數逐漸趨向于與馬赫數無關,即高超聲速馬赫數無關理論。另外,俯仰力矩系數作為其中最難計算準確的一個氣動系數,在AFPE和實驗測量值及其他計算工具、CFD工具對比中可以發現,AFPE計算所得的俯仰力矩系數精度是滿足初步設計時氣動力估算的需求的。同時,AFPE雖然預測所得的升力系數等存在一定的誤差,但是預測所得的氣動力系數變化趨勢和實驗測量結果很吻合,說明快速計算平臺的設計思路和方法是正確的。

圖9為快速估算平臺計算的HL-20在Ma4.5時的氣動力系數云圖。可以看出在不同迎角時,壓力系數的分布存在著很大的差異,相比于負迎角,大的正迎角能顯著增大機體下部和翼面下部的壓力系數分布,這一點也可從圖6中的流線分布得出相同的結論。另外,壓力系數分布變化最明顯的區域是座艙頭部區域,相比于小迎角狀態,大迎角時座艙頭部壓力系數明顯減小,因為這時座艙頭部已處于機體頭部的遮擋之下,受氣流壓力逐漸減弱。從氣動系數云圖和曲線來看,AFPE能夠很好地從物理機理上快速估算高超聲速飛行器在較寬馬赫數范圍內的氣動受力情況,滿足設計初衷和基本要求。

圖9 Ma 4.5壓力系數云圖

4 協同仿真

對于某基于RBCC動力的水平起飛/水平著陸的亞軌道可重復使用運載器一級構型[25-26],利用本文開發的AFPE氣動力估算平臺進行超聲速和高超聲速氣動力估算,并利用計算的氣動數據進行飛行器彈道仿真。假設飛行器只在縱向平面內運動,則其彈道方程如(17)式[27]所示,式中V為速度;α為迎角;θ為彈道傾角;ωy為俯仰角速度;My為俯仰力矩,Jy為繞y軸的轉動慣量;x,y為飛行器空間坐標;?為俯仰角;m為飛行器質量;ms為燃料消耗引起的飛行器總質量變化。

(17)

圖10為仿真計算所得的飛行器飛行動力學結果。由圖10中可以看出,飛行器在開始后一段時間內有加速運動,馬赫數逐漸升高,但是隨著高度降低,阻力逐漸增大,飛行器逐漸減速。由于本文在仿真時并沒有使用復雜控制制導方法,因而在大約t=160 s以后飛行器俯仰角速率開始發散,飛行器姿態逐漸失去控制。因此,要使AFPE計算所得的氣動數據適用于此飛行器,需要進行控制系統設計,然而這并不影響通過此仿真驗證AFPE和彈道計算平臺的協同仿真能力。

圖10 RBCC動力的飛行器動力學協同仿真結果

5 結 論

本文根據目前公開的高超聲速氣動力快速估算工具的不足和工程實際需求,基于在壓力/摩擦阻力理論/工程估算方法的基礎上提出的組合模型數據庫,使用流線追蹤方法建立了高超聲速飛行器氣動力快速估算平臺AFPE,通過使用該平臺與CFD計算數據、實驗數據和其他計算方法的對比驗證，并在實際高超聲速飛行器彈道計算中的應用,得出以下結論:

1) AFPE能夠快速且準確地估算較寬馬赫數范圍內的高超聲速飛行器氣動受力情況。相比于大多數的估算工具,AFPE不僅能夠計算壓力系數分布,還能夠使用組合摩阻計算模型計算摩擦阻力分布,因而相比于傳統的快速估算工具,AFPE計算可靠性和精度會有較大提升。

2) 與實驗數據、CFD數據以及其他計算工具結果的對比認為AFPE在氣動力系數發展趨勢上預測很吻合,但平均誤差相比于經典的氣動力估算工具和CFD略大,即使如此,AFPE依然能夠以其較強的適應范圍和計算時間彌補誤差的不足。

3) 通過與彈道仿真平臺的協同仿真計算,驗證了本文開發的高超聲速氣動力快速估算平臺AFPE的協同仿真能力,同理也可以與構型優化分析平臺協同,為研究人員和工程技術人員在高超聲速飛行器初步方案遴選、性能評估以及設計優化階段提供強有力的工具保證。

4) AFPE目前處于初期階段,需要在以后的研究中逐步提高其穩定性和應用范圍,例如加入壁面熱流的計算模塊,這些將進一步提升AFPE的適用性。