非均質(zhì)推進(jìn)劑燃燒演化模擬:算法與驗(yàn)證

張晏榕,倉(cāng) 宇,楊 碩,趙 瑜,李坤霖,王利坡

(1.上海交通大學(xué) 密西根學(xué)院,上海 200240;2.上海航天化工應(yīng)用研究所,上海 201109)

引 言

固體燃料燃燒問(wèn)題對(duì)于基礎(chǔ)研究和工程應(yīng)用都十分重要。該問(wèn)題的復(fù)雜性體現(xiàn)在多物理過(guò)程的耦合,包括氣相燃燒、固相熱傳導(dǎo)、固相結(jié)構(gòu)的表征以及氣-固界面的質(zhì)量能量與動(dòng)量交換。

由于固體燃料使用安全、能量密度高并且成本較低而被固體火箭發(fā)動(dòng)機(jī)廣泛使用[1],其中氧化物顆粒高氯酸胺(AP)與聚合黏合劑端羥基聚丁二烯(HTPB)是一種常見(jiàn)的組合。在自維持燃燒開(kāi)始后,由氣相化學(xué)反應(yīng)產(chǎn)生的熱量導(dǎo)致固相燃料的分解,進(jìn)而為氣相反應(yīng)提供反應(yīng)物。這種雙向耦合決定了界面的溫度以及界面退移速率等條件。

固體燃料燃燒過(guò)程涉及多尺度,如何有效模擬不同尺度相關(guān)的子過(guò)程引起了該領(lǐng)域?qū)W者的廣泛關(guān)注。與傳統(tǒng)的實(shí)驗(yàn)試錯(cuò)方法相比,高質(zhì)量的數(shù)值模擬對(duì)于設(shè)計(jì)及預(yù)測(cè)燃燒行為更為先進(jìn)和可靠。早期研究通常是利用一維模型或一些統(tǒng)計(jì)平均信息。Beckstead等[2]發(fā)展了基于統(tǒng)計(jì)平均氧化物顆粒的AP混合物燃燒模型。而后Beckstead等[3]展示了混合推進(jìn)劑火焰的總體結(jié)構(gòu),包括預(yù)混單質(zhì)火焰、基本擴(kuò)散火焰以及最終擴(kuò)散火焰等不同部分,并指出不同火焰之間的相互作用對(duì)于理解燃燒過(guò)程的重要性,并且也討論了壓強(qiáng)對(duì)于火焰結(jié)構(gòu)的影響。Cohen和Strand[4]結(jié)合改進(jìn)的AP單質(zhì)模型考慮了黏合劑部分的能量守恒關(guān)系,并對(duì)黏合劑的退移速率模型進(jìn)行修正。為捕捉更多的物理細(xì)節(jié),高維真實(shí)的仿真計(jì)算不可或缺。Miccio[5]提出一個(gè)使用五步化學(xué)反應(yīng)的二維數(shù)值計(jì)算模型,其中推進(jìn)劑的幾何特征由表面拓?fù)鋻呙铔Q定。Hegab等[6]基于周期性邊界條件的三明治構(gòu)型做了數(shù)值計(jì)算并討論不同參數(shù)對(duì)火焰?zhèn)鞑ニ俣鹊挠绊?。Jackson和Buckmaster[7]首次對(duì)AP和HTPB的混合燃料進(jìn)行了三維模擬,實(shí)現(xiàn)了材料微觀結(jié)構(gòu)的表征以及氣固兩相的共軛傳熱等真實(shí)物理現(xiàn)象。此外,Jackson和Buckmaster還考慮氣相與固相的強(qiáng)耦合以解決界面的非穩(wěn)態(tài)且不均勻的退移過(guò)程[8]。Massa等采用類(lèi)似的方法通過(guò)映射貼體網(wǎng)格到規(guī)則的笛卡爾網(wǎng)格來(lái)解決三維界面不規(guī)則的問(wèn)題[9-10]。 Buckmaster研究非均質(zhì)推進(jìn)劑燃燒關(guān)于壓強(qiáng)擾動(dòng)的響應(yīng)[11]。對(duì)推進(jìn)劑微尺度結(jié)構(gòu)的有效刻畫(huà)使得結(jié)構(gòu)非均質(zhì)性對(duì)燃燒過(guò)程的定量影響成為可能[12]。更全面的非均質(zhì)推進(jìn)劑燃燒建模綜述討論可以參考文獻(xiàn)[13]。

