王革， 李冬冬， 韓萬(wàn)之， 張瑩

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院, 黑龍江 哈爾濱 150001)

基于水平集函數(shù)的固體火箭發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)模擬方法研究

王革， 李冬冬， 韓萬(wàn)之， 張瑩

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院, 黑龍江 哈爾濱 150001)

針對(duì)固體火箭發(fā)動(dòng)機(jī)裝藥燃面退移下的瞬態(tài)內(nèi)流場(chǎng)模擬方法進(jìn)行研究，采用水平集方法捕獲發(fā)動(dòng)機(jī)工作過程中的裝藥燃面，多孔介質(zhì)模型約束裝藥區(qū)域的流動(dòng)特性，建立了一種可以準(zhǔn)確模擬發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)瞬態(tài)流動(dòng)的水平集和多孔介質(zhì)耦合方法(LSPM)。采用該方法對(duì)有壁面退移的圓管通道流動(dòng)問題進(jìn)行了計(jì)算，對(duì)短內(nèi)燃管形裝藥和復(fù)雜翼柱形裝藥發(fā)動(dòng)機(jī)進(jìn)行了模擬。研究結(jié)果表明：LSPM計(jì)算結(jié)果與動(dòng)網(wǎng)格方法吻合較好，可以較好地處理有界面退移的加質(zhì)流動(dòng)問題；LSPM壓力和燃面的計(jì)算結(jié)果與零維內(nèi)彈道結(jié)果基本一致，可以準(zhǔn)確計(jì)算裝藥燃面退移下的發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)。

兵器科學(xué)與技術(shù)； 裝藥燃面退移； 瞬態(tài)內(nèi)流場(chǎng)； 水平集方法； 多孔介質(zhì)模型

0 引言

隨著火箭發(fā)動(dòng)機(jī)性能的不斷提高和裝藥設(shè)計(jì)的復(fù)雜化程度不斷增加，火箭發(fā)動(dòng)機(jī)瞬態(tài)流場(chǎng)的仿真模擬占據(jù)著越來(lái)越重要的地位。

目前，國(guó)內(nèi)外在裝藥燃面退移計(jì)算方法上研究比較多，主要有通用坐標(biāo)法[1-4]、實(shí)體造型法[5-8]、動(dòng)網(wǎng)格方法[9-12]、界面追蹤法[13-16]和最小距離函數(shù)(MDF)法[17-20]。其中實(shí)體造型法主要依托計(jì)算機(jī)輔助設(shè)計(jì)(CAD)軟件平臺(tái)，進(jìn)行二次開發(fā)，可以得到任意時(shí)刻燃面的大小，獲得與藥柱體積相關(guān)的信息，可視化顯示藥柱退移過程；動(dòng)網(wǎng)格方法易于實(shí)現(xiàn)界面退移和流動(dòng)的耦合計(jì)算，但處理復(fù)雜幾何邊界的流動(dòng)問題較為困難；界面追蹤法將發(fā)動(dòng)機(jī)燃面視為流場(chǎng)中的自由邊界，通過一定手段追蹤或者捕獲燃面位置，然后求得燃面面積和裝藥體積，其中用的最多的是水平集(level-set)方法；MDF法可以處理燃面計(jì)算中普遍存在的通用性和穩(wěn)定性問題。

裝藥燃面退移和瞬態(tài)流場(chǎng)耦合計(jì)算方法的研究還不夠成熟。Michael等[17]、Willcox等[18]利用MDF方法開展了三維裝藥燃面退移問題的相關(guān)研究，發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)計(jì)算采用一維單向可壓縮流方程，在此基礎(chǔ)上進(jìn)行了固體火箭發(fā)動(dòng)機(jī)內(nèi)彈道仿真計(jì)算。賀征等[21]采用動(dòng)網(wǎng)格方法對(duì)星形裝藥燃面退移及內(nèi)流場(chǎng)進(jìn)行了一體化計(jì)算，計(jì)算結(jié)果和準(zhǔn)穩(wěn)態(tài)計(jì)算結(jié)果以及試驗(yàn)數(shù)據(jù)的對(duì)比表明，動(dòng)網(wǎng)格方法可以得到更好的結(jié)果。李強(qiáng)等[12]使用“面偏移”方法進(jìn)行了藥柱燃面退移問題研究，配合網(wǎng)格光順和網(wǎng)格重構(gòu)技術(shù)處理流體域網(wǎng)格的變形和扭曲，使用經(jīng)驗(yàn)公式確定裝藥燃面的當(dāng)?shù)赝艘扑俣?。使用該方法模擬的三維藥柱燃面退移下的內(nèi)流場(chǎng)結(jié)果和試驗(yàn)結(jié)果吻合比較好。 Gueyffier等[22]使用level-set和傅里葉變換等多種復(fù)體網(wǎng)格適應(yīng)技術(shù)避免網(wǎng)格的重新劃分，將這套技術(shù)嵌入計(jì)算可壓縮、多組分、湍流的模擬程序中，用來(lái)計(jì)算復(fù)雜藥型的內(nèi)流場(chǎng)，并且得到了比較好的結(jié)果。

