肖輝鵬,王 浩,陶如意

(1.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094;2.中國(guó)船舶重工集團(tuán)公司750試驗(yàn)場(chǎng),云南 昆明 650051)

近年來(lái),隨著計(jì)算機(jī)硬件的快速發(fā)展與多相流模型及其數(shù)值算法的日臻成熟,計(jì)算內(nèi)彈道學(xué)為探究武器發(fā)射系統(tǒng)燃燒室內(nèi)流場(chǎng)變化規(guī)律提供了重要研究手段。計(jì)算內(nèi)彈道從最初的零維集總參數(shù)法發(fā)展到雙一維兩相流計(jì)算內(nèi)彈道,甚至到現(xiàn)在的二維或者三維計(jì)算內(nèi)彈道兩相流技術(shù),在這個(gè)過(guò)程中內(nèi)彈道兩相流數(shù)值計(jì)算取得了諸多的研究成果[1-4]。

目前,在應(yīng)用內(nèi)彈道兩相流算法針對(duì)不同的武器發(fā)射系統(tǒng)進(jìn)行內(nèi)彈道仿真分析時(shí),計(jì)算內(nèi)彈道兩相流算法往往是獨(dú)立編程實(shí)現(xiàn)的。獨(dú)立編程實(shí)現(xiàn)內(nèi)彈道兩相流程序存在如下問(wèn)題:①仿真周期長(zhǎng),獨(dú)立編寫的內(nèi)彈道兩相流程序都是針對(duì)每一個(gè)問(wèn)題專門編寫一個(gè)程序,因此針對(duì)工程問(wèn)題仿真周期比較長(zhǎng);②代碼復(fù)用率低,代碼一般采用面向過(guò)程的編程思想,而且針對(duì)不同結(jié)構(gòu)燃燒室開(kāi)發(fā)不同的仿真程序;③不具備復(fù)雜計(jì)算域的求解能力,在代碼開(kāi)發(fā)前期往往需要將模型簡(jiǎn)化為一維或者二維模型,因此所開(kāi)發(fā)的程序不能對(duì)同結(jié)構(gòu)燃燒室進(jìn)行三維內(nèi)彈道仿真。

為了克服上述問(wèn)題,本文擬在開(kāi)源計(jì)算流體軟件OpenFOAM框架上開(kāi)發(fā)三維內(nèi)彈道兩相流求解器,減少內(nèi)彈道仿真周期,提高代碼復(fù)用率以及該程序復(fù)雜計(jì)算域求解的能力。

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

內(nèi)彈道兩相流基本假設(shè)詳見(jiàn)文獻(xiàn)[1],采用雙歐拉描述的三維內(nèi)彈道兩相流數(shù)學(xué)物理模型如下所示。

1)氣相質(zhì)量方程。

(1)

2)固相質(zhì)量方程。

(2)

3)氣相動(dòng)量方程。

(3)

式中:p為當(dāng)?shù)貧庀鄩毫?Fdrag為相間阻力。

4)固相動(dòng)量方程。

(4)

式中:Rp為顆粒間的應(yīng)力。

5)氣相能量方程。

(5)

式中:hg為火藥燃?xì)獗瓤傡?即單位質(zhì)量火藥氣體的總焓;Wg為火藥氣體對(duì)火藥顆粒所做的功;Qg為相間熱傳遞熱量。

相關(guān)輔助方程,限于篇幅本文不再列出,詳細(xì)可見(jiàn)文獻(xiàn)[5]。

2 壓力修正算法

壓力修正算法是以迭代的方式求解非線性耦合的動(dòng)量方程與壓力方程的數(shù)值方法,到目前為止,已經(jīng)將壓力修正算法成功地應(yīng)用于多相流全聲速流動(dòng)問(wèn)題。本文的三維內(nèi)彈道兩相流求解器主體是在開(kāi)源計(jì)算流體軟件OpenFOAM中的twoPhaseEuler求解器框架上實(shí)現(xiàn)的。twoPhaseEuler求解器是采用兩相流壓力修正算法實(shí)現(xiàn)的,該求解算法主要以學(xué)者Rushe和Weller所做的工作為理論依據(jù)[6-7]。

2.1 控制方程的離散