本研究提出一種基于重疊網(wǎng)格的非均質(zhì)固體推進(jìn)劑燃燒三維模擬框架。該計(jì)算框架避免整體求解包含氣固兩相的強(qiáng)耦合方程組,采用分離式的數(shù)值求解框架,具體而言,標(biāo)準(zhǔn)笛卡爾網(wǎng)格用于離散氣相以保證計(jì)算穩(wěn)定性,非結(jié)構(gòu)化四面體網(wǎng)格對(duì)固相進(jìn)行離散以更加精準(zhǔn)地刻畫(huà)由隨機(jī)填充算法生成AP粒子結(jié)構(gòu)。兩套網(wǎng)格在計(jì)算過(guò)程中都是固定的,無(wú)需網(wǎng)格重建,所以有利于計(jì)算效率和數(shù)值準(zhǔn)確性。另外水平集方法因其在處理不規(guī)則幾何上的靈活性而被采用來(lái)有效地追蹤氣-固界面的演化退移。對(duì)于浸沒(méi)在笛卡爾網(wǎng)格中的邊界單元,可以通過(guò)施加源項(xiàng)以滿(mǎn)足氣相的質(zhì)量、組分和能量的守恒條件。該套計(jì)算框架所有功能的集成基于開(kāi)源平臺(tái)OpenFOAM[1]構(gòu)建。從數(shù)值上考慮,算法實(shí)現(xiàn)清晰便捷,對(duì)于不同邊界條件都可以兼容處理,可為進(jìn)一步利用數(shù)值模擬研究該類(lèi)問(wèn)題提供基礎(chǔ)。

1 控制方程和氣-固耦合條件

1.1 氣相

三維低馬赫數(shù)帶化學(xué)反應(yīng)的變密度流動(dòng)由以下連續(xù)性方程、動(dòng)量方程、組分方程以及能量方程來(lái)描述,下標(biāo)“g”和“c”將分別表示氣相和固相:

(1)

(2)

(3)

(4)

λg=1.08×10-4T+0.0133

(5)

Prandtl數(shù)和Lewis數(shù)假定為1,采用如下理想氣體狀態(tài)方程:

p0=ρgRT∑YiMWi

(6)

式中:背壓p0設(shè)定為均一定值(由于低馬赫數(shù)流動(dòng)條件);理想氣體常數(shù)R=8.314J/(mol·K);MWi表示i組分的摩爾質(zhì)量。

1.2 固相

固相的異質(zhì)形態(tài)是復(fù)合推進(jìn)劑燃燒中具有挑戰(zhàn)性的問(wèn)題,退移速度和燃料的分解過(guò)程這些燃燒過(guò)程中的物理現(xiàn)象包括很大程度上由固體結(jié)構(gòu)的異質(zhì)性決定。模擬這種微尺度結(jié)構(gòu)通常需要隨機(jī)填充算法[14-15]。固相中的溫度方程如下:

(7)

式中:密度ρc、比熱容Cp和熱導(dǎo)率λc根據(jù)該點(diǎn)處于AP內(nèi)還是HTPB內(nèi)進(jìn)行取值。

1.3 氣固界面的邊界條件

固相氣化為氣相提供了質(zhì)量、動(dòng)量和能量的源項(xiàng);與此同時(shí),高溫的氣體也作用于氣固界面進(jìn)而影響固相內(nèi)溫度。這種耦合的數(shù)值精度對(duì)整個(gè)燃燒演化過(guò)程影響很大。

(8)

φt-rb|φ|=0

(9)

根據(jù)熱裂解關(guān)系[16-17]:

(10)

