部分裂解煤油的旋轉爆震發動機數值模擬

王 丹，周晨初，陳宏玉，嚴 宇，洪 流

(西安航天動力研究所液體火箭發動機技術重點實驗室，陜西 西安 710100)

0 引言

旋轉爆震發動機(rotating detonation engine，RDE)是一種利用爆震波在燃燒室內連續旋轉傳播，并產生穩定推力的發動機[1-5]。圖 1為典型旋轉爆震發動機的燃燒室結構和爆震波三維傳播結構。

圖1 旋轉爆震發動機燃燒室結構Fig.1 Combustion chamber structure of RDE

由于旋轉爆震發動機具有結構簡單，在亞聲速和超聲速入流條件下均可工作，只需要一次起爆的優勢，具備廣闊的應用前景，美國、俄羅斯、法國、波蘭和日本等國家已經開展了大量的研究工作[6-13]。在數值仿真方面，Zhdan等針對RDE結構特點對其流場進行了簡化，采用兩步化學反應模型進行了二維數值模擬[14]。Hishida采用加密網格對RDE開展了二維數值模擬，精確捕捉到了爆震波面上的三波點結構，數值胞格結果表明，受側向膨脹影響，靠近入口處的胞格尺寸較小，而靠斜激波處的胞格尺寸較大[15]。Davidenk等采用氫氧基元反應模型對RDE流場結構進行了二維分析，并對比研究了燃燒室長度對發動機性能的影響[16]。Eude等針對RDE開展了三維數值模擬，研究結果表明受曲率影響，爆震波面后燃燒室內、外壁處壓力不等，此壓力梯度在流場內引起了橫向擾動[17]。

國內方面，國防科技大學劉世杰等人采用二維數值模擬，驗證了RDE的推力矢量調節能力，采用三維數值模擬，對比研究了燃燒室長度、寬度、尾噴管等幾何參數對發動機性能的影響[18]。北京大學王健平等使用兩步化學反應模型對RDE進行了二維數值模擬研究。數值計算獲得了同軸圓管腔中間層曲面上連續爆震的多個循環過程，分析了燃料入射、提前燃燒、爆震波結構和波傳播速度等幾個關鍵問題[19]。南京理工大學馬虎等[20]采用 9 組分 19 步的基元反應模型，對以氫氣、氧氣、氮氣為反應混合物的RDE進行了二維模擬，研究結果表明，采用分段填充方法可以有效地形成沿同一方向傳播的爆震波。

國內外采用氫氣作為燃料的旋轉爆震發動機燃燒室研究取得了眾多進展，但是氫氣存在貯存困難、價格昂貴、危險性高等缺點，采用常溫且無毒的煤油作為燃料將極大地拓展旋轉爆震發動機的應用領域，且能夠實現低成本和長時間貯存。然而，采用液態煤油作為燃料存在起爆困難的問題，本文以煤油—空氣作為推進劑的旋轉爆震發動機為研究對象，提出一種采用預先加熱模式使得液態煤油發生部分裂解，剩余被加熱的煤油及裂解產物再與空氣混合進行點火起爆的方案。通過裂解產生氫氣、乙烯等高活性成分，降低煤油起爆的難度。采用數值模擬的手段開展部分裂解煤油旋轉爆震發動機的仿真分析，分析流場內燃料組分分布對爆震波傳播方向的影響，以及煤油裂解率對發動機性能的影響。

1 數值方法及計算模型

采用部分裂解方式實現點火起爆的煤油—空氣旋轉爆震發動機基本原理為:將發動機設計為兩部分，第一部分為加熱裂解部分，第二部分為爆震燃燒部分。其中第一部分通過加熱的方式使得液態煤油發生汽化，其中小部分在高溫作用下產生裂解，裂解過程將使大分子長碳鏈碳氫化合物變為小分子短鏈碳氫化合物，從而提高了燃料的化學活性。對煤油的加熱方式有多種，其一可采用電加熱模式，即通過將電加熱絲盤繞于煤油管路之上，在傳輸過程中實現汽化和熱裂解。其二可采用預先燃燒的方式，取全部煤油和少量空氣進行富燃燃燒，以此過程作為煤油加熱、裂解的途徑。然而，采用預先燃燒方式所獲得的裂解煤油組分與單純電加熱時組分不同，因此僅針對電加熱模式的產物進行分析。

