周述光，國義軍，賀立新，劉驍

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

再入彈頭三維非對稱燒蝕外形模擬

周述光1,*，國義軍2，賀立新2，劉驍2

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

遠程戰(zhàn)略導(dǎo)彈再入大氣層，熱環(huán)境十分嚴酷，燒蝕量很大，外形呈現(xiàn)非對稱性。燒蝕變形、質(zhì)量引射和表面燒蝕花紋反作用于氣動力、氣動熱及表面溫度場。它們之間的相互作用十分復(fù)雜，影響落點精度。依據(jù)等效轉(zhuǎn)換，將有迎角、側(cè)滑角的繞流轉(zhuǎn)換成總迎角繞流。對于不同子午面，將其變換到縱向坐標軸沿等效流動方向的氣流坐標系，通過等效迎角化為零迎角繞流。采用Monte-Carlo統(tǒng)計方法，生成偽隨機數(shù)，對彈頭表面粗糙度隨機抽樣，表征不同子午線上轉(zhuǎn)捩初始位置的隨機分布。改進了零迎角燒蝕外形計算方法，拓展了其應(yīng)用范圍，并通過與能反映轉(zhuǎn)捩點移動方向的松耦合計算策略相結(jié)合，采用數(shù)值求解燒蝕外形方程，成功模擬了小傾角再入過程中端頭非對稱燒蝕外形。

再入彈頭；非對稱體；熱防護；氣動熱；燒蝕；碳/碳材料；松耦合模擬

傳統(tǒng)遠程導(dǎo)彈彈頭以零迎角慣性再入，再入時間短，熱流高，燒蝕嚴重，燒蝕外形基本上呈現(xiàn)對稱性。近年來，新型遠程戰(zhàn)略導(dǎo)彈為提高機動性，采用持續(xù)有迎角再入彈道，再入時間長，盡管峰值熱流比傳統(tǒng)彈頭小很多，但總熱量更大，燒蝕外形為非對稱。燒蝕外形/氣動力/氣動熱/溫度場/燒蝕等緊密耦合，相互影響，不確定因素多，預(yù)測精度差。為了提高飛行安全性能，通常加大防熱系統(tǒng)的安全裕度，使得防熱結(jié)構(gòu)重量增加。為此，諸如碳/碳等燒蝕型防熱材料的燒蝕外形變形規(guī)律得到了大量研究。利用電弧加熱設(shè)備、高頻等離子風洞、燃氣流設(shè)備等，結(jié)合對飛行回收物的分析和理論研究，對彈頭燒蝕外形有了基本認識。1970年，Thyson等預(yù)測了層流和湍流中理論外形的變化[1]，并比較了地面試驗結(jié)果，意識到粗糙度的重要性。隨后，Baker也預(yù)測了零入射角和零迎角下的燒蝕[2]，利用地面試驗設(shè)備研究了有迎角下樟腦制成的球錐試樣燒蝕特性，并推斷是物體表面粗糙度引起了湍流熱流的增加。轉(zhuǎn)捩依據(jù)粗糙度的高低先后發(fā)生，熱流也隨之發(fā)生變化。不同流態(tài)下的燒蝕速率相差甚大，燒蝕外形與邊界層的流動狀態(tài)密切相關(guān)。例如，球錐外形的三向編織碳/碳端頭，依據(jù)彈頭再入飛行過程中邊界層流動狀態(tài)的變化，一般會出現(xiàn)三類6種燒蝕外形[2]：層流燒蝕外形(A、B)、轉(zhuǎn)捩燒蝕外形(C、D)和湍流燒蝕外形(E、F)，如圖1所示[3],RD、RE、RF分別為外形D、E、F的極徑。

在再入彈頭熱防護計算研究的早期階段，忽視端頭外形的小不對稱性，例如，20世紀70年代，Crowell研究了在彈頭零迎角再入過程中，碳/碳材料的端頭燒蝕外形變化規(guī)律[4]。Chin研究了端頭的非對稱燒蝕快速計算方法，發(fā)展了適用于軸對稱外形的三維燒蝕計算程序[5]，方程基于球面、可變形、移動坐標系，計算表明石墨端頭在10°迎角下燒蝕外形呈現(xiàn)顯著的非對稱性。之后，Dirling[6]根據(jù)端頭回收體(Nosetip Recovery Vehicle, NRV，見圖2)項目的飛行試驗實物[7]，針對彈頭非對稱外形，采用材料表面粗糙度的隨機分布來解釋周向邊界層轉(zhuǎn)捩不一致性，并引入了Monte-Carlo分析法。Hall和Nowlan對控制方程采取有限差分-積分法和統(tǒng)計散布分析的近似技術(shù)[8]，分析了小不對稱端頭的氣動特性。