控制方程的離散包括質(zhì)量方程的離散、動(dòng)量方程的離散、能量方程的離散以及壓力方程的離散。通過(guò)質(zhì)量方程離散可以得到相分?jǐn)?shù)方程,通過(guò)動(dòng)量方程離散可以得到速度預(yù)測(cè)方程。上述方程的離散及輔助方程具體細(xì)節(jié)參考文獻(xiàn)[5]所做的工作。本文著重討論火藥著火模型數(shù)值計(jì)算方法和火藥顆粒特征尺寸數(shù)值計(jì)算方法,該數(shù)值模型是實(shí)現(xiàn)內(nèi)彈道兩相流求解器的關(guān)鍵。

2.1.1 火藥顆粒特征尺寸數(shù)值計(jì)算方法

內(nèi)彈道兩相流數(shù)值問(wèn)題,常常需要計(jì)算與火藥顆粒特征尺寸有關(guān)的物理量,例如:計(jì)算雷諾數(shù)、放熱系數(shù)時(shí)都需要預(yù)先得到火藥顆粒的等效直徑;計(jì)算火藥生成量需要預(yù)先得到火藥顆粒表面積。

由于火藥顆粒燃燒采用等面燃燒的假設(shè),因此只要計(jì)算得到火藥已燃相對(duì)厚度,其他的火藥特征參數(shù)就可以導(dǎo)出。本文采用火藥已燃相對(duì)厚度作為待求解的火藥特征參數(shù),其他火藥特征參數(shù)可由火藥已燃相對(duì)厚度導(dǎo)出。假設(shè)空間任意一點(diǎn)的火藥已燃相對(duì)厚度為位于該點(diǎn)控制體內(nèi)的所有火藥顆粒已燃厚度的平均值,為了能在歐拉坐標(biāo)系下捕捉火藥已燃相對(duì)厚度,根據(jù)文獻(xiàn)[8]計(jì)算火藥特征等效直徑思想,建立已燃相對(duì)厚度的輸運(yùn)方程:

(6)

2.1.2 火藥著火模型數(shù)值計(jì)算方法

火藥表面溫度著火模型是采用拉格朗日坐標(biāo)系描述的,為了使火藥著火模型在內(nèi)彈道雙歐拉兩相假設(shè)下實(shí)現(xiàn),將火藥著火方程轉(zhuǎn)化為歐拉描述:

(7)

式中:ap為火藥導(dǎo)溫系數(shù),λp為火藥導(dǎo)熱系數(shù),δt為時(shí)間增量。對(duì)于單相不可壓縮流體而言,速度v的散度為0,即·v=0,將其帶入方程(7)中可得:

(8)

火藥表面溫度方程采用全顯式求解,根據(jù)線性疊加原理,火藥表面溫度升高為由輸運(yùn)引起的溫度升高與源項(xiàng)引起的溫度升高之和,火藥表面溫度方程可寫成如下形式:

(9)

(10)

2.2 內(nèi)彈道兩相流求解器的實(shí)現(xiàn)

內(nèi)彈道兩相流求解的實(shí)現(xiàn)需在開(kāi)源計(jì)算流體OpenFOAM框架中加入與火藥相關(guān)的物理模型,及對(duì)兩相流twoPhaseEuler求解器的修改,主要包括在相應(yīng)的位置對(duì)源項(xiàng)、相間阻力項(xiàng)、顆粒間應(yīng)力項(xiàng)的處理,以及火藥表面溫度計(jì)算與火藥顆粒特征參數(shù)計(jì)算。

2.2.1 相關(guān)物理模型的實(shí)現(xiàn)

火藥氣體模型:在OpenFOAM中沒(méi)有火藥氣體,火藥氣體模型主要提供火藥氣體可壓縮性的計(jì)算以及氣體密度的計(jì)算。

火藥顆粒間阻力模型:火藥顆粒阻力模型主要提供系數(shù)Ad的計(jì)算,因?yàn)橄嚅g阻力項(xiàng)的離散采用半隱藏離散,因此需要計(jì)算系數(shù)Ad,該系數(shù)定義為

(11)

式中:Cf為火藥顆粒床摩擦系數(shù)。

火藥顆粒應(yīng)力模型:該模型提供火藥顆粒之間應(yīng)力的計(jì)算,火藥顆粒間應(yīng)力采用顯式離散。