本文在綜合國(guó)內(nèi)外研究的基礎(chǔ)上，提出了基于level-set方法和多孔介質(zhì)模型的裝藥燃面退移下的固體火箭發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)模擬方法，將該方法稱為水平集和多孔介質(zhì)耦合方法(LSPM)，并通過必要的算例，驗(yàn)證了該方法在發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)計(jì)算中的準(zhǔn)確性和處理復(fù)雜三維裝藥的能力。

1 計(jì)算方案和計(jì)算方法

1.1 裝藥燃面退移模擬方法

采用level-set方法和MDF方法組合技術(shù)計(jì)算燃面退移。level-set方法求解計(jì)算所需的初場(chǎng)和裝藥燃面退移速度，求解后可以得到當(dāng)前燃面位置以及流體區(qū)域和裝藥區(qū)域的劃分；MDF方法的初始化需要燃面位置和流體與固體(簡(jiǎn)稱流固)區(qū)域劃分信息，變換后可以為level-set方程的求解提供所需的初場(chǎng)和燃面退移速度場(chǎng)，這兩種方法互相配合，整個(gè)過程如圖1所示。

圖1 燃面退移計(jì)算方案Fig.1 Calculation scheme of grain burning surface regression

通過對(duì)level-set求解后獲得的水平集函數(shù)φ?qǐng)龅奶幚?，可以獲得裝藥燃面的位置、燃面面積、網(wǎng)格內(nèi)裝藥的體積分?jǐn)?shù)等信息，具體求解過程可以參照文獻(xiàn)[23]。

1.2 燃面退移下的瞬態(tài)流場(chǎng)計(jì)算模型

1.2.1 物理問題分析與簡(jiǎn)化

裝藥燃面退移是推進(jìn)劑燃燒的結(jié)果，而推進(jìn)劑燃燒是一個(gè)復(fù)雜的過程，從推進(jìn)劑燃燒的角度去模擬燃面退移下的瞬態(tài)內(nèi)流場(chǎng)在時(shí)間和空間尺度上都是比較困難的。相對(duì)的，不直接從燃燒機(jī)理出發(fā)，而使用學(xué)者們總結(jié)的燃面退移規(guī)律，更容易實(shí)現(xiàn)燃面退移下的發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)模擬。

為了簡(jiǎn)化計(jì)算，做如下假設(shè)：

1)裝藥組分及其物理化學(xué)性質(zhì)是均勻的；

2)整個(gè)裝藥燃燒表面同時(shí)燃燒；

3)燃燒表面上的各點(diǎn)都以當(dāng)?shù)厝妓傺胤ň€向裝藥內(nèi)部推進(jìn)；

4)燃?xì)鉃榻M分單一的常物性完全氣體；

5)不考慮侵蝕燃燒的影響：目前大都采用一維的侵蝕燃燒模型來(lái)處理侵蝕燃燒問題，但這種模型一般僅限于工程計(jì)算，對(duì)于二維、三維的流場(chǎng)計(jì)算并不適用。因此，暫不考慮侵蝕燃燒的影響；

1.2.2 控制方程

忽略重力及燃燒室內(nèi)輻射傳熱，其流動(dòng)控制方程為

(1)

level-set方程為

(2)

狀態(tài)方程為

p=ρRgT，

(3)

式中：Rg為氣體常數(shù)。

湍流方程采用Realizablek-ε模型，詳見文獻(xiàn)[24]。

1.2.3 燃面加質(zhì)加熱

假設(shè)某一體積為ΔVc的網(wǎng)格，在時(shí)間Δt范圍內(nèi)存在燃面，如圖2所示，網(wǎng)格內(nèi)燃面面積為Sc，燃?xì)饷芏葹棣裧，裝藥密度為ρp，裝藥體積分?jǐn)?shù)為vof，分析可得該網(wǎng)格中的質(zhì)量源項(xiàng)形式為

(4)

圖2 網(wǎng)格燃面示意圖Fig.2 Schematic diagram of burning surface in mesh

使用單相流模擬裝藥燃面退移下的內(nèi)流場(chǎng)，采用多孔介質(zhì)模型固定裝藥區(qū)域流速來(lái)表征裝藥的固體特性。在裝藥區(qū)域事先存在填充于該區(qū)域的燃?xì)猓虼藢?duì)質(zhì)量源項(xiàng)進(jìn)行如下修正：

(5)

能量源項(xiàng)的形式為

(6)

式中：cp為燃?xì)舛▔罕葻?；Ttot為燃?xì)饪倻兀籘r為參考溫度，取為298.15 K.

1.2.4 裝藥區(qū)域流速的限定

通過多孔介質(zhì)模型添加人工黏性阻力和慣性阻力，限制計(jì)算區(qū)域內(nèi)裝藥區(qū)域流動(dòng)速度為0 m/s. 當(dāng)裝藥區(qū)域和流體區(qū)域間存在比較大的壓差時(shí)，也允許裝藥區(qū)域流動(dòng)緩慢進(jìn)行，以平衡裝藥區(qū)域和流體區(qū)域的壓差，使燃面退移時(shí)流場(chǎng)壓力不會(huì)因?yàn)橥蛔兌l(fā)散。

假設(shè)多孔介質(zhì)是各向同性的，則多孔介質(zhì)模型的應(yīng)用可以表示為在動(dòng)量方程中添加如下源項(xiàng)：