圖1 彈頭典型燒蝕外形
Fig.1 Classical ablation shapes of missile nosetip

圖2 整體石墨NRV再入飛行返回端頭燒蝕外形
Fig.2Ablative shape of nosetip of recovered bulk graphite NRV flight test

國內(nèi)從20世紀70年代也開始了再入彈頭燒蝕防熱的工程計算方法研究，在莊逢甘院士、張涵信院士等推動下，姜貴慶等開展了燒蝕防熱機理研究及其工程計算研究[9-11]，先后建立了零迎角燒蝕計算方法、有迎角解耦計算方法等，并基于炮風洞開展了粗糙度對轉(zhuǎn)捩、熱交換等影響的試驗研究。20世紀80年代以來，對非對稱燒蝕問題進行了持續(xù)研究，其中又以在風洞模擬環(huán)境下測量氣動力、熱、燒蝕的研究為主[12-15]。常用的燒蝕防熱快速算法大部分為對稱燒蝕[16-17]，在非對稱燒蝕外形方面的研究并不多見，袁湘江等[18]根據(jù)Dirling[19]介紹的方法，采取解耦計算方法，大致復(fù)現(xiàn)了零迎角狀態(tài)NRV的小不對稱燒蝕整體外形，模擬的溝槽不夠明顯，這可能與計算方法和姿態(tài)角選取有一定關(guān)系。由于端頭外形對升阻比和穩(wěn)定性影響很大，時至今日仍有報道，最近，陳自發(fā)等針對有迎角或滾轉(zhuǎn)角的情況，計算得出碳基端頭迎風面燒蝕量大于背風面的結(jié)論[20]。

在所有燒蝕外形計算方法中，解耦計算方法是最常見、最簡單的，即首先采用初始外形沿彈道計算所有時刻的熱環(huán)境，然后根據(jù)熱環(huán)境再計算燒蝕外形，不考慮燒蝕外形變化對熱環(huán)境的影響，可簡單表述為“初始外形→熱環(huán)境→燒蝕外形”。燒蝕量小的載人飛船返回艙、戰(zhàn)術(shù)導(dǎo)彈等常采取解耦計算[4-5,21]。而燒蝕防熱型的戰(zhàn)略彈頭，由于燒蝕量大，外形的重要性凸現(xiàn)出來，不能再采取解耦計算。燒蝕對氣動力影響很大，體現(xiàn)在3個方面：燒蝕外形變化，會明顯改變飛行器的氣動性能；燒蝕過程中的質(zhì)量引射對邊界層狀態(tài)的影響也會反映到氣動性能上；燒蝕形成的小不對稱與質(zhì)量、慣性的分布不對稱性，對運動產(chǎn)生較大影響，反過來影響氣動性能。文獻[10]指出氣動力、氣動熱、彈道運動的耦合計算內(nèi)容包括：飛行器運動特性、外形計算、氣動特性、燒蝕計算。彈頭再入過程中，如果考慮它們之間的作用，即為所謂(緊)耦合計算。由于涉及面廣，難度很大，目前還未見到這方面的研究報道。

總的來說，國內(nèi)對帶迎角、側(cè)滑角的飛行狀態(tài)下的燒蝕模擬技術(shù)欠缺；如何合理地模擬彈頭表面粗糙度影響值得研究；為了提高模擬精準度，還需要考慮每一時刻的端頭外形對氣動熱環(huán)境的作用。

為此，本文燒蝕外形計算中，不考慮燒蝕外形變化時氣動特性的改變對運動特性(彈道參數(shù))的影響，沿著彈道，根據(jù)上一時刻的燒蝕外形，計算當前時刻熱環(huán)境參數(shù)(熱流、恢復(fù)焓、溫度、邊界層外緣速度、壓力等)，然后根據(jù)燒蝕方法，刻畫出這一時刻的外形，作為下一時刻熱環(huán)境的輸入依據(jù)，即所謂燒蝕外形松耦合計算，可簡單表示為“熱環(huán)境?燒蝕外形”。本文在零迎角、松耦合計算方法驗證算例的基礎(chǔ)上[22]，研究了由于粗糙度沿周向非均勻分布對燒蝕外形產(chǎn)生的影響，建立了帶迎角、側(cè)滑角的松耦合計算方法，實現(xiàn)了非對稱燒蝕外形預(yù)測。

