肖志平，韓　波，張　亮，職世君，邢國強

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

?

嵌銀絲端燃裝藥耦合傳熱的內(nèi)彈道數(shù)值分析

肖志平，韓波，張亮，職世君，邢國強

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

針對內(nèi)埋金屬絲端燃裝藥燃燒過程中出現(xiàn)壓強峰值的問題，耦合金屬絲區(qū)域的傳熱過程，建立了內(nèi)埋金屬絲端燃裝藥燃速計算模型。通過對有限元軟件，進行二次開發(fā)實現(xiàn)了內(nèi)埋金屬絲裝藥傳熱及燃面退移過程的模擬。結果顯示，沿金屬絲被加熱推進劑的溫度、燃速、增速比受燃燒室壓強影響較大。在此基礎上，進行了嵌單根銀絲端燃裝藥內(nèi)彈道仿真，得到了壓強，燃面面積，沿金屬絲燃速隨時間變化規(guī)律。結果表明，計算模型可有效模擬內(nèi)埋金屬絲端燃裝藥燃燒過程，與實驗數(shù)據(jù)吻合較好，從理論上解釋了燃燒過程中出現(xiàn)壓強峰值的實驗現(xiàn)象。

嵌金屬絲裝藥；二次開發(fā)；內(nèi)彈道；傳熱；燃速

0　引言

端燃藥柱通過內(nèi)埋金屬絲，可增加燃燒面積，達到提高總壓的目的。該過程的基本原理是在裝藥燃燒過程中，金屬絲受到高溫燃氣作用，溫度急劇升高，從而對金屬絲附近的裝藥進行了加熱，提高了裝藥的溫度。金屬絲附近藥柱溫度提高，導致其燃速較其他各處要快，在金屬絲處將形成錐面，隨著錐面的增大，燃燒面積也不斷增大，整個燃燒室的壓強也相應提高。

試驗顯示，嵌金屬絲后，由于燃面面積增大，壓強會顯著增加，并在初始時刻出現(xiàn)壓強峰[1]，之后壓強達到平衡。

國外King M 等[2]用數(shù)值方法較系統(tǒng)地研究了金屬絲的材質(zhì)、直徑以及裝藥的性質(zhì)、間隙等對嵌金屬絲端燃藥柱燃速的影響。Gossant B等[3]計算了嵌金屬絲端燃裝藥的增速比。國內(nèi)張有為等探討了降低嵌金屬絲端燃藥柱固體火箭發(fā)動機點火初始推力峰值的方法[4]，并先后進行一系列嵌金屬絲端燃裝藥的特性研究[5-6]。關于該方面的數(shù)值計算，研究人員都開展了大量仿真工作[7-9]。但由于問題的復雜性，現(xiàn)階段對嵌金屬絲端燃裝藥的燃面變化過程還缺乏深入的研究，上述文獻多采用簡化模型，將燃面等效為母線為直線的錐面，計算得到的內(nèi)彈道曲線并沒有捕捉到初始壓強峰值，與試驗不符。而實際過程中，燃面形貌在不斷發(fā)生變化，壓強會出現(xiàn)較高峰值，最后達到穩(wěn)態(tài)時，燃燒錐面的母線才是直線。

本文通過對有限元軟件進行二次開發(fā)耦合嵌金屬絲端燃裝藥的傳熱和燃面變化過程，通過數(shù)值仿真，得到了壓強-時間變化曲線，從理論上解釋了初始壓強峰值出現(xiàn)的原因。

1　數(shù)學模型

內(nèi)埋銀絲端燃裝藥內(nèi)彈道過程中，存在著內(nèi)彈道微分方程與傳熱方程的相互耦合，計算模型示意圖見圖1。

圖1　嵌銀絲裝藥模型示意圖

1.1傳熱方程

熱傳導微分方程[10-11]表達形式為

(1)

式中T、ρ、c、ki、Q分別表示溫度、密度、比熱容、熱導率、單位體積熱生成率。

燃氣對金屬絲的熱傳導邊界條件采用對流換熱邊界，熱流密度[2]采用下式：

(2)

式中ρprop、rw、cp、Tgas、Tw分別為推進劑密度、沿金屬絲推進劑燃速、推進劑定壓比熱容、火焰溫度、金屬絲表面溫度。

在計算過程中，隨著裝藥燃面的不斷推移；包裹金屬絲的推進劑燃燒完后，金屬絲會暴露出來，從而形成新的與燃氣相互作用的邊界；另外，還要考慮金屬絲在燃氣中溫度達到熔點后的熔化過程。

