高超聲速導(dǎo)彈多場耦合仿真

徐世南，吳催生

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

0 引 言

新一代高超聲速導(dǎo)彈采用超燃沖壓發(fā)動機等新型動力，從而實現(xiàn)短時間內(nèi)打擊長距離軍事目標(biāo)[1]。為了保證燃料供給充分，又不增大飛行阻力，彈體結(jié)構(gòu)需要在減小彈徑的同時，增加發(fā)動機艙長度，因此大長細比成為高超聲速導(dǎo)彈結(jié)構(gòu)形式的一個發(fā)展方向。此結(jié)構(gòu)形式保證導(dǎo)彈實現(xiàn)高速度、超機動性和廣泛作戰(zhàn)空間的同時，也增大了結(jié)構(gòu)的柔性，導(dǎo)致結(jié)構(gòu)在氣動載荷作用下會產(chǎn)生更大的變形，氣動與結(jié)構(gòu)間的耦合效應(yīng)顯著增強[2-5]。

由于飛行速度快，高超聲速導(dǎo)彈面臨更加嚴峻的熱環(huán)境，結(jié)構(gòu)剛度的改變造成結(jié)構(gòu)的承載能力降低，飛行器內(nèi)部溫度的增高致使元器件等部件損害率加大，這也是目前高超聲速飛行器面臨的共性問題。面對熱防護設(shè)計和結(jié)構(gòu)設(shè)計帶來的挑戰(zhàn)，眾多試驗項目以此為背景相繼進行：北京航天長征飛行器研究所采用主動式熱疏導(dǎo)技術(shù)對金屬尖化前緣模型進行試驗研究[6];吳大方等[7]對某飛行器翼面進行了熱/振聯(lián)合試驗;成都飛機設(shè)計研究所對飛行器熱環(huán)境下的結(jié)構(gòu)固有振動特性進行了試驗[8];中國空氣動力研究所與發(fā)展中心分別進行了MF-1航天模型飛行試驗以及激波風(fēng)洞試驗[9-10]。但是試驗要求技術(shù)難度大，經(jīng)費高，采用數(shù)值仿真對飛行器氣動力、結(jié)構(gòu)溫度等進行有效的預(yù)測，可為高超聲速導(dǎo)彈熱防護設(shè)計與結(jié)構(gòu)設(shè)計提供指導(dǎo)。

由于耦合效應(yīng)不可忽視，對高超聲速導(dǎo)彈流場與結(jié)構(gòu)溫度仿真研究需采用多場耦合方法。目前關(guān)于耦合數(shù)值計算方法，可以分為兩大類：一體化求解方法和分區(qū)求解方法。一體化求解方法也被稱為直接耦合求解方法，是將不同子系統(tǒng)的控制方程放在一個系統(tǒng)方程組中求解[11]。黃唐等[12]進行了二維流固熱一體化數(shù)值模擬，以二維圓柱繞流問題為例，流場和熱響應(yīng)部分分別采用總變差消格式和迦遼金方法進行方程離散。季衛(wèi)棟等[13]將流場與結(jié)構(gòu)控制方程寫成統(tǒng)一形式，對整個物理場進行統(tǒng)一的迎風(fēng)格式有限體積方法離散，實現(xiàn)了二維圓管模型的一體化數(shù)值模擬。雖然一體化求解中流體與結(jié)構(gòu)的解變量完全耦合，求解精度更高，但由于流體和結(jié)構(gòu)的控制方程具有完全不同的物理特性和數(shù)學(xué)性質(zhì)，目前僅在簡單的二維問題展開研究。分區(qū)求解方法是將系統(tǒng)分為流體和固體兩個獨立的模塊[14]。劉健等[15]、周印佳等[16]與黃杰等[17]均采用分區(qū)耦合方法，對高超聲速二維圓管進行數(shù)值模擬，計算結(jié)果與試驗值相吻合。夏剛等[18]將新型上風(fēng)格式與結(jié)構(gòu)傳熱的Galerkin方法相結(jié)合，實現(xiàn)鈍體高超聲速飛行器的分區(qū)耦合計算。張勝濤等[19]采用自適應(yīng)耦合時間步長和混合差值策略實現(xiàn)了高超聲速飛行器進氣道前緣結(jié)構(gòu)的分區(qū)耦合仿真分析。肖軍等[20]通過在子代步上氣動與結(jié)構(gòu)方程交替求解實現(xiàn)了一種全隱式分區(qū)耦合算法。李凱倫等[21]基于有限元數(shù)值分析方法建立了一種多場分區(qū)耦合模型，對高超聲速環(huán)境中功能梯度薄板熱氣動彈性進行了研究。梁杰等[22]基于一種分區(qū)求解策略，采用直接模擬蒙特卡洛方法，對某飛行器氣動力、氣動熱進行了數(shù)值模擬。