1 燒蝕外形模擬技術(shù)

1.1 坐標變換

氣動力和氣動熱參數(shù)計算均在坐標原點位于彈頭頭部的坐標系中進行，而燒蝕導(dǎo)致頭部物面向后退縮，坐標原點隨之后移，每一時刻需要重新確定物面坐標值。此外，頭部頂點還徑向偏離原始物體頂點。對于此兩類情況，一個較好的辦法是采用瞬時坐標系，通過坐標變換，建立新舊坐標系之間的聯(lián)系(圖3)。設(shè)Oxyz為初始時刻的固連直角坐標系，坐標原點O位于彈頭頂點。定義O*x*y*z*為當前計算時刻的瞬時直角坐標系，坐標原點O*始終位于彈頭燒蝕外形的前部頂點處，各坐標軸的方向與初始坐標系一致。子午面分布不變。瞬時坐標系中的某一物面點P，在原始坐標系中為(xP，yP，zP)，子午角為φ0。瞬時坐標系O*x*y*z*與初始坐標系Oxyz之間的關(guān)系為

(1)

式中：δx、δy、δz為兩個坐標系的坐標原點之間的距離Δ分別在x、y、z軸上的投影。

圖3 初始坐標和瞬時坐標示意圖
Fig.3Schematic of original and instantaneous coordinates

旋成體非常適合采用球坐標系(R,θ,φ)來描述，設(shè)初始時刻的球坐標原點N′位于x軸上，某一時刻移原點至N，R為原點到物面的極徑；θ和φ分別為球心角和子午角。則初始直角坐標系與球坐標系之間的關(guān)系為

(2)

由燒蝕計算可確定各子午面新的R值，并由式(2)計算物面上各點的直角坐標x、y、z，然后尋找整個物面上的最小x坐標點，從而確定δx、δy、δz。最后由式(1)確定新坐標(x*，y*，z*)，開始下一時刻的流場參數(shù)計算。

1.2 燒蝕外形數(shù)值求解

燒蝕外形嚴格說來是一個三維問題，但對于實際飛行的彈頭而言，燒蝕外形沿周向的變化率遠遠小于縱向變化率。因此，通常當作二維問題處理，僅在每個子午面內(nèi)解二維燒蝕外形方程。球坐標系下，描述燒蝕外形變化歷程的控制方程為[4，23]

(3)

(4)

記式(4)右邊為函數(shù)f，它是Rθ、R、θ、φ、來流馬赫數(shù)Ma∞的函數(shù)，即

則式(4)的空間半離散方程為

對應(yīng)的修正方程為

(Rt)j=f(Rθj+ε,Rj,θj,Ma∞,φ)

(5)

式中：ε為一小量。式(5)展開得

(6)

2 姿態(tài)角等效處理

(7)

由此，可將有迎角和側(cè)滑角的繞流問題轉(zhuǎn)化為單純總迎角繞流問題。

對所研究的子午面，將其變換到縱向坐標軸沿等效流動方向(與x坐標軸成迎角αEB)的氣流

圖4 等價變化坐標示意圖
Fig.4Schematic of equivalence transformation of coordinate

坐標系軸線，可以進一步化為零迎角繞流問題。運用等效體方法，可得到相對于等效零迎角旋成體繞流的物面等價傾角，即

(8)

式中：n為物面外法線方向的單位矢量。在球坐標系中，定義物面方程F為

Fr,φ,x=r-fx,φ=0

(9)

式中：r為物面上的點到軸線(x軸)的距離。則物面外法線單位向量為

(10)

式中：ex、er、eφ為球坐標系中的3個單位向量。由式(11)定義物面角(切向和周向)，即

(11)

式中：-π/2≤θb≤π/2；-π/2≤δφ≤π/2。由式(11)可導(dǎo)出

(12)

將式(11)和式(12)代入式(10)，可得

n=-exsinθb+ercosθbcosδφ-eφcosθbsinδφ

可記為

n=-isinθb+jcosθbcosφ-δφ+

kcosθbsinφ-δφ

(13)

式中：i、j、k為空間直角坐標系中的3個單位向量。

速度矢量在瞬時坐標系中的表達式為

V∞=V∞icosβcosα+jcosβsinα+ksinβ

(14)

將式(13)、式(14)代入式(8)，可得