1.2內(nèi)彈道計算方程

內(nèi)彈道計算微分方程[12]為

(3)

式中pc、t、R、T0、Vc、S、r、ρp、CD、At分別表示燃燒室總壓、時間、氣體常數(shù)、總溫、燃燒室自由容積、燃面面積、燃速、裝藥密度、流量系數(shù)、噴管喉部截面積。

裝藥燃速計算公式：

(4)

燃速系數(shù)a(T)采用King M提出的改進模型進行計算[2]。

(5)

式中T為推進劑溫度。

金屬絲對推進劑進行加熱，相當于提高了附近推進劑初溫T，燃速系數(shù)隨推進劑溫度是單調(diào)遞增的，溫度越高，燃速系數(shù)越大，燃速也相應的增大。通過實驗可確定式(5)中的相關系數(shù)。

由于燃速隨壓強和裝藥溫度的變化規(guī)律可通過實驗來確定，因此不需計算火焰對裝藥的傳熱過程，避免了引入更多的不確定因素。

實驗表明，增速比隨金屬絲直徑先增大、后減小，即有獲得最大增速效果的直徑[2]。仿真計算表明，若直接取燃面處金屬絲與裝藥接觸點上的溫度，或者接觸面軸向增量燃燒厚度上的平均溫度來計算燃速系數(shù)，會得到金屬絲越細、增速比(沿金屬絲被加熱裝藥燃速和基礎燃速比值)越大的結論，不符合試驗規(guī)律。通過實驗和仿真交互驗算 ，本文修正了King M給出的公式[2]，采用下式作為金屬絲附近裝藥溫度：

(6)

式中α為推進劑的熱擴散系數(shù)；rwire為金屬絲半徑；r′為沿金屬絲半徑方向距離。

2　定壓下金屬絲對推進劑傳熱過程仿真

大量的實驗數(shù)據(jù)顯示，嵌金屬絲端燃裝藥，在出現(xiàn)壓強峰后，最終壓強會達到穩(wěn)定，說明在其他參數(shù)如金屬絲直徑、裝藥物性都確定的情況下，對于定壓強，沿金屬絲裝藥燃速會達到一個穩(wěn)定的數(shù)值，因為如果沿金屬絲裝藥燃速沒有達到穩(wěn)定，則燃面和壓強都不可能達到穩(wěn)定平衡態(tài)。選擇不同定壓情況，進行金屬絲對推進劑傳熱過程仿真，探究沿金屬絲推進劑燃速、增速比與壓強的關系，為內(nèi)彈道計算和分析建立基礎。

2.1仿真計算模型

仿真計算模型：金屬絲采用銀絲，直徑為0.36 mm，熔點為1 235 K；裝藥初始溫度298 K，初溫燃速系數(shù)a=1.272 8×10-4，壓強指數(shù)為n=0.35?；鹧鏈囟? 600 K，燃燒室壓強取定值，氣體常數(shù)R=287.43 J/(mol·K)。相關材料參數(shù)如表1。

由于銀絲對燃面上裝藥溫度的影響范圍很小，徑向影響范圍在0.1 mm量級。為減小計算量，傳熱計算模型取金屬絲附近0.3 mm半徑范圍的裝藥，并假設在該極小范圍內(nèi)燃面是平的，沒有傾角。定壓下傳熱計算與整個裝藥燃面的變化過程是無關的，只與金屬絲附近裝藥推移有關系。因此，只需要建立耦合金屬絲附近裝藥推移過程的傳熱分析，就可計算沿金屬絲燃速，傳熱分析模型如圖2所示。熱傳導分析可計算出溫度分布，按照式(6)計算沿金屬絲裝藥溫度，并由此得到裝藥燃速。

表1　計算模型物性參數(shù)

圖2　傳熱仿真模型示意圖

2.2仿真方法

通過python編程對有限元軟件進行二次開發(fā)，實現(xiàn)金屬絲對裝藥的傳熱仿真計算，程序計算步驟如下：

(1)建立包含金屬絲和裝藥的參數(shù)化傳熱計算模型。

(2)在Δt時間內(nèi)進行傳熱仿真計算。

(3)程序提取結果文件中藥柱燃面溫度Ti，計算裝藥燃速系數(shù)和燃速ri，燃面推移距離Δd=riΔti。確定金屬絲到達熔點的單元。

(4)通過腳本參數(shù)化方法改變計算模型，等同于去除沿金屬絲裝藥燃燒的單元和金屬絲熔斷單元，并將溫度場插值到新構型下，作為此構型下的初始溫度場，更新傳熱邊界。

