基于動網(wǎng)格的裝藥燃燒的燃面退移仿真

苗海玉，劉少偉，朱柏銀

（空軍工程大學(xué)防空反導(dǎo)學(xué)院，西安 710051）

0 引言

垂直彈射技術(shù)［1］是近些年來發(fā)展起來的一門新型的發(fā)射技術(shù)，而彈射器燃燒室裝藥的燃燒過程對導(dǎo)彈彈射過程能否正常順利進行起著極為重要的作用。在彈射器工作過程中，燃燒室裝藥形成的燃?xì)饬鲌鍪鞘謴?fù)雜的，裝藥的燃面在燃燒過程中不斷地退移，流場疊加形成了更加復(fù)雜的幾何流場結(jié)構(gòu)。由于在計算過程的不斷推進中，流場的計算區(qū)域是不斷變化的，就需要準(zhǔn)確地跟蹤變化的燃面。采用傳統(tǒng)的內(nèi)流場計算方法［3-5］形成流場的計算邊界是固定的，給定的燃燒燃面與實際情況并不吻合，因而傳統(tǒng)計算方法無法準(zhǔn)確地描述復(fù)雜的非定常燃?xì)饬鲌龅募?xì)節(jié)。

在Fluent中運用移動網(wǎng)格技術(shù)［6-8］對裝藥的燃面進行跟蹤，通過提取裝藥表面的壓力并在裝藥燃速方程中計算得出單位時間內(nèi)裝藥燃燒的厚度，能夠精確地控制每一個燃面點的運動，從而精確地跟蹤運動的燃面，實現(xiàn)網(wǎng)格在計算過程中自動更新。目前國內(nèi)外針對固體火箭發(fā)動機流場的研究已經(jīng)有很多，固體裝藥在彈射器中應(yīng)用較多，但針對裝藥在燃燒室內(nèi)的燃?xì)饬鲌鲞M行實際的仿真較少。文中采用動網(wǎng)格技術(shù)，通過UDF編譯［9］，采用Fluent中的隱式、非定常的基于壓力的耦合式求解器，結(jié)合Realizable k-ε湍流模型，對彈射器燃燒室的燃?xì)饬鲌鲞M行一體化計算，研究移動燃面下的瞬態(tài)流場，并分析了燃燒室的壓力分布和裝藥附近的激波傳遞情況，實現(xiàn)了裝藥復(fù)雜燃燒化學(xué)反應(yīng)的數(shù)值仿真，計算結(jié)果能夠很好地描述彈射器燃燒室的內(nèi)彈道性能變化規(guī)律。

1 數(shù)學(xué)物理模型

1.1 基本假設(shè)

根據(jù)所研究的固體裝藥結(jié)構(gòu)特點，其采用的是單孔管狀裝藥，在彈射器中經(jīng)常被采用。建立燃燒室裝藥三維非定常模型，采用Realizable k-ε兩方程湍流模型；不考慮裝藥的點火過程，將點火壓強作為燃燒室裝藥燃燒的初始條件處理；考慮裝藥表面氣相加質(zhì)層的燃?xì)饧淤|(zhì)，而不引入具體的化學(xué)反應(yīng)。為方便計算，對燃燒室的燃燒過程作出了如下簡化：1）所有裝藥的燃燒規(guī)律［10］服從幾何燃燒定律的假設(shè)，即裝藥是按照平行層或同心層的規(guī)律逐層進行燃燒的情況下進行的；2）所有裝藥的表面同時被點燃，燃燒氣體服從理想氣體狀態(tài)方程；3）所有裝藥的燃燒滿足給定的燃燒定律，且裝藥是在平均壓力下進行燃燒的；4）裝藥在燃燒過程中不考慮侵蝕效應(yīng)；5）裝藥在燃燒室內(nèi)的燃燒過程看成是絕熱過程，高壓室壁面密封性能良好，無漏氣，是絕熱的。

1.2 控制方程

在研究燃燒室中燃?xì)饬鲃拥淖兓?guī)律時，燃?xì)饬鲃犹幱谕牧鳡顟B(tài)，因此，建立的控制方程包括質(zhì)量守恒方程、能量守恒方程以及湍流控制方程。

質(zhì)量守恒方程