θbEQ=arcsin(cosαcosβsinθb-sinαcosβcosθb·

cosφ-δφ-sinβcosθbsinφ-δφ)

(15)

因為φ=0°為背風面，φ=180°為迎風面，所以式(15)中的α必須反號，即-α。定義等價迎角為

αEB=θb-θbEQ

(16)

至此，建立了任意子午面的等效迎角與實際飛行迎角和側(cè)滑角的關(guān)系。通過以上變換，將有迎角和側(cè)滑角的繞流問題化為等價零迎角繞流問題，這樣就可以運用以往成熟的零迎角情況下的計算方法，求得表面壓力分布、激波形狀以及表面熱流分布等流場參數(shù)。

3 熱環(huán)境和燒蝕計算方法

3.1 邊界層外緣參數(shù)計算

球錐外形彈頭的外部流場考慮端頭燒蝕形狀的變化，表面壓力分布按照端頭部位和錐身部位分別計算。為了計算激波層中熵梯度效應(yīng)對熱流的影響，首先計算每個時刻燒蝕外形的頭部激波形狀，獲得實際激波角和激波脫體距離，再計算激波后參數(shù)，根據(jù)這些參數(shù)對邊界層外緣焓進行修正，利用均值法，獲得邊界層外緣平均速度。

3.2 表面熱流計算及粗糙壁熱增量修正

對于不可穿透物形表面的層流、湍流熱流，采用有效長度法計算；考慮熵梯度效應(yīng)，應(yīng)用平均值法；粗糙壁熱增量即燒蝕表面粗糙度對熱流的影響，采用粗糙壁熱流放大因子進行修正，即

qj=Fjqjh0,u0,ρ0,μ0j=l,t

(17)

式中：Fj為粗糙壁熱增因子；q為熱流；h為焓值；u為邊界層外緣速度；ρ為密度；μ為黏性系數(shù)。上標“0”表示外部主流區(qū)；下標l、t表示層流和湍流狀態(tài)。

轉(zhuǎn)捩區(qū)熱流用間歇因子Γ加權(quán)平均，即

qtr=(1-Γ)ql+Γqt

(18)

轉(zhuǎn)捩區(qū)間歇因子Γ的值在0(層流)～1(湍流)之間變化，在不同雷諾數(shù)范圍內(nèi)有不同的表達式，即

(19)

式中：Re1為轉(zhuǎn)捩開始雷諾數(shù)；Re2為平均雷諾數(shù)；Re3為結(jié)束雷諾數(shù)；Reθ為當?shù)乩字Z數(shù)。

3.3 燒蝕速率計算方法

端頭的物面法向線燒蝕速率由材料的有效焓計算，將物面視為準穩(wěn)定的加熱面，對于碳/碳復(fù)合材料，有

(20)

4 粗糙度隨機分布模擬

彈體表面粗糙度在燒蝕外形計算中扮演著重要的角色，它不僅影響邊界層轉(zhuǎn)捩，也影響表面熱流，決定了燒蝕速率的快慢，使得端頭外形出現(xiàn)深淺不一的溝槽，逐漸呈現(xiàn)小不對稱性。通常情況下，表面粗糙度分布是隨機的，另外還有一些小不對稱因素也是隨機分布的，因此，燒蝕外形也具有隨機的小不對稱性。

kj=kL+σKXj

(21)

式中：kL為測量試樣得到的名義粗糙度高度；σK為粗糙度測量值的名義標準偏差；Xj為偽隨機數(shù)列，由標準的偽隨機數(shù)列程序產(chǎn)生，周期足夠長，以保障其隨機性。在任意的一組隨機抽樣值中，選出其最大值Kπ1，最小值Kπ2，作為對應(yīng)子午線φ=φ*和子午線φ=π+φ*的粗糙度值，通過這兩個值，確定中間子午線上的Kπ值，即

(22)

式中：φ*為最大粗糙度值所在的子午角。

5 計算結(jié)果分析

5.1 結(jié)果驗證

彈頭外形參考NRV的設(shè)計[6]，為球-錐旋成體構(gòu)型，端頭半徑為31.75 mm，底部半徑為151.13 mm，總長為1 168.4 mm，半錐角為6°。根據(jù)Dirling的建議[6]，粗糙度因子取為15.5 μm，均方差為8.1 μm。綜合文獻給出的部分彈道參數(shù)[6-7]，見表1，并將彈頭姿態(tài)角沿彈道全程設(shè)為0。驗證計算結(jié)果見圖5(縱橫比非1∶1，下面所有的外形、輪廓圖均如此)，給出了端頭外形。可知，模擬出端頭上多道深淺不一的燒蝕溝槽，與回收物燒蝕外形類似。端頭駐點后退量約為7.6 mm，落在Otey和English預(yù)測的7.4～8.9 mm范圍內(nèi)[6]，而對飛行試驗回收的端頭進行測量，獲得的后退量約為8.5 mm。考慮到彈道參數(shù)不全，粗糙度影響等原因，可以認為本次驗證計算結(jié)果合理，程序可靠。