(7)

式中：D和C分別為黏性阻力系數(shù)和慣性阻力系數(shù)。為了能夠同時(shí)模擬流體區(qū)域和裝藥區(qū)域，D和C需要滿足(見圖3)：1)對(duì)于裝藥區(qū)域，D和C為極大的數(shù)(如1010)，限定裝藥區(qū)域流動(dòng)速度為0 m/s；2)流體區(qū)域D和C的取值為0，確保對(duì)流體區(qū)域無(wú)任何影響；3)界面處系數(shù)能夠很好過渡，以免影響燃面的加質(zhì)和計(jì)算效率。

圖3 計(jì)算域阻力系數(shù)分布示意圖Fig.3 Resistance coefficient in computational domain

假設(shè)D和C是裝藥體積分?jǐn)?shù)vof的函數(shù)，將裝藥區(qū)域阻力和流體區(qū)域阻力類比于兩電阻并聯(lián)，假定裝藥提供的阻力系數(shù)極大，取值為Rs(如1010)，流體區(qū)域阻力系數(shù)較小，取值為Rf(如102)，結(jié)合網(wǎng)格內(nèi)裝藥體積分?jǐn)?shù)，則網(wǎng)格內(nèi)的阻力系數(shù)Dm可以表示為

(8)

(8)式在裝藥區(qū)域(vof=1.0)取值為Rs，在燃?xì)鈪^(qū)域(vof=0.0)取值為Rf，這樣的表達(dá)式雖然能夠比較好地使阻力系數(shù)從Rs過渡到Rf，但是不滿足條件(2)式，即流體區(qū)域阻力系數(shù)不為0，為此對(duì)(8)式乘一系數(shù)vof，即為

(9)

其曲線如圖4所示，C取值和D一致。

圖4 黏性阻力系數(shù)D分布曲線Fig.4 Distribution curve of viscous resistance coefficients

1.3 燃面退移與流場(chǎng)計(jì)算的耦合

采用level-set方法和多孔介質(zhì)模型組合技術(shù)模擬燃面退移下的發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)。其中l(wèi)evel-set方法計(jì)算裝藥燃面退移，為多孔介質(zhì)模型提供裝藥燃面面積、裝藥體積分?jǐn)?shù)等信息用于處理裝藥表面的加質(zhì)和加熱，以模擬發(fā)動(dòng)機(jī)的瞬態(tài)內(nèi)流場(chǎng)；利用發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)模擬獲得的燃燒室壓力，通過經(jīng)驗(yàn)公式可以得到裝藥燃速，為燃面退移計(jì)算提供依據(jù)。具體耦合算法如圖5所示。

圖5 燃面退移瞬態(tài)內(nèi)流場(chǎng)計(jì)算方案Fig.5 Calculation scheme of transient internal flow field with burning surface regression

2 界面退移下的瞬態(tài)流場(chǎng)計(jì)算方法驗(yàn)證

為檢驗(yàn)LSPM是否能夠用來(lái)模擬復(fù)雜裝藥發(fā)動(dòng)機(jī)燃燒室內(nèi)流場(chǎng)，對(duì)以下3種工況進(jìn)行計(jì)算分析：

1)壁面退移圓管通道流動(dòng)計(jì)算。對(duì)比LSPM和動(dòng)網(wǎng)格方法的計(jì)算結(jié)果，判斷LSPM是否可以準(zhǔn)確、合理地模擬壁面退移加質(zhì)流動(dòng)；

2)短內(nèi)燃管形發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)計(jì)算。將LSPM計(jì)算結(jié)果與零維內(nèi)彈道進(jìn)行對(duì)比，判斷LSPM是否能夠用來(lái)計(jì)算發(fā)動(dòng)機(jī)內(nèi)彈道；

3)復(fù)雜翼柱形裝藥發(fā)動(dòng)機(jī)瞬態(tài)流場(chǎng)計(jì)算。進(jìn)一步驗(yàn)證LPSM方法計(jì)算復(fù)雜裝藥發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)的能力和準(zhǔn)確性。

2.1 壁面退移圓管通道流動(dòng)

2.1.1 物理模型

如圖6所示，有一半徑為R(t)的圓柱形流場(chǎng)，圓柱面為加質(zhì)面，加質(zhì)面退移速度大小為vr，氣體注入速度大小為vi.

圖6 壁面退移通道流計(jì)算物理模型Fig.6 Physical model of channel flow with wall regression

計(jì)算中作如下假設(shè)：

1)氣體是不可壓縮的流體，不考慮換熱；

2)忽略質(zhì)量力；

3)氣體是常物性的，黏性系數(shù)為常數(shù)；

4)界面加質(zhì)是均勻的；

5)不考慮旋轉(zhuǎn)效應(yīng)，速度在周向上分量為0 m/s；

6)流動(dòng)為層流流動(dòng)。

根據(jù)動(dòng)網(wǎng)格方法和LSPM的不同，建立如下計(jì)算模型，如圖7所示。動(dòng)網(wǎng)格方法計(jì)算中計(jì)算域隨著加質(zhì)面的退移自動(dòng)擴(kuò)大，只需要給出流體區(qū)域初始幾何，如圖7(a)所示；LSPM是一種在固定網(wǎng)格下計(jì)算界面退移的方法，流體區(qū)域和裝藥區(qū)域需要根據(jù)實(shí)際問題預(yù)先確定，如圖7(b)所示。