2.2.2 內(nèi)彈道兩相流程序求解流程

內(nèi)彈道兩相流程序計(jì)算流程主要求解過(guò)程如下:

①開(kāi)始PIMPLE壓力修正算法循環(huán)。

②通過(guò)固相的相分?jǐn)?shù)方程求解固相相分?jǐn)?shù),用得到的固相相分?jǐn)?shù)φp更新火藥氣體的相分?jǐn)?shù)(空隙率),氣體相分?jǐn)?shù)通過(guò)關(guān)系式φg=1-φp得到。

③更新與相間阻力有關(guān)的系數(shù)Ad以及火藥顆粒間的應(yīng)力Rp。

④用前一次迭代計(jì)算得到的速度構(gòu)建動(dòng)量方程。

⑤求解氣相能量方程,用計(jì)算得到的火藥氣體能量更新火藥氣體的溫度場(chǎng)、密度場(chǎng)。

⑥壓力校正循環(huán)。

a. 分別用固相通量預(yù)測(cè)方程與氣相通量預(yù)測(cè)方程計(jì)算預(yù)測(cè)通量;

b. 構(gòu)造并求解壓力方程;

c. 修正固相通量方程與氣相通量方程;

d. 重構(gòu)固相速度與氣相速度。

⑦計(jì)算與火藥顆粒有關(guān)的參數(shù)。

a. 求解火藥表面溫度輸運(yùn)方程,計(jì)算火藥表面溫度;

b. 求解火藥已燃相對(duì)厚度輸運(yùn)方程,計(jì)算火藥已燃相對(duì)厚度;

c. 更新火藥顆粒表面積,更新火藥體積、火藥等效直徑、火藥燃燒線速度;

d. 求解火藥生成的質(zhì)量。

圖1給出了內(nèi)彈道兩相流程序計(jì)算流程圖。圖中,n(Correctors),n(OutCorrectors)為OpenFOAM計(jì)算手冊(cè)中規(guī)定的標(biāo)準(zhǔn)參數(shù),分別表示壓力修正循環(huán)需要迭代的次數(shù)和PIMPLE循環(huán)需要迭代的次數(shù)。

圖1 內(nèi)彈道兩相流程序計(jì)算流程

3 算例驗(yàn)證

本文從三方面對(duì)基于OpenFOAM實(shí)現(xiàn)的三維內(nèi)彈道兩相流程序進(jìn)行了驗(yàn)證:一是壓力波捕捉能力驗(yàn)證,二是全局守恒性驗(yàn)證,三是內(nèi)彈道試驗(yàn)驗(yàn)證。

3.1 壓力波捕捉能力驗(yàn)證

內(nèi)彈道兩相流流動(dòng)問(wèn)題常常伴隨壓力波的傳播,精確捕捉壓力波傳播過(guò)程中波形與位置是實(shí)現(xiàn)內(nèi)彈道兩相流程序的重點(diǎn)[4]。為了捕捉壓力波的傳播,內(nèi)彈道兩相流程序需具備高分辨壓力波捕捉能力。為了驗(yàn)證基于OpenFOAM實(shí)現(xiàn)的內(nèi)彈道兩相流程序的壓力捕捉能力,本文擬以“經(jīng)典算例”一維激波管問(wèn)題為例,將一維激波管問(wèn)題的精確解與內(nèi)彈道兩相流程序得到的數(shù)值解進(jìn)行對(duì)比。

激波管的結(jié)構(gòu)示意圖如圖2所示,激波管總長(zhǎng)為10 m,左端高壓區(qū)域長(zhǎng)度為5 m,右端低壓區(qū)域長(zhǎng)度為5 m,低壓區(qū)域與高壓區(qū)域的初始條件如表1所示。

圖2 激波管示意圖

表1 激波管初始條件

區(qū)域p/kPaT/Kρ/(kg·m-3)高壓區(qū)域1003001.158 623 0低壓區(qū)域103000.115 862 3

圖3分別給出了t=7 ms時(shí),一維激波管問(wèn)題OpenFOAM計(jì)算得到的數(shù)值解與解析解的壓力p和速度vx曲線的對(duì)比圖。從圖中可以得到,精確解與OpenFOAM數(shù)值解無(wú)論是在激波管的膨脹區(qū)域還是間斷區(qū)域都符合較好,說(shuō)明基于OpenFOAM的內(nèi)彈道兩相流程序具有足夠的精度,對(duì)壓力波捕捉具有高分辨率。

