基于樣條曲線描述的超聲速噴管型面優化設計

吳盛豪, 廖達雄, 陳吉明, 陳欽, 裴海濤

(中國空氣動力研究與發展中心 設備設計與測試技術研究所, 四川 綿陽 621000)

超聲速噴管的目的是保證超聲速風洞試驗段獲得設計馬赫數，并使試驗段馬赫數分布均勻、保持平直氣流。在傳統二維噴管設計方法中，一般基于特征線原理設計噴管位流型線[1]，在一定的邊界層假設條件下，對噴管的位流型面進行邊界層修正[2]，得到噴管型面曲線。隨著計算流體力學和計算機技術的快速發展，分析設計的方法(design by analysis)被引入到噴管設計中[3-4]，廣泛應用于火箭發動機噴管和高超聲速飛行器尾噴管優化設計[5-6]，并取得良好的應用效果。超聲速風洞噴管與火箭發動機和超聲速飛行器尾噴管設計存在一定相似性，但對流場品質要求更高，擬將基于CFD的搜索算法應用到超聲速風洞噴管設計中，以期獲得更優的試驗段流場品質。在實際應用噴管型面優化設計時，面臨2個重要問題：噴管型面的參數化和優化目標的選取。噴管型面參數化的方法較多，高超聲速飛行器的尾噴管可采用Akima三次曲線作為參數化的方法[5]，將燃燒室出口高度及噴管下板長度和偏轉角作為設計變量，火箭發動機噴管可采用2段三次樣條曲線進行構造[7],由于本文研究的超聲速噴管對流場分布的均勻性要求較高，因此對型面的要求也更為精細，擬采用多段Spline樣條曲線對噴管進行參數化；優化目標的選取直接關系到優化設計結果的質量，考慮到多段樣條曲線參數化的噴管，流場品質不確定度較高，擬選擇試驗段軸向馬赫數與設計值的均方根誤差作為優化目標，即可保證噴管達到設計馬赫數的同時具有較優的流場分布。

1 優化問題的描述

1.1 優化設計對象

本文研究對象是中國空氣動力研究與發展中心的0.6 m連續式跨超聲速風洞的二維半柔壁型式低超聲速噴管(見圖1)。

圖1 0.6 m連續式風洞噴管結構示意圖

風洞噴管膨脹段位流型面采用具有連續曲率的噴管型面設計方法——Sivells方法得到,利用經驗估算的方法進行邊界層修正[2]，即近似認為噴管超聲速段型面上的邊界層順氣流方向線性增長。該風洞在低超聲速運行過程中，采用可調柔性板模擬噴管型面，試驗段為四壁實壁，為消除氣流流動過程中沿壁面邊界層逐漸增厚對試驗段流場品質的影響，試驗段的上下壁板擴開角可在(-0.5°,1°)范圍內實現連續調節。

1.2 優化問題的定義

風洞流場校測試驗中，利用(1)式計算風洞模型區內沿軸向馬赫數與均值的標準差(σM)作為表征流場均勻性的指標,據此本文選擇(2)式計算得出的試驗段實際軸向馬赫數與設計值的均方根誤差(ReM)作為優化目標。

(1)

(2)

式中,Mi為試驗段內軸向測點馬赫數,Mdesign為噴管設計馬赫數,本文選擇M1.5為設計對象。據此可提出如下優化設計問題:

優化目標: minReM

設計變量:h={h0,h1,h2…hn,θ}T

式中,h為設計變量,其中h0為喉道相對中心軸線高度,其后h1…hn為插值點高度;θ為噴管出口傾角亦為試驗段壁板擴開角;bl,bu為設計變量的上下限,Ah≤b限定了分變量h1～hn按照遞增的順序排列。

優化問題中模型區內ReM直接反映流場波動情況,當ReM→0可判定模型區內斜激波或膨脹波基本被消除,流場均勻性達到要求;同時由于ReM是試驗段實際馬赫數與設計的誤差,所以當其趨于0時,可判定試驗段馬赫數達到設計值。

1.3 CFD計算驗證

優化設計過程中涉及到對候選噴管型面進行氣動性能評估,選擇CFD的方法進行評估。對風洞高速段建立計算模型,計算域包含收縮段、噴管段和試驗段,利用對稱條件建立四分之一模型。噴管選擇M1.5型面噴管靜調值,試驗段壁板擴開角設為風洞調試過程中選擇的0.1°,流動條件為風洞常規運行工況:總壓100 kPa,常溫,干燥空氣介質。選擇kw-sst湍流模型求解,給定高速段進出口壓力作為邊界條件,入口湍流度按試驗測量值0.4%給定。