(5)回到步驟(2)，不斷地循環(huán)計算，直到第i+1步的沿銀絲推進劑燃速ri+1與第i步燃速ri差值很小時，則認為燃速達到平衡，結束傳熱仿真計算。

2.3數(shù)值計算結果

通過熱傳導和金屬絲處裝藥燃面推移仿真計算過程，得到了不同壓強下沿金屬絲裝藥溫度和燃速隨時間變化規(guī)律。

圖3顯示了沿金屬絲裝藥溫度-時間變化規(guī)律。裝藥初始溫度為298 K，由于高溫燃氣提高了金屬絲的溫度，金屬絲向附近裝藥進行加熱，沿金屬絲裝藥溫度迅速上升。伴隨著傳熱過程，裝藥燃燒，燃面不斷的推移，且沿金屬絲推進劑溫度越高，燃速越大，推移越快。裝藥燃面推移速度的加快，會使得溫度提升變緩，且當推移速度達到一定程度時，燃面上被加熱的裝藥溫度就無法繼續(xù)上升，最終將達到穩(wěn)定的平衡態(tài)。壓強從1～5 MPa變化時，沿金屬絲裝藥的平衡態(tài)溫度從466 K降低到440 K。

沿金屬絲裝藥燃速rw主要受裝藥溫度和壓強控制，裝藥溫度越高，燃速系數(shù)越大，燃速也越大，由于溫度隨時間不斷增加，所以燃速也隨時間單調(diào)遞增。但裝藥燃速和溫度又相互影響，溫度達到平衡時，燃速也同時達到平衡態(tài)，見圖4中曲線。壓強從1～5 MPa變化時，沿金屬絲裝藥燃速從42.1 mm/s增大到52.9 mm/s。

圖3　裝藥溫度-時間關系曲線

圖4　燃速rw-時間關系曲線

通過多個定壓情況下的仿真分析，可計算出沿金屬絲裝藥燃面上平衡態(tài)溫度、燃速與壓強的關系。圖5、圖6中曲線分別顯示了沿金屬絲裝藥燃面上平衡態(tài)溫度隨壓強單調(diào)遞減，沿金屬絲推進劑燃速隨壓強單調(diào)遞增。壓強越高，推移越快，相當于縮短了金屬絲對裝藥加熱的時間，則其平衡態(tài)溫度將有所降低。沿金屬絲燃速由于受燃速系數(shù)(燃速系數(shù)減小)和壓強(壓強增大)兩個因素影響，仍是增大的。

圖5　平衡溫度-壓強曲線

圖6　沿金屬絲裝藥燃速-壓強曲線

圖7中，增速比隨壓強單調(diào)遞減。增速比滿足：

(7)

從式(7)可見，增速比等于燃速系數(shù)的比值，基礎燃速的燃速系數(shù)是不變的，而沿金屬絲裝藥溫度隨壓強遞減，則燃速系數(shù)比值隨壓強的增大減小。所以，增速比隨壓強增大也是減小的。

(8)

由此可得，在燃面變化過程中，金屬絲附近錐角與增速比成反比；壓強增大時，增速比減小，金屬絲處燃面錐角變平緩；壓強降低時，增速比增大，金屬絲處燃面錐角變陡峭。該結論在下章中解釋燃面變化過程中，將起到至關重要的作用。

圖7　增速比-壓強曲線

3　嵌金屬絲端燃裝藥內(nèi)彈道計算

嵌金屬絲端燃裝藥在燃面推移過程中會出現(xiàn)初始壓強峰，試驗曲線如圖8所示。

圖8　試驗壓強-時間曲線

試驗顯示，壓強峰是在點火后數(shù)秒時才出現(xiàn)的，且嵌金屬絲裝藥燃面為錐形，流量較小，壓強峰不可能由侵蝕燃燒造成。同時，恒定增速比模型(保持燃面錐角不變)數(shù)值仿真結果顯示，壓強隨時間增大到最大值后平衡，沒有出現(xiàn)達到峰值后回落的現(xiàn)象，因此，增速比為恒定數(shù)值的假定也是不能成立的。

采用python語言編程耦合傳熱過程和內(nèi)彈道計算，通過傳熱計算得到沿金屬絲裝藥燃速、增速比以及金屬絲處燃面錐角，傳入到內(nèi)彈道計算程序代碼中，進行燃面推移，提取燃面推移后燃面面積，并得到相應的壓強，然后回代到傳熱計算程序中，再次計算沿金屬絲燃速，通過增量耦合計算能夠得到壓強-時間曲線。仿真計算結果與試驗結果吻合較好。