式中，ρ為密度，t為時間，u、v、w 是速度矢量在 x、y、z方向上的分量。

能量守恒方程

式中，T為溫度，k為傳熱系數(shù)，cp為比熱容，ST為粘性耗散項。

湍流控制方程

1.3 計算網(wǎng)格劃分

所采用燃燒室模型中裝藥、燃燒室的結(jié)構(gòu)都是圓柱體，且在裝藥燃燒減小的過程中，不斷地有內(nèi)部網(wǎng)格移動變形和新網(wǎng)格的生成更新，所以在劃分體網(wǎng)格時，把燃燒室劃分為適應(yīng)能力強的非結(jié)構(gòu)化網(wǎng)格［11］，主要以四面體網(wǎng)格形式劃分，在適當(dāng)?shù)奈恢每梢园骟w、錐體和楔形單元。燃燒室是裝藥燃燒產(chǎn)生大量燃?xì)獾牡胤剑瑢τ嬎阋蟾?，所以非結(jié)構(gòu)化網(wǎng)格的劃分較密，燃燒室內(nèi)部空腔為主裝藥燃燒的流場計算區(qū)域，共有143 144個四面體非結(jié)構(gòu)化網(wǎng)格單元。由于管狀裝藥內(nèi)部和中心無網(wǎng)格，所以將管狀裝藥簡化為柱狀模型。燃燒室的計算模型與網(wǎng)格如圖1所示。

圖1 燃燒室計算模型與網(wǎng)格

1.4 邊界條件

1.5 燃面退移與動網(wǎng)格技術(shù)

動態(tài)網(wǎng)格的控制方程為：

式中，u為流體的流速矢量，ug為邊界網(wǎng)格的移動速度矢量，φ為待求變量，Γ為擴散系數(shù)，Sφ為源項，?V為有限容積的邊界。

裝藥燃燒的流場仿真中，需要分別編寫裝藥加質(zhì)燃燒的UDF程序和裝藥燃面退移的UDF程序。在裝藥加質(zhì)燃燒的過程中，通過模型指定宏中的DEFINE_SOURCE宏來定義裝藥燃燒運輸方程的源項，包括質(zhì)量源項、能量源項和動量源項，從而實現(xiàn)裝藥的加質(zhì)過程。在裝藥燃燒燃面退移的過程中，需要通過動網(wǎng)格模型DEFINE宏中的DEFINE_GRID_MOTION宏定義燃面的變化規(guī)律，從而準(zhǔn)確地控制裝藥燃面的移動。在網(wǎng)格節(jié)點的運動過程中，采用彈性光順法和網(wǎng)格重構(gòu)法實現(xiàn)網(wǎng)格的實時更新［15-16］。

1.6 計算流程

圖2為流場計算程序流程圖，主程序主要有流場計算模塊和動網(wǎng)格更新模塊，首先進行初始化，然后求解質(zhì)量方程、動量方程和能量方程，選擇相應(yīng)的湍流模型進行求解，如果收斂，則進行下一時間步的計算；如果不收斂，繼續(xù)進行迭代計算。動網(wǎng)格更新模塊的作用是通過調(diào)用UDF實現(xiàn)網(wǎng)格的更新，計算并更新下一時刻燃面的位置，然后進行下一時間步的仿真計算，直到預(yù)定的仿真時間結(jié)束計算。

2 計算結(jié)果及分析

圖2 流場計算程序流程圖

由于燃燒室是一個具有軸對稱性的圓柱體狀的燃燒室，在對燃燒室進行研究時只需要分析其中的各個方向的截面就能說明問題。圖3為在裝藥燃燒剛開始時燃燒室內(nèi)部的速度矢量分布圖，從圖3中可以看出在裝藥表面的氣體微團的運動方向，裝藥開始燃燒后，裝藥表面不斷有高溫燃?xì)鈴难b藥表面析出，燃?xì)獾倪\動方向是指向燃燒室四周的，隨著裝藥不斷地燃燒，流向燃燒室空間的燃?xì)庠黾虞^快。

圖3 5 ms時速度矢量分布圖

圖4 5 ms時壓力和溫度分布

