翁春生 倪曉冬 續(xù)晗

摘 要：旋轉(zhuǎn)爆轟發(fā)動機因其較高的熱循環(huán)效率在航空航天領(lǐng)域備受關(guān)注。 為了研究氣固混合燃料在空氣氛圍中的旋轉(zhuǎn)爆轟特性， 在圓柱坐標(biāo)系中建立了氣固混合相燃料爆轟理論模型， 并基于三維守恒元與求解元方法， 對圓盤形燃燒室內(nèi)的爆轟過程進行三維數(shù)值模擬。 計算獲得了氣固混合相旋轉(zhuǎn)爆轟波穩(wěn)定傳播時的流場結(jié)構(gòu)， 分析了爆轟流場的熱力學(xué)參數(shù)分布特征、 化學(xué)反應(yīng)區(qū)分布以及旋轉(zhuǎn)爆轟波后的波系分布特點。 研究結(jié)果表明， 微米級煤粉在氫氣輔助作用下可快速反應(yīng)并支持旋轉(zhuǎn)爆轟波的穩(wěn)定傳播， 但因盤形燃燒室的扁平結(jié)構(gòu)特征和非預(yù)混噴注方式， 導(dǎo)致爆轟波波陣面呈現(xiàn)為不規(guī)則的曲面。 煤粉和氫氣射流穿透空氣層的能力差異導(dǎo)致氣固混合燃料的化學(xué)反應(yīng)區(qū)發(fā)生分離， 最終以雙反應(yīng)區(qū)的爆轟組織形式支持旋轉(zhuǎn)爆轟波穩(wěn)定傳播。 旋轉(zhuǎn)爆轟波不規(guī)則的曲面結(jié)構(gòu)使其在扁平的燃燒室中發(fā)生多次反射， 進而在爆轟波后有規(guī)律地出現(xiàn)多道反射激波。

關(guān)鍵詞：? 氣固混合相爆轟； 爆轟波； 粉末燃料； 煤粉； 非預(yù)混； 數(shù)值模擬

中圖分類號：TJ760； V231

文獻標(biāo)識碼： A

文章編號：? 1673-5048（2024）02-0060-11

DOI： 10.12132/ISSN.1673-5048.2024.0036

0 引? 言

旋轉(zhuǎn)爆轟發(fā)動機（Rotating detonation engine， RDE）是一種利用爆轟波在環(huán)形［1-2］、 空桶形［3-5］或圓盤形［6-10］等具有圓柱面約束的燃燒室中連續(xù)旋轉(zhuǎn)傳播所產(chǎn)生的高溫高壓氣體產(chǎn)物由發(fā)動機出口高速排出， 進而產(chǎn)生推力的新概念發(fā)動機。 與定壓燃燒相比， 爆轟過程接近于定容燃燒， 因此RDE相較于傳統(tǒng)發(fā)動機具有更高的熱循環(huán)效率。 此外， RDE憑借其一次點火便可連續(xù)穩(wěn)定運行的工作模式， 以及推力穩(wěn)定、 推重比大、 寬范圍的流入速度等優(yōu)勢［11］受到世界各國的廣泛關(guān)注。 在軍機、 導(dǎo)彈、 超高聲速飛行器等領(lǐng)域具有廣闊的應(yīng)用前景。

目前， 旋轉(zhuǎn)爆轟推進技術(shù)的開發(fā)應(yīng)用主要以氣體和液體燃料為主， 而隨著軍事應(yīng)用場景的不斷擴大， 以含能材料等固體推進劑為燃料的旋轉(zhuǎn)爆轟發(fā)動機的研制被提上日程。 相較于氣體和液體燃料， 固體燃料具有儲存方便、 能量密度高、 抗沖擊性能好等優(yōu)點。 將固體燃料制成粉末狀應(yīng)用于RDE， 使其繼承固體燃料優(yōu)點的同時， 又兼?zhèn)淞藲怏w、 液體燃料燃燒效率高、 流量可調(diào)等優(yōu)點。 因此， 近幾年作為航天及軍事強國的俄羅斯和美國， 先后開展了固體粉末旋轉(zhuǎn)爆轟發(fā)動機的實驗研究。 俄羅斯的Bykovskii等人［12-16］率先在不同規(guī)格大小的圓盤形燃燒室中對各類型煤粉的旋轉(zhuǎn)爆轟過程進行了相關(guān)實驗研究， 在不嚴(yán)格控制燃料流量的情況下驗證了各類型煤粉末旋轉(zhuǎn)爆轟的可行性。 美國的Dunn等人［17-19］隨后對煤粉在環(huán)形燃燒室中的旋轉(zhuǎn)爆轟過程進行了實驗研究， 并通過理論及數(shù)值研究進一步分析了煤粉爆轟的可行性。 國內(nèi)方面， 與Dunn等人同期， Xu等人［6-8］在控制燃料流量穩(wěn)定的情況下于直徑150 mm的盤形燃燒室中對煤粉、 鋁粉等燃料的爆轟特性進行了實驗研究， 并對爆轟波波頭高度、 強度、 傳播速度以及發(fā)動機推力和燃料比沖進行了系統(tǒng)的分析。

粉末旋轉(zhuǎn)爆轟發(fā)動機作為近幾年剛開始研究的課題， 研究初期， 由于固體燃料的旋轉(zhuǎn)爆轟特性尚不明確， 許多關(guān)鍵技術(shù)需要探索與驗證。 因此， 選用一種取材方便、 成本較低、 安全可靠的實驗材料非常重要。 煤炭， 作為生活中常見的固體燃料， 與旋轉(zhuǎn)爆轟中經(jīng)常用到的氣態(tài)甲烷、 液態(tài)煤油等同為碳?xì)浠剂希? 在世界范圍內(nèi)分布廣、 儲量大、 開采成本低， 且燃燒產(chǎn)物以無害的CO2和H2O為主， 因此成為現(xiàn)階段粉末旋轉(zhuǎn)爆轟發(fā)動機的理想燃料。 但煤的化學(xué)反應(yīng)活性較低， 研究其旋轉(zhuǎn)爆轟特性需借助于氫氣等活性氣體的輔助。 因此， 上述研究工作中關(guān)于煤粉的爆轟特性研究均是在氫氣等活性氣體輔助的情況下進行的。 然而， 因為研究時間較短， 且實驗測試手段的局限， 煤粉-氫氣-空氣旋轉(zhuǎn)爆轟燃燒室中大量的物理參數(shù)無法觀測和記錄， 因此對于燃燒室中的流場參數(shù)分布及變化規(guī)律缺乏直觀的認(rèn)識和理解， 使得粉末旋轉(zhuǎn)爆轟的研究留下巨大的空白。 為此， 需借助于數(shù)值模擬的方式進一步展開研究。 Zhu等人［20-21］對環(huán)形燃燒室中煤粉、 鋁粉的旋轉(zhuǎn)爆轟過程進行了二維模擬研究， 圍繞爆轟波前顆粒三角區(qū)與空氣三角區(qū)不完全重疊的現(xiàn)象進行了一系列研究。 然而， 現(xiàn)實中的旋轉(zhuǎn)爆轟燃燒室均為三維非預(yù)混噴注結(jié)構(gòu)， 煤粉-氫氣-空氣旋轉(zhuǎn)爆轟波的組織結(jié)構(gòu)也呈現(xiàn)出三維的結(jié)構(gòu)特點。 為了更好地了解煤粉-氫氣-空氣旋轉(zhuǎn)爆轟波的傳播特性， 且對前期實驗研究工作［6， 8］進行補充， 本文建立了三維煤粉-氫氣-空氣混合相爆轟理論模型， 并基于時空守恒元與求解元（The Space-Time Conservation Element and Solution Element， CE/SE）方法［22-27］， 模擬了盤形燃燒室中非預(yù)混噴注條件下煤粉-氫氣-空氣旋轉(zhuǎn)爆轟波的傳播過程， 進一步分析了煤粉氫氣混合燃料爆轟波的三維結(jié)構(gòu)特征、 化學(xué)反應(yīng)區(qū)分布特征以及熱力學(xué)參數(shù)、 波系結(jié)構(gòu)在不同方向上的分布特點。