圖7 初始計(jì)算幾何示意圖Fig.7 Illustration of initial geometry

2.1.2 計(jì)算工況

參照文獻(xiàn)[25]，進(jìn)行如下無(wú)量綱化處理：

邊界條件處理如下：

1)前端面：無(wú)滑移壁面，即ur=0 m/s，uz=0 m/s；

2)加質(zhì)面：在r=R(t)處有ur=-vi，uz=0 m/s，其中vi根據(jù)Rei具體設(shè)置；

3)軸線：對(duì)稱邊界條件，在r=0 m處有?uz/?r=0，ur=0 m/s；

4)出口：壓力出口，壓力大小為101 325.0 Pa;

5)LSPM中裝藥區(qū)域：使用多孔介質(zhì)模型對(duì)z方向和r方向分別添加源項(xiàng)：

(10)

(11)

針對(duì)以下不同工況，分別以動(dòng)網(wǎng)格方法和LSPM計(jì)算壁面退移加質(zhì)流場(chǎng)：1)Rei=10.0，Rer=10.0；2)Rei=1 000.0，Rer=10.0；3)Rei=10.0，Rer=100.0；4)Rei=1 000.0，Rer=100.0.

2.1.3 計(jì)算結(jié)果對(duì)比

1) 流體區(qū)域變化。分別提取動(dòng)網(wǎng)格方法和LSPM計(jì)算中界面退移距離e=0.005 m和e=0.073 m時(shí)的計(jì)算域，如圖8和圖9所示，圖8、圖9中深色區(qū)域?yàn)榱黧w區(qū)域。從圖8和圖9中可以看出：①在動(dòng)網(wǎng)格計(jì)算方法下，隨著加質(zhì)面的不斷退移，計(jì)算網(wǎng)格不斷增加，計(jì)算區(qū)域(流體區(qū)域)不斷擴(kuò)大；②在LSPM計(jì)算方法下，隨著加質(zhì)面的不斷退移，計(jì)算區(qū)域和網(wǎng)格雖然都是不變的，但是多孔裝藥區(qū)域不斷變小，流體區(qū)域不斷擴(kuò)大?？偟膩?lái)說，盡管方法不同，但是所得界面退移計(jì)算結(jié)果是一致的。

圖8 動(dòng)網(wǎng)格方法計(jì)算域變化圖Fig.8 Change of computational domain in dynamic mesh

圖9 LSPM計(jì)算域變化圖Fig.9 Change of computational domain in LSPM method

2)流線。提取各工況下的流線，如圖10所示。從圖10中可以看出，盡管計(jì)算方法不同，但是在相同的邊界條件下LSPM和動(dòng)網(wǎng)格方法計(jì)算流場(chǎng)得到的流線是一致的。

圖10 流線對(duì)比圖Fig.10 Comparison of stream lines

3) 速度分布。根據(jù)文獻(xiàn)[25]可知在此邊界條件下，流道的Uz/Um無(wú)量綱分布與軸向位置無(wú)關(guān)(除頭部壁面附近)，提取Z=5.0 處的軸向速度數(shù)據(jù)，變換求得Uz/Um，繪制曲線，如圖11所示。從圖11中可以發(fā)現(xiàn)，LSPM和動(dòng)網(wǎng)格方法計(jì)算得到的速度分布是基本一致的。

圖11 無(wú)量綱軸向壓力分布對(duì)比圖Fig.11 Comparison of dimensionless pressure distribution along axial direction

對(duì)壁面退移圓管通道流動(dòng)計(jì)算結(jié)果的分析表明，本文所建立的LSPM在計(jì)算界面退移下的流場(chǎng)時(shí)是合理可靠的。

2.2 短內(nèi)燃管形發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)

2.2.1 物理模型

圖12為計(jì)算中采用的幾何模型，其具體參數(shù)如表1、表2和表3所示，a為燃速系數(shù)，np為燃速壓力指數(shù)，At、Ae分別為噴管喉部面積和出口面積。

圖12 短內(nèi)燃管形發(fā)動(dòng)機(jī)計(jì)算模型圖Fig.12 Calculation model of SRM with short tubular internal combustion grain

參數(shù)數(shù)值ρp/(kg·m-3)1800.0a0.015np0.3

表2 燃?xì)鈪?shù)

考慮到實(shí)際發(fā)動(dòng)機(jī)中點(diǎn)火器的工作和推進(jìn)劑的點(diǎn)火壓力，給定如下初場(chǎng)：1)初始溫度3 532 K；2)初始?jí)毫Γ喝S2.5 MPa，二維4.0 MPa；3)根據(jù)燃面位置使用逐點(diǎn)迭代獲得最小距離函數(shù)Φe；4)初始水平集函數(shù)：

表3 噴管尺寸

(12)

式中：L取值為平均網(wǎng)格尺寸的4～5倍。

主要邊界條件(見圖13)為：1)出口：壓力出口，背壓為101 325.0 Pa，溫度為300.0 K；2)壁面：無(wú)滑移壁面條件；3)裝藥區(qū)域：使用多孔介質(zhì)模型在動(dòng)量方程中添加如(7)式所示的動(dòng)量源項(xiàng)；4)加質(zhì)面：在含有燃面的網(wǎng)格中添加如(5)式所示質(zhì)量源項(xiàng)。