式中:TAP,s、TB,s分別表示界面上AP和黏合劑HTPB的溫度。在具體實(shí)施中,只有比較大的AP顆粒可以被建模求解,過(guò)小AP顆粒被均一地加入到黏合劑中。在這種情況下,根據(jù)Chen等的分析,混合材料的退移速率可以表示為[18]:

(11)

式中:α表示作為均質(zhì)加入黏合劑中AP的體積分?jǐn)?shù)。本研究中相關(guān)的參數(shù)見(jiàn)表1。

表1 計(jì)算所取參數(shù)值[8,10,16]Table 1 The parameter values used in calculation[8,10,16]

氣固界面需要滿(mǎn)足如下的溫度連續(xù)熱流相等的邊界條件,

Tg=Tc

(12)

(13)

1.4 化學(xué)反應(yīng)

為提高計(jì)算效率并且不影響過(guò)多的真實(shí)燃燒物理現(xiàn)象,本研究采用了已有文獻(xiàn)的簡(jiǎn)化全局化學(xué)動(dòng)力學(xué)[2,13],具體包含界面和氣相反應(yīng)的機(jī)理[9]。其中界面氣化反應(yīng)包括:

如果化學(xué)反應(yīng)放熱,此處Qs,1(2)為正值。由于實(shí)際界面反應(yīng)異常復(fù)雜(物理、化學(xué)參數(shù)的確定及凝聚相影響等),該兩步機(jī)理是對(duì)各種因素的總體體現(xiàn)。

熱分解定律可以衡量氣化的質(zhì)量流量

(14)

式中:Ac為反應(yīng)常數(shù);ρc為固相密度;Ec為反應(yīng)活化能;c為AP/HTPB;Ts為界面溫度。

氣相反應(yīng)包括:

式中:βs為化學(xué)計(jì)量參數(shù),類(lèi)似地,βs在放熱反應(yīng)中為正值。相應(yīng)的化學(xué)反應(yīng)速率為:

(15)

式中:Di(i取1,2,3)為反應(yīng)常數(shù)。這3個(gè)氣相反應(yīng)對(duì)應(yīng)了3個(gè)不同的火焰[9]: R1對(duì)應(yīng)單質(zhì)推進(jìn)劑火焰,R2對(duì)應(yīng)基本擴(kuò)散火焰,R3對(duì)應(yīng)最終擴(kuò)散火焰。

2 數(shù)值計(jì)算方法

在本研究中氣相和固相的計(jì)算采用了兩套獨(dú)立重疊的網(wǎng)格。對(duì)于氣相使用均一的笛卡爾網(wǎng)格,固相采用非結(jié)構(gòu)的四面體網(wǎng)格以刻畫(huà)非均質(zhì)結(jié)構(gòu)。采用Lubachevsky-Stillinger算法生成雙球混合結(jié)構(gòu)[18-19]。圖1展示的填充結(jié)構(gòu)是基于Miller研究中縮寫(xiě)為M24的SD-III-24[21]算例,紅球和藍(lán)球分別代表了直徑為197.6μm和49.4μm的AP粒子,小顆粒被均勻添加到黏合劑中。

圖1 M24填充模型示意圖Fig.1 Schematic diagram of M24 filling model

2.1 氣相和固相

氣相使用投影法求解壓強(qiáng)泊松方程,耦合處理連續(xù)性方程和動(dòng)量方程[20-21]?；瘜W(xué)反應(yīng)放熱帶來(lái)的效果采用Pierce提出的迭代法以實(shí)現(xiàn)[24]。具體算法流程如下:

使用Φ=(Yi,T)表示標(biāo)量,在時(shí)間步從tn推進(jìn)到tn+1過(guò)程中首先進(jìn)入序號(hào)為l的內(nèi)循環(huán)。

第一步:在時(shí)間步tn中的每一個(gè)內(nèi)循環(huán)步l使用以下初始猜測(cè)值:

(ρui)0,l=0=(ρui)n

(16)

(ρΦ)0,l=0=(ρΦ)n

(17)

ρ0=2ρn-ρn-1

(18)

(19)

第四步:根據(jù)更新后的密度ρ1,l更新Φ,表示如下:

Φ1,l=(ρΦ)1,l/ρ1,l

(20)

第六步:通過(guò)解泊松方程對(duì)壓強(qiáng)進(jìn)行修正:

(21)

第七步:更新速度和壓強(qiáng):

(22)

(23)

p1,l=p0,l+δp

(24)

重新進(jìn)行第二步,重復(fù)內(nèi)循環(huán)直至達(dá)到預(yù)設(shè)步數(shù)。通常在時(shí)間步從tn推進(jìn)到tn+1中采用3次內(nèi)循環(huán)l=0, 1, 2。這些循環(huán)完成后,對(duì)化學(xué)反應(yīng)求解以更新組分和能量源項(xiàng),之后時(shí)間步從tn推進(jìn)到tn+1。在以上過(guò)程中,空間離散采用基于體心定義的有限體積法,對(duì)流項(xiàng)使用多維限制器(MLP)[23-24],方程中其他項(xiàng)使用中心差分。

由于氣相固相網(wǎng)格重疊,該算法需要對(duì)浸沒(méi)在固相的氣相網(wǎng)格添加限制,具體如下:

ρg=const,ρgΦ=const

(25)

(26)

對(duì)于固相控制方程[式(7)]的空間離散選用有限體積法,時(shí)間上采用后向歐拉法。氣相的時(shí)間推進(jìn)使用龍格庫(kù)塔法[27]。

2.2 界面條件

氣固界面追蹤對(duì)于固體推進(jìn)劑燃燒尤為重要。圖2網(wǎng)格展示浸沒(méi)邊界法[3]中的網(wǎng)格區(qū)分。粗線表示氣固界面,在粗線以下為固相,以上為氣相。如果一個(gè)網(wǎng)格的體心在氣相/固相中,那么就定義這個(gè)網(wǎng)格為氣相/固相網(wǎng)格。賦予源項(xiàng)的邊界網(wǎng)格定義在氣固相之間,并且體心處于氣相內(nèi)的網(wǎng)格。網(wǎng)格中藍(lán)色、棕色和紅色分別代表氣相、固相和浸沒(méi)被賦予源項(xiàng)的邊界網(wǎng)格。

(27)

式中:ρs為固相密度;V為氣相網(wǎng)格的體積。

圖3 一個(gè)源項(xiàng)網(wǎng)格中通過(guò)兩個(gè)連續(xù)時(shí)間步水平集劃分的界面計(jì)算Fig.3 Interface calculation of in one source term grid identified by two consecutive time steps level

(28)

數(shù)值上,水平集界面追蹤通過(guò)顯示正系數(shù)法求解式(8)[28]。對(duì)于大量非結(jié)構(gòu)化網(wǎng)格,本研究采用快速并行求解器,包括窄帶快速迭代法[29]和補(bǔ)丁窄帶法[30]以分別計(jì)算式(9)中的距離變換及界面演化。

燃燒界面上的溫度十分關(guān)鍵,本研究通過(guò)非迭代的方式對(duì)式(13)進(jìn)行離散求解[31]。在每個(gè)時(shí)間步,代表氣固界面的等值面φ等于0被顯式插值構(gòu)建,在固相的非結(jié)構(gòu)化網(wǎng)格中進(jìn)行這種幾何操作會(huì)生成三角形和四邊形單元,如圖4所示,每個(gè)單元會(huì)被分配到由各自臨近格點(diǎn)組成的在氣固相之間的控制體積內(nèi)?；诠?jié)點(diǎn)的格林高斯公式由于其在高斜網(wǎng)格情況下具有較好的魯棒性被用來(lái)離散溫度梯度。界面溫度通過(guò)對(duì)包括每個(gè)界面單元的線性稀疏矩陣進(jìn)行非迭代求解獲得。這套求解方式被證明擁有二階精度[31]。

圖4 為計(jì)算界面溫度構(gòu)建的控制體積Fig.4 Control volumes constructed for interface temperature calculation

3 結(jié)果與討論

