基于OpenFOAM的三維內(nèi)彈道兩相流程序開發(fā)及其應(yīng)用

肖輝鵬,王 浩,陶如意

(1.南京理工大學 能源與動力工程學院,江蘇 南京 210094;2.中國船舶重工集團公司750試驗場,云南 昆明 650051)

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

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

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

1 數(shù)學物理模型

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

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

(1)

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

(2)

3)氣相動量方程。

(3)

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

4)固相動量方程。

(4)

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

5)氣相能量方程。

(5)

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

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

2 壓力修正算法

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

2.1 控制方程的離散

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

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

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

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

(6)

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

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

(7)

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

(8)

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

(9)

(10)

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

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

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

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

火藥顆粒間阻力模型:火藥顆粒阻力模型主要提供系數(shù)Ad的計算,因為相間阻力項的離散采用半隱藏離散,因此需要計算系數(shù)Ad,該系數(shù)定義為

(11)

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

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

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

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

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

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

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

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

⑤求解氣相能量方程,用計算得到的火藥氣體能量更新火藥氣體的溫度場、密度場。

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

a. 分別用固相通量預測方程與氣相通量預測方程計算預測通量;

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

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

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

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

a. 求解火藥表面溫度輸運方程,計算火藥表面溫度;

b. 求解火藥已燃相對厚度輸運方程,計算火藥已燃相對厚度;

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

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

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

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

3 算例驗證

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

3.1 壓力波捕捉能力驗證

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

激波管的結(jié)構(gòu)示意圖如圖2所示,激波管總長為10 m,左端高壓區(qū)域長度為5 m,右端低壓區(qū)域長度為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時,一維激波管問題OpenFOAM計算得到的數(shù)值解與解析解的壓力p和速度vx曲線的對比圖。從圖中可以得到,精確解與OpenFOAM數(shù)值解無論是在激波管的膨脹區(qū)域還是間斷區(qū)域都符合較好,說明基于OpenFOAM的內(nèi)彈道兩相流程序具有足夠的精度,對壓力波捕捉具有高分辨率。

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

3.2 全局守恒性驗證

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

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

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

表2給出了2種方法計算得到的物理量,表中,te表示火藥燃燒結(jié)束時間,me表示火藥燃燒結(jié)束時燃氣總質(zhì)量,Ee表示火藥燃燒結(jié)束時燃氣總能量。密閉空間中生成的火藥燃氣總質(zhì)量結(jié)果分別為46.33 g與46.12 g,誤差為0.45%,可以認為OpenFOAM實現(xiàn)的三維內(nèi)彈道兩相流程序滿足全局質(zhì)量守恒性。密閉空間生成火藥燃氣總能量值分別為47 725 J與47 710 J,誤差為0.03%,可以認為火藥燃燒釋放的能量等于火藥燃氣增加的總能,這驗證了OpenFOAM實現(xiàn)的三維內(nèi)彈道兩項流程序滿足全局能量守恒性。

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

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

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

3.3 試驗驗證

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

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

3.3.1 仿真模型的建立

本文將全部火藥燃氣可流動區(qū)域作為求解計算域。為了減少計算量,建立拋撒系統(tǒng)燃氣可流動區(qū)域的軸對稱模型,計算域網(wǎng)格如圖7所示。

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

3.3.2 仿真與試驗結(jié)果的對比

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

圖8 測壓孔位置壓力-時間對比曲線

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

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

3.3.3 點傳火過程特性分析

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

圖9 不同時刻的火藥著火程度云圖

圖9(a)～9(c)說明點火基座小孔破膜后,高溫火藥燃氣與火藥顆粒從點火基座上的小孔流入傳火管中,一部分高溫點火藥燃氣與已燃火藥顆粒從傳火孔流入中心管,另一部分繼續(xù)沿著傳火管軸向流動。由于火藥燃氣到達每個傳火孔的時間具有延遲性,因此靠近首孔位置的拋撒藥首先被引燃,然后依次將靠近下一個傳火孔處的拋撒藥引燃,說明點火過程具有延遲性,延遲性的大小取決于傳火管內(nèi)點火藥燃氣的流動速度。

由圖9(d)～9(f)可見,首孔位置附近區(qū)域的拋撒藥已經(jīng)全部達到著火溫度,此處的拋撒藥燃燒產(chǎn)生高溫、高壓燃氣,高壓的燃氣將拋撒藥向四周推散,拋撒藥床與中心管內(nèi)壁的空間屬于“自由空間”,已燃的拋撒藥與拋撒藥燃氣主要向自由空間流動,而部分拋撒藥與拋撒藥燃氣向未燃燒的拋撒藥床區(qū)域推進。由于中心管管壁附近存在自由空間,流入此處的已燃拋撒藥與高溫燃氣快速地向中心管底部推進,從而加速引燃靠近管壁附近的拋撒藥床。因此,在設(shè)計中心炸管式拋撒系統(tǒng)裝藥結(jié)構(gòu)時,拋撒藥與管壁之間要預留一定自由空間,以利于提高拋撒藥引燃的速度。

從圖9(g)～9(h)可以看出,末端傳火孔處的拋撒藥已燃區(qū)域要大于前面幾個傳火孔處的拋撒藥已燃區(qū)域,其原因是,在傳火管內(nèi)高速流動燃氣到達末端后,撞擊壁面后經(jīng)過反射,使末端傳火孔處的壓力以及能量迅速增大,導致末端傳火孔處的拋撒藥燃燒區(qū)域要大于前面幾個傳火孔附近拋撒藥已燃區(qū)域。在2.69～3.5 ms時間段內(nèi),在末尾11個傳火孔處的拋撒藥床已燃區(qū)域是從上至下擴大的,說明該處拋撒藥床主要是被管壁附近的拋撒藥燃氣所引燃,彌補了傳火能量在此處能量的不足。

4 結(jié)論

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

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

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

③內(nèi)彈道兩相流求解器具有全局守恒性:火藥燃氣減少的總質(zhì)量等于火藥燃氣增加的總質(zhì)量,火藥燃燒產(chǎn)生的總能量等于全部內(nèi)能和動能的總和。

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