1 混合相爆轟理論模型與計算方法

1.1 混合相爆轟理論模型

實踐中， 圓盤形RDE與離心式壓氣機和徑流渦輪表現(xiàn)出良好的匹配特性， 同時也具有便于可視化觀測和布置監(jiān)測點的優(yōu)勢， 因此對處于初始研究階段的粉末旋轉(zhuǎn)爆轟， 其多數(shù)實驗研究工作是在盤形燃燒室中進行的。 為了銜接并進一步完善前期研究工作， 本文的數(shù)值研究工作也圍繞盤形燃燒室開展。 為此， 在考慮了圓盤形燃燒室主要結(jié)構(gòu)特征的基礎(chǔ)上， 建立了扁平的煤粉-氫氣-空氣非預(yù)混圓盤形計算域， 如圖1所示。 計算域上下端面分別為燃燒室上壁面和下壁面； 外圓柱面為噴注面板， 噴注面板上陣列分布有燃料及空氣噴注入口， 位于上方的入口陣列為燃料入口， 位于下方的入口陣列為空氣入口， 其余部分為固壁（為便于描述， 噴注面板的固壁部分稱為噴注面）； 內(nèi)圓柱面為出口。 該計算域出口半徑為R1、 噴注面半徑為R2、 厚度為L。

1.1.1 基本假設(shè)

由于煤粉-氫氣-空氣旋轉(zhuǎn)爆轟過程既包含氣相燃料的爆轟， 也包含固相粉末燃料的爆轟， 整個過程極其復(fù)雜， 為了方便研究與計算， 需對該過程進行簡化。

在開展粉末爆轟特性研究的工作中， 通常以微米或亞微米級的超細(xì)粉為研究對象。? 由于顆粒尺寸極小且在粉末發(fā)動機流量供給條件下顆粒數(shù)目極大， 因此， 本文數(shù)值模型主要用于描述大量超細(xì)固體顆粒在氣體中的宏觀運動以及氣固兩相之間的相互作用過程。 旋轉(zhuǎn)爆轟發(fā)動機工作時， 固相燃料顆粒群與氣相燃料充分混合后經(jīng)小孔噴入燃燒室， 整個過程顆粒群懸浮于氣相環(huán)境， 且顆粒數(shù)量足夠大， 因此在宏觀上可以使用連續(xù)體假設(shè)來將固相顆粒群近似為“準(zhǔn)流體”， 而單個顆粒的各物理特征則均勻一致。 數(shù)值計算時， 氣相與固相顆粒群可被處理為共享空間區(qū)域的雙流體， 且用流體動力學(xué)方程來描述這兩相的運動， 而質(zhì)量、 動量與能量的源項則被用來描述兩相流體之間的熱-機械-化學(xué)相互作用。 由于固相的真實密度遠(yuǎn)大于氣相， RDE中固相燃料的體積分?jǐn)?shù)始終不超過總體的1%， 因此， 固相體積的變化對流場的影響十分有限。 此外， 爆轟過程是瞬態(tài)變化過程， 爆轟波的傳播速度通常達(dá)到幾千米/秒， 因此燃燒室流場與壁面的傳熱作用對爆轟波的傳播影響較小， 分析爆轟流場特征時可忽略其與壁面的熱交換過程。

目前， 大量的旋轉(zhuǎn)爆轟發(fā)動機數(shù)值研究工作基于Euler方程開展。 Bruce等人［28］指出， 粘性對非定常的激波運動過程幾乎沒有影響。 林玲等人［26］和Oran等人［29］則通過對比Euler方程和Navier-Stokes方程計算得到的爆轟流場， 發(fā)現(xiàn)粘性對爆轟波的傳播及結(jié)構(gòu)特征也沒有明顯影響。 因此， Euler方程完全滿足對RDE中旋轉(zhuǎn)爆轟波傳播過程及流場結(jié)構(gòu)特征研究的需求。

基于上述分析和論證， 本文在使用Euler方程的基礎(chǔ)上總結(jié)出如下幾條假設(shè)： （1）爆轟過程中燃燒室壁面絕熱； （2）粉末燃料視為固體顆粒群， 顆粒及其反應(yīng)后產(chǎn)生的氣體在每個網(wǎng)格內(nèi)與氣體氧化劑瞬間混合均勻； （3）顆粒在整個爆轟過程中均保持球狀， 且固體顆粒間互不影響， 顆粒內(nèi)部溫度均勻分布； （4）忽略氣固混合物中顆粒體積對流場的影響。

1.1.2 氣固混合相爆轟控制方程

基于圓盤形燃燒室的結(jié)構(gòu)特征與前文假設(shè)， 在柱坐標(biāo)系中建立如下控制方程組：

Ut+Fr+Grθ+Hz=R－Fr （1）

式中： U為守恒項； F， G， H為對流項； R為源項。

U=ρgYO2ρgYH2ρgYCO2ρgYH2OρgYN2ρgugρgvgρgwgρgEgρsρsusρsvsρswsρsEs ；

F=ρgugYO2ρgugYH2ρgugYCO2ρgugYH2OρgugYN2ρgu2g+pρgugvgρgugwg（ρgEg+p）ugρsusρsu2sρsusvsρsuswsusρsEs；

G=ρgvgYO2ρgvgYH2ρgvgYCO2ρgvgYH2OρgvgYN2ρgugvgρgv2g+pρgvgwg（ρgEg+p）vgρsvsρsusvsρsv2sρsvswsvsρsEs； H=ρgwgYO2ρgwgYH2ρgwgYCO2ρgwgYH2OρgwgYN2ρgugwgρgvgwgρgw2g+p（ρgEg+p）wgρswsρsuswsρsvswsρsw2swsρsEs；

R=－8m·s/3－16ω·－2ω·11m·s/318ω·0m·sus－Frm·svs－Fθm·sws－Fz－Qconv－（Frus+Fθvs+Fzws）+m·s（Es+p/ρs）－m·s－m·sus+Fr－m·svs+Fθ－m·sws+FzQconv+（Frus+Fθvs+Fzws）－m·s（Es+p/ρs）。 其中： 下標(biāo)g， s分別表示氣相和固相； 變量ρ， u， v， w， E， p分別表示密度、 徑向速度、 周向速度、 軸向速度、 總能和氣相壓力； YO2， YH2， YCO2， YH2O， YN2分別表示氣相中氧氣、 氫氣、 二氧化碳、 水蒸汽和氮氣的質(zhì)量分?jǐn)?shù)， 且YO2+YH2+YCO2+YH2O+YN2=1； m·s為固相燃料質(zhì)量消耗速率， ω·為氫氣反應(yīng)摩爾濃度消耗速率； Fr， Fθ， Fz分別為固體顆粒群在徑向、 周向和軸向上受到的氣相作用力； Qconv表示氣相與固體顆粒群之間的對流換熱。