如圖1所示,本研究選用Miller實(shí)驗(yàn)中基于AP的M21和M24推進(jìn)劑[21]以驗(yàn)證本套算法。與實(shí)驗(yàn)的壓強(qiáng)范圍一致,選取從0.69MPa到20.7MPa的壓強(qiáng)以研究基于仿真計(jì)算的燃速。

對(duì)于氣相,計(jì)算域尺寸在x,y,z的3個(gè)方向上分別為500、500、1000μm。使用均一大小50×50×100的笛卡爾網(wǎng)格,每個(gè)網(wǎng)格的尺寸為10μm。速度邊界條件在計(jì)算域底部被設(shè)置為無(wú)滑移,頂部選用零梯度,側(cè)面選取滑移條件。壓強(qiáng)在出口設(shè)置為定值條件。對(duì)于標(biāo)量(組分和溫度),所有邊界條件均為零梯度。

對(duì)于固相,本研究在尺寸為500μm的計(jì)算域中生成了13.5萬(wàn)個(gè)非結(jié)構(gòu)網(wǎng)格單元。溫度在底部被設(shè)置為定溫300K,側(cè)面為絕熱條件。頂部邊界隨著水平集捕捉的氣固界面進(jìn)行演化并由氣固兩相間的耦合條件確定。

如圖5(a)所示壓強(qiáng)6.89MPa條件下時(shí)間演化到0.01s的M24構(gòu)型溫度場(chǎng)。由于非均勻的局部退移速率,與真實(shí)物理過(guò)程一致,氣固界面變得扭曲。在氣相場(chǎng)中,最高溫度達(dá)到約3200K。從流場(chǎng)中可以發(fā)現(xiàn)溫度與速度的相互作用,圖5(b)為相同設(shè)置下M24構(gòu)型的溫度云圖及流線圖。為了更清楚地展示等值面,只顯示出在x方向上一半的等值面。在氣固界面附近,速度方向在局部與界面垂直。界面的不規(guī)則性也導(dǎo)致了速度方向和大小的強(qiáng)烈變化。遠(yuǎn)離界面,流場(chǎng)受界面的影響逐漸變?nèi)?速度方向近乎垂直向上。在界面附近主要是單質(zhì)推進(jìn)劑火焰,溫度低于遠(yuǎn)處流場(chǎng)中的最終擴(kuò)散火焰。

圖5 6.89MPa下M24構(gòu)型在0.01s時(shí)的三維示意圖Fig.5 3D schematic diagram M24 configuration at 6.89MPa and 0.01s

圖6和圖7單獨(dú)顯示M24構(gòu)型界面及固相區(qū)域的計(jì)算結(jié)構(gòu)。其中圖6為同一時(shí)刻下不同壓強(qiáng)的情況,圖7為高壓下(13.8MPa)不同時(shí)刻的結(jié)果。觀察發(fā)現(xiàn),溫度在AP與黏合劑交界附近更高。在燃燒過(guò)程中,界面會(huì)發(fā)生扭曲變形。并且熱擴(kuò)散層非常薄,較高的界面溫度只能擴(kuò)散到固相內(nèi)很有限的距離中。在低壓情況下(2.07MPa),界面溫度較低,整體反應(yīng)相對(duì)較為緩慢,界面扭曲不明顯。中壓情況(4.83和6.89MPa)界面扭曲較高壓(13.8和20.7MPa)更為劇烈,這表明較高的壓強(qiáng)會(huì)削弱界面上退移速度的差異,產(chǎn)生一個(gè)相對(duì)更平整的界面。與此同時(shí),可以發(fā)現(xiàn)界面溫度隨壓強(qiáng)增長(zhǎng)而提高的趨勢(shì)。

圖6 不同壓強(qiáng)下M24構(gòu)型在0.004s時(shí)的燃燒界面Fig.6 Burning interface of M24 configuration under different pressures at 0.004s

圖7 13.8MPa下不同時(shí)刻M24構(gòu)型的燃燒界面Fig.7 Burning interface of M24 configuration at 13.8MPa under different instants