3.1計算模型

內(nèi)彈道仿真計算模型：推進劑半徑為60 mm，初溫298 K，初溫燃速系數(shù)a=1.272 8×10-4，壓強指數(shù)為n，中心處嵌一根0.36 mm銀絲，燃氣溫度3 600 K，喉部半徑12 mm，氣體常數(shù)R=287.43 J/(mol·K)。

模型計算重點考慮金屬絲對內(nèi)彈道燃面變化的影響，為簡化模型，計算從端燃(無金屬絲時)平衡態(tài)開始。

通過python控制有限元軟件cae中的點、線的生成和移動，可實現(xiàn)燃面的自動推移過程，燃面推移方法如圖9所示。由于金屬絲處錐角會不斷地變化，因此只能用折線近似模擬連續(xù)的燃面。金屬絲處錐角決定了燃面的變化趨勢，離開金屬絲的燃面按平行推移。計算過程中在增量足夠小、點數(shù)很多的情況下，能夠達到較高的精度。

圖9　燃面推移方法示意圖

3.2內(nèi)彈道仿真計算結果

采用增量計算方法，每兩個燃面之間時間間隔為0.03 s，由于計算燃面過于密集，為便于顯示，每隔10步輸出一次結果，燃面變化過程如圖10所示。

圖10　燃面變化圖

由燃面變化圖可看出，燃面從開始的小錐面逐漸擴展，面積不斷增大，到邊界時，形狀為一個外凸的錐面，隨著燃面的變化，最終燃面演化為母線為直線的圓錐面，并達到平衡態(tài)。平衡態(tài)增速比為2.03，錐面角度約為60 °。

燃面面積隨時間變化規(guī)律見圖11，由圖11中曲線可看到，燃面面積在0.4 s之前基本沒有變化，僅有微弱的增加，隨后出現(xiàn)急劇上升段，且在4.92 s達到峰值后回落，在7 s左右出現(xiàn)一個極小值，之后微弱上升，并達到平衡。

壓強隨燃面變化，與燃面具有相同的規(guī)律，見圖12。4.92 s峰值壓強為4.01 MPa，平衡態(tài)壓強為3.34 MPa，若不嵌金屬絲，平衡壓強只有1.12 MPa，嵌金屬絲后整個壓強增加了3倍。對比圖8試驗曲線和圖12中仿真得到的壓強-時間曲線，可看到數(shù)值計算結果與試驗結果壓強變化規(guī)律具有很好的一致性。

圖11　燃面面積-時間曲線

圖12　壓強-時間曲線

金屬絲處燃面錐角是由沿金屬絲燃速與基礎燃速的比值即增速比決定的。增速比-時間關系曲線見圖13。增速比在初期很短的時間內(nèi)增大到最大值，而后下降，在4.92 s時達到一個極小值，最終達到平衡態(tài)。說明沿金屬絲錐角開始由平緩急劇變尖，但隨著壓強增大，角度開始變大、變平，最后形成一個固定角度。

為說明燃面變化規(guī)律，選取了幾個具有代表意義的燃面，作出了燃面變化示意圖，如圖14所示。

圖13　增速比-時間關系曲線

圖14　燃面變化示意圖

初期錐面形成過程：金屬絲對其附近裝藥進行加熱，使得沿金屬絲裝藥燃速較快，在金屬絲處將形成錐角，隨著平行推移過程，錐面將不斷的擴大，燃面面積增大，壓強也隨之提高。壓強的增大，導致增速比rw/r降低，即金屬絲處錐面角度會變平緩。整個錐面母線是外凸的， 見圖14中初期燃面。

最大燃面形成過程：由于錐面向邊界推移和擴展，所以整個燃面面積不斷增大，最后錐面擴展到邊界附近時，燃面面積和壓強都將達到最大值，見圖14中最大燃面。此時，燃面特點是金屬絲處燃面平緩，而邊界處燃面陡峭。

燃面減小過程：錐面到達邊界后，由于燃面是外凸的，隨著燃面推移，金屬絲附近區(qū)域平緩的燃面會擴大，而邊界處斜度較大的燃面會在平行推移過程中迅速縮小并消失， 造成整個燃面面積迅速下降，壓強也隨之降低，最后燃面和壓強都達到極小值，見圖14中減小燃面。