能量方程中， 氣相和固相總能分別滿足如下關(guān)系：

Eg=h－pρ+12（u2g+v2g+w2g） （2）

Es=cvTs+12（u2s+v2s+w2s） （3）

式中： h為氣相的比焓； cv為固相的比熱容； Ts為固相溫度。

本文所用固相燃料為煤粉， 鑒于前期所用煤粉固定碳含量為95%以上， 故在計算時主要考慮碳顆粒群的反應(yīng)。 目前的研究中， 碳顆粒的燃燒主要受到表面動力學(xué)或物質(zhì)擴散的限制。 擴散限制的反應(yīng)速率依賴于氧化劑的濃度， 而動力學(xué)限制的反應(yīng)速率則取決于顆粒表面處的反應(yīng)速率。 因此， 鑒于Salvadori等人［19］的分析與驗證， 碳顆粒群的總反應(yīng)速率被定義為

m·s=4Ncπr2cpO21/ks+1/kd（4）

式中： Nc為顆粒數(shù)密度； rc為顆粒半徑； pO2為氧氣分壓； ks和kd分別代表煤粉顆粒的動力學(xué)反應(yīng)速率和擴散反應(yīng)速率（單位為kg·m-2·s-1·Pa-1）。

ks=0.86exp－1.495×108RgTs（5）

kd=9.72DrcRg（Tg+Ts） （6）

式中： Tg為氣相溫度； Rg為氣體常數(shù)； D為擴散系數(shù)。

D=4.24×10－7T 3/2g（1/MO2+1/MCO2）1/2p（V1/3O2+V1/3CO2）2（7）

式中： MO2， MCO2分別為氧氣和二氧化碳的摩爾質(zhì)量； VO2， VCO2分別為氧氣和二氧化碳的摩爾體積。

本文所用氣相燃料為氫氣， 其主要作用為流化并輔助微米級煤粉爆轟， 因此， 這里采用一步反應(yīng)［20］來描述其化學(xué)反應(yīng)過程， 相應(yīng)的反應(yīng)速率由Arrhenius公式計算：

ω·=A［H2］m［O2］nexp（－EaRgTg）（8）

式中： A為化學(xué)反應(yīng)指前因子， 取值為9.87×108； ［H2］， ［O2］分別為氫氣、 氧氣的摩爾濃度； m， n為反應(yīng)級數(shù)， 這里均取值為1； Ea為活化能， 取值為5×107J/kmol。

氣固兩相流場中， 氣體與固體顆粒間存在拖曳力［24］， 結(jié)合粉末顆粒的濃度， 可以在宏觀上得到氣相在各方向上對固相的作用力為

Fr=12Ncπr2cCDρg|Vg－Vs|（ug－us）

Fθ=12Ncπr2cCDρg|Vg－Vs|vg－vs）

Fz=12Ncπr2cCDρg|Vg－Vs|（wg－ws）? （9）

式中： Vg， Vs分別為氣相和固相的速度矢量； CD為阻力系數(shù)。

CD與氣固兩相之間的相對雷諾數(shù)Re存在如下關(guān)系［24］：

CD=24Re1+16Re2/3（10）

Re=2ρgrc|Vg－Vs|μg（11）

式中： μg為氣體粘度。

兩相流動過程中， 由于溫差和速度差等因素的影響， 氣相與固相之間通過對流換熱傳遞的熱量為［24］

Qconv=2πrcNcλgNu（Tg－Ts） （12）

式中： λg為氣體的導(dǎo)熱系數(shù)； Nu為努塞爾數(shù)， 且

Nu=2+0.459Re0.55Pr0.33（13）

此處， 普朗特數(shù)Pr=0.72。

1.2 計算方法

爆轟過程的一個主要特征是有激波等強間斷的傳播。 為此， 爆轟過程的數(shù)值模擬對強間斷的捕捉能力具有較高的要求。 故本文采用適用于強間斷問題求解的時空守恒元與求解元方法（CE/SE方法）來對上述控制方程組進行求解。

CE/SE方法于1995年被Chang正式提出以來已發(fā)展30年左右， 它是基于空間通量與時間通量的守恒性原理推導(dǎo)而來， 目前被大量應(yīng)用于爆轟等強間斷問題的數(shù)值計算［1， 23-25］。 該方法作為一種全新的守恒型方程計算格式， 與傳統(tǒng)數(shù)值方法（如有限差分法和有限體積法）不同， 它將時間視為第四個維度， 與空間中的三個維度進行統(tǒng)一劃分處理。 在求解計算時， 該方法將整個時間-空間區(qū)域劃分為若干個求解元。 在每個求解元內(nèi)， 假設(shè)流場的變量是連續(xù)的， 并可用泰勒級數(shù)展開， 穿過相鄰求解元的邊界， 流場的變量可以是不連續(xù)的。 對于每個網(wǎng)格點對應(yīng)的守恒元， 空間-時間密度矢量的積分通量是守恒的。

CE/SE方法計算精度高， 且無需采用算子分裂即可保證局部和整體的通量守恒， 無需Riemann分解即可實現(xiàn)激波捕捉， 在保證守恒性的同時， 能夠準(zhǔn)確地獲取流場的激波和其他接觸間斷現(xiàn)象。

1.2.1 計算格式

在采用CE/SE方法對控制方程組進行求解時， 首先將整個計算域劃分為交替的六面體網(wǎng)格， 網(wǎng)格點（i， j， k）位于六面體的中心位置， 且i， j， k分別是柱坐標(biāo)系中的r， θ， z方向的網(wǎng)格點數(shù)。 通常， 在一個完整的時間步Δt內(nèi)， 間隔分布的網(wǎng)格點分兩步計算完成， 即前半個時間步求解一半網(wǎng)格點， 后半個時間步計算另一半網(wǎng)格點。 三維CE/SE方法守恒元與求解元的定義涉及到四維結(jié)構(gòu)的構(gòu)建， 其計算格式的推導(dǎo)過程詳見文獻［27］， 這里直接給出其在柱坐標(biāo)系下的推導(dǎo)結(jié)果：

（U）li， j， k－Δt2（R－F/r）li， j， k=16 U+Δr4Ur－Δt4（Gθ/r+Hz）+? 3ΔtΔrF+Δr6Fr+Δt4Ftl－1/2i－1/2，? j，? k+

U－Δr4Ur－Δt4（Gθ/r+Hz）－3ΔtΔrF－Δr6Fr+Δt4Ftl－1/2i+1/2，? j，? k+U+Δθ4Uθ－Δt4（Fr+Hz）+

3ΔtΔθG+Δθ6Gθ/r+Δt4Gtl－1/2i， j－1/2， k+U－Δθ4Uθ－Δt4（Fr+Hz）－ 3ΔtΔθG－Δθ6Gθ/r+Δt4Gtl－1/2i， j+1/2， k+U+Δz4Uz－Δt4（Fr+Gθ/r）+3ΔtΔzH+Δz6Hz+Δt4Htl－1/2i，? j，? k－1/2+

U－Δz4Uz－Δt4（Fr+Gθ/r）－

3ΔtΔzH－Δz6Hz+Δt4Htl－1/2i，? j，? k+1/2 （14）

式中： l為時間點。 Fr， Gθ， Hz， Ft， Gt， Ht與Ur， Uθ， Uz之間存在如下關(guān)系：

Fr=Fr=FU Ur=AUr

Gθ=Gθ=GU Uθ=BUθ