目前關(guān)于高超聲速多場耦合的研究主要在耦合方法的實現(xiàn)上，但是順利完成高超聲速導(dǎo)彈的設(shè)計，不僅需要有效的多場耦合仿真技術(shù)去實現(xiàn)氣動熱力學(xué)環(huán)境預(yù)示，還需對各場之間耦合關(guān)系定量化分析，否則會造成設(shè)計不當(dāng)而引發(fā)導(dǎo)彈高速飛行時結(jié)構(gòu)破壞等不利后果。通過研究高超聲速導(dǎo)彈溫度場、壓力場與結(jié)構(gòu)變形之間的相互關(guān)系，分析耦合效應(yīng)對其影響，對導(dǎo)彈熱防護以及結(jié)構(gòu)設(shè)計具有指導(dǎo)意義。

本文基于分區(qū)耦合方法，以大長細比導(dǎo)彈為研究對象，實現(xiàn)其流/固/熱數(shù)值模擬。對導(dǎo)彈不同攻角下的結(jié)構(gòu)變形和溫度場、壓力場進行計算分析，并對導(dǎo)彈氣動力、氣動熱與結(jié)構(gòu)變形間的耦合效應(yīng)進行定量化研究。

1 耦合分析策略和控制方程

1.1 耦合方法

耦合模型如圖1所示。在流體介質(zhì)內(nèi)部，求解統(tǒng)一的流體控制方程得到熱流qf和壓力pf，并通過流-固耦合界面向固體介質(zhì)提供熱載荷qf和力載荷pf。在固體介質(zhì)內(nèi)部，求解熱傳導(dǎo)控制方程和熱彈性力學(xué)控制方程得到結(jié)構(gòu)壁面溫度Ts和位移us，并通過流-固耦合界面向流體介質(zhì)提供溫度條件Ts和結(jié)構(gòu)變形條件us。

基于多場耦合模型，采用分區(qū)求解方法建立如圖2所示的緊耦合分析策略，流體與固體區(qū)域求解器均為瞬態(tài)求解，求解器所需數(shù)據(jù)通過耦合界面數(shù)據(jù)不斷交換提供，耦合分析策略實施流程描述如下：

1)首先根據(jù)結(jié)構(gòu)的初始值，即結(jié)構(gòu)壁面初始溫度Ts與位移us，作為瞬態(tài)耦合分析的初始條件，通過耦合界面?zhèn)鬟f給流體域。

2)從流體域內(nèi)進行瞬態(tài)流場計算，得到熱流qf和壓力(流體應(yīng)力)pf通過耦合界面?zhèn)鬟f給固體域。

3)固體域內(nèi)進行瞬態(tài)求解得到結(jié)構(gòu)壁面溫度Ts和位移us通過耦合界面?zhèn)鬟f給流體域。

4)檢查流體應(yīng)力和結(jié)構(gòu)位移是否滿足收斂標(biāo)準(zhǔn)，如果不滿足，則返回第2步繼續(xù)進行迭代求解，直到滿足收斂標(biāo)準(zhǔn)或者達到設(shè)置的迭代次數(shù)；如果滿足收斂標(biāo)準(zhǔn)，則t0時間步的計算結(jié)束，輸出流體和結(jié)構(gòu)的結(jié)果，進行t1時間步的耦合計算。