燃面達到平衡過程：當壓強達到極小值，增速比增大，金屬絲處燃面錐角變陡，燃面面積開始緩慢增加，最終燃面在不斷調(diào)整過程中逐漸變直，形成母線為直線的錐面，見圖14中平衡燃面。隨后的推移過程中，燃面面積不變。此時，燃速、增速比、壓強均達到穩(wěn)定平衡態(tài)。

4　結論

(1)在金屬絲材質(zhì)，直徑和裝藥性質(zhì)，燃燒室總溫等參數(shù)確定的情況下，沿金屬絲裝藥傳熱計算顯示：壓強增大時，沿金屬絲裝藥燃速也隨之增大，金屬絲附近被加熱的裝藥溫度相對降低，增速比減小，金屬絲處燃面錐角變平緩。

(2)內(nèi)彈道仿真計算結果顯示，嵌金屬絲端燃裝藥初始壓強峰是由于銀絲對裝藥增速效果隨燃燒室壓強發(fā)生變化，即增速比和金屬絲處燃面錐角變化引起的，燃面最大時，錐面是外凸的，燃面平衡時，壓強穩(wěn)定，燃面是一個母線為直線的錐面。

(3)本文所建立的嵌銀絲端燃裝藥耦合傳熱過程的內(nèi)彈道計算方法，成功地模擬了嵌金屬絲端燃裝藥點火起始階段的壓強峰形成原因，為端燃裝藥點火過程優(yōu)化提供了較準確的模型基礎和計算方法。

[1]張煒,朱慧,夏智勛,等.嵌銀絲貧氧推進劑燃燒特性[J].國防科技大學學報,2001,23(6):17-19.

[2]King M.Analytical modeling of effects of wires on solid motor ballistics[R].AIAA 89-2784.

[3]Gossant B,Godfroy F,Robert P H.Theoretical calculus of burning rate ratio in grains with embedded metal wires[R].AIAA 88-3255.

[4]張有為,王曉宏.嵌金屬絲端燃藥柱固體發(fā)動機水下點火數(shù)值模擬研究[J].固體火箭技術,2007,30(1):30-38.

[5]張有為,王曉宏.利用嵌金屬絲藥柱調(diào)節(jié)固體火箭發(fā)動機推力的計算研究[J].兵工學報,2007,28(10):1218-1222.

[6]張有為,王曉宏,楊舉賢.嵌金屬絲端燃藥柱燃燒過程的數(shù)值研究[J].推進技術,2007,28(1):4-8.

[7]肖志斌,王家鑫,王繼.嵌金屬絲端燃裝藥方案絕熱層設計方法[J].上海航天,2010(2):61-64.

[8]于寧,楊濤,劉巍,等.嵌銀絲端燃藥柱燃氣發(fā)生器理論性能分析[J].彈箭與制導學報,2009,29(1):180-183.

[9]肖志平,張澤遠,邢國強.星型藥柱固化降溫瞬態(tài)熱力耦合分析[J].航空兵器,2014(1):40-43.

[10]王勖成.有限單元法[M].北京:清華大學出版社,2003.

[11]楊世銘,陶文銓.傳熱學(4版)[M].北京:高等教育出版社,2006.

[12]董師顏,張兆良.固體火箭發(fā)動機原理[M].北京:北京理工大學出版社,1996.

(編輯：崔賢彬)

Numerical study on the interior ballistics of grain embedded with silver-wire coupled convective heat transfer process

XIAO Zhi-ping，HAN Bo，ZHANG Liang，ZHI Shi-jun，XING Guo-qiang

(Air-to-Air Missile Research Institute，Luoyang471000，China)

The study focuses on the mechanism of the appearance of initial pressure peak when using the embedded wire on the end burning grain.The convective heat transfer process was modeled and the computational model for the velocity of burning was created.The heat transfer process and the motion of the burning grain were simulated via the secondary development of the finite element software.The result indicates that the temperature,the burning velocity and augmentation ratio of the grain along the wire is strongly affected by the pressure in the combustion chamber.Based on this finding,the interior ballistics simulation of the grain with only one embedded wire was simulated ,thereby the pressure ,the burning area and the burning velocity along the wire versus time were obtained.The good agreement with the experimental data indicates that the computational model can be used properly to simulate the burning process of the grain with wires. In addition,the appearance of the initial pressure peak during the burning process was explained theoretically.

embedded wire grain；secondary development；interior ballistics；heat transfer；burning rate

2015-04-28；

2015-06-15。

國家自然科學基金(11402283)。

肖志平(1986—)，男，工程師，研究方向為發(fā)動機設計。E-mail：zhiping120@126.com

V435

A

1006-2793(2016)02-0188-06

10.7673/j.issn.1006-2793.2016.02.007