圖13 邊界條件設(shè)置示意圖Fig.13 Setting of boundary conditions

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

1) 燃面面積和燃燒室壓力。短內(nèi)燃管形裝藥發(fā)動(dòng)機(jī)在尺寸上長(zhǎng)徑比比較小，燃?xì)庠谕ǖ纼?nèi)的流動(dòng)速度較低，沒有顯著的軸向分布，比較適合用零維內(nèi)彈道理論來(lái)計(jì)算燃燒室壓力。

首先給出有關(guān)零維內(nèi)彈道的解，零維內(nèi)彈道假設(shè)如下：

①燃燒室內(nèi)壓力均勻一致，不計(jì)因燃?xì)饬鲃?dòng)而造成的壓力下降；

②裝藥燃面上各點(diǎn)的燃速均勻一致，不計(jì)侵蝕燃燒；

③燃燒產(chǎn)物是常物性的單一氣體，服從完全氣體狀態(tài)方程；

④噴管流動(dòng)是準(zhǔn)定常的；

⑤燃燒室內(nèi)無(wú)熱損失。

根據(jù)以上假設(shè)，按以下過程計(jì)算零維內(nèi)彈道曲線：

②計(jì)算當(dāng)燃層厚度為iΔe1(i=0,1,2,…)時(shí)對(duì)應(yīng)的燃面面積Si(i=0,1,2,…)；

⑥計(jì)算燃燒各段肉厚的燃燒時(shí)間Δti=Δe1/r；

提取短內(nèi)燃管形裝藥發(fā)動(dòng)機(jī)各時(shí)刻的燃面面積和燃燒室壓力，與零維內(nèi)彈道進(jìn)行對(duì)比，如圖14和圖15所示，可知：LSPM在二維和三維情況下計(jì)算得到的燃面面積變化與建模軟件Pro/E輸出的面積變化基本一致；壓力曲線的初始段比零維內(nèi)彈道要高，其他位置與零位內(nèi)彈道基本一致。分析認(rèn)為產(chǎn)生這種現(xiàn)象的主要原因是：發(fā)動(dòng)機(jī)內(nèi)穩(wěn)定流動(dòng)的建立需要一定的時(shí)間，在此之前，由燃燒表面進(jìn)入流場(chǎng)的燃?xì)赓|(zhì)量流率大于由噴管流出的質(zhì)量流率，使得燃燒室壓力暫時(shí)高于穩(wěn)定值(零維內(nèi)彈道理論)。流動(dòng)達(dá)到穩(wěn)定值的過程受初場(chǎng)和發(fā)動(dòng)機(jī)初始自由容積的影響，這也與發(fā)動(dòng)機(jī)的實(shí)際工作過程是相符的。

圖14 裝藥燃面面積對(duì)比圖Fig.14 Comparison of burning surface areas

圖15 燃燒室壓力對(duì)比圖Fig.15 Comparison of combustion chamber pressures

2) 裝藥區(qū)域速度分布。提取裝藥區(qū)域的速度分布，如圖16所示。從圖16中可以看出，裝藥區(qū)域流速非常小，流動(dòng)速度都在10-3～10-2m/s量級(jí)，整個(gè)裝藥區(qū)域流動(dòng)性較小，表明多孔介質(zhì)模型可以比較好地表征裝藥區(qū)域的固體特性。

圖16 裝藥區(qū)域速度分布Fig.16 Distribution of velocity in grain domain

綜上所述，燃面面積的對(duì)比分析表明了LSPM中燃面計(jì)算方法的有效性和準(zhǔn)確性；裝藥區(qū)域流速的分布表明通過多孔介質(zhì)模型限制裝藥區(qū)域流速是有效的，可以確保裝藥區(qū)域的固體特征；內(nèi)彈道對(duì)比分析表明該方法中流動(dòng)求解方法同樣可靠?？偟膩?lái)說，所建立的LSPM在計(jì)算燃面退移下的發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)問題時(shí)，可以獲得合理的發(fā)動(dòng)機(jī)二維或三維內(nèi)彈道參數(shù)分布。

2.3 復(fù)雜翼柱形發(fā)動(dòng)機(jī)瞬態(tài)流場(chǎng)

對(duì)翼柱形裝藥發(fā)動(dòng)機(jī)進(jìn)行計(jì)算分析，并將結(jié)果與零維內(nèi)彈道結(jié)果進(jìn)行對(duì)比，以進(jìn)一步檢驗(yàn)LSPM計(jì)算復(fù)雜三維裝藥燃面退移下的發(fā)動(dòng)機(jī)瞬態(tài)內(nèi)流場(chǎng)的能力。

2.3.1 物理模型