圖3 t=7 ms時(shí)激波管問(wèn)題數(shù)值解與精確解的對(duì)比

3.2 全局守恒性驗(yàn)證

分別采用基于OpenFOAM的內(nèi)彈道兩相流程序與經(jīng)典內(nèi)彈道集總參數(shù)法仿真密閉空間內(nèi)火藥顆粒被引燃的問(wèn)題,驗(yàn)證密閉空間總質(zhì)量守恒性以及密閉空間的總能量守恒性,最后對(duì)內(nèi)彈道兩相流求解器與集總參數(shù)法計(jì)算得到的內(nèi)彈道物理量曲線進(jìn)行對(duì)比。

在10 cm×10 cm×10 cm密閉空間內(nèi)裝填46.12 g的2#小粒黑火藥。采用基于OpenFOAM實(shí)現(xiàn)的三維內(nèi)彈道兩相流程序模擬上述裝填條件下火藥在密閉空間內(nèi)燃燒過(guò)程,整體計(jì)算域網(wǎng)格示意圖如圖4(a)所示。網(wǎng)格總數(shù)為8 000,每個(gè)正方形單元網(wǎng)格邊長(zhǎng)為0.5 mm,2#小粒黑火藥均勻分布在計(jì)算域中心的4 cm×4 cm×4 cm正方形區(qū)域內(nèi),如圖4(b)所示。

圖4 計(jì)算域網(wǎng)格

表2給出了2種方法計(jì)算得到的物理量,表中,te表示火藥燃燒結(jié)束時(shí)間,me表示火藥燃燒結(jié)束時(shí)燃?xì)饪傎|(zhì)量,Ee表示火藥燃燒結(jié)束時(shí)燃?xì)饪偰芰?。密閉空間中生成的火藥燃?xì)饪傎|(zhì)量結(jié)果分別為46.33 g與46.12 g,誤差為0.45%,可以認(rèn)為OpenFOAM實(shí)現(xiàn)的三維內(nèi)彈道兩相流程序滿足全局質(zhì)量守恒性。密閉空間生成火藥燃?xì)饪偰芰恐捣謩e為47 725 J與47 710 J,誤差為0.03%,可以認(rèn)為火藥燃燒釋放的能量等于火藥燃?xì)庠黾拥目偰?這驗(yàn)證了OpenFOAM實(shí)現(xiàn)的三維內(nèi)彈道兩項(xiàng)流程序滿足全局能量守恒性。

表2 計(jì)算結(jié)果的比較

圖5為OpenFOAM內(nèi)彈道兩相流與集總參數(shù)法計(jì)算得到的已燃相對(duì)厚度曲線與壓力曲線對(duì)比圖,其中OpenFOAM計(jì)算得到的已燃相對(duì)厚度與壓力是取所有單元的平均值。由圖可知,2種方法計(jì)算得到的對(duì)應(yīng)物理量值幾乎近似相等,其原因如下:集總參數(shù)法將火藥的三維密閉燃燒問(wèn)題轉(zhuǎn)換為零維問(wèn)題,所有物理量取其平均值,不考慮壓力波的傳播;而OpenFOAM內(nèi)彈道兩相流方法,雖然考慮了壓力波在三維密閉空間的傳播,但是火藥全部被點(diǎn)燃的初始條件、裝藥條件以及密閉空間形狀,共同決定了空間壓力擾動(dòng)能迅速傳到三維空間每個(gè)計(jì)算單元所在的位置。

圖5 OpenFOAM內(nèi)彈道兩相流與集總參數(shù)計(jì)算結(jié)果對(duì)比曲線

3.3 試驗(yàn)驗(yàn)證

圖6為中心炸管式拋撒系統(tǒng)結(jié)構(gòu)示意圖。本節(jié)擬采用基于OpenFOAM的內(nèi)彈道兩相流程序模擬拋撒機(jī)構(gòu)定容階段內(nèi)彈道過(guò)程,將測(cè)壓孔處試驗(yàn)測(cè)得的壓力曲線與數(shù)值計(jì)算得到的壓力曲線進(jìn)行對(duì)比,分析拋撒機(jī)構(gòu)的點(diǎn)傳火過(guò)程。