Hz=Hz=HU Uz=CUz （15）

Ft=A（R－F/r－AUr－BUθ/r－CUz）

Gt=B（R－F/r－AUr－BUθ/r－CUz）

Ht=C（R－F/r－AUr－BUθ/r－CUz）? （16）

其中， A=FU， B=GU和C=HU均為雅克比系數(shù)矩陣。 因此， 在使用該計算格式時只需再提供Ur， Uθ和Uz的算法即可。 因旋轉(zhuǎn)爆轟波等強間斷的存在， 本文在計算Ur， Uθ和Uz時宜采用加權(quán)函數(shù)法求解。 為此， 定義

（U±r）li，? j， k=±［（U′）li±1/2， j， k－（U）li， j， k］Δr/2（U±θ）li， j， k=±［（U′）li， j±1/2， k－（U）li， j， k］Δθ/2（U±z）li， j， k=±［（U′）li， j， k±1/2－（U）li， j， k］Δz/2? （17）

式中： （U′）li±1/2， j， k=U+Δt2Utl－1/2i±1/2， j， k；

（U′）li， j±1/2， k=U+Δt2Utl－1/2i， j±1/2， k；

（U′）li， j， k±1/2=U+Δt2Utl－1/2i， j， k±1/2。

進一步定義U（κ）r， U（κ）θ和U（κ）z（κ = 1， 2， 3， …， 8）分別對應(yīng)于（U±r）li， j， k， （U±θ）li， j， k和（U±z）li， j， k的8種組合， 即

U（1）r=（U+r）li， j， k， U（1）θ=（U+θ）li， j， k， U（1）z=（U+z）li， j， k；

U（2）r=（U+r）li， j， k， U（2）θ=（U+θ）li， j， k， U（2）z=（U－z）li， j， k；

U（3）r=（U+r）li， j， k， U（3）θ=（U－θ）li， j， k， U（3）z=（U+z）li， j， k；

U（4）r=（U+r）li， j， k， U（4）θ=（U－θ）li， j， k， U（4）z=（U－z）li， j， k；

U（5）r=（U－r）li， j， k， U（5）θ=（U+θ）li， j， k， U（5）z=（U+z）li， j， k；

U（6）r=（U－r）li， j， k， U（6）θ=（U+θ）li， j， k， U（6）z=（U－z）li， j， k；

U（7）r=（U－r）li， j， k， U（7）θ=（U－θ）li， j， k， U（7）z=（U+z）li， j， k；

U（8）r=（U－r）li， j， k， U（8）θ=（U－θ）li， j， k， U（8）z=（U－z）li， j， k。

基于此， Ur， Uθ和Uz可通過如下加權(quán)函數(shù)求解：

（Ur）li， j， k=0， 若ψ（ζ）=0 （ζ=1， 2， 3，…， 8）∑8κ=1［（W（κ））αU（κ）r］/∑8κ=1（W（κ））α， 其他

（Uθ）li， j， k=0， 若ψ（ζ）=0 （ζ=1， 2， 3， …， 8）

∑8κ=1［（W（κ））αU（κ）θ］/∑8κ=1（W（κ））α， 其他

（Uz）li， j， k=0， 若ψ（ζ）=0 （ζ=1， 2， 3， …， 8）

∑8κ=1［（W（κ））αU（κ）z］/∑8κ=1（W（κ））α， 其他? （18）

式中： α是一個可調(diào)常數(shù)， 通常取為1或者2；? W（κ）=∏8ζ=1， ζ≠κψ（ζ） ， ψ（ζ）=（U（ζ）r）2+（U（ζ）θ）2+（U（ζ）z）2。

1.2.2 時間步長及源項的處理

數(shù)值計算過程中， CE/SE方法時間步長Δt的選取通常應(yīng)滿足穩(wěn)定性條件：

Δt=minεΔr|u|+c， εrΔθ|v|+c， εΔz|w|+c （19）

式中： c為聲速； ε為小于1的數(shù)， 弱波的情況下可以取大一些， 而對于強波則可以取小一些。

由于爆轟過程中， 化學(xué)反應(yīng)的特征時間遠(yuǎn)小于物理特征時間， 因此其源項是剛性的。 這里采用四階龍格庫塔法對剛性源項進行求解， 且求解的時間步長按如下關(guān)系選?。?/p>

Δtr－k=Δt/Nr－k （20）

式中： 迭代次數(shù)Nr－k通常取5 ～ 20。

本文剛性問題的處理步驟為： 在不考慮源項的情況下， 采用CE/SE方法先求解（U）li， j， k， 然后將所得結(jié)果作為初值， 再對常微分方程組dU/dt=R－F/r進行求解。

1.3 初始條件與邊界條件

初始條件： 初始時刻燃燒室內(nèi)填充有常溫常壓下總當(dāng)量比為1.0的煤粉-氫氣-空氣可燃混合物， 其中氫氣/空氣當(dāng)量比為0.7， 煤粉/空氣當(dāng)量比為0.3。 圖2所示為點火區(qū)域分布圖， 為快速獲得穩(wěn)定傳播的旋轉(zhuǎn)爆轟波， 在圖2中角度為15°的紅色區(qū)域設(shè)置高溫高壓高速的周向氣流作為點火條件， 并在指定時間內(nèi)將紅色區(qū)域左側(cè)的計算域縫合面設(shè)置為固壁。

邊界條件： 如圖1所示， 上壁面、 下壁面和噴注面為固壁邊界； 噴注面上層入口陣列為煤粉-氫氣混合燃料的壓力入口邊界， 下層入口陣列為空氣的壓力入口邊界； 作為出口的內(nèi)圓柱面為無反射自由邊界［3］； 圓周縫合處在初始爆轟波傳播270°后由初始的固壁邊界轉(zhuǎn)變?yōu)橹芷谶吔纭?對于壓力入口邊界， 參數(shù)的計算與設(shè)置可參考文獻［10］。 同時， 為便于燃料與空氣的混合， 并增加燃料在燃燒室中的停留時間， 參照前期實驗研究工作［6-8］， 燃料噴孔處煤粉-氫氣混合燃料的噴注速度向空氣噴孔一側(cè)偏轉(zhuǎn)， 而空氣噴孔處氣流的噴注速度則向圓周切線方向發(fā)生偏轉(zhuǎn)。 對于出口無反射自由邊界， 參數(shù)的具體計算與設(shè)置可參考文獻［30］， 即出流為亞聲速時， 采用外展出流邊界條件， 而出流為超聲速時， 采用自由出流邊界條件。

2 計算結(jié)果與分析

網(wǎng)格大小的選擇既要考慮計算精度， 同時也需要兼顧計算成本。 在進行數(shù)值模擬之前， 首先進行網(wǎng)格無關(guān)性驗證。 本文分別在Mesh-1（1 040×124×44）， Mesh-2 （720×90×32）和Mesh-3（560×56×20） 三種網(wǎng)格規(guī)模下進行爆轟流場的數(shù)值計算， 得到各網(wǎng)格尺寸條件下同一時刻燃燒室外圓壓力沿周向的分布情況如圖3所示。 可知以上三種網(wǎng)格均能有效捕捉到爆轟波的強間斷面， 且各網(wǎng)格規(guī)模下爆轟波壓力差異均不超過2%。 因此， 本文選取720×90×32網(wǎng)格開展數(shù)值研究。