圖17為翼柱形裝藥的結(jié)構(gòu)簡(jiǎn)圖，計(jì)算所使用的翼柱形裝藥幾何參數(shù)、推進(jìn)劑參數(shù)、噴管參數(shù)和燃?xì)鈪?shù)分別如表4、表5、表6和表7所示，幾何結(jié)構(gòu)如圖18所示。圖17、表4和表6中，L為藥柱總長(zhǎng)度，Ls為裝藥翼柱段長(zhǎng)度，d為藥柱外徑，Rc為藥柱半徑，N為翼個(gè)數(shù)，R1為翼柱段內(nèi)孔最小半徑；R2為圓弧半徑，L2為藥柱軸線到翼柱段圓弧的距離；dt為噴管喉部直徑，de為噴管出口直徑，α為翼的半角。

圖17 翼柱型裝藥幾何尺寸Fig.17 Geometric dimensions of slotted tube grain

參數(shù)數(shù)值d/mm150.0L/mm66.7N4R1/mm15.0R2/mm5.0L2/mm20.0Ls/mm30.0

表5 推進(jìn)劑參數(shù)表

表6 噴管幾何參數(shù)表

表7 燃?xì)鈪?shù)表

圖18 翼柱型裝藥計(jì)算模型圖Fig.18 Calculation model of SRM with slotted tube grain

初始條件參照文獻(xiàn)[26]，初始?jí)毫o定為1.2 MPa，初始溫度為3 200 K. 邊界條件如圖19所示，具體細(xì)節(jié)與短內(nèi)燃管形的邊界條件相同。

圖19 翼柱型裝藥邊界條件設(shè)置示意圖Fig.19 Boundary condition of slotted tube grain

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

對(duì)燃面退移下的瞬態(tài)計(jì)算結(jié)果進(jìn)行分析，提取燃燒室壓力和裝藥燃面的變化，檢驗(yàn)LSPM處理復(fù)雜裝藥發(fā)動(dòng)機(jī)瞬態(tài)流動(dòng)的能力及其準(zhǔn)確性。

圖20 翼柱型裝藥發(fā)動(dòng)機(jī)燃面退移示意圖Fig.20 Burning surface regression of slotted tube grain

1) 燃面變化對(duì)比。提取裝藥燃面退移過程中各個(gè)時(shí)刻的燃面位置，如圖20所示，從中可以看到翼柱形裝藥燃面退移的整個(gè)過程。圖21為燃面面積變化及其與建模軟件Pro/E輸出的面積值的對(duì)比情況。

圖21 裝藥燃面面積變化曲線Fig.21 Change of burning surface area

2) 內(nèi)彈道變化對(duì)比。圖22給出了裝藥燃面退移下不同時(shí)刻的流體區(qū)域壓力分布，提取燃燒室頭部壓力并與零維內(nèi)彈道理論解做對(duì)比，如圖23所示。

圖22 流體區(qū)域壓力分布Fig.22 Distribution of pressure in fluid region

圖23 燃燒室壓力變化Fig.23 Change of combustion chamber pressure

通過對(duì)圖22和圖23中數(shù)據(jù)的分析可以得出：

1)LSPM計(jì)算得到的燃面面積和建模軟件Pro/E計(jì)算得到的面積基本吻合。Pro/E得到的面積變化趨勢(shì)略提前于LSPM. 主要原因是Pro/E計(jì)算中未考慮燃燒室壓力上升段，壓力一直處于較高水平，燃面退移比考慮了初始低壓的LSPM要快。

2)LSPM計(jì)算得到的燃燒室頭部壓力從初始給定的1.2 MPa經(jīng)過約0.1 s達(dá)到穩(wěn)定工作壓力，這與隨后的最大壓力滯后時(shí)間0.1 s基本吻合； LSPM計(jì)算得到的壓力要大于零維內(nèi)彈道理論所獲得的值，分析認(rèn)為是由于LSPM計(jì)算中考慮的三維流動(dòng)造成實(shí)際的噴管喉部面積要小于幾何喉部(零維內(nèi)彈道理論采用此值)造成的。

3)圖23中壓力變化與2.2節(jié)中壓力表現(xiàn)出的現(xiàn)象有所不同，分析認(rèn)為是較大的自由容積和較小的初始?jí)毫?dǎo)致達(dá)到穩(wěn)定工作狀態(tài)的時(shí)間變長(zhǎng)，從而也避免了壓力在初始的極短時(shí)間內(nèi)突然增加到大于穩(wěn)定工作壓力值的情況。

4)圖22為L(zhǎng)SPM計(jì)算得到的不同時(shí)刻下流體區(qū)域壓力分布，可以直觀地看到整個(gè)流動(dòng)區(qū)域的壓力分布以及燃面退移下壓力分布的改變，這類信息的獲得可以為分析發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)特征提供幫助。

對(duì)裝藥燃面退移過程中燃面及燃面面積、內(nèi)彈道壓力與零維理論解的對(duì)比分析表明，LSPM可以較好地模擬三維復(fù)雜裝藥燃面退移下的瞬態(tài)內(nèi)流場(chǎng)。不僅可以獲得一般設(shè)計(jì)所需的內(nèi)彈道曲線，還可以提供有關(guān)流動(dòng)細(xì)節(jié)的信息。

3 結(jié)論

本文采用level-set方法和多孔介質(zhì)模型，實(shí)現(xiàn)了復(fù)雜裝藥燃面退移和瞬態(tài)內(nèi)流場(chǎng)的耦合計(jì)算，通過對(duì)不同界面退移下瞬態(tài)流場(chǎng)的計(jì)算和分析，得到以下結(jié)論：