從圖4可以看出，裝藥已經(jīng)加質(zhì)燃燒成功，實現(xiàn)燃?xì)獾馁|(zhì)量、動量、能量向燃燒室的注入。在裝藥表面壓力升高，但是由于裝藥表面燃?xì)馑俣容^大，導(dǎo)致裝藥表面的壓強比其他位置低，說明裝藥開始燃燒。裝藥燃燒釋放熱量使燃燒室溫度升高，由于裝藥固相面要吸收大量熱量用于裝藥的軟化和融化導(dǎo)致靠近裝藥燃燒面處溫度較低。說明能量已成功地注入裝藥的表面，裝藥開始加質(zhì)燃燒。

圖5 燃燒室網(wǎng)格隨時間變化

圖5為燃燒室網(wǎng)格隨時間變化，從圖5可以看出，燃燒室裝藥的燃燒基本為等速燃燒，隨著裝藥的燃燒，燃燒室的網(wǎng)格在不斷地更新，裝藥藥柱在不斷地減小。

圖6為裝藥燃燒過程截面壓力變化，燃?xì)庠谌紵覂?nèi)的傳播規(guī)律，從圖6可以看出燃?xì)饬餮刂b藥通道傳播，使得燃燒室的壓強逐漸升高，同時燃?xì)饬髯矒羧紵冶诿嫘纬杉げê团蛎洸?，燃燒室裝藥的頭部壓強與尾部壓強變化較大。這是由于裝藥燃燒所產(chǎn)生的高溫高速燃?xì)夂芸斓刈⑷胙b藥通道中，在裝藥通道軸向上的速度分量為零，當(dāng)高溫高速燃?xì)膺M入主流區(qū)后，沿著燃燒室壁面網(wǎng)燃燒室兩端流動，迅速與主流燃?xì)獍l(fā)生碰撞并混合在一起，獲得一定的軸向速度，主流燃?xì)庥捎谑艿阶⑷氲母邷馗咚偃細(xì)獾呐鲎埠突旌希沟没旌虾蟮娜細(xì)饧铀龠\動，燃燒室的壓強迅速升高。

圖7為燃燒室裝藥燃燒過程中裝藥藥柱垂直于軸向某一截面上的流場速度變化。由圖7（a）可看出，裝藥燃燒初期，生成的高溫高速燃?xì)獠粩嘞蜓b藥四周擴散，在裝藥周圍形成一個環(huán)形腔，在靠近藥柱的部分，速度云圖呈現(xiàn)出近似于環(huán)形的形狀。燃?xì)庠谶\動過程中撞擊到燃燒室壁面后，發(fā)生折轉(zhuǎn)。

圖6 裝藥燃燒過程截面壓力變化

圖7 燃燒室裝藥藥柱軸向Z=80 mm截面速度變化圖

由圖7可看出，在裝藥燃燒的時間段內(nèi)，隨著裝藥的不斷加質(zhì)燃燒，藥柱的燃燒速度減小，燃面不斷地發(fā)生退移。對比圖7（b）～圖7（d）發(fā)現(xiàn)，在計算過程中，隨著燃面的不斷退移，裝藥截面的燃?xì)饬魉俨粩嗉哟?，這是因為隨著燃燒室內(nèi)壓強的增大使得燃?xì)馍伤俾首兇蟆慕孛娴牧鲌鏊俣仍茍D的形狀來看，在靠近裝藥藥柱周圍壁面處，壓力云圖基本呈圓形，但在遠(yuǎn)離藥柱的地方，燃?xì)獾牧鲃犹幱诓环€(wěn)定狀態(tài)，燃?xì)庾矒舾邏菏冶诿婧蟀l(fā)生反彈和折轉(zhuǎn)，形成一個渦流區(qū)，有的渦流區(qū)速度變化很大，有的渦流區(qū)速度變化很小。

由圖7（e）可知，100 ms之后裝藥藥柱比初始狀態(tài)已經(jīng)明顯減小，藥柱快接近燃完，只剩下少許裝藥繼續(xù)燃燒，從速度的分布云圖可以看出，二維速度云圖分布越來越不明顯，說明燃?xì)饬鲗θ紵冶诿娴臎_擊作用逐漸減弱，已無明顯渦流區(qū)，此時高壓室空間的燃?xì)饬魉俣确植驾^均勻。

圖8 燃燒室上下表面溫度變化曲線

圖9 燃燒室內(nèi)彈道壓力變化曲線