表1 化學反應組分Tab.1 Chemical reaction group

第二部分的爆震燃燒由第一部分所產生的部分裂解混合物與主燃燒室的空氣摻混，通過起爆管點火起爆。由于有第一部分加熱裂解過程的存在，進入第二部分的燃料活性明顯提升，大幅降低了點火起爆的難度。與此同時，發生部分裂解的混合物與液態煤油相比，組分更加復雜，化學反應過程也更加復雜，需要通過數值模擬的方式著重對第二部分的爆震燃燒過程進行深入分析。將部分裂解煤油作為燃料邊界條件，空氣為氧化劑，模擬燃燒室的起爆過程。利用Fluent求解器，基于理想氣體假設。采用隱式格式求解二維歐拉控制方程，忽略黏性、熱傳導和擴散等輸運效應。用三階MUSCL格式對對流項進行離散，時間項采用具有二階精度的四步龍格—庫塔法。計算采用壓力基、非穩態計算，化學反應模型為有限速率模型，反應速率常數采用Arrhenius公式計算。

由于煤油成分中的長鏈化合物在高溫下不穩定，很容易裂解為碳鏈更短的烷烴、烯烴等，詳細的化學反應機理十分復雜，產物有上百種，NASA Langley提出了17組分，30步反應的簡化機理[21]，在此機理的基礎上，選取主要裂解產物：H2、CH4、C2H4、C3H6，其余為不參與化學反應的N2，忽略占比小于1%的其他產物。以上5種裂解產物加上被加熱汽化的C12H23，共6種成分的混合物作為燃料，空氣(21%O2+79%N2)作為氧化劑進行燃燒，化學反應采用多組分模型，各組分的化學反應均采用總包模型，化學反應模型來源于CHEMKIN軟件[22]，組分列于表1，表中PEF(pre-exponential factor)為指前因子，AE(activation energy)為活化能。

計算針對的燃燒室為環形燃燒室，由于其厚度相對于其周長較小，計算中忽略其厚度的影響(見圖1)，將其沿其母線拉直后得到如圖 2所示的二維計算區域。計算域為尺寸為100 mm×260 mm 的矩形，對應圖 1將圓柱沿x方向剪開，沿y方向拉直(圓柱圓形截面周長為260 mm，圓柱x方向長度L為100 mm)。

計算域的上下邊界為周期性邊界，模擬爆震波在周向不間斷的旋轉過程。計算域的左邊界為入口邊界(入口1和入口2)，為部分裂解煤油，劃分為兩個入口是為了分析初始爆震波的傳播過程，實際上兩入口的物理意義相同，都代表環腔頭部環縫向燃燒室不間斷供應的燃料(對應圖 1左側的環縫)，初始過程中，在入口區域初始化一個高溫高壓區域作為點火起爆源，該區域為0≤x≤10,40≤y≤50。計算域右邊界為壓力出口邊界。

圖2 計算模型Fig.2 Simulation model

2 網格劃分及模型校驗

采用ICEM軟件繪制二維結構化網格，為了保證對爆震波的準確捕捉，估算爆震波波速在2 000 m/s附近，計算時間步長10-7s，每個時間步爆震波通過的距離為0.2 mm，因此將網格尺寸選擇為0.1 mm，保證對爆震波的可靠捕捉，網格數量為260萬。

為了驗證計算模型的可靠性，進行典型爆震管起爆過程計算。爆震管初始結構為直徑30 mm，長度600 mm，采用二維軸對稱模型，取一半進行計算。初始計算時，在爆震管區域填充常溫的、壓力為0.1 MPa的煤油—空氣當量混合比預混燃氣，外界區域填充常溫常壓的空氣。在爆震管入口側設計一個矩形的點火起爆源，位置在x：0～5 mm，y:0～15 mm，點火能量為1 000 J，保證充足的點火能量使得發生爆震。起爆源在0時刻開始工作，持續時間10-5s。計算結果列于圖3。仿真模型能夠準確模擬點火起爆過程，模擬獲得的爆震波波速與理論C—J速度計算值誤差不超過15%，計算模型可靠。

