基于OpenFOAM的三維內彈道兩相流程序開發及其應用

肖輝鵬,王 浩,陶如意

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

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

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

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

1 數學物理模型

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

1)氣相質量方程。

(1)

2)固相質量方程。

(2)

3)氣相動量方程。

(3)

式中:p為當地氣相壓力,Fdrag為相間阻力。

4)固相動量方程。

(4)

式中:Rp為顆粒間的應力。

5)氣相能量方程。

(5)

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

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

2 壓力修正算法

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

2.1 控制方程的離散

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

2.1.1 火藥顆粒特征尺寸數值計算方法

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

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

(6)

2.1.2 火藥著火模型數值計算方法

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

(7)

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

(8)

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

(9)

(10)

2.2 內彈道兩相流求解器的實現

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

2.2.1 相關物理模型的實現

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

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

(11)

式中:Cf為火藥顆粒床摩擦系數。

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

2.2.2 內彈道兩相流程序求解流程

內彈道兩相流程序計算流程主要求解過程如下:

①開始PIMPLE壓力修正算法循環。

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

③更新與相間阻力有關的系數Ad以及火藥顆粒間的應力Rp。

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

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

⑥壓力校正循環。

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

b. 構造并求解壓力方程;

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

d. 重構固相速度與氣相速度。

⑦計算與火藥顆粒有關的參數。

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

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

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

d. 求解火藥生成的質量。

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

圖1 內彈道兩相流程序計算流程

3 算例驗證

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

3.1 壓力波捕捉能力驗證

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

激波管的結構示意圖如圖2所示,激波管總長為10 m,左端高壓區域長度為5 m,右端低壓區域長度為5 m,低壓區域與高壓區域的初始條件如表1所示。

圖2 激波管示意圖

表1 激波管初始條件

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

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

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

3.2 全局守恒性驗證

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

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

圖4 計算域網格

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

表2 計算結果的比較

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

圖5 OpenFOAM內彈道兩相流與集總參數計算結果對比曲線

3.3 試驗驗證

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

圖6 中心炸管式拋撒系統結構示意圖

3.3.1 仿真模型的建立

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

圖7 求解計算域網格示意圖

3.3.2 仿真與試驗結果的對比

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

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

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

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

3.3.3 點傳火過程特性分析

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

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

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

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

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

4 結論

本文給出了三維內彈道兩相流的數值方法實現過程,并基于OpenFOAM開發了三維內彈道兩相流求解器,最后進行了求解器的正確性驗證,得出以下結論:

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

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

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

④基于OpenFOAM的內彈道兩相流求解器能夠準確地模擬中心炸管式拋撒系統實際點傳火過程,計算程序能為類似結構的拋撒裝置結構優化提供良好的理論依據。