徐世南，吳催生

(中國(guó)空空導(dǎo)彈研究院，河南洛陽 471000)

0 引言

目前世界各國(guó)都大力投入高超聲速導(dǎo)彈的研制，試圖搶占高超聲速實(shí)戰(zhàn)化先機(jī)，以其為應(yīng)用背景的驗(yàn)證項(xiàng)目相繼進(jìn)行。但也面臨試驗(yàn)要求技術(shù)難度大、經(jīng)費(fèi)高等困難。采用計(jì)算機(jī)仿真對(duì)高超聲速導(dǎo)彈進(jìn)行流場(chǎng)、結(jié)構(gòu)溫度等準(zhǔn)確預(yù)測(cè)，可以對(duì)其熱防護(hù)與結(jié)構(gòu)設(shè)計(jì)提供指導(dǎo)[1-3]。

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

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

1 控制方程和耦合策略

1.1 控制方程

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

(1)

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

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

(2)

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

由于結(jié)構(gòu)溫度與變形存在相互耦合關(guān)系，利用有限單元法對(duì)其進(jìn)行統(tǒng)一求解?；陟o氣動(dòng)彈性，將變分原理應(yīng)用于熱傳導(dǎo)控制方程、結(jié)構(gòu)控制方程以及它們相應(yīng)的邊界條件，可得到如下形式的熱-結(jié)構(gòu)有限元矩陣方程：

(3)

式中：K為結(jié)構(gòu)剛度矩陣；KT為熱傳導(dǎo)矩陣；KuT為熱彈性剛度矩陣；u為位移向量；T為溫度向量；F為力載荷向量；QT為熱載荷向量。

1.2 耦合策略

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

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

2)在流體區(qū)域進(jìn)行氣動(dòng)熱力學(xué)仿真計(jì)算，得到熱流qf和壁面壓力pf；

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

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

5)對(duì)流體區(qū)域和固體區(qū)域計(jì)算結(jié)果進(jìn)行收斂檢查，檢查流體壓力pf和結(jié)構(gòu)溫度Ts是否滿足收斂標(biāo)準(zhǔn)，如果不收斂，繼續(xù)在t0時(shí)間步進(jìn)行迭代計(jì)算，直到滿足收斂標(biāo)準(zhǔn)或達(dá)到設(shè)置的迭代次數(shù)，如果滿足收斂結(jié)果，輸出此時(shí)間步計(jì)算結(jié)果并進(jìn)入t0+Δt時(shí)間步的迭代計(jì)算；

6)不斷循環(huán)該過程，直至設(shè)置的時(shí)間終止值，完成多場(chǎng)耦合數(shù)值計(jì)算。

2 仿真校驗(yàn)

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

圖2 二維計(jì)算網(wǎng)格

3 計(jì)算結(jié)果分析

3.1 仿真模型

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

圖3 壁面歸一化壓力分布與實(shí)驗(yàn)結(jié)果對(duì)比

圖4 駐點(diǎn)溫度隨時(shí)間變化

由于溫度場(chǎng)、壓力場(chǎng)與結(jié)構(gòu)變形之間存在耦合效應(yīng)，將影響因素分為以下3種：

1)熱因素，即仿真計(jì)算時(shí)結(jié)構(gòu)變形僅由熱變形造成；

2)氣動(dòng)力因素，即仿真計(jì)算時(shí)結(jié)構(gòu)變形僅由氣動(dòng)力造成；

3)無變形因素，即仿真計(jì)算時(shí)結(jié)構(gòu)無變形。

表1 材料參數(shù)

3.2 結(jié)構(gòu)變形分析

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

圖5 仿真模型

圖6 結(jié)構(gòu)變形圖

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

3.3 溫度與壓力分析

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

圖8 溫度與壓力曲線

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

3.4 耦合算法適用性分析

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

4 結(jié)論

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

1) 導(dǎo)彈彈體結(jié)構(gòu)發(fā)生軸向拉伸和徑向膨脹，且氣動(dòng)加熱對(duì)結(jié)構(gòu)變形影響較大。

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

3) 結(jié)構(gòu)變形對(duì)壓力和溫度仿真影響結(jié)果較小。

圖9 耦合效應(yīng)對(duì)溫度與壓力影響

圖10 不同長(zhǎng)細(xì)比仿真值對(duì)比

4) 對(duì)導(dǎo)彈進(jìn)行氣動(dòng)熱和氣動(dòng)力分析時(shí)，可不考慮耦合效應(yīng)以提高仿真計(jì)算效率。