表2 邊界條件Tab.2 Boundary condition

圖3 爆震管仿真Fig.3 Detonation tube simulation

3 起爆過程計算及分析

隨著裂解率的升高，混合物中小分子烯烴的比例不斷提升，化學活性增強，更加容易起爆。考慮到盤管電熱絲加熱能力以及煤油在管路中的停留時間，無法使得煤油發生充分的裂解，因此從工程設計角度考慮，著重針對裂解率20%以下的工況進行計算，研究爆震波的形成和發展規律，獲得可靠的煤油爆震發動機起爆條件。裂解組分含量由CEA(chemical equilibrum with applications)軟件計算獲得。

表2列出了計算的邊界條件，二維計算忽略燃料與氧化劑的摻混過程，進入燃燒室的工作介質是部分裂解煤油與空氣完全摻混的混合氣，入口1及入口2為速度入口，需要給定組分、溫度及初始壓力。根據煤油發生裂解的起始溫度，選擇850 K作為燃料混合物溫度。速度邊界通過自定義udf給定。計算方法為：當邊界處壓力p大于預混燃料噴射總壓p0=0.35 MPa時，流動被阻塞，速度為0，臨界噴射壓力

當p0>p>pcr時，速度通過等熵膨脹計算，即

(2)對動漫行業典型人物進行訪談。典型人物訪談側重對其個人成長與奮斗經歷的了解與認識，并與自己的職業生涯規劃相聯系，撰寫2000字左右的訪談報告。

當p

式中：R為預混燃料氣體常數；T0為總溫；γ為比熱比。

上下邊界采用周期性邊界，模擬發動機的旋轉過程。點火源初始化為溫度為3 000 K，壓力2 MPa。出口為空氣，0.1 MPa壓力邊界，溫度300 K。在旋轉爆震試驗過程中，為了保證爆震波單向傳播，往往采用切向入射的方式。而在二維數值仿真計算中，難以模擬切向入射的過程，高溫高壓的點火區域將會向外產生球形爆震波，向上下兩個方向傳播，當兩道爆震波相遇時，局部壓力高于噴射壓力，入口噴射速度為0，預混燃氣無法進入燃燒室，此時缺乏維持爆震的預混燃料，可能導致火焰熄滅，因此要求在點火處產生的爆震波最好僅沿一個方向傳播。因此在計算中通過分區初始化來模擬該過程，初始時刻在0≤x≤10，50≤y≤260區域填充850 K的預混燃料，其余部分則填充常溫空氣(模擬部分裂解煤油從頭部環縫進入一段時間后充填于燃燒室頭部區域的狀態)，在0≤x≤10，0≤y≤40區域填充空氣而不填充燃料，是為了人為控制爆震波向上傳播(下方無燃料，僅能傳播壓力波而難以形成爆震波)，選擇850 K是對裂解煤油的一個溫度預估，為了對比裂解率的影響而將溫度保持恒定，方便對比分析。圖 4為計算初始化示意圖。

圖5 算例1起爆過程Fig.5 Initiation process of test 1

圖4 初始條件Fig.4 Initial condition

首先針對算例1(煤油裂解率5%)進行分析。圖 5為裂解率5%工況下壓力及C12H23質量分數隨時間的變化情況，其中C12H23的質量分數代表了主要燃料的占比。點火時間為0時刻，25 μs時，形成了正向的爆震波(y軸正方向為正向)，反方向僅有壓力波。50 μs時，由于入口燃料的加入，反向也形成了一道爆震波，隨后兩個爆震波相遇對撞，對撞后的透射激波沿原方向繼續傳播。同時計算區域下方還產生了一個較弱的向上傳播的壓力波，150 μs時向上傳播的壓力波與向下傳播的透射激波對撞，而此刻計算區域上半部分填入的燃料較多，二次對撞后透射的正向壓力波得到了燃料的補充，迅速發展，形成了穩定的爆震波，而其余壓力波由于未得到充分的燃料補充，逐漸耗散，最終在流場內只剩下了正向傳播的一道爆震波，240 μs后該爆震波一直保持穩定自持傳播，460 μs時經過調整過程和3個完整的旋轉爆震周期，爆震波狀態未發生明顯改變。