1) LSPM實(shí)現(xiàn)了界面退移下的加質(zhì)管流流動(dòng)的模擬，其結(jié)果與傳統(tǒng)的動(dòng)網(wǎng)格方法計(jì)算得到的結(jié)果一致，表明LSPM可以較好地模擬界面退移下的加質(zhì)流場(chǎng)。

2) LSPM實(shí)現(xiàn)了短內(nèi)燃管形裝藥和翼柱形裝藥燃面退移下的瞬態(tài)內(nèi)流場(chǎng)模擬，體現(xiàn)了該方法處理三維復(fù)雜結(jié)構(gòu)裝藥的能力；得到的燃面面積和壓力變化與零維內(nèi)彈道結(jié)果一致，體現(xiàn)了該方法的準(zhǔn)確性；獲得的三維流動(dòng)參數(shù)分布還可以為分析發(fā)動(dòng)機(jī)工作過程中的流動(dòng)特征提供參考。

目前，本文提出的LSPM基本可以實(shí)現(xiàn)不同結(jié)構(gòu)裝藥的發(fā)動(dòng)機(jī)燃面退移和流場(chǎng)流動(dòng)的耦合計(jì)算，但仍有一些問題需要進(jìn)一步研究和解決。例如，侵蝕燃燒、裝藥熱及結(jié)構(gòu)與流動(dòng)的耦合求解方法，從推進(jìn)劑燃燒的角度考慮燃速的計(jì)算等。

References)

[1] Peterson E G, Nielsen G C. Generalized coordinate grain design and internal ballistics evaluation program[C]∥Solid Propulsion Conference. Atlantic City, NJ, US:AIAA, 1968: 490.

[2] Coats D E, Levine J N, Nickerson G R. A computer program for the prediction of solid propellant rocket motor performance, volume 1[R]. US: Air Force Rocket Propulsion Laboratory, 1975.

[3] 侯曉, 蹇澤群. 三維藥柱燃面的通用積分計(jì)算法[J]. 固體火箭技術(shù), 1993, 9(3):1-6. HOU Xiao, JIAN Ze-qun. General integration calculation for burning surface of three dimensional grain[J]. Journal of Solid Rocket Technology, 1993, 9(3):1-6. (in Chinese)

[4] 鮑福廷, 李逢春. 固體發(fā)動(dòng)機(jī)裝藥CAD[J]. 固體火箭技術(shù), 1994, 9(3):1-7. BAO Fu-ting, LI Feng-chun. Computer aided design of propellant grains for solid rocket motors[J]. Journal of Solid Rocket Technology, 1994, 9(3):1-7. (in Chinese)

[5] 田維平, 王錕. 固體發(fā)動(dòng)機(jī)藥柱CAD及燃燒模擬分析[J]. 固體火箭技術(shù), 1993,8(3):14-22. TIAN Wei-ping, WANG Kun. Computer aided design and burning simulation analysis of solid rocket motor grains[J]. Journal of Solid Rocket Technology, 1993,8(3):14-22. (in Chinese)

[6] 顏仙榮, 于勝春，王慶官. 實(shí)體造型技術(shù)與固體發(fā)動(dòng)機(jī)裝藥燃面計(jì)算[J]. 固體火箭技術(shù), 2003, 26(2):20-22. YAN Xian-rong, YU Sheng-chun, WANG Qing-guan. The solid modeling technology and the grain-burning-area calculation[J]. Journal of Solid Rocket Technology, 2003, 26(2): 20-22. (in Chinese)

[7] 劉偉. 基于Pro/E和Qt平臺(tái)的固體火箭發(fā)動(dòng)機(jī)內(nèi)彈道性能計(jì)算[D]. 哈爾濱:哈爾濱工程大學(xué), 2014. LIU Wei. Performance computing of solid propellant rocket motor[D]. Harbin: Harbin Engineering University, 2014. (in Chinese)

[8] 肖志平, 韓波. 嵌銀絲端燃裝藥耦合傳熱的內(nèi)彈道數(shù)值分析[J]. 固體火箭技術(shù), 2016, 39(2):188-193. XIAO Zhi-ping, HAN Bo. Numerical study on the interior ballistics of grain embedded with silver-wire coupled convective heat transfer process[J]. Journal of Solid Rocket Technology, 2016, 39(2):188-193. (in Chinese)

[9] Breton P L, Ribéreau D, Godfroy F. SRM performance analysis by coupling bodimensional surface burnback and pressure field computations[C]∥34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Cleveland, OH, US：AIAA/ASME/SAE/ASEE, 1998.

[10] 沈偉, 邢耀國(guó). 基于非結(jié)構(gòu)網(wǎng)格的燃面推進(jìn)算法[J].固體火箭技術(shù), 2005, 28(3):176-179. SHEN We, XING Yao-guo, A simulation method of burning surface regression based on unstructured mesh[J]. Journal of Solid Rocket Technology, 2005, 28(3):176-179. (in Chinese)

[11] Jiao X M. Face offsetting: a unified approach for explicit moving interfaces[J]. Journal of Computational Physics, 2007, 220(2):612-625.

[12] Li Q, He G Q. Coupled simulation of fluid flow and propellant burning surface regression in a solid rocket motor[J]. Computers & Fluids, 2014, 93(8):146-152.