圖6 中心炸管式拋撒系統(tǒng)結(jié)構(gòu)示意圖

3.3.1 仿真模型的建立

本文將全部火藥燃?xì)饪闪鲃?dòng)區(qū)域作為求解計(jì)算域。為了減少計(jì)算量,建立拋撒系統(tǒng)燃?xì)饪闪鲃?dòng)區(qū)域的軸對(duì)稱模型,計(jì)算域網(wǎng)格如圖7所示。

圖7 求解計(jì)算域網(wǎng)格示意圖

3.3.2 仿真與試驗(yàn)結(jié)果的對(duì)比

圖8給出了測(cè)壓孔位置測(cè)試結(jié)果與計(jì)算結(jié)果的壓力-時(shí)間對(duì)比曲線。

圖8 測(cè)壓孔位置壓力-時(shí)間對(duì)比曲線

由圖8可見(jiàn),在0～2 ms時(shí)測(cè)壓孔位置的壓力曲線一直基本保持水平。在中心炸管式拋撒系統(tǒng)進(jìn)行點(diǎn)傳火過(guò)程中,當(dāng)點(diǎn)火具達(dá)到破膜壓力后,點(diǎn)火具內(nèi)火藥燃?xì)鈴狞c(diǎn)火基座的小孔噴射到傳火管中,傳火管中高溫高壓火藥燃?xì)庖徊糠謴膫骰鹂字辛鞒?進(jìn)入燃燒室并引燃拋撒藥,另一部分在傳火管內(nèi)沿著軸向流動(dòng);在t=2 ms時(shí),從I區(qū)壓力-時(shí)間曲線局部放大圖中可以看出,測(cè)壓孔處的壓力開(kāi)始上升,說(shuō)明點(diǎn)傳火過(guò)程已經(jīng)從中心管頂端位置進(jìn)行到中心管底端位置;在t=5 ms時(shí),中心管內(nèi)的拋撒藥進(jìn)入全部燃燒階段,拋撒藥迅速燃燒,中心管內(nèi)的壓力呈指數(shù)增長(zhǎng)的趨勢(shì)迅速增加,直至測(cè)壓孔壓力達(dá)到37.6 MPa時(shí)中心管炸裂。

通過(guò)比較測(cè)壓孔位置的測(cè)試壓力曲線和計(jì)算壓力曲線,可以看出計(jì)算結(jié)果與試驗(yàn)結(jié)果符合較好,說(shuō)明基于OpenFOAM的內(nèi)彈道兩相流程序能夠準(zhǔn)確地模擬中心炸管系統(tǒng)實(shí)際點(diǎn)傳火過(guò)程,計(jì)算程序能為類似結(jié)構(gòu)燃?xì)獍l(fā)生器優(yōu)化提供良好的理論依據(jù)。

3.3.3 點(diǎn)傳火過(guò)程特性分析

圖9給出了不同時(shí)刻中心炸管式拋撒系統(tǒng)燃燒室內(nèi)火藥著火程度云圖,圖中,T*為無(wú)量綱溫度。云圖中主要有3種顏色分布,紅色區(qū)域表示該位置的火藥顆粒已經(jīng)達(dá)到著火溫度,對(duì)于已經(jīng)達(dá)到著火溫度(大于600 K)的火藥顆粒在云圖中用T*=1表示;灰白色區(qū)域表示火藥顆粒表面溫度為常溫(300 K),在云圖取T*=0.5;藍(lán)色區(qū)域沒(méi)有火藥顆粒的存在,在編程處理時(shí),只有存在火藥顆粒的區(qū)域火藥顆粒表面溫度才有值,對(duì)于沒(méi)有火藥顆粒存在的區(qū)域火藥表面溫度的值取T*=0。

圖9 不同時(shí)刻的火藥著火程度云圖