為了監測爆震波的形成和發展規律，在計算域中取了7個點來分析該點溫度、壓力等參數的變化情況，各點的坐標如表3所示，圖 2中示意畫出了監控點P1～P7的位置。

對于爆震發動機的性能評估，最重要的參數就是發動機比沖和爆震頻率，由于二維計算與實際發動機的情況略有差異，采用P1點壓力的變化曲線來反映爆震頻率，因為爆震波集中在入口側傳播，當調整階段結束后可認為每次在P1點形成一個高壓就意味著爆震波經過一次，即為一個周期。

表3 監控點坐標分布Tab.3 Coordinate distribution of monitoring points 單位:mm

圖6和圖7分別列出了入口附近沿y方向間隔50 mm均布的P1、P2、P3三點的壓力及y方向速度的變化情況。從圖中可以看出，不考慮200 μs前爆震波的形成和調整過程，僅分析之后較穩定的3個周期。可以看出，沿y軸方向，不同點的壓力峰值有所差異，最低2 MPa，最高3 MPa以上，且P1壓力小于P2，P2壓力小于P3，具有一定的規律性。各點y方向的氣流速度峰值基本保持一致，基本在1 000 m/s附近。可見，爆震波在環形燃燒室中傳播過程中，流場內氣流的傳播速度基本保持不變，壓力則呈現規律性變化，存在一個由弱變強再由強變弱的過程，可見爆震波在傳播過程中受到入射預混氣速度及壓力的影響、波后燃燒產物分布的影響，壓力并非時刻保持穩定，存在起伏現象。用P1、P2、P3三點出現壓力峰的時間間隔來估算爆震波傳播的平均速度。例如，P1和P2點y方向距離為50 mm，相繼出現壓力峰的時間間隔為25.5 μs，因此爆震波從P1點傳播到P2點過程中的平均速度為距離與時間的比值，即1 940 m/s。經計算，不同周期不同監測位置爆震波平均傳播速度基本相同。

圖6 算例1監控點P1、P2、P3的壓力變化Fig.6 Pressure history in monitor point P1,P2 and P3 of test 1

圖8列出了沿x軸方向間隔20 mm均布的P4、P5、P6三點壓力變化情況。從圖2中可以看出，雖然P4、P5、P6與P1在同一條直線上，但是其壓力峰值明顯偏低，這是由于穩定后的旋轉爆震波是由爆震波和斜激波組合而成的，具有側向膨脹的性質，爆震波峰面寬度僅有20 mm左右。從穩定后的壓力分布可以看出，P4壓力大于P5，P5壓力大于P6，可見隨著向出口的發展，斜激波逐漸膨脹，壓力降低。

圖7 算例1監控點P1、P2 、P3的y方向速度變化Fig.7 Velocity history of y direction in monitor point P1,P2 and P3 of test 1

圖8 算例1監控點P4、P5、P6的壓力變化Fig.8 Pressure history in monitor point P4,P5 and P6 of test 1

4 裂解率對發動機性能的影響

以上分析均針對算例1(裂解率5%)，由于裂解率不同將影響發動機的燃料組分，為了分析組分對發動機性能的影響，增加算例2(裂解率10%)和算例3(裂解率20%)2種工況開展計算。3種工況僅燃料組分不同，其余計算條件相同。

算例2在50 μs時，兩道面對面傳播的壓力波中間區域發生緩燃，圖 9為此時刻的壓力和組分云圖，與圖 5不同，裂解率升高后預混燃料活性更高，發生了期望外的緩燃，將消耗掉一定量的預混燃料，不利于穩定爆震波的快速形成。壓力波對撞后形成了一些峰值較小的分散的壓力波，在120 μs時甚至看不出存在明顯的壓力波，然而在170 μs時其中一個壓力波得到了燃料的支持并不斷增強，最終形成了正向傳播的爆震波。