圖8為0.08 s燃燒室上下表面溫度變化曲線，可以看出，此時燃燒室上下表面溫度分布出現(xiàn)差別，上表面的溫度較之下表面的溫度要高，壁面最高溫度可達(dá)2 400 K，這是由于產(chǎn)生的高溫燃?xì)馔ㄟ^前燃燒室空間向下傳播，燃燒室下表面的溫度要高于上表面，隨著燃?xì)馀c下表面發(fā)生碰撞和折轉(zhuǎn)，高溫高速的燃?xì)忾_始往回流向上表面，上表面的溫度開始上升，直到燃燒室上下表面的溫度接近。由于高溫高速的燃?xì)饬鞑粩鄾_擊燃燒室的壁面，因此，必須對燃燒室的壁面進行有效的熱防護。

圖9為燃燒室內(nèi)彈道壓力變化曲線，由圖9可知，數(shù)值計算結(jié)果和理論計算結(jié)果吻合較好，10 ms之后裝藥藥柱開始減小，燃燒室壓強迅速增加，比初始狀態(tài)明顯增大。從圖中可以看出在105 ms時藥柱快接近燃完，只剩下少許裝藥繼續(xù)燃燒，燃燒室空間的燃?xì)饬鲃于呌诜€(wěn)定，燃燒室壓強仍然在緩慢增加，之后隨著裝藥停止燃燒，燃燒室的燃?xì)獠辉倭鲃?，燃?xì)廒呌诜€(wěn)定。

圖10 不同時刻燃面瞬態(tài)退移速率曲線

圖10為不同時刻燃面瞬態(tài)退移速率曲線，可以發(fā)現(xiàn)，在裝藥開始燃燒后的一段時間里燃面退移速率有較大變化，不同時刻的變化趨勢基本相同，都是在增加，這是由于截面的燃?xì)饬魉僖恢痹龃?，裝藥在減面燃燒的過程中，在裝藥外部空間同一個截面區(qū)域，存在不同的速度分布，由于燃?xì)馔ǖ赖闹饾u擴張，燃?xì)獠粩嗯蛎?，速度變大，裝藥藥柱燃燒生成的氣體不斷地充滿燃燒室的自由容積，使燃燒室內(nèi)壓強逐漸達(dá)到平衡壓強的穩(wěn)定工作狀態(tài)。

綜上所述，裝藥在燃燒的過程中，同一截面上不同時刻的氣流速度分布在不斷變化，隨著燃面的不斷退移，燃面不斷減小，呈現(xiàn)出減面燃燒的特點，隨著燃燒室空間燃?xì)獾牟粩嘣黾?，裝藥的燃面減小使得燃燒室燃?xì)饬鲃拥目臻g擴大，在裝藥外部空間的氣流速度逐漸穩(wěn)定。

3 結(jié)論

本文針對彈射器燃燒室管狀裝藥在燃燒過程中的變化特性，通過對燃燒室裝藥的非穩(wěn)態(tài)燃燒過程進行三維數(shù)值仿真，得出如下結(jié)論：1）利用動網(wǎng)格技術(shù)能夠較好地仿真管狀裝藥在非穩(wěn)態(tài)燃燒過程中燃面隨著時間的變化，真實地再現(xiàn)燃燒室裝藥的燃燒燃面退移過程，得到彈射器燃燒室的溫度、壓強等內(nèi)彈道參數(shù)在燃面退移的過程中的變化規(guī)律；2）隨著燃面的退移，燃燒室燃?xì)馔ǖ赖淖兇螅济嫱艘频乃俾手饾u增加，燃?xì)獠粩嗟嘏蛎洠瑢ρb藥周圍的流場結(jié)構(gòu)產(chǎn)生較大影響；3）在裝藥燃燒的過程中，燃燒室的壓強逐漸升高，燃燒室上表面和下表面存在一定的溫度差，燃?xì)饬髯矒羧紵冶诿嫘纬杉げǎ?）隨著裝藥燃燒逐漸停止，燃面退移速率在單位時間內(nèi)的變化減小，整個過程的內(nèi)彈道性能曲線與指標(biāo)要求相吻合，仿真結(jié)果對燃燒室的性能研究和結(jié)構(gòu)設(shè)計提供了參考。