表1 彈道參數(shù)Table 1 Parameters of trajectory

圖5 驗證計算端頭燒蝕外形
Fig.5Verification simulation of ablation shape of nosetip

5.2 姿態(tài)角與端頭燒蝕外形非對稱性

不同迎角和側(cè)滑角會引起迎風面和背風面的熱流不同，燒蝕速率也就出現(xiàn)差異，導(dǎo)致燒蝕外形的不對稱性。這里將粗糙度方差設(shè)置為0，即消除粗糙度不均帶來的影響。圖6給出了同一飛行條件下有/無姿態(tài)角的端頭燒蝕外形，dH表示當前時刻外形與原始外形之間變化量；圖7給出了其中3個狀態(tài)下縱向(z=0)剖面輪廓圖，橫坐標為軸向坐標x，縱坐標為剖面上輪廓的周向直徑r。端頭頂點后退量均約為7.4～7.5 mm，輪廓線的上半部分為背風面，下半部分為迎風面，不對稱性隨總迎角的增大而增強。圖8給出了Λ=15°的姿態(tài)角彈道、不同高度下的端頭后退量?？芍?，有姿態(tài)角時，迎風面燒蝕量大于背風面，呈現(xiàn)非對稱燒蝕外形；此外，在30 km以上，燒蝕量較小，25 km以下，燒蝕量增大。當迎角和側(cè)滑角相同時，來流方向在子午角45°平面里，如圖9所示(α=β=10.6°，Λ=15.0°)，這意味φ=45°的子午面輪廓不對稱性最顯著，φ=135°的子午面輪廓為對稱性，縱向?qū)ΨQ剖面(φ=90°)和水平剖面(φ=0°)的輪廓一致，均為程度較輕的小對稱性。

圖6 端頭燒蝕小不對稱外形
Fig.6 Small asymmetric ablation shape of nosetip

圖7 端頭燒蝕小不對稱縱向剖面輪廓圖
Fig.7Outlines of small asymmetric ablation shape of nosetip at z=0 section

5.3 粗糙度與端頭燒蝕外形非對稱性

在粗糙度研究中，模型姿態(tài)角取為0°。由于熱流會明顯影響燒蝕外形，為了便于觀察，這里選用一條環(huán)境嚴酷的彈道進行模擬。由圖10可知，燒蝕的不對稱性是隨機分布的，其中RL為層流粗糙度因子，RT為湍流粗糙度因子，MVL為層流粗糙度因子分布均方差，MVT為湍流粗糙度因子分布均方差。無論是層流粗糙度還是湍流粗糙度，隨著粗糙度的降低，燒蝕外形的不對稱性程度減輕。與湍流狀態(tài)相比較，層流狀態(tài)下的粗糙度因子對熱流作用更顯著。這可能源自層流流動對擾動更敏感，更容易向轉(zhuǎn)捩發(fā)展，熱流變化相對大些，燒蝕速率差異較大；湍流抗干擾能力較強，流動穩(wěn)定性更好，熱流變化也相對緩和，燒蝕速率較為均衡。

5.4 耦合-解耦計算結(jié)果分析