圖9 算例2初始過程Fig.9 Initial process of test2

算例3在50 μs時，兩道面對面傳播的壓力波中間區域同樣發生緩燃，對撞后90 μs形成了一個比較明顯向下傳播的壓力波，但是由于預混燃料并未填充進去，壓力又衰減了。在120 μs到240 μs的時間段內流場里并沒有明顯的壓力波，280 μs時產生了一個向下傳播的壓力波，隨著燃料的不斷填充，最終形成了穩定的反向傳播的爆震波。圖10為該工況穩定工作時的壓力及組分云圖。

圖10 算例3穩定傳播結果Fig.10 Stable propagation result of test 3

對比了不同組分的爆震情況可以看出，燃料組分對爆震波的形成和發展過程有非常重要的影響，3種工況除了燃料組分外其余條件完全相同，但是其傳播方向有所差異，形成穩定爆震波的時間也不同。圖 11列出了P1點壓力變化情況。從圖11中可以看出，算例1和算例2第一個穩定的爆震波形成于200 μs，而算例3的第一個穩定的爆震波則形成于400 μs，中間有很長一段時間并未形成爆震波(相當于穩定工作爆震周期的1.5倍)。這是由于裂解率提升后預混燃料活性提升，容易在燃燒室產生緩燃，消耗預混燃料，不利于穩定爆震波的形成，流場內的微小波系需要較長的調整時間，等待足夠的預混燃料填入，方能形成穩定的爆震波，導致調整階段時間延長。調整階段過后3種工況的爆震波頻率基本相同，壓力峰值基本不變。

圖11 不同裂解率工況P1的壓力變化Fig.11 Pressure history in P1 under different crack rate

挑選一個周期來對比燃料組分對發動機比沖的影響，為了方便比較，將3種裂解率放在同一橫坐標起點進行對比。圖12列出了3種裂解率工況在穩定工作后某1個周期P1點的壓力變化情況，從圖中可以看出，3種裂解率情況下，1個周期的時間均為134 μs，旋轉爆震波頻率為7 500 Hz左右，且P1點的壓力峰值均為2.3 MPa左右。

為了比較裂解率對性能的影響，通過作用在出口平面的流體性質來計算比沖。比沖由出口區域上產生的推力與流過的流量之比來確定，具體由下式計算

圖12 一個周期內P1點壓力變化情況(x軸平移為同一零時刻)Fig.12 Pressure history in P1 during one period (translate x axis to the same zero time)

計算中分子分母均在1個周期內對時間進行積分，可以獲得3種工況在1個周期內的平均比沖(基于燃料、氧化劑及稀釋氣體混合物的比沖)。同時為了排除計算誤差及爆震波傳播穩定性的影響，計算了3個周期內的平均比沖，計算結果列于表4。

表4 比沖計算結果Tab.4 Impulse simulation result

可以看出，1個周期內，裂解率5%和10%時，比沖非常接近，然而裂解率20%時，比沖略有提升；3個周期內，不同裂解率下比沖基本相同。可見，裂解率對發動機穩定工況的平均比沖無明顯影響。

5 結論

1)在爆震波形成的初始階段，具有一定的不確定性，該不確定性與流場的初始狀態有緊密的聯系，最終穩定爆震波的傳播方向取決于調整階段燃燒場的波系變化及組分分布。

2)燃料組分對爆震波的形成和發展過程有非常重要的影響，裂解率提升導致預混燃料活性提升，容易在燃燒室產生緩燃從而消耗預混燃料，流場內的微小波系需要較長的調整時間等待足夠的預混燃料填入，方能形成穩定的爆震波，導致形成穩定爆震波的時間延后。

3)調整階段過后在不同裂解率預混燃料的工況下爆震波頻率基本相同，壓力及速度峰值也基本不變。

4)煤油裂解率20%以下時，改變裂解率對發動機穩定工作后的性能無明顯影響。