圖8為M24構(gòu)型和M21構(gòu)型壓強(qiáng)從0.69MPa至20.7MPa下的計(jì)算平均退移速率與實(shí)驗(yàn)值及Bojko等采用Rocfire軟件包計(jì)算結(jié)果的對(duì)比[21,4]。

圖8 M24構(gòu)型和M21構(gòu)型不同壓強(qiáng)下界面平均退移速率的計(jì)算結(jié)果與對(duì)比Fig.8 Calculation and comparison of the mean regression rate of M24 configuration and M21 configuration, respectively under differen pressures

計(jì)算結(jié)果給出了合理的燃速大小并與實(shí)驗(yàn)趨勢(shì)吻合較好。一些可能導(dǎo)致差異的原因包括:數(shù)值誤差和模型誤差,具體來(lái)說(shuō)比如近界面網(wǎng)格分辨率可能不足以解析界面附近很小的溫度流場(chǎng)結(jié)構(gòu),以及化學(xué)動(dòng)力學(xué)的選擇和固體結(jié)構(gòu)刻畫(huà)中的誤差。

圖9展示了M24構(gòu)型在0.007s時(shí)分別在3.45MPa和20.7MPa的情況下,由式(15)定義的氣相R1、R2及R3的化學(xué)反應(yīng)速率在y方向中點(diǎn)切面上的云圖。如之前介紹,反應(yīng)R1、R2及R3分別代表了單質(zhì)推進(jìn)劑火焰,基本擴(kuò)散火焰和最終擴(kuò)散火焰?？梢园l(fā)現(xiàn),在兩種壓強(qiáng)工況下單質(zhì)推進(jìn)劑火焰和初級(jí)擴(kuò)散火焰都非常接近界面。在高壓下,單質(zhì)推進(jìn)劑火焰和初級(jí)擴(kuò)散火焰的化學(xué)反應(yīng)速率峰值更高,說(shuō)明兩種火焰在高壓情況都變得更為劇烈。對(duì)于最終擴(kuò)散火焰可以發(fā)現(xiàn)其反應(yīng)速率在高壓下更高并且火焰整體更加接近于燃面。

圖9 M24構(gòu)型在3.45MPa和20.7MPa的化學(xué)反應(yīng)速率Fig.9 Chemical reaction rates of M24 configuration at 3.45MPa and 20.7MPa

需要說(shuō)明的是對(duì)于當(dāng)前這個(gè)典型的多尺度問(wèn)題,由于火焰離界面的距離非常近(大約幾微米),在高壓下,該距離將會(huì)更小,因此計(jì)算的分辨率很難滿(mǎn)足要求,單質(zhì)推進(jìn)劑火焰及初級(jí)擴(kuò)散火焰這些尺度更小的結(jié)構(gòu)可能沒(méi)有很好地被解析。本研究重點(diǎn)強(qiáng)調(diào)該套新提出的計(jì)算框架的合理性及初步定性定量的驗(yàn)證,更高分辨率的計(jì)算和對(duì)結(jié)果的深入分析將會(huì)在進(jìn)一步的研究中進(jìn)行。

4 結(jié) 論

(1)本研究計(jì)算中氣相計(jì)算采用均一大小的笛卡爾網(wǎng)格保證計(jì)算的穩(wěn)定性,固相使用非結(jié)構(gòu)化網(wǎng)格以方便刻畫(huà)非均質(zhì)粒子結(jié)構(gòu)。兩套固定網(wǎng)格無(wú)需隨時(shí)間步重建,因此數(shù)值效率和準(zhǔn)確性較高。

(2)非均勻退移的氣固界面由水平集方法和浸沒(méi)邊界法追蹤,質(zhì)量、動(dòng)量和能量的守恒關(guān)系通過(guò)對(duì)氣相內(nèi)的浸沒(méi)邊界網(wǎng)格施加源項(xiàng)以保證。

(3)針對(duì)兩種構(gòu)型在不同壓強(qiáng)條件下的燃燒過(guò)程進(jìn)行的討論,該套計(jì)算框架得到的燃速與實(shí)驗(yàn)及已有計(jì)算結(jié)果接近,證明該數(shù)值框架的可行性與較廣的適用性。