本文對以氫氣-煤粉氣固混合物為燃料、 空氣為氧化劑的盤形旋轉(zhuǎn)爆轟發(fā)動機進行三維流場數(shù)值計算。 前期的實驗研究在文獻［8］中已經(jīng)證實， 在半徑R2=75 mm， 厚度L=16 mm的圓盤形燃燒室中， 由260 g/s的空氣、 5.3 g/s的氫氣和6.7 g/s的微米級多孔球狀無煙煤粉末組成的可燃混合物， 經(jīng)高能點火可實現(xiàn)單個旋轉(zhuǎn)爆轟波的穩(wěn)定傳播。 為進一步探究實驗現(xiàn)象的本質(zhì)， 同時驗證數(shù)值計算的準(zhǔn)確性， 本文數(shù)值計算均按上述條件進行設(shè)置， 根據(jù)實驗中所用粉末材料的粒徑檢測結(jié)果， 計算中的球狀顆粒粒徑設(shè)為5 μm。

圖4所示為旋轉(zhuǎn)爆轟波在盤形燃燒室穩(wěn)定傳播階段， 距離噴注面35 mm（即距燃燒室中軸線R=40 mm）處下壁面監(jiān)測點的實驗與計算壓力時程曲線。 可以看出， 兩曲線的波形及幅值大小基本一致， 證明文中所建計算模型可準(zhǔn)確模擬非預(yù)混噴注條件下煤粉-氫氣-空氣在盤形燃燒室中的旋轉(zhuǎn)爆轟過程。

2.1 氣固混合相旋轉(zhuǎn)爆轟結(jié)構(gòu)特征

從高溫高壓點火開始， 經(jīng)過數(shù)十圈的傳播， 盤形燃燒室中最終形成一個沿逆時針傳播的旋轉(zhuǎn)爆轟波。 圖5所示為下壁面上緊貼噴注面（距燃燒室中軸線R=75 mm）處監(jiān)測點壓力及爆轟波傳播頻率時程曲線， 可以看到該測點處所測壓力峰值周期性振蕩且平均值為5.54 MPa。 顯然， 該值遠(yuǎn)大于圖4監(jiān)測點（R=40 mm）處壓力峰值。 這主要是因為該點處于爆轟波波頭可覆蓋范圍， 且位于兩壁面交匯處， 壁面對激波的反射作用使其壓力成倍增加， 而圖4監(jiān)測點距離噴注面35 mm， 距爆轟波波頭相對較遠(yuǎn)， 所測壓力峰值為斜激波壓力， 壓力峰值相對較低。 然而， 各測點處壓力峰值大小的差異并不影響爆轟波傳播頻率的獲取。 通過求取壓力時程曲線上爆轟波或斜激波兩相鄰峰值時間差的倒數(shù)， 可近似得到爆轟波的瞬時傳播頻率。 圖5顯示， 該爆轟波的平均傳播頻率為4.211 kHz， 若以測點所在圓周（R=75 mm）計算， 則該旋轉(zhuǎn)爆轟波的傳播速度為1 984 m/s。 圖6所示為8.25 ms時刻煤粉-氫氣混合相燃料在空氣中穩(wěn)定爆轟時盤形燃燒室的壓力及溫度分布圖， 可以看到爆轟波在噴注面的約束作用下沿著圓盤外緣穩(wěn)定傳播。

圖6（a） 顯示， 非預(yù)混噴注條件下， 旋轉(zhuǎn)爆轟波波陣面因壁面約束和非預(yù)混噴注等原因呈現(xiàn)為不規(guī)則的曲面， 從噴注面方向看爆轟波陣面主要由上下兩個弧面交匯而成， 交匯點位于盤形燃燒室下半層， 即靠近空氣噴孔一側(cè)。 噴注面上兩弧面交匯處壓力最大， 噴注面下半層整體壓力明顯高于上半層。 這一方面是因為激波交匯導(dǎo)致壓力升高， 另一方面是交匯點位于噴注流量較大的空氣一側(cè)， 該區(qū)域密度較大， 相應(yīng)地， 爆轟波經(jīng)過時強度較大。 此外， 爆轟波沿圓周傳播過程中， 爆轟波與上壁面之間存在一層很薄的已燃?xì)怏w層， 由于該薄層無燃料及氧化劑的噴注， 因此爆轟波在該薄層的燃燒產(chǎn)物中產(chǎn)生一透射激波隨爆轟波波陣面一起向前傳播。 而沿徑向方向， 不規(guī)則的爆轟波波陣面在爆轟波與燃燒室出口之間又透射產(chǎn)生一系列強度較弱的激波。

圖6（b）顯示， 爆轟波前的區(qū)域在低溫燃料及空氣的高速噴注作用下， 燃燒室外緣靠近噴注面板的中間區(qū)域形成低溫燃料填充區(qū)域， 而在靠近上下壁面的噴注盲區(qū)， 因上一循環(huán)留下的高溫產(chǎn)物駐留而保持相對較高的溫度。 噴注面附近， 爆轟波掃過后， 在稀疏膨脹波的作用下， 燃燒室壓力迅速下降， 噴孔在壓差的作用下向燃燒室中噴注低溫燃料及空氣， 使得波后噴孔附近混合氣體的溫度相較于貼近上下壁面的噴注盲區(qū)明顯下降。 從噴注面方向看， 旋轉(zhuǎn)爆轟波附近的溫度場呈現(xiàn)出上下兩個“V”型高溫區(qū)分布的特點。

2.2 氣固混合相爆轟化學(xué)反應(yīng)分布特征

為了進一步探究煤粉-氫氣-空氣旋轉(zhuǎn)爆轟波的組織形式， 圖7給出了煤粉-氫氣-空氣旋轉(zhuǎn)爆轟過程中各燃料組分的反應(yīng)放熱速率透視圖。 從圖中可以清晰地看到， 混合燃料反應(yīng)放熱速率的空間分布特征與旋轉(zhuǎn)爆轟波的結(jié)構(gòu)高度一致， 表明此處混合燃料（煤粉-氫氣）與空氣反應(yīng)放熱支持了旋轉(zhuǎn)爆轟波的持續(xù)傳播。

與實驗噴注情形一致， 混合均勻的氣固混合相燃料經(jīng)燃料噴孔陣列向空氣的噴注區(qū)域噴注。 但因為氣相燃料與固相粉末燃料在物理性質(zhì)上的巨大差異（主要表現(xiàn)為固體顆粒相較于氣體具有更大的密度與慣性）， 混合燃料經(jīng)空氣的“過濾”作用后在燃燒室中會呈現(xiàn)出分層分布的特點。 進而導(dǎo)致各燃料的反應(yīng)放熱區(qū)域也呈現(xiàn)出分層分布的特點。 因此， 從圖7中各燃料的放熱速率空間分布來看， 煤粉-氫氣混合相燃料旋轉(zhuǎn)爆轟波在燃料與空氣分開噴注的情形下依靠雙反應(yīng)區(qū)的爆轟組織形式自持傳播。 其中， 煤粉與空氣中的氧氣主要在燃燒室貼近下壁面的區(qū)域反應(yīng)放熱以支持旋轉(zhuǎn)爆轟波的傳播， 而氫氣與氧氣發(fā)生反應(yīng)的區(qū)域則主要位于煤粉反應(yīng)區(qū)上方。 圖7中煤粉化學(xué)反應(yīng)區(qū)分布與圖6中爆轟波處壓力、 溫度空間分布的高度耦合， 證明了煤粉參與爆轟并支持混合相旋轉(zhuǎn)爆轟的連續(xù)傳播。