在耦合-解耦計算策略影響研究中，模型姿態(tài)角取為0°，粗糙度均方差設(shè)置為0，燒蝕外形將具有對稱性。圖11給出了沿同一條彈道計算， 耦合計算與非耦合計算的結(jié)果。由圖可知，耦合計算結(jié)果顯示為湍流燒蝕外形，頂點后退量約為28.4 mm；而解耦計算結(jié)果顯示為轉(zhuǎn)捩燒蝕外形， 頂點后退量約為27.4 mm。耦合情況下，流場參數(shù)根據(jù)燒蝕外形計算，開始時轉(zhuǎn)捩點位于聲速點附近，轉(zhuǎn)捩后的熱流增加，燒蝕量變大，隨著外形變化，轉(zhuǎn)捩點向端頭頂點移動，致使端頭越來越尖，趨于雙錐外形。而解耦情況下，采用原始外形計算熱流，對于保持不變的外形，轉(zhuǎn)捩位置基本固定，轉(zhuǎn)捩點后熱流高，燒蝕量大；轉(zhuǎn)捩點前為層流，燒蝕量小。與之前的非耦合計算結(jié)果比較[22]，通過優(yōu)化解耦計算方法，消除了以前的非物理波動，模擬結(jié)果不再為凹陷外形。而耦合計算考慮了“熱環(huán)境?燒蝕外形”的雙向作用，模擬更精細，結(jié)果更符合真實情況。此外，差分離散依據(jù)情況自動選擇隱式或顯示格式，與無波動、無自由參數(shù)、耗散的(Non-oscillatory containing No free parameters and Dissipative，NND)格式相比較[22]，對燒蝕外形仿真的模擬程度二者相當，但前者效率更高。

圖8 端頭不同子午剖面的輪廓圖
Fig.8 Outline of nosetip at different meridian sections

圖9 迎角與側(cè)滑角相等時氣流所在子午面
Fig.9Meridian section of flow direction with angle of attack being equal to sideslip angle

圖10 層流和湍流粗糙度條件下的隨機燒蝕外形
Fig.10 Asymmetric ablation shapes with different roughs in laminar flow and turbulent flow

圖11 耦合-解耦計算燒蝕外形
Fig.11Ablation shapes with coupled/decoupled simulation

6 結(jié) 論

本文通過坐標變換、外形方程數(shù)值求解等外形重構(gòu)技術(shù)，建立了三維燒蝕外形計算方法。通過定義總迎角和等效迎角的處理方法，將零迎角燒蝕方法拓展到非零迎角、側(cè)滑角飛行條件下，應(yīng)用于不同姿態(tài)角下高超聲速飛行器的熱環(huán)境模擬問題，適用于旋成體彈頭小傾角再入燒蝕外形計算。

將Monte-Carlo隨機抽樣技術(shù)應(yīng)用于彈頭物面的層流和湍流粗糙度，結(jié)合粗糙壁邊界層轉(zhuǎn)捩及熱流增益效應(yīng)研究，模擬了彈頭燒蝕外形的非對稱溝槽。降低粗糙度可以減輕小不對稱程度，因此，通過提高彈頭表面致密性和表面光潔度將有助于燒蝕外形的對稱性和連續(xù)性。

將氣動熱、燒蝕和三維外形變化進行雙向耦合計算，采用此策略避免了解耦計算中轉(zhuǎn)捩點位置固定，而是隨著燒蝕外形的不斷變化，逐漸向端頭頂點移動，符合實際情況，因此提高了對彈頭再入燒蝕過程的模擬程度。

在現(xiàn)有零迎角燒蝕外形計算方法的基礎(chǔ)上，考慮飛行姿態(tài)角和表面粗糙度兩個重要因素，結(jié)合松耦合計算，建立了一套新的燒蝕計算方法，成功地對燒蝕非對稱外形進行了預(yù)示，這將為精細的氣動特性計算提供重要的外形參數(shù)。

[1] THYSON N, NEURINGER J, PALLONE A. Nose tipe change predictions during atmospheric re-entry[C]∥AIAA 5th Thermophysics Conference. Reston, VA: AIAA, 1970.

[2] BAKER R L. Low temperature ablator nosetip shape change at angle of attack[C]∥AIAA 10th Aerospace Sciences Meeting. Reston, VA: AIAA, 1972.

[3] 李素循. 燒蝕形狀的地面實驗?zāi)M[C]∥航空航天工業(yè)部科學技術(shù)司論文集, 1990：379-384.

LI S X. Wind tunnel tests simulate ablation shapes[C]. Proceedings of Department of Science and Technology，Ministry of Aerospace Industry, 1990：379-384 (in Chinese).

[4] CROWELL P G. Finite difference schemes for the solution of the nosetip shape change equations: SAMSO-TR-76-8 [R]. Philadelphia: General Electric Company, 1976.

[5] CHIN J H. Shape change and conduction for nosetips at angle of attack[C]∥7th Fluid and Plasma Dynamics Conference. Reston, VA: AIAA, 1974.

[6] DIRLING R B. Asymmetric nosetip shape change during atmospheric entry[C]∥AIAA 12th Thermophysics Conference. Reston, VA: AIAA, 1977.

[7] OTEY G R, ENGLISH E A. High-βre-entry vehicle recover[J]. Spacecraft，1977, 14(5): 290-293.