5)該過程不斷重復(fù)直到計算時間結(jié)束。

1.2 流體控制方程

流體控制方程采用積分形式的N-S方程，其形式如下[23]：

(1)

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

1.3 熱-結(jié)構(gòu)控制方程

結(jié)構(gòu)傳熱方程計算采用Fourier定律，不考慮內(nèi)熱源，并假設(shè)無輻射，基于Fourier定律的三維瞬態(tài)熱傳導(dǎo)控制微分方程為Fourier定律的三維瞬態(tài)熱傳導(dǎo)控制微分方程[23]：

(2)

式中：ρ為密度，c為比熱容，T為溫度，t為時間，λ為導(dǎo)熱系數(shù)。

結(jié)構(gòu)方程為熱彈性力學(xué)控制方程，采用剛度矩陣求解法，假設(shè)結(jié)構(gòu)產(chǎn)生的靜變形是一個緩慢的過程，不考慮慣性力，其方程為[23]：

(3)

K是熱結(jié)構(gòu)剛度矩陣，熱-結(jié)構(gòu)控制方程通過K實現(xiàn)耦合，us是結(jié)構(gòu)位移向量，F(xiàn)s是作用在結(jié)構(gòu)上的力矩陣。溫度的改變會造成材料特性改變，變化后的結(jié)構(gòu)剛度矩陣為KT；溫度還會使結(jié)構(gòu)內(nèi)部產(chǎn)生熱應(yīng)力，結(jié)構(gòu)間的相互約束使其受熱后膨脹受到限制引發(fā)預(yù)拉壓應(yīng)力，這種熱應(yīng)力產(chǎn)生的附加剛度為KΔT。

2 仿真校驗

經(jīng)典圓管繞流試驗作為算例[24]。試驗所用圓管為不銹鋼，圓管內(nèi)、外半徑分別為25.4 mm和 38.1 mm，熱力學(xué)參數(shù)：導(dǎo)熱系數(shù)16.27 W/(m·℃)，比熱容502.48 J/(kg·℃)，密度8030 kg/m3，彈性模量1.2×1013Pa，熱膨脹系數(shù)1.68×10-5(1/℃)，泊松比0.3。圓管內(nèi)壁設(shè)為等溫壁，溫度值與初始環(huán)境溫度一致，為21.4 ℃。來流參數(shù)：溫度-31.5 ℃，壓力648 Pa，速度Ma6.47。圖3為二維計算網(wǎng)格，流場與圓管網(wǎng)格交界面為耦合界面，流固通過耦合界面進行數(shù)據(jù)傳輸實現(xiàn)多場耦合。圖4和圖5分別為計算初始時圓管表面壓力和熱流分布與文獻[24]的試驗結(jié)果對比，圖中p0，q0分別為駐點壓強與熱流，θ為物面到圓心的連線與X軸的夾角。可以看出壓p強和熱流q分布與文獻[24]符合得很好。圖6為駐點2 s時刻的溫度值與文獻的對比，對比結(jié)果較好。由此，校驗了此耦合方法的有效性。

3 計算結(jié)果與分析

3.1 計算模型

導(dǎo)彈幾何模型如圖7所示，僅考慮彈體結(jié)構(gòu)，簡化彈翼、舵面和內(nèi)部元器件等。導(dǎo)彈結(jié)構(gòu)材料采用鈦合金，材料參數(shù)如表1所示。在彈身前段迎風(fēng)和背風(fēng)位置分別布置監(jiān)測點p1，p2。導(dǎo)彈頭部頂端為坐標(biāo)原點。