在0.6 m連續式風洞中,M1.5低超聲速運行時,利用單點總壓探頭(如圖2所示)移測的方式對試驗段第一菱形區和模型區內軸向馬赫數分布進行測量。其中單點總壓探頭直接測得波后總壓,再利用(3)式通過迭代求解出測點馬赫數。

圖2 核心流波后單點總壓探頭

(3)

圖3給出沿試驗段中心軸向馬赫數分布的試驗和數值模擬結果,由于噴管出口和試驗段入口存在折角,造成模型區內馬赫數軸向分布出現劇烈波動,試驗得到模型區(距試驗段入口1 000～1 600 mm范圍)流場指標:馬赫數標準偏差σM=0.017;軸向馬赫數梯度dM/dx=0.060。數值計算得到模型區內流場指標:馬赫數標準偏差σM=0.012;軸向馬赫數梯度dM/dx=0.049。

圖3 試驗段軸向馬赫數分布

對比數值模擬結果和試驗結果,試驗段沿軸向馬赫數分布規律吻合較好,σM和dM/dx指標計算誤差在允許范圍內。同時圖4給出該狀態下,CFD計算得到的試驗段縱切面上密度梯度云圖,從圖中可以看出由于噴管出口和試驗段壁板傾角不同,在連接處出現一道弱斜激波,經過壁板反射之后,造成試驗段內的馬赫數波動,與試驗結果相同。

圖4 試驗段縱切面密度梯度云圖(入口坐標:0)

上述CFD仿真計算可獲得較為精確的試驗段馬赫數分布,對激波和膨脹波的捕捉也較為可靠,因此可將該CFD計算方法用于對候選噴管型面的氣動性能評估。

2 噴管型面曲線參數化

2.1 噴管型面構造

喉道前的收縮段曲線通過(4)式描述的雙三次曲線給出,式中:xm為兩曲線前后連接點;D2收縮段出口截面半徑,亦為喉道高度;D1收縮段進口截面半徑;D軸向距離為x處的截面半徑。

(4)

對于喉道后的擴張段,撓性噴管利用撓性板的彈性曲線模擬理論氣動型面,其曲率變化一定是連續的,因此選擇具有二階連續導數的三次Spline樣條曲線擬合噴管氣動型面可保證氣動型面具有連續曲率。如圖5所示，P0,P1,…,Pn,Pout為插值節點,P0喉道位置,Pout為噴管出口位置。

圖5 噴管型面曲線構造過程示意圖

為求解三次樣條插值函數的系數矩陣,給出端點處的第一邊界條件即一階導數,喉道處型面傾角為0,噴管出口處型面傾角取為壁板擴開角θ,據此可給出(5)式所示的第一邊界條件。

(5)

2.2 插值點的選擇

所選插值點應能夠精確擬合噴管型面,而增加插值點數目導致設計變量增多,影響求解的效率,因此需對所選插值點的個數和位置進行研究,以找出最佳方案。

以0.6 m連續式跨超聲速風洞原M1.5噴管設計型面為參照對象,將擬合型面與原設計型面曲線最大高度差δmax作為評估,找出最佳插值點數目和分布方案,使得δmax≤0.05 mm(噴管型面靜調誤差)。

利用噴管出口高度和擴張段長度分別對型面縱坐標和橫坐標做無量綱處理,利用搜索算法以minδmax為目標,搜尋不同型面控制點數目(3～8)下的,最優控制點分布方案。從表1中數據看出,隨著控制點數目的增加,擬合型面與原型面誤差逐漸降低,當控制點數目為6時,按照搜尋到的最佳節點分布方案,擬合誤差為0.01 mm，達到0.05 mm的誤差要求。

表1 不同分布插值點的型線擬合誤差

3 基于高斯過程模型的重啟優化算法

在應用優化算法進行噴管型面設計時,面臨2個問題:優化算法的全局特性和氣動性能CFD計算的效率問題:為提高全局特性,本文構造梯度優化算法的重啟機制;同時利用高斯過程模型構建設計變量與氣動性能指標之間的關系,以減少CFD評估次數。

3.1 重啟搜索算法

基于梯度的最優化算法具有求解次數少、收斂速度快的優點。但是這類方法都屬于局部最優化方法,能否收斂到全局最優點,取決于搜索的起始點。為提高梯度搜索的全局特性,借用NSGA算法中提出的擁擠距離概念,構建梯度算法的重啟機制[8]。可在由某初始點x0執行梯度搜索到達局部最優點之后,重新產生新的起始點xrs搜索。為了提高算法的效率,則需要使xrs位于沒有探索過的區域。這個條件可以用最大化擁擠距離來定義。假設已評估的設計點集合X?D,新的設計點x∈D,定義擁擠距離為如下的函數:

(6)

3.2 高斯過程模型介紹

(7)

3.3 算法框架優化準則

基于高斯過程替代模型的優化算法框架如圖6所示。

1.procedure NOZZLE OPTIMIZATION(lb,ub)2. Input:nmax,Ω3. D←[bl,bu]4. neval←05. D←?6. K←?7. n←dim(D)8. n0←11n-19. d0←{d10,…,dn00}?D10. Evaluation:κ←{κ(d0,Ω)}11. neval←neval+n012. D←D∪d0,K←K∪κ13. Modelling:Φ:DK←DK14. while(neval

圖6 基于替代模型的優化算法框架

D為設計空間(設為n維),Ω為流動條件,κ為所選的氣動性能參數。首先從設計空間中利用拉丁超空間試驗設計方法(Latin hyperspace design, LHS)選擇11n-1個初始采樣點,并評估其氣動性能。根據已評估的設計點參數和氣動性能建立高斯過程模型;利用高斯過程替代模型來探索設計空間,選擇需要評估的候選設計點;評估候選設計點的氣動性能。然后重復建模-選擇候選點-評估候選點的過程,直到得到滿意的設計結果或者優化資源消耗完畢。在這個算法中,實際的氣動性能評估計算只發生在圖6中的第10步和16步,優化過程則基于替代模型完成(第15步)。

3.4 算法性能驗證

針對經典2維算例,限定計算目標函數的值100次,考察算法的性能,其優化結果如表2所示。在這6個函數中,Sphere和Ellipsoid是單模態函數,其余均為多模態函數。函數最小值都為0。從表2中可以看到,除模態數量極為龐大的Griewank函數和Rastrigin函數之外,均可以得到1×10-4量級的優化結果。

表2 基于高斯過程模型的重啟優化算法性能

4 噴管優化設計結果

針對0.6 m連續式跨聲速風洞,M1.5型面、運行總壓100 kPa的工況,對1.2節中描述的優化問題開展單目標優化。優化設計過程中允許的總氣動評估次數為1 000次,梯度算法最大迭代次數2 000,重啟點個數89,初始采樣點個數81,同時CFD進行氣動評估的過程中利用前次結果為后次計算賦初值,以加快收斂。氣動評估次數達到600次后,分別在第534次和第561次得到較優的2個設計結果,考慮到設計變量多達8個,因此求解效率相對較高。

表3給出了優化設計結果,結合圖7給出的噴管擴張段型面曲線,可以看出,2個優化設計結果型面曲線差別較小,同時流場品質指標也較為接近。

與原設計型面對比,優化設計型面喉道高度由243.55 mm變為252.9 mm有效解決了原噴管型面馬赫數偏高的問題,提高了噴管設計型面馬赫數的精準度。

圖7 噴管擴張段型面曲線對比

圖8 試驗段軸向馬赫數分布

優化設計結果計算步數變量設計結果/mmh0h1h2h3h4h5h6θ流場品質σMdM/dxorigin243.55246.43271.61275.83294.00297.45298.740.100.0120.049 result 1534252.98256.70274.53277.99297.19298.63299.630.110.001 30.001 result 2561252.93256.41274.65278.03297.12298.66299.620.120.001 90.003

圖8給出優化設計得到的試驗段軸向馬赫數分布,結合表3的流場指標,可以看出優化設計得到的試驗段內馬赫數分布均勻性較好,流場均方根偏差由0.012降低到0.001～0.002左右,試驗段模型區內軸向馬赫數梯度在0.001～0.003范圍內,滿足設計指標的要求。由于本將試驗段與噴管進行一體化設計,即設計變量之一的噴管出口角度θ和試驗段壁板擴開角相等,改善了噴管出口和試驗段入口曲率不連續的問題,試驗段內的弱壓縮波和弱膨脹波較原型面減弱明顯,從圖9優化設計前后密度梯度分布云圖亦可看出。

圖9 優化前后試驗段密度梯度分布云圖

5 結 論

基于本文提出的噴管優化設計框架可有效地獲得流場品質較高的超聲速噴管型面:選擇具有二階連續導數的三次spline樣條曲線和最優插值點分布方案可較為精確地描述噴管型面;將試驗段實際軸向馬赫數與設計值的均方根誤差(ReM)作為優化目標,可在確保設計馬赫數準確度的同時保證試驗段流場品質;基于高斯過程模型的重啟優化算法,改善了優化問題全局特性的同時較大地提高了求解效率。

基于本文提出的優化框架,均方根偏差降低到0.001,仍有提升的空間,下一步可考慮將算法框架擴展到包含多目標優化、多設計變量如噴管長度、最大膨脹角等參數,以進一步提升噴管性能。