[8] HALL D W, NOWLAN D T. Aerodynamics of ballistic re-entry vehicles with asymmetric nosetips[J]. Spacecraft, 1978, 15(1): 55-61.

[9] 姜貴慶, 劉連元. 高速氣流傳熱與燒蝕熱防護[M]. 北京: 國防工業(yè)出版社, 2003.

JIANG G Q, LIU L Y. Heat transfor of hypersonic gas and ablation thermal protection[M]. Beijing: National Defense Industry Press, 2003(in Chinese).

[10] 黃志澄. 航天空氣動力學[M]. 北京:中國宇航出版社, 1994.

HUAGN Z C. Aerospace aerodynamics[M]. Beijing: China Astronautic Publishing House, 1994 (in Chinese).

[11] 薛成位, 王喜軍, 張文杰, 等. 戰(zhàn)略彈頭氣動攻關(guān)論文集[M]. 綿陽: 中國空氣動力研究與發(fā)展中心, 1990.

XUE C W, WANG X J, ZHANG W J, et al. Proceedings of strategic missile's aerodynamics tackled[M]. Mianyang: China Aerodynamics Research and Development Center, 1990 (in Chinese).

[12] 蔡金獅. 再入體滾轉(zhuǎn)共振與小不對稱氣動力[J]. 空氣動力學學報, 1985(3): 72-81.

CAI J S. Roll resonance and small asymmetrical aerodynamics of reentry body[J]. Acta Aerodynamica Sinica, 1985(3): 72-81(in Chinese).

[13] 黃興中, 蔡慧娟. 小不對稱體滾動力矩風洞實驗研究[J]. 空氣動力學學報, 1985(2): 69-73.

HUANG X Z, CAI H J. Roll moment research of small asymmetrical body with wind tunnel test[J]. Acta Aerodynamica Sinica, 1985(2): 69-73(in Chinese).

[14] 白葵, 馮明溪, 付光明. 小不對稱再入體滾轉(zhuǎn)氣動力測量技術(shù)[J]. 流體力學實驗與測量, 2002, 16(3): 63-67.

BAI K, FENG M X, FU G M. Experimental technique for rolling aerodynamic of slight asymmetric re-entry body[J]. Experiments and Measurements in Fluid Mechanics, 2002, 16(3): 63-67 (in Chinese).

[15] 趙俊波, 付增良, 梁斌, 等. 再入彈頭小不對稱俯仰氣動特性測量技術(shù)研究[J]. 實驗流體力學, 2015, 29(5): 55-59.

ZHAO J B, FU Z L, LIANG B, et al. Research on the wind tunnel test techniques for micro-pitching-aerodynamics of re-entry body[J]. Journal of Experiments in Fluid Mechanics, 2015, 29(5): 55-59(in Chinese).

[16] 張毅, 苗育紅, 周江華. 燒/侵蝕耦合效應(yīng)及其對再入體落點的影響[J]. 飛行力學, 1997, 15(2): 91-96.

ZHANG Y, MIAO Y H, ZHOU J H. The influence of the coupled ablation/erosion effect on reentry-body impact-point[J]. Flight Dynamics, 1997, 15(2): 91-96(in Chinese).

[17] 陸海波, 陳偉芳, 袁雷, 等. 再入體碳基材料燒蝕特性分析與工程計算[J]. 彈道學報, 2008, 2(3): 107-110.

LU H B, CHEN W F, YUAN L, et al. Ablation analysis and engineering calculation on carbon based material of reentry vehicles[J]. Journal of Ballistics, 2008, 2(3): 107-110(in Chinese).

[18] 袁湘江, 楊茂昭, 潘梅林. 小不對稱燒蝕外形的典型化處理方法[J]. 空氣動力學學報, 1997, 15(3): 386-392.

YUAN X J, YANG M Z, PAN M L. Typical processing method of slight asymmetric ablative shape[J]. Acta Aerodynamica Sinica, 1997, 15(3): 386-392(in Chinese).

[19] DIRLING R B. Performance technology program(PTP-SII) Vol. VI. A statistical model of nosetip shape change for reentry vehicles: BMO-TR-80-53 [R]. Irvine: Science Application, Inc., 1980.

[20] 陳自發(fā), 張曉晨, 王振峰, 等. 高超聲速飛行器碳基頭錐燒蝕外形計算[J]. 航空學報, 2016, 37(S1): S38-S45.