計算模型與網(wǎng)格示意圖分別如圖8、圖9所示。仿真計算模型流場區(qū)域分為外流場和內(nèi)流場，外流場為導(dǎo)彈外部氣流，內(nèi)流場為導(dǎo)彈結(jié)構(gòu)內(nèi)部氣體。外流場對稱面采用對稱邊界條件，其它面采用外部流邊界條件，在此邊界條件上設(shè)置壓力、速度和溫度載荷，載荷值分別為5529 Pa，Ma5，-56 ℃。固體結(jié)構(gòu)場為彈體結(jié)構(gòu)，在彈體尾部采用固定邊界條件。分別在外流場與彈體外壁面交界處，內(nèi)流場與彈體內(nèi)壁面交界設(shè)置外場與內(nèi)場耦合邊界條件，通過耦合邊界實現(xiàn)分區(qū)耦合的數(shù)據(jù)傳輸。計算攻角分別為0°，5°，10°，15°，20°，仿真時間40 s，初始內(nèi)、外流場和結(jié)構(gòu)溫度假設(shè)為20 ℃，壓力值5529 Pa。

彈性模量/GPa熱膨脹系數(shù)/(℃-1·10-6)熱傳導(dǎo)系數(shù)/(W·m-1·℃-1)比熱容/(J·kg-1·℃-1)75～1099.20～11.066.8～11.8610～702

為研究導(dǎo)彈氣動力、氣動熱與結(jié)構(gòu)變形之間的相互關(guān)系，將計算工況分為耦合計算(考慮耦合效應(yīng)影響)、僅考慮熱因素(結(jié)構(gòu)僅有熱變形)的計算、僅考慮氣動力因素(結(jié)構(gòu)僅有氣動變形)的計算和無結(jié)構(gòu)變形因素計算四種情況。

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

導(dǎo)彈結(jié)構(gòu)在X，Y，Z三個方向發(fā)生變形，以攻角20°的位移云圖為例，如圖10所示。由圖10可知，彈體結(jié)構(gòu)發(fā)生拉伸、膨脹和彎曲變形。導(dǎo)彈為軸對稱結(jié)構(gòu)且仿真計算未考慮側(cè)偏角，彈體結(jié)構(gòu)在Y方向相對軸線僅發(fā)生膨脹而無彎曲變形，所以主要研究彈體結(jié)構(gòu)軸向(X方向)和橫向(Z方向)變形。

彈體結(jié)構(gòu)沿導(dǎo)彈軸向的變形量為Δx=(Δx1+Δx2)/(2L)×100%，其中L為導(dǎo)彈全長，Δx1為結(jié)構(gòu)迎風(fēng)側(cè)沿X方向的位移改變量，Δx2為結(jié)構(gòu)背風(fēng)側(cè)沿X方向的位移改變量。彈體結(jié)構(gòu)軸向變形曲線如圖11所示，橫坐標(biāo)為彈體結(jié)構(gòu)所處位置，x=0為導(dǎo)彈頭部頂端位置。

由圖11可知，導(dǎo)彈高速飛行時產(chǎn)生氣動加熱，使結(jié)構(gòu)產(chǎn)生軸向的拉伸，頭部氣動加熱最劇烈，拉伸量大。隨著攻角增加，軸向變形量減小，因結(jié)構(gòu)沿軸向拉伸時，撓度隨攻角的增加而增加，相應(yīng)沿軸向變形量減小。

彈體結(jié)構(gòu)沿導(dǎo)彈橫向的變形量為Δz=(Δz1+Δz2)/(2R)×100%，其中R為導(dǎo)彈半徑，Δz1為結(jié)構(gòu)迎風(fēng)側(cè)沿Z方向的位移改變量，Δz2為結(jié)構(gòu)背風(fēng)側(cè)沿Z方向的位移改變量。其結(jié)構(gòu)橫向變形曲線如圖12所示。

由圖12可知，0°攻角時溫度和氣動力均布，彈體結(jié)構(gòu)未發(fā)生軸向變形，隨著攻角增大，結(jié)構(gòu)發(fā)生彎曲變形，且頭部變形量大。氣動力/熱共同作用使結(jié)構(gòu)產(chǎn)生彎曲變形：彈體末端采用固定邊界條件，類似于懸臂梁模型，梁承受Z方向的氣動力使結(jié)構(gòu)發(fā)生橫向彎曲；同時氣動加熱作用使結(jié)構(gòu)發(fā)生軸向拉伸，結(jié)構(gòu)迎風(fēng)側(cè)溫度高，拉伸量大，結(jié)構(gòu)背風(fēng)側(cè)溫度低，拉伸量小，造成彈體結(jié)構(gòu)沿背風(fēng)側(cè)位置發(fā)生彎曲；氣動力和氣動加熱在頭部作用明顯，彎曲變形量大。