2.3 流場三維分布特征

煤粉-氫氣混合燃料與空氣在軸向上分開噴注的方式， 使得燃料與氧化劑的摻混在軸向上存在較大差異， 進而導(dǎo)致流場剖面在軸向上發(fā)生變化， 而軸向上燃料噴孔向空氣噴孔一側(cè)傾斜噴注的方式也使得最佳摻混區(qū)域向下壁面偏移。 因此， 根據(jù)上一節(jié)中燃燒室壓力溫度在爆轟波處的軸向分布特征， 在軸向上選取若干切面及特征點展開分析與討論。

圖8所示為盤形燃燒室軸向不同切面上的壓力與溫度分布圖。 因噴注摻混區(qū)域整體向下壁面偏移， 徑向噴注速度的不均勻分布導(dǎo)致旋轉(zhuǎn)爆轟波在上、 下壁面處的激波前沿在徑向偏轉(zhuǎn)趨勢上存在明顯差異。 中心切面處因受上下波陣面交匯的影響， 此處波陣面在徑向上發(fā)生明顯轉(zhuǎn)折， 且轉(zhuǎn)折處因激波匯聚而壓力最大。 與中心切面不同， 上下壁面處的壓力最大區(qū)域緊貼噴注面， 這主要是圓柱噴注面的收斂壓縮作用導(dǎo)致的結(jié)果。 因燃料空氣摻混區(qū)域向下壁面偏移， 所以爆轟波經(jīng)過時靠近下壁面處化學(xué)反應(yīng)更為劇烈， 可以看到， 爆轟波經(jīng)過時下壁面處壓力梯度較大， 上壁面處壓力梯度較小。

圖8（b）中各切面的溫度分布顯示， 盤形燃燒室中的高溫燃?xì)獬市郎u狀流出燃燒室。 這一方面是因為噴注空氣的切向分速度引起的旋轉(zhuǎn)流動； 另一方面是旋轉(zhuǎn)爆轟波傳播過程中， 沿周向高速運動的斜激波壓縮氣體， 增加了氣體的切向速度分量。 從軸向中心切面的溫度分布可以看到， 低溫燃料與空氣的噴注摻混過程在徑向上形成一環(huán)形的新鮮燃料預(yù)混區(qū)， 該區(qū)域的穩(wěn)定存在保證了盤形燃燒室中旋轉(zhuǎn)爆轟波的穩(wěn)定傳播。 因空氣噴孔所在切面距離下壁面較近， 在高速的低溫空氣對流和新鮮的煤粉顆粒群吸熱作用下， 下壁面上靠近噴注面的區(qū)域與整個切面其他區(qū)域相比， 在爆轟波到來之前也存在一溫度較低的環(huán)形帶。 爆轟波掃過時， 經(jīng)過一段時間預(yù)熱的煤粉與附近氧氣快速反應(yīng)釋放熱量， 支持旋轉(zhuǎn)爆轟波的傳播。 相比之下， 燃燒室上半?yún)^(qū)因?qū)儆谛迈r來流噴注的盲區(qū)， 且距離噴孔相對較遠(yuǎn)， 在上壁面處并未形成明顯的低溫環(huán)形帶。

圖9所示為徑向坐標(biāo)R=75 mm時， 下壁面、 中心切面和上壁面處監(jiān)測點的壓力與溫度時程曲線。 在摻混區(qū)向燃燒室下半?yún)^(qū)偏移以及圓柱形噴注面對流場的收斂壓縮作用下，? 下壁面處監(jiān)測點平均壓力峰值最高，? 達(dá)到5.54 MPa。 其次是距離摻混區(qū)較近的中心切面處測點， 平均壓力峰值為2.42 MPa； 相比之下， 離摻混區(qū)最遠(yuǎn)的上壁面處， 監(jiān)測點平均壓力峰值最低， 為2.18 MPa。 因處于噴注死角， 上壁面貼近噴注面處監(jiān)測點的溫度值基本維持在1? 500 K以上，? 且在爆轟波傳播至對應(yīng)相位角時， 平均溫度峰值可達(dá)到3 063 K。 與之相比， 距離噴注區(qū)域更近的下壁面處， 貼近噴注面的位置雖位于噴注盲區(qū)， 但強烈的對流作用使得監(jiān)測點溫度在每一個旋轉(zhuǎn)周期內(nèi)都有一段時間維持在1 000 K以下。 噴注區(qū)域的中心切面處， 因新鮮低溫氧化劑與燃料的噴注， 監(jiān)測點溫度在一半以上的時間內(nèi)維持在500 K以下， 僅在爆轟波掃過時溫度短暫超過1 500 K， 平均溫度峰值可達(dá)到2 464 K。

盤形燃燒室扁平的結(jié)構(gòu)特征， 使其在軸向上對爆轟波的傳播產(chǎn)生了明顯的約束。 為了進一步探究上下壁面對爆轟波傳播的影響以及爆轟波后壓力場的分布特點， 圖10給出了距圓盤形燃燒室中心軸線不同距離的圓周切面壓力分布展開圖。 從各切面的壓力分布間斷來看， 因非預(yù)混噴注條件下爆轟波波陣面為不規(guī)則的曲面結(jié)構(gòu)， 使得其與上下壁面之間形成不同的傾斜角度， 進一步導(dǎo)致了爆轟波波后區(qū)域激波在上下壁面處的交替反射。 圓周切面的半徑越小， 距離出口越近， 距噴注面越遠(yuǎn)， 壓力場受到爆轟波以及噴注面收斂壓縮作用的影響越小。 因此， 緊貼噴注面的切面處因反射激波而產(chǎn)生的壓力變化最為明顯。 而隨著圓周切面半徑的減小， 同時考慮到側(cè)向膨脹作用的影響， 因爆轟波傳播引起的斜激波及波后一系列反射激波逐漸減弱， 切面上壓力變化逐漸減緩。

從圖8中壓力和溫度沿徑向的分布來看， 新鮮燃料及氧化劑的噴注混合區(qū)域高度有限， 導(dǎo)致爆轟波在徑向上的高度有所限制。 隨著與噴注面距離的增加， 燃燒室中的強壓縮擾動由爆轟波轉(zhuǎn)變?yōu)樾奔げǎ?其中， 斜激波以“弓形激波系”的形式呈現(xiàn)。 圖8中各軸向圓盤切面的壓力分布顯示， 因爆轟波陣面在不同切面上形狀有所差異， 爆轟波與燃燒室出口之間產(chǎn)生了一系列強度較弱的弓形激波。 這些弓形激波均匯聚于不規(guī)則的爆轟波陣面， 而發(fā)散于燃燒室出口。 因燃燒室為圓盤形， 激波沿圓周切向傳播相同的距離， 距燃燒室中軸線越近的位置相位差越大， 因此由爆轟引起的“弓形激波系”前沿在相位上始終領(lǐng)先于爆轟波前沿。 這些規(guī)律在圖例取值范圍較小的圖10中得到進一步體現(xiàn)。 從圖10中可以看到， 隨著與噴注面距離的增加， 爆轟波波陣面受非預(yù)混噴注的影響變小， 波陣面逐漸趨于平整。 并且從R=65 mm切面處開始， 各弓形激波之間的相位差隨徑向坐標(biāo)R的減小快速增加。