CHEN Z F, ZHANG X C, WANG Z F, et al. Hypersonic aircraft’s carbon-based nose ablation shape calculation[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(S1): S38-S45(in Chinese).

[21] 黃振中. 燒蝕端頭的瞬變外形及內(nèi)部溫度分布[J]. 空氣動力學學報, 1981(1): 53-65.

HUANG Z Z. Transient changing shape and temperature distribution of ablative nosetip[J]. Acta Aerodynamica Sinica, 1981(1):53-65(in Chinese).

[22] 國義軍, 童福林, 桂業(yè)偉. 燒蝕外形方程差分計算方法研究(II:耦合計算)[J].空氣動力學學報, 2010, 28(4):441-445.

GUO Y J, TONG F L, GUI Y W. Finite difference schemes for solution of the nosetip shape change equation(Part II, coupling calculation)[J]. Acta Aerodynamica Sinica, 2010, 28(4):441-445 (in Chinese).

[23] 張志成. 高超聲速氣動熱和熱防護[M]. 北京: 國防工業(yè)出版社, 2003.

ZHANG Z C. Aerothermal heating and thermal protection of spacecraft at hypersonic conditions[M]. Beijing: National Defend Industry Press, 2003(in Chinese).

[24] HALL D W. The three-dimensional shock and pressure (3DSAP) approximate flow field technique: SAMSO-TR-77-C145 [R]. Philadelphia: General Electric Company, 1977.

Simulationof3Dasymmetricablationshapeofreentrymissile

ZHOUShuguang1,*,GUOYijun2,HELixin2,LIUXiao2

1.StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China2.ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

Strategicmissilesareheatedbyhypersonicgaswhiletheyarereenterintotheearth’satmosphereatsmallattitudeangles.Thethermalenvironmentissoharshthatthelossofsurfacematerialfromspacecraftthroughevaporationormeltingcausedbyfrictionwiththeatmosphereisaconsiderableamountofablation,therebyformingavisibleasymmetricshape.Shapechange,matterinjecttoboundarylayerandanumberofdifferenttypesofpatternsbyablationwillaffecttheaerodynamics,aerothermodynamicsandtemperaturedistributioninthematters.Thermalenvironment,ablation,andtheshapearethecauseandeffectforeachother,andthereforemaketheimpactscattertospreadout.Accordingtoequivalenttransformanalysis,theflowfieldrelativetothevehicleatsomeanglesofattackandsideslipcanbeequaltothatatatotalangleofincidence.Theverticalaxisatameridianplaneistransformedtotheflowdirection,whereforetheflowfliedisatthezeroangleofattack.RandomsamplingofsurfaceroughsisconductedwithMonte-Carlostatisticalmeasuretosimulatetherandomdistributionofthetransitionregionsatdifferentcircumferentialpositions.Theroughsizesarecreatedbytherandomizer.Thecalculationprocedurehasbeenimprovedwiththesetechniques,andthelooselycoupledsimulationthatcanexhibitthemovementofthetransitionpointtowardsthevertexofthenosetipisalsoemployed.Thenumericalcalculationmethodisalsousedtosolethenonlinearpartialdifferentialequationfortheablativeshape.Theirregularshapeofthenosetipofballisticmissilesatsmallincidenceangleinreentryintotheearth’satmospherecanbesuccessfullysimulatedwiththemethodsappliedproposed.

reentrymissile;asymmetricbody;thermalprotection;aerothermal;ablation;carbon/carbonmaterial;looselycoupledsimulation

2017-05-08；

2017-05-31；

2017-06-20；Publishedonline2017-06-291323

URL：http://hkxb.buaa.edu.cn/CN/html/20171211.html

NationalKeyResearchandDevelopmentProgramofChina(2016YFA0401200)

.E-mailzhoushuguang111@163.com

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121397

2017-05-08；退修日期2017-05-31；錄用日期2017-06-20；網(wǎng)絡(luò)出版時間2017-06-291323

http://hkxb.buaa.edu.cn/CN/html/20171211.html

國家重點研發(fā)計劃(2016YFA0401200)

.E-mailzhoushuguang111@163.com

周述光，國義軍，賀立新，等．再入彈頭三維非對稱燒蝕外形模擬J．航空學報,2017,38(12):121397.ZHOUSG,GUOYJ,HELX,etal.Simulationof3DasymmetricablationshapeofreentrymissileJ.ActaAeronauticaetAstronauticaSinica,2017,38(12):121397.

V11

A

1000-6893(2017)12-121397-12

李明敏)