[13] Sethian J A. Theory, algorithms, and applications of level set methods for propagating interfaces[J]. Acta Numerica, 1996, 5(5): 309-395.

[14] 秦飛. 固體火箭發(fā)動(dòng)機(jī)復(fù)雜裝藥燃面算法研究[D]. 西安：西北工業(yè)大學(xué), 2003. QIN Fei. Method research for burning surface calculation of solid rocket motor with complicated grain[D]. Xi’an: Northwestern Polytechnical University, 2003. (in Chinese)

[15] Wang D H, Fei Y, Hu F, et al. An integrated framework for solid rocket motor grain design optimization[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2014, 228(7): 1156-1170.

[16] Yao H H, Chiang C H. Tracking methods to study the surface regression of the solid-propellant grain[J]. International Journal of Engineering and Technology Innovation, 2014, 4(4):213-222.

[17] Michael M A, Brewster Q M, Tang K C, et al. Solid rocket motor internal ballistics simulation using three-dimensional grain burnback[J]. Journal of Propulsion and Power, 2007, 23(3):575-584.

[18] Willcox M A, Brewster Q M, Tang K C, et al. Solid propellant grain design and burnback simulation using a minimum distance function[J]. Journal of Propulsion and Power, 2007, 23(2):465-475.

[19] 馬長(zhǎng)禮. 固體火箭發(fā)動(dòng)機(jī)MDF燃面計(jì)算方法研究[D]. 長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué), 2007. MA Chang-li. Research of MDF burning surface calculation method for solid rocket motor[D]. Changsha: National University of Defense Technology, 2007.(in Chinese)

[20] 熊文波, 劉宇. 基于單元法的三維裝藥通用燃面計(jì)算[J]. 航空學(xué)報(bào), 2009, 30(7):1176-1180. XIONG Wen-bo, LIU Yu. Generalized burning surface calculation of three dimensional propellant based on element method[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(7):1176-1180. (in Chinese)

[21] 賀征, 郜冶. 動(dòng)態(tài)網(wǎng)格在固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)計(jì)算中的應(yīng)用研究[J]. 彈箭與制導(dǎo)技術(shù), 2008,28(1):164-166. HE Zheng, GAO Ye. The study of dynamic grids application to calculate the inner flow field of solid rocket motor[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2008,28(1):164-166. (in Chinese)

[22] Gueyffier D, Roux F X. High-order computation of burning propellant surface and simulation of fluid flow in solid rocket chamber[C]∥50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. Cleveland, OH, US: AIAA/ASME/SAE/ASEE, 2014:1-16.

[23] 王革, 韓萬(wàn)之, 李冬冬, 等. 基于水平集方法和最小距離函數(shù)法的復(fù)雜裝藥燃面退移問題研究[J]. 兵工學(xué)報(bào), 2017, 33(2):280-290. WANG Ge, HAN Wan-zhi, LI Dong-dong, et al. Research on grain burning surface regression based on level-set method and minimum distance function[J]. Acta Armamentarii, 2017, 33(2):280-290.(in Chinese)

[24] 吳德銘, 郜冶. 實(shí)用計(jì)算流體力學(xué)基礎(chǔ)[M].哈爾濱:哈爾濱工程大學(xué)出版社, 2008. WU De-ming, GAO Ye. Practical computational fluid mechanics[M]. Harbin:Harbin Engineering University Press, 2008. (in Chinese)

[25] 王革, 方磊. 固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)解析解與數(shù)值解[J]. 推進(jìn)技術(shù), 2007, 28(3):225-229. WANG Ge, FANG Lei. Comparison of analytical and numerical solutions of mean flow fields in solid rocket motors[J]. Journal of Propulsion Technology, 2007, 28(3): 225-229. (in Chinese)

[26] Puskulcu G, Ulas A. 3-D grain burnback analysis of solid propellant rocket motors: part 2-modeling and simulations[J]. Aerospace Science and Technology, 2008, 12:585-591.

Simulation of Transient Internal Flow Field of Solid Rocket Motor Based on Level-set Method

WANG Ge, LI Dong-dong, HAN Wan-zhi, ZHANG Ying

(College of Architecture and Aerospace Engineering, Harbin Engineering University, Harbin 150001, Heilongjiang, China)

A method considering combustion surface regression is developed to simulate the transient flow field of solid rocket motor. The method combining the level-set method with the porous media model can accurately simulate the transient flow field. In this method, the level-set method is used to capture the fluid-solid interface, and the porous media model is used to constrain the flow characteristics of grain region. The transient flows in the circular tube channel with interface regression and the solid rocket motor with short tube and complex wing-shaped charge are simulated to demonstrate the simulation accuracy of the proposed method. The calculated results agree well with the results obtained from dynamic grid method and zero-dimensional internal trajectory theory. It shows that the method can accurately simulate the transient flow with combustion surface regression.

ordnance science and technology; combustion surface regression; transient internal flow field; level-set method; porous media model

2017-02-08

王革( 1966—)，男，教授，博士生導(dǎo)師。E-mail: wangge@hrbeu.edu.cn

V435+.12

A

1000-1093(2017)08-1520-12

10.3969/j.issn.1000-1093.2017.08.009