圖11所示為下壁面上“弓形激波系”的主要分布區(qū)域內(nèi)， 同一相位角度不同徑向位置的壓力時程曲線。 從圖中各監(jiān)測點捕捉的第一道弓形激波壓力擾動的上升沿開始時間點， 可以看出“弓形激波系”率先覆蓋距燃燒室中軸線較近的位置。 而對以上各時間點作差， 從R=40 mm處開始， 到R=65 mm， R每增加5 mm， 測點捕捉到壓力上升沿的時間間隔Δt依次為4.2 μs， 4.2 μs， 4.1 μs， 4.1 μs和4.2 μs。 因此， 在角速度基本不變的情況下， “弓形激波系”前沿的相位角隨徑向坐標(biāo)R的增大線性減小。 結(jié)合2.1節(jié)中旋轉(zhuǎn)爆轟波的傳播頻率， 可以計算得到徑向坐標(biāo)每間隔5 mm， “弓形激波系”前沿相位差Δθ=2πf·Δt= 0.108 ～ 0.111 rad。 因弓形激波強度較弱， 前文提到的因上下壁面反射而產(chǎn)生的反射激波在該區(qū)域內(nèi)也相對較弱， 圖11中各反射激波峰值壓力的大小對比再次證明， 隨著徑向坐標(biāo)R的減小， 反射激波逐漸減弱， 在R小于55 mm時， 弓形激波后無明顯的反射激波存在。

3 結(jié)? 論

為給固體燃料旋轉(zhuǎn)爆轟發(fā)動機的研制提供技術(shù)支持與驗證， 在國內(nèi)外關(guān)于煤粉-氫氣-空氣旋轉(zhuǎn)爆轟實驗研究持續(xù)開展的背景下， 本文基于雙流體模型， 在圓柱坐標(biāo)系下建立了三維氣固混合相燃料爆轟理論模型， 采用CE/SE方法計算求解了盤形燃燒室中煤粉-氫氣-空氣混合相旋轉(zhuǎn)爆轟流場； 證實了微米級煤粉可支持旋轉(zhuǎn)爆轟波的傳播， 同時分析總結(jié)了粉末旋轉(zhuǎn)爆轟的三維分布特征。 計算結(jié)果表明：

（1） 煤粉-氫氣-空氣旋轉(zhuǎn)爆轟波波陣面因壁面約束和非預(yù)混噴注等原因呈現(xiàn)為不規(guī)則的曲面， 從噴注面方向看， 波陣面主要由上下兩個弧面交匯而成， 噴注面上兩弧面交匯處壓力最大。 因噴注方式采用非預(yù)混小孔噴注， 噴注面與上下壁面之間形成噴注盲區(qū)， 旋轉(zhuǎn)爆轟波附近的溫度場呈現(xiàn)出上下兩個“V”型高溫區(qū)分布的特點。

（2） 因煤粉顆粒與氫氣穿透空氣區(qū)域的能力存在較大差異， 煤粉氫氣混合流在穿越空氣噴注區(qū)域時分層分布， 煤粉-氫氣混合相燃料旋轉(zhuǎn)爆轟波在非預(yù)混噴注的情形下依靠雙反應(yīng)區(qū)的爆轟組織形式自持傳播。 其中， 煤粉與空氣中的氧氣主要在燃燒室貼近下壁面的區(qū)域反應(yīng)放熱以支持旋轉(zhuǎn)爆轟波的傳播， 而氫氣與氧氣發(fā)生反應(yīng)的區(qū)域則主要位于煤粉反應(yīng)區(qū)上方。

（3） 因非預(yù)混噴注條件下爆轟波波陣面為不規(guī)則的曲面結(jié)構(gòu)， 使得其與上下壁面之間形成不同的傾斜角度， 進一步導(dǎo)致了爆轟波波后區(qū)域出現(xiàn)多重反射激波。

（4） 盤形燃燒室中爆轟波與燃燒室出口之間產(chǎn)生了一系列強度較弱的弓形激波， 這些弓形激波均匯聚于不規(guī)則的爆轟波陣面， 發(fā)散于燃燒室出口， 形成“弓形激波系”； 而該“弓形激波系”前沿的相位角坐標(biāo)隨燃燒室中徑向坐標(biāo)的增加線性減小。

參考文獻：

［1］ 李寶星， 翁春生. 氣體與液體兩相連續(xù)旋轉(zhuǎn)爆轟發(fā)動機爆轟波傳播特性三維數(shù)值模擬研究［J］. 兵工學(xué)報， 2017， 38（7）： 1358-1367.

Li Baoxing， Weng Chunsheng. Three-Dimensional Numerical Simulation on the Propagation Characteristics of Detonation Wave in Gas-Liquid Two-Phase Continuous Rotating Detonation Engine［J］. Acta Armamentarii， 2017， 38（7）： 1358-1367.（in Chinese）

［2］ Meng H L， Xiao Q， Feng W K， et al. Air-Breathing Rotating Detonation Fueled by Liquid Kerosene in Cavity-Based Annular Combustor［J］. Aerospace Science and Technology， 2022， 122： 107407.

［3］ Tang X M， Wang J P， Shao Y T. Three-Dimensional Numerical Investigations of the Rotating Detonation Engine with a Hollow Combustor［J］. Combustion and Flame， 2015， 162（4）： 997-1008.

［4］ Lin W， Zhou J， Liu S J， et al. An Experimental Study on CH4/O2 Continuously Rotating Detonation Wave in a Hollow Combustion Chamber［J］. Experimental Thermal and Fluid Science， 2015， 62： 122-130.

［5］ Liu S J， Huang S Y， Peng H Y， et al. Characteristics of Methane-Air Continuous Rotating Detonation Wave in Hollow Chambers with Different Diameters［J］. Acta Astronautica， 2021， 183： 1-10.

［6］ Xu H， Ni X D， Su X J， et al. Experimental Investigation on the Application of the Coal Powder as Fuel in a Rotating Detonation Combustor［J］. Applied Thermal Engineering， 2022， 213： 118642.

［7］ 續(xù)晗， 羅永晨， 倪曉冬， 等. 鋁粉燃料連續(xù)旋轉(zhuǎn)爆轟發(fā)動機工作特性［J］. 兵工學(xué)報， 2022， 43（5）： 1046-1053.

Xu Han， Luo Yongchen， Ni Xiaodong， et al. Operating Characte-ristics of Aluminum Powder Rotating Detonation Engine［J］. Acta Armamentarii， 2022， 43（5）： 1046-1053.（in Chinese）

［8］ Ni X D， Xu H， Su X J， et al. Effects of Different Physical Properties of Anthracite Powder Fuel on Detonation Characteristics of a Rotating Detonation Engine［J］. Physics of Fluids， 2023， 35（5）： 053325.

［9］ 夏鎮(zhèn)娟， 張義寧， 馬虎， 等. 點火位置對圓盤結(jié)構(gòu)下旋轉(zhuǎn)爆震波起爆過程的影響［J］. 推進技術(shù)， 2021， 42（4）： 805-814.

Xia Zhenjuan， Zhang Yining， Ma Hu， et al. Effects of Ignition Position on Initiation Process of Rotating Detonation Wave in Plane-Radial Structure［J］. Journal of Propulsion Technology， 2021， 42（4）： 805-814.（in Chinese）

［10］ 夏鎮(zhèn)娟， 武曉松， 馬虎， 等. 圓盤結(jié)構(gòu)下旋轉(zhuǎn)爆震波的二維數(shù)值研究［J］. 推進技術(shù)， 2017， 38（6）： 1409-1418.

