高超聲速導彈流場與結構溫度場耦合數值分析*

徐世南，吳催生

(中國空空導彈研究院，河南洛陽 471000)

0 引言

目前世界各國都大力投入高超聲速導彈的研制，試圖搶占高超聲速實戰(zhàn)化先機，以其為應用背景的驗證項目相繼進行。但也面臨試驗要求技術難度大、經費高等困難。采用計算機仿真對高超聲速導彈進行流場、結構溫度等準確預測，可以對其熱防護與結構設計提供指導[1-3]。

對于高超聲速飛行器熱、力等載荷環(huán)境的數值模擬，傳統(tǒng)方法忽略耦合效應的影響，即未考慮結構變形與溫度、壓力之間的耦合關系。仲繼澤等[4]發(fā)展了一種基于流-固單向耦合的方法，得到了機翼的力載荷。韓玉琪等[5]通過定常流動換熱分析計算飛行器渦輪盤的熱、力載荷，傳遞給盤體得到其應力場。此方法雖然計算效率高，但當耦合效應明顯時，無法精確預示熱力學環(huán)境。隨著仿真技術的發(fā)展，考慮耦合效應影響的求解方法得到了廣泛開展，保證計算結果更加精確。周印佳等[6]采用分區(qū)求解方法，將計算區(qū)域分為流體區(qū)域和固體區(qū)域，在耦合交界面進行數據傳輸，完成了高超聲速流動與結構的耦合分析。由于此方法耦合子物理場在一個時間步內只進行一次數據交換，時間精度上僅能達到一階，安效民等[7]在流固之間引入內迭代，使分區(qū)耦合方法在時間方向上達到了二階時間精度。肖軍等[8-9]將氣動和結構動力方程各自構造子迭代求解，也獲得了每一物理時間步的高精度結果。黃唐等[10]建立和流固之間的聯(lián)系，實現(xiàn)了流場、熱、結構一體化數值模擬。

高超聲速導彈處于一個復雜的多場耦合環(huán)境，采用合適的計算方法對其進行仿真分析，能夠更有效的得到導彈熱/力載荷。文中基于分區(qū)耦合方法，建立高超聲速導彈平飛狀態(tài)的仿真模型，對其流場和結構溫度場進行仿真分析，并研究溫度場、壓力場和結構變形之間的關系。

1 控制方程和耦合策略

1.1 控制方程

在進行流-固-熱耦合計算時，各自域內遵循基本的守恒原則，其中流體控制方程采用積分形式的N-S方程[11]：

(1)

式中：Q為守恒變量；G為無黏通量；GV為粘性通量；t為物理時間；?Ω為某一固定區(qū)域Ω的邊界；dS為面積微元；n為控制邊界法向單位矢量。

對式(1)按有限體積法進行空間離散可得：

(2)

式中：VI為第I個控制體單元；NF為包圍第I個單元的所有面數；ΔSm為第m個表面的法向量。

由于結構溫度與變形存在相互耦合關系，利用有限單元法對其進行統(tǒng)一求解。基于靜氣動彈性，將變分原理應用于熱傳導控制方程、結構控制方程以及它們相應的邊界條件，可得到如下形式的熱-結構有限元矩陣方程：

(3)

式中：K為結構剛度矩陣；KT為熱傳導矩陣；KuT為熱彈性剛度矩陣；u為位移向量；T為溫度向量；F為力載荷向量；QT為熱載荷向量。

1.2 耦合策略

采用分區(qū)求解方法，將計算區(qū)域分為流體區(qū)域和固體區(qū)域。在流體區(qū)域內求解統(tǒng)一的流體控制方程得到熱流qf和壓力pf，并將其通過耦合界面?zhèn)鬟f給固體區(qū)域；在固體區(qū)域內，求解熱-結構控制方程得到結構壁面溫度Ts和位移us，并將其通過耦合界面?zhèn)鬟f給流體區(qū)域。基于ADINA軟件，完成對高超聲速導彈氣動/熱/結構的多場耦合分析，具體實現(xiàn)形式如圖1所示。

1)將初始設置的結構壁面溫度Ts和結構位移us作為仿真分析初始條件通過耦合界面?zhèn)鬟f給流體區(qū)域；

2)在流體區(qū)域進行氣動熱力學仿真計算，得到熱流qf和壁面壓力pf；

3)熱流qf和流體壁面壓力pf通過耦合界面?zhèn)鬟f給固體區(qū)域；

4)在固體區(qū)域進行熱-結構仿真計算，得到結構壁面溫度Ts和結構位移us；

5)對流體區(qū)域和固體區(qū)域計算結果進行收斂檢查，檢查流體壓力pf和結構溫度Ts是否滿足收斂標準，如果不收斂，繼續(xù)在t0時間步進行迭代計算，直到滿足收斂標準或達到設置的迭代次數，如果滿足收斂結果，輸出此時間步計算結果并進入t0+Δt時間步的迭代計算；

6)不斷循環(huán)該過程，直至設置的時間終止值，完成多場耦合數值計算。

2 仿真校驗

為驗證耦合計算方法有效性，以經典圓管繞流實驗作為算例[12]。實驗來流參數見圖2，圓管內半徑25.4 mm，外半徑38.1 mm。圓管材料為不銹鋼，材料熱力學參數可參見文獻[12-14]。圓管內壁設為等溫壁，溫度值為294.4 K。圖3為計算初始時圓管表面壓力與文獻[12]實驗結果對比，圖中p0為駐點壓強，θ為物面到圓心的連線與x軸的夾角。圖4為計算的駐點2 s的溫度變化，仿真值與文獻結果相符。驗證了耦合仿真方法的有效性。

圖2 二維計算網格

3 計算結果分析

3.1 仿真模型

采用某無翼/舵布局高超聲速導彈，不考慮內部元器件與輻射效應；導彈頭部為陶瓷材料，彈身為鈦合金材料，具體材料參數見表1。導彈仿真模型如圖5所示，長細比L/D為20，L為導彈全長，D為導彈彈徑。仿真邊界條件設置為：在AB、BC、CD邊界施加速度、壓力和溫度載荷；導彈結構為固體場，外場和內場為流體場，流體場與固體場交界面為耦合界面。仿真參數為：來流壓力、溫度、速度分別為7 494 Pa、216.5 K、5Ma，攻角為0°，初始溫度293 K，初始壓力7 494 Pa。采用瞬態(tài)計算，并基于仿真終止時間40 s進行研究。

圖3 壁面歸一化壓力分布與實驗結果對比

圖4 駐點溫度隨時間變化

由于溫度場、壓力場與結構變形之間存在耦合效應，將影響因素分為以下3種：

1)熱因素，即仿真計算時結構變形僅由熱變形造成；

2)氣動力因素，即仿真計算時結構變形僅由氣動力造成；

3)無變形因素，即仿真計算時結構無變形。

表1 材料參數

3.2 結構變形分析

研究結構變形，定義軸向變形量為ΔL/L×100%，式中ΔL(假設拉伸為正值，壓縮為負值)為結構沿X軸方向的位移量；定義徑向變形量為ΔD/D×100%，式中ΔD(假設膨脹為正值，收縮為負值)為結構沿Y軸方向的位移量。橫坐標進行歸一化處理，坐標數值為：x/L，導彈頭部最前端位置為坐標原點，x表示導彈軸向位置，L為導彈全長。其變形量如圖6所示，結構沿X軸發(fā)生軸向拉伸，沿Y軸發(fā)生徑向膨脹。

圖5 仿真模型

圖6 結構變形圖

以頭部末端位置為例，研究氣動熱和氣動力對結構變形的影響，如圖7所示，在導彈平飛狀態(tài)，熱因素使結構發(fā)生拉伸和膨脹，氣動力因素與之相反，其中結構變形主要由熱因素引起。

3.3 溫度與壓力分析

導彈外壁溫度與壓力分布如圖8所示，頭部溫度與壓力載荷值大，越靠近頭部位置，來流發(fā)生的摩擦和壓縮越劇烈，動能轉化成的熱能越多，速度的變化梯度也越大。

圖8 溫度與壓力曲線

以頭部末端位置為例，研究耦合效應對溫度和壓力的影響，如圖9所示，在導彈平飛狀態(tài)，耦合效應對溫度和壓力的影響較小，結構變形量不大，不足以改變流場與結構溫度場分布。

3.4 耦合算法適用性分析

雖然考慮耦合效應，采用雙向耦合算法仿真結果更精確，但是計算效率比不考慮耦合效應的低，在實際應用中，仿真值在許用誤差范圍內即可。通過改變彈身長度研究不同長細比下耦合效應對導彈氣動熱、氣動力載荷的影響。以頭部末端位置為例，如圖10所示，圖中ΔT、ΔP分別為考慮耦合效應與不考慮耦合效應的溫度差和壓力差，P為考慮耦合效應的壓力值。由圖10可知，導彈在不同長細比下，耦合效應對仿真計算結果影響均較小，實際分析時可忽略耦合效應以提高計算效率。

4 結論

基于高超聲速導彈平飛狀態(tài)分析得到以下結論：

1) 導彈彈體結構發(fā)生軸向拉伸和徑向膨脹，且氣動加熱對結構變形影響較大。

2) 頭部位置氣動熱與氣動力環(huán)境較為嚴酷。

3) 結構變形對壓力和溫度仿真影響結果較小。

圖9 耦合效應對溫度與壓力影響

圖10 不同長細比仿真值對比

4) 對導彈進行氣動熱和氣動力分析時，可不考慮耦合效應以提高仿真計算效率。