圖9(a)～9(c)說(shuō)明點(diǎn)火基座小孔破膜后,高溫火藥燃?xì)馀c火藥顆粒從點(diǎn)火基座上的小孔流入傳火管中,一部分高溫點(diǎn)火藥燃?xì)馀c已燃火藥顆粒從傳火孔流入中心管,另一部分繼續(xù)沿著傳火管軸向流動(dòng)。由于火藥燃?xì)獾竭_(dá)每個(gè)傳火孔的時(shí)間具有延遲性,因此靠近首孔位置的拋撒藥首先被引燃,然后依次將靠近下一個(gè)傳火孔處的拋撒藥引燃,說(shuō)明點(diǎn)火過(guò)程具有延遲性,延遲性的大小取決于傳火管內(nèi)點(diǎn)火藥燃?xì)獾牧鲃?dòng)速度。

由圖9(d)～9(f)可見(jiàn),首孔位置附近區(qū)域的拋撒藥已經(jīng)全部達(dá)到著火溫度,此處的拋撒藥燃燒產(chǎn)生高溫、高壓燃?xì)?高壓的燃?xì)鈱伻鏊幭蛩闹芡粕?拋撒藥床與中心管內(nèi)壁的空間屬于“自由空間”,已燃的拋撒藥與拋撒藥燃?xì)庵饕蜃杂煽臻g流動(dòng),而部分拋撒藥與拋撒藥燃?xì)庀蛭慈紵膾伻鏊幋矃^(qū)域推進(jìn)。由于中心管管壁附近存在自由空間,流入此處的已燃拋撒藥與高溫燃?xì)饪焖俚叵蛑行墓艿撞客七M(jìn),從而加速引燃靠近管壁附近的拋撒藥床。因此,在設(shè)計(jì)中心炸管式拋撒系統(tǒng)裝藥結(jié)構(gòu)時(shí),拋撒藥與管壁之間要預(yù)留一定自由空間,以利于提高拋撒藥引燃的速度。

從圖9(g)～9(h)可以看出,末端傳火孔處的拋撒藥已燃區(qū)域要大于前面幾個(gè)傳火孔處的拋撒藥已燃區(qū)域,其原因是,在傳火管內(nèi)高速流動(dòng)燃?xì)獾竭_(dá)末端后,撞擊壁面后經(jīng)過(guò)反射,使末端傳火孔處的壓力以及能量迅速增大,導(dǎo)致末端傳火孔處的拋撒藥燃燒區(qū)域要大于前面幾個(gè)傳火孔附近拋撒藥已燃區(qū)域。在2.69～3.5 ms時(shí)間段內(nèi),在末尾11個(gè)傳火孔處的拋撒藥床已燃區(qū)域是從上至下擴(kuò)大的,說(shuō)明該處拋撒藥床主要是被管壁附近的拋撒藥燃?xì)馑?彌補(bǔ)了傳火能量在此處能量的不足。

4 結(jié)論

本文給出了三維內(nèi)彈道兩相流的數(shù)值方法實(shí)現(xiàn)過(guò)程,并基于OpenFOAM開(kāi)發(fā)了三維內(nèi)彈道兩相流求解器,最后進(jìn)行了求解器的正確性驗(yàn)證,得出以下結(jié)論:

①該求解器可以針對(duì)不同結(jié)構(gòu)燃燒室內(nèi)彈道兩相流問(wèn)題進(jìn)行仿真分析,減少了對(duì)工程問(wèn)題仿真分析周期,借助OpenFOAM的前處理能力,所開(kāi)發(fā)的內(nèi)彈道兩相流求解器可以分析復(fù)雜的三維計(jì)算域火藥燃燒流動(dòng)問(wèn)題。

②內(nèi)彈道兩相流求解器具有足夠的精度,能夠精確捕捉壓力波傳播過(guò)程中的波形與位置,無(wú)需加入人工黏性項(xiàng),對(duì)壓力波捕捉具有高分辨率。

③內(nèi)彈道兩相流求解器具有全局守恒性:火藥燃?xì)鉁p少的總質(zhì)量等于火藥燃?xì)庠黾拥目傎|(zhì)量,火藥燃燒產(chǎn)生的總能量等于全部?jī)?nèi)能和動(dòng)能的總和。

④基于OpenFOAM的內(nèi)彈道兩相流求解器能夠準(zhǔn)確地模擬中心炸管式拋撒系統(tǒng)實(shí)際點(diǎn)傳火過(guò)程,計(jì)算程序能為類似結(jié)構(gòu)的拋撒裝置結(jié)構(gòu)優(yōu)化提供良好的理論依據(jù)。