Xia Zhenjuan， Wu Xiaosong， Ma Hu， et al. Two-Dimensional Numerical Simulation of Rotating Detonation Wave in Plane-Radial Structure［J］. Journal of Propulsion Technology， 2017， 38（6）： 1409-1418.（in Chinese）

［11］ Ma J Z， Luan M Y， Xia Z J， et al. Recent Progress， Development Trends， and Consideration of Continuous Detonation Engines［J］. AIAA Journal， 2020， 58（12）： 4976-5035.

［12］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Continuous and Pulsed Detonation of a Coal-Air Mixture［J］. Doklady Phy-sics， 2010， 55（3）： 142-144.

［13］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Detonation Combustion of Coal［J］. Combustion， Explosion， and Shock Waves， 2012， 48（2）： 203-208.

［14］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Continuous Spin Detonation of a Coal-Air Mixture in a Flow-Type Plane-Radial Combustor［J］. Combustion， Explosion， and Shock Waves， 2013， 49（6）： 705-711.

［15］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Detonation Burning of Anthracite and Lignite Particles in a Flow-Type Radial Combustor［J］. Combustion， Explosion， and Shock Waves， 2016， 52（6）： 703-712.

［16］ Bykovskii F A， Vedernikov E F， Zholobov Y A. Detonation Combustion of Lignite with Titanium Dioxide and Water Additives in Air［J］. Combustion， Explosion， and Shock Waves， 2017， 53（4）： 453-460.

［17］ Dunn I B， Malik V， Flores W， et al. Experimental and Theoretical Analysis of Carbon Driven Detonation Waves in a Heteroge-neously Premixed Rotating Detonation Engine［J］. Fuel， 2021， 302： 121128.

［18］ Dunn I B， Flores W， Morales A， et al. Carbon-Based Multi-Phase Rotating Detonation Engine［J］. Journal of Energy Resources Technology， 2022， 144（4）： 042101.

［19］ Salvadori M， Dunn I B， Sosa J， et al. Numerical Investigation of Shock-Induced Combustion of Coal-H2-Air Mixtures in a Unwrapped Non-Premixed Detonation Channel［C］∥AIAA Scitech 2020 Forum， 2020.

［20］ Zhu W C， Wang Y H， Wang J P. Flow Field of a Rotating Detonation Engine Fueled by Carbon［J］. Physics of Fluids， 2022， 34（7）： 073311.

［21］ 李世全， 楊帆， 王宇輝， 等. 當(dāng)量比對鋁粉/空氣旋轉(zhuǎn)爆轟發(fā)動機流場影響的數(shù)值模擬［J/OL］. 航空動力學(xué)報， doi： 10.13224/j.cnki.jasp. 20210560.

Li Shiquan， Yang Fan， Wang Yuhui， et al. Numerical Simulations on Effect of Equivalence Ratio on Flow Field in Aluminum/Air Rotating Detonation Engines［J/OL］. Journal of Aerospace Power， doi： org/10.13224/j.cnki.jasp.20210560. （in Chinese）

［22］ Chang S C. The Method of Space-Time Conservation Element and Solution Element-A New Approach for Solving the Naver-Stokes and Euler Equations［J］. Journal of computational physics， 1995， 119（2）： 295-324.

［23］ Wang G， Zhang D L， Liu K X， et al. An Improved CE/SE Scheme for Numerical Simulation of Gaseous and Two-Phase Detonations［J］. Computers and Fluids， 2010， 39（1）： 168-177.

［24］ Zhang Z J， Wen C Y， Liu Y F， et al. Application of CE/SE Method to Gas-Particle Two-Phase Detonations under an Eulerian-Lagrangian Framework［J］. Journal of Computational Physics， 2019， 394： 18-40.

［25］ 馬丹花， 翁春生. 二維守恒元和求解元方法在兩相爆轟流場計算中的應(yīng)用［J］. 燃燒科學(xué)與技術(shù)， 2010， 16（1）： 85-91.

Ma Danhua， Weng Chunsheng. Application of Two-Dimensional CE/SE Method Calculation of Two-Phase Detonation Flow Field［J］. Journal of Combustion Science and Technology， 2010， 16（1）： 85-91.（in Chinese）

［26］ 林玲， 翁春生. 等離子體射流點火對燃燒轉(zhuǎn)爆轟影響的二維數(shù)值計算［J］. 兵工學(xué)報， 2014， 35（9）： 1428-1435.

Lin Ling， Weng Chunsheng. Two-Dimensional Numerical Calculation for the Influence of Plasma Jet Ignition on Deflagration-to-Deto-nation Transition［J］. Acta Armamentarii， 2014， 35（9）： 1428-1435.（in Chinese）

［27］ 翁春生， 王浩. 計算內(nèi)彈道學(xué)［M］. 北京： 國防工業(yè)出版社， 2006： 290-306.

Weng Chunsheng， Wang Hao. Computational Interior Ballistics［M］. Beijing： National Defense Industry Press， 2006： 290-306.（in Chinese）

［28］ Bruce P J K， Babinsky H. Unsteady Shock Wave Dynamics［J］. Journal of Fluid Mechanics， 2008， 603： 463-473.

［29］ Oran E S， Weber J W， Stefaniw E I， et al. A Numerical Study of a Two-Dimensional H2-O2-Ar Detonation Using a Detailed Chemical Reaction Model［J］. Combustion and Flame， 1998， 113（1/2）： 147–163.

［30］ 王健平， 姚松柏. 連續(xù)爆轟發(fā)動機原理與技術(shù)［M］. 北京： 科學(xué)出版社， 2018： 68-72.

Wang Jianping， Yao Songbai. Principle And Technology Of Continuous Detonation Engine［M］. Beijing： Science Press， 2018： 68-72.（in Chinese）

Numerical Study on the Two-Phase Detonation Flow Field of

Coal Powder-Hydrogen in a Rotating Detonation Engine

Weng Chunsheng， Ni Xiaodong*， Xu Han

（National Key Laboratory of Transient Physics， Nanjing University of Science and Technology， Nanjing 210094， China）

Abstract： Rotating detonation engines （RDEs） have attracted significant attention in the aerospace field due to their high thermal cycle efficiency. This study explores the characteristics of rotating detonation in gas-solid mixed fuel within an air atmosphere. A theoretical detonation model for gas-solid mixed-phase fuel is developed in a cylindrical coordinate system. Using the three-dimensional conservation element and solution element （CE/SE） method， numerical simulations of the detonation process in a disc combustor are conducted in three dimensions. The flow field structure during the stable propagation of gas-solid mixed-phase rotating detonation waves is calculated. Subsequently， the analysis is conducted on the distribution characteristics of thermodynamic parameters in the detonation flow field， the distribution of chemical reaction zone， and the characteristics of wave system following the rotating detonation wave. The research results demonstrate that micron-sized coal powder can undergo a rapid reaction and support the stable propagation of rotating detonation waves with the assistance of hydrogen. However， due to the flat structural characteristics of the disc combustor and the non-premixed injection mode， the wave front of detonation wave appears as an irregular surface. The difference in the ability of coal powder and hydrogen jet to penetrate the air layer leads to the separation of the chemical reaction zone of the gas-solid mixed fuel， ultimately supporting the stable propagation of rotating detonation wave in the form of a dual reaction zone. The irregular surface structure of the rotating detonation wave leads to multiple reflections in the flat combustor， consequently causing the regular appearance of multiple reflected shock waves after the detonation wave.

Key words： gas-solid mixed-phase detonation； detonation wave； powder fuel； coal powder； non-premixed； numerical simulation