研究氣動加熱與氣動力對導(dǎo)彈結(jié)構(gòu)變形的影響，以彈身前段位置為例，變形量由p1和p2測點數(shù)據(jù)計算而得，其中對結(jié)構(gòu)軸向變形的影響見圖13，對結(jié)構(gòu)橫向變形的影響見圖14。

由圖13可知，彈體結(jié)構(gòu)的軸向變形主要由氣動加熱引起；氣動力使彈體結(jié)構(gòu)產(chǎn)生軸向壓縮但影響很小，可忽略不計。由圖14可知，彈體結(jié)構(gòu)在橫向的彎曲變形由氣動熱和氣動力共同引起，且導(dǎo)彈采用大長細比結(jié)構(gòu)時，結(jié)構(gòu)柔度增加，氣動力造成結(jié)構(gòu)軸向的彎曲變形影響更大。

3.3 溫度分析

導(dǎo)彈溫度云圖如圖15所示。在彈體軸向位置，來流在頭部位置發(fā)生摩擦和壓縮最劇烈，動能轉(zhuǎn)化為熱能最多，氣動加熱劇烈，且隨著攻角增大，溫度最大值點從頭部最前端往頭部迎風(fēng)側(cè)后側(cè)偏移。攻角為0°時溫度均布，隨著攻角增大，迎風(fēng)側(cè)與背風(fēng)側(cè)位置溫度差也隨之增大。

研究耦合效應(yīng)對溫度仿真結(jié)果影響，以彈身前段迎風(fēng)側(cè)p1監(jiān)測點為例，如圖16所示。當(dāng)攻角較小時，耦合效應(yīng)對對溫度仿真結(jié)果影響不大；隨著攻角增大，耦合效應(yīng)對仿真結(jié)果的影響增大，因為導(dǎo)彈大攻角飛行時，彈體結(jié)構(gòu)因氣動熱和氣動力發(fā)生彎曲變形，攻角由α增加為α+Δα，且此Δα較大造成彈體結(jié)構(gòu)的溫度分布發(fā)生改變，加劇了彈體結(jié)構(gòu)的氣動加熱。

3.4 壓力分析

導(dǎo)彈壓力場云圖如圖17所示，分布規(guī)律與溫度場分布規(guī)律相似，頭部氣動力環(huán)境嚴酷。

研究耦合效應(yīng)對壓力仿真結(jié)果影響，以彈身前段迎風(fēng)側(cè)p1監(jiān)測點為例，如圖18所示。由圖18可知，當(dāng)攻角較小時，耦合效應(yīng)對壓力仿真結(jié)果影響不大；隨著攻角增大，耦合效應(yīng)對仿真結(jié)果的影響增大，彈體結(jié)構(gòu)的彎曲變形造成攻角α改變，增加一個不可忽略的增量Δα，使導(dǎo)彈的壓力分布產(chǎn)生改變，氣動力環(huán)境更加惡劣。

4 結(jié) 論

1)大長細比導(dǎo)彈在高速飛行時，隨著攻角增大，結(jié)構(gòu)發(fā)生軸向拉伸與橫向彎曲，其中軸向拉伸主要由氣動加熱引起，橫向彎曲主要由氣動力引起。

2)導(dǎo)彈在頭部溫度和壓力載荷大，且隨著攻角的增大，迎風(fēng)位置溫度與壓力載荷增大。

3) 對于大長細比導(dǎo)彈，以小攻角飛行時耦合效應(yīng)對溫度和壓力仿真結(jié)果影響不大；以大攻角飛行時耦合效應(yīng)明顯，造成導(dǎo)彈氣動熱力學(xué)環(huán)境嚴酷，導(dǎo)彈設(shè)計進行熱力學(xué)環(huán)境預(yù)示時需要考慮此效應(yīng)。