溫度擾動對ODW結構影響的數值模擬

陳楠， Sudip BHATTRAI， 唐豪，*

(1. 南京航空航天大學 能源與動力學院， 南京 210016；2. 尼泊爾工程學院 機械工程部， 帕坦 44700)

爆震燃燒是一種效率極高的燃燒形式，其傳播速度可達每秒千米級，具有很廣闊的應用前景。目前實現爆震推進的方式主要有3種：脈沖爆震發動機、斜爆震發動機以及旋轉爆震發動機。其中工作方式最為簡單的是斜爆震發動機，其基本原理為以一道駐定于燃燒室內的斜激波(OSW)來誘發超聲速來流混氣燃燒，從而得到穩定自持的斜爆震波(ODW)，本文所選用計算域正是一種簡易的ODW燃燒室。

20世紀80年代后期，對駐定ODW的數值模擬工作得到了快速發展，由于計算能力的限制，多是采用較為簡單的模型進行模擬。Fujiwara等[1-2]用迎風顯式 TVD 格式對駐定ODW進行二維和三維數值模擬；Cambier等[3]用帶化學反應的 TVD 格式對燃燒室的駐定ODW問題進行了數值計算；Dennis[4]和Yungster[5]等都用不同的方法對駐定ODW問題進行研究，并且將計算出的ODW波形與實驗得到的照片進行比較。隨著計算方法與運算能力的大幅提高，數值模擬已經成為主要的研究手段。Powers和Gonthier[6]采用二步不可逆反應模擬了弱驅斜爆震結構。Pimentel等[7]采用詳細反應動力學模型研究了OSW與斜爆震轉變區域的不同結構，以及斜劈面角度變化產生的Taylor波對爆震結構的影響，得到的3種不同類型的轉換形式:①OSW角逐漸地增大，平滑轉換；②ODW形成，OSW角陡峭轉換；③ODW形成，OSW角平滑轉換。Berlyand等[8]研究了爆震波在斜面上的反射出現的穩態與非穩態2種現象。國內學者對ODW的研究起步較晚，但近年來發展迅速，Teng等[9]研究了不同過驅度下OSW向ODW轉變的模式特征，得到了誘導區轉變形式的臨界過驅度，并且指出OSW和ODW面的夾角可以用作判斷激波向爆震波轉化的有效參考[10]，隨后的研究模擬了活化能和斜楔面對波面胞格結構的影響[11-12]，發現存在2種橫波結構：單向橫波和雙向橫波，并探討了橫波由單向到雙向的轉變規律，在最新發表的文章中[13]，運用CVC(Constant Volume Combustion)理論研究了來流動力學參數對誘導區特性的影響。Gui等[14]研究了ODW面胞格結構，同樣發現單向橫波和雙向橫波結構，并深入分析了這種微觀結構的周期性碰撞現象。Liu等[15]運用R-H(Rankine-Hugoniot)方法分析了ODW向上游傳播的誘因是誘導區尾部的壓力超過ODW的脫體壓力。此外Liu等[16]利用提高部分區域反應進程變量的方法觸發ODW，得到了存在著單向橫波的完全耦合ODW和存在著雙向橫波的部分耦合ODW，并且進一步驗證了橫波對ODW的穩定自持有重要影響。同時Liu等[16]在爆震波不同區域引入局部溫度擾動以檢驗ODW的穩定性，發現穩定傳播的ODW在傳播過程中穩定性會不斷得到增強，但其并未深入討論溫度擾動的傳播形式。鑒于此，本文可以被視為對其工作的一項延伸和拓展。

雖然前人對ODW的研究已經進行了大量的工作，但其內部結構的發展規律并不完全清楚，并且絕大部份數值研究都在均勻來流條件下進行。發動機實際工作中，由于燃料摻混不均勻、蒸發不完善、進氣道出口壓力波動等因素，都可能在燃燒室入口處產生來流物性參數的變化。這些擾動都會對ODW結構產生很大的影響，其中又以溫度擾動最為顯著。本文出發點在于，通過對燃燒室內的駐定ODW施加來自上游的溫度擾動，觀察擾動在燃燒室內的傳播特點及擾動對ODW結構的影響。本文側重于ODW大致結構的變化情況，并不對爆震波內部微小結構(胞格結構等)進行分析。

1 計算區域與網格劃分方法

圖1 計算區域及ODW結構示意圖Fig.1 Schematic of ODW structure and computational domain

計算區域為一個210 mm×120 mm的燃燒室，如圖1所示,RCJ ODW為反射CJ ODW。反應物為混合均勻的氫氣和空氣，當量比為1，進口來流馬赫數Ma=5.95，溫度T=796 K，壓強P=35.04 kPa，斜楔角為16°，圖1給出了計算域區域及ODW的基本結構，其總體架構與Ghorbanian和Sterling[17]所描述的架構相近。左邊界為進口端，進口端邊界上來流參數設定為初始值不變。

為保證數值結果的收斂，本文選取了3種網格劃分的方法，分別為756×400、854×600和1 132×800，3種網格下的庫朗特數分別取0.2、0.4和0.8，最大時間步長取5×10-9s。圖2給出了t=0.09 ms時，3種網格下的溫度分布云圖，能夠看出，3種網格下都得到了能夠穩定自持的ODW，隨著網格密度加大，ODW的微觀結構逐漸顯現，圖2(b)中自ODW面向下游產生橫波，在圖2(c)這種橫波變得更加清晰。

圖3顯示了斜楔面上的壓力變化情況，可看出3種網格下CJ ODW后的壓力在理論壓力(PCJ)附近上下波動，這種波動是自斜爆震產生的橫波作用的結果，隨著網格密度的增大，數值波動也更為劇烈。由圖3還可看出，854×600和1 132×800下壓力的變化趨勢基本吻合，而756×400下的壓力偏離較大，為進一步比較3種網格下的數值結果，取平行于斜楔面一直線(y=4 mm)，得到溫度和OH-質量濃度沿該直線的分布情況，如圖4所示。從圖4可看出，854×600和1 132×800下數值結果的變化趨勢已基本吻合，數值波動也十分接近。

本文目的在于研究溫度擾動引起的ODW結構的變化，并不要求網格能夠精細描述ODW的微觀結構。綜上所述，854×600下能夠得到完整的ODW結構，數值結果與1 132×800下的數值結果接近，并且波后壓力能夠符合理論預測，故是合適的網格劃分。

圖2 t=0.09 ms時3種網格劃分下溫度分布云圖Fig.2 Temperature distribution contours of three kinds of mesh density at t=0.09 ms

圖3 t=0.09 ms時3種網格劃分下沿斜楔面的壓力分布(Ma=5.95，T=796 K，P=35.04 kPa)Fig.3 Pressure distribution along ramp at t=0.09 ms in three kinds of mesh density (Ma=5.95，T=796 K and P=35.04 kPa)

圖4 t=0.09 ms時3種網格劃分下y=4 mm上的溫度和OH-質量濃度分布(Ma=5.95，T=796 K，P=35.04 kPa)Fig.4 Distribution of temperature and OH- mass density at t=0.09 ms measured at y=4 mm in three kinds of mesh density (Ma=5.95，T=796 K and P=35.04 kPa)

2 控制方程與數值方法

與bhattraI和Hao[18]相同，本文選取了二維非穩態無黏可壓歐拉方程來構建數學模型。采用開源CFD軟件OpenFOAM進行運算。控制方程如下：

連續方程

(1)

動量方程

(2)

能量方程

(3)

組份方程

(4)

理想氣體狀態方程

(5)

k=ATβexp(-Ea/(RT))

(6)

式中：Ea為各分步反應的活化能;A和β為分步反應的2個變量因子，它們決定各分步反應的頻率因子ATβ。

本文采用的化學動力模型為氫氣/空氣9組份[O2,H,OH,O,H2,H2O,N2,HO2,H2O2] 19步可逆反應模型[19]，該反應模型在模擬燃燒過程的同時，能夠得到的反應中各原子團的變化情況，這可以與實驗結果進行直接對比，其可靠性和有效性已被大量研究所驗證[20-22]，因此在爆震燃燒研究中被廣泛應用。具體反應機制如表1所示，反應速率及摩爾生成速率由CHEMKIN計算而得。

表1 19步氫氣/空氣反應模型[19]Table 1 19-step H2-Air reaction model[19]

注：第3體效率:a-fH2O=12.0,fH2=2.5；b-fH2O=12.0,fH2=0.73;c-fH2O=14.0,fH2=1.3。

本文數值方法采用AUSM+通量離散格式[18]。每步迭代中，采用黎曼外推法從各網格值推導各網格中心點的左右面通量值，在標量梯度的計算中，采用壓力梯度項用以限制數據變化范圍，化學反應采用解耦的方式進行處理。進口邊界處來流參數保持不變，混氣為均勻摻混的氫氣和空氣，當量比為1，氮氣質量分數為0.743 5。出口邊界處采用零梯度條件，上下壁面處采用有滑移的絕熱邊界條件。

3 結果與討論

3.1 ODW的基本結構

當t=0.1 ms時，燃燒室內部形成了能夠穩定自持的ODW結構，如圖5所示，從流場的密度云圖能夠看出完整的ODW結構。氣流撞擊楔面產生OSW,激波誘導混氣反應形成一段誘導區，在誘導區末端混氣開始劇烈反應生成強烈的燃燒波，燃燒波不斷疊加形成ODW，Ghorbanian和Sterling[17]通過數值模擬證明對于在過驅狀態下，誘導區末端形成的爆震波為CJ ODW，圖5中CJ爆震波的角度為36.8°，與CJ理論值(βCJ)37.4°很接近。OSW與CJ ODW匯聚形成駐定的ODW，3種波相交于一點，該點被稱為三波點，自三波點處向下游延伸產生一條滑移線，在滑移線上下兩側壓力相等，但溫度、密度和速度等其他參數會發生突躍變化。由于三波點是高溫高壓的區域，通常CJ ODW在三波點處會發生反射，形成RCJ ODW。

誘導區后的反應區被滑移線分為2個區域，滑移線以上到爆震波間為爆震區(detonation zone)，其以下到壁面之間則為爆燃區(deflagration zone)。CJ ODW與其后形成的膨脹波(expansion waves)耦合，形成向誘導區彎曲的低溫區域。RCJ ODW在反應區內經壁面和ODW反射傳播并不斷衰減，Choi等[23]指出，RCJ ODW的傳播可能與ODW面的胞格結構變化有關，也可能影響爆震區內的橫波傳播。自ODW波面形成的橫波向下游傳播，導致了反應區內微觀結構的不斷變化，Dupré等[24]指出這種不穩定性對爆震波的自持傳播起到本質作用。

圖5 t=0.1 ms時密度分布云圖Fig.5 Density distribution contours at t=0.1 ms

3.2 溫降擾動在燃燒室內的傳播

對于幾何結構不可變的燃燒室，進口邊界處任一參數動可能會引起其他物性參數的瞬態波動。在t=0.1 ms時，在流場進口邊界引入一擾動，即溫度突然下降100 K(溫降擾動)，為保證來流參數物理上的合理性，考慮到燃燒室的實際工作情況，保持進口來流速度不變，壓力、密度和來流馬赫數等進口參數做出相應調整。擾動在進口邊界均勻分布并保持不變，觀察溫度擾動在燃燒室內傳播的過程。R-H分析指出，在弱過驅爆震區間內ODW能夠穩定自持，擾動傳播過程中ODW波結構必然會經歷調整，圖6顯示了自擾動傳入至ODW重新調整到穩定狀態的完整過程。

將擾動前圖6(a)和擾動流出后流場圖6(f)對比，易發現爆震區結構發生了劇烈變化。圖7給出了擾動前后，三波點附近的密度變化情況。由圖7(a)可看出，擾動前爆震區結構變化較平滑，在接近三波點的區域，爆震區胞格結構不明顯，直到流場下游，才逐漸顯現。然而擾動傳播過后，爆震區內結構不再平滑，胞格結構在結構新三波點處不遠便十分清晰，火焰波陣面呈現明顯周期性波動，并從波陣面生成面向下游的多層橫波結構，如圖7(b)所示。

圖6 溫降擾動傳入至ODW重新調整到穩定狀態的時序圖Fig.6 Sequence charts for variation of ODW structure from introducing temperature drop disturbance to retrieving stabilization

圖7 三波點附近擾動前后的密度分布灰度圖Fig.7 Grayscale images for density nearby triple point before and after disturbance

為了研究擾動在誘導區、爆燃區和爆震區中擾動的傳播形式，在3個區域分別選取了3個觀測點，如圖6(f)中①、②、③所示。由圖6可看出擾動在反應區中的傳波形態復雜，大致可分為3個階段：弓形CJ ODW及雙三波點的形成；新三波點的形成；激波、膨脹波和弱壓縮波的組合傳播。

3.2.1 第1階段

圖6(a)和(b)中，在誘導區內擾動引起的膨脹波穿過誘導區，隨后與CJ ODW作用使其發生彎曲，形成一道弓形CJ ODW，同時形成新的三波點，為與第2階段的新三波點相區別，這里稱第1階段的新三波結構為擾動三波點結構。圖8為此階段的速度局部放大云圖，可看出擾動三波點沿新生成的CJ ODW向下游延伸，同時不斷改變原結構中CJ ODW的形態。

自擾動三波點生成一道近似垂直于壁面的激波，其下游側生成一組膨脹波，這種激波和膨脹波的組合隨擾動三波點向下游傳播。在第1階段，誘導區后段形成了獨特的雙三波點結構。最終，擾動三波點與原波結構中的三波點匯聚形成新波的三波點。至此誘導區便完成了從原三波至新的三波結構的轉變。

圖8 第1階段局部放大速度云圖Fig.8 Partially amplified velocity contour for the first period

3.2.2 第2階段

在該階段，第1階段的雙三波點結構完成了匯聚，新的ODW結構逐漸形成。然而數值模擬發現，新形成的結構存在較大的不穩定性，物性參數出現劇烈的變化，如圖9所示。當t=0.134 ms時新三波點形成，此后新三波點處的溫度、壓力和速度較擾動前都出現了劇烈的數值波動，甚至會產生局部熱點。這表明在溫降擾動下，雖然ODW波結構能夠保持結構的完整性不被破壞，但內部微觀結構的不穩定性卻被放大了，這與擾動傳出后胞格結構更明顯的現象一致，這一點后文會進一步闡述。

圖9 新三波點處溫度、速度和壓力隨時間變化曲線Fig.9 Profiles for temperature, velocity and pressure nearby new triple point changing with time

3.2.3 第3階段

在該階段，ODW的結構變化復雜，存在著不同種類波之間的相互作用，圖10顯示了擾動引起的激波、膨脹波及弱壓縮波在反應區內的傳播過程的簡化圖，其時序與圖6(c)～(f)對應。

擾動在進入反應區后，上游新的ODW結構逐漸形成。在第1階段形成的射向壁面的激波變為弓形激波，激波下游自避免生成一段膨脹區，激波下游氣體經膨脹波加速，為平衡上下游兩側的氣體狀態，自壁面生成一道向下游傳播的弱壓縮波，如圖10(a)所示。

圖10 溫降擾動引起的激波、膨脹波及弱壓縮波在反應區內傳播的時序簡圖Fig.10 Sequence sketches for propagation of shock wave, expansion waves and weak compression wave in reaction zone caused by temperature drop disturbance

圖11 ΔT=-100 K 時3個區域內觀測點的Ma變化Fig.11 Variation of Ma in sampling points from three zones at ΔT=-100 K

在壁面中端，弓形激波在壁面發生了明顯的馬赫反射，在壁面處生成一道馬赫桿，但這種馬赫反射現象持續時間不長，隨著激波向下游傳播，馬赫桿逐漸變短，反射點向壁面逼近，直至在壁面中后段，馬赫反射退化為普通的規則反射，如圖10(b)和(c)所示。由圖10(d)可以看出，在壁面末段，反射激波幾乎消失，只以一道近乎垂直于壁面的激波向下游傳播，直至擾動傳出燃燒室，在此階段，激波強度變得很弱，不像此前能夠引發ODW面和滑移線的強烈的卷曲。

圖11為誘導區、爆燃區和爆震區中3個觀測點的Ma的變化情況，在爆燃區③處，膨脹波在激波下游傳播，Ma存在上升、下降再上升的過程，而在爆震區②處，Ma先下降后上升，這幾乎與前者完全相反。這種現象的原因在于，激波下游的膨脹波在爆燃區能夠保持足夠強度，但在進入爆震區后迅速減弱，以致并不足以引起觀測點②處強烈的參數變化，因此②處Ma的劇烈下降反應的是膨脹波上游的激波的傳播，在激波經過②處后，由于在新的ODW結構中，②處的Ma遠高此前的最低值，必然存在一組膨脹波以使激波上游的氣體壓力降低，從而使壓力相匹配。在觀測點③處，情況則完全不同，貼近壁面處的膨脹波較強，最先引起③處Ma的明顯上升，當上游激波經過③處時引起了Ma的驟降，最后激波上游的膨脹波又將此處的氣體平衡至新的爆燃區所匹配的狀態。

在分析爆震區和爆燃區內Ma變化截然相反的現象時需要結合定性和定量的方法。由前文可知，激波、膨脹波和弱壓縮波在擾動傳播過程中起著不同程度的作用，其中激波和膨脹波的作用能在圖11中明顯觀察到，而弱壓縮波的作用并不明顯，這與觀測點的選取位置有關。此外，爆震區和爆燃區對擾動的響應也不同，達到穩定狀態后，Ma在爆震區波動明顯比爆燃區劇烈，這與前文新三波點處的劇烈數值波動一致，側面體現了不穩定性對爆震波穩定傳播起著本質作用。

3.3 溫升擾動在燃燒室內的傳播

3.2.3節著重研究溫度降低情況下的擾動傳播情況，為了研究溫度上升時的情況，在t=0.1 ms時刻又引入了另一溫升擾動，即溫度上升100 K，選取了相同位置的觀測點，分析方法與3.2.3節相同。

溫升擾動同樣引發了激波、膨脹波和弱壓縮波的復雜組合，圖12(a)為t=0.14 ms時，流場內的Ma分布云圖，圖12(b)為其簡化圖，可看出3種波自爆震波面射出，彼此相互間隔，并且沿ODW面向下游傳播，這與前文引入的溫降擾動沿壁面傳播的情形完全不同。圖13中，爆震區②處Ma數先下降后上升，隨后又有一個下降、上升的過程，這與圖12(b)中3種波相互間隔一一對應。

圖13給出了3個區域內Ma的變化情況。對比圖13和圖10可知， 2種溫度擾動下，爆震區(②處)內的Ma的分布形式大體相反，而在爆燃區(③處)內則完全相同，造成這種結果的主要原因是在溫升擾動下，弱壓縮波的強度明顯較溫降擾動下強，如圖13局部放大圖所示，弱壓縮波穿過滑移線同時影響了爆震區和爆燃區，從而引發Ma在2個區域內都出現了小幅度的下降。此外在溫升擾動的影響下，爆震區內呈現的不穩定性也較爆燃區強烈，這與3.2.3節一致。

圖12 t=0.14 ms時Ma分布云圖和溫升擾動引起的激波、膨脹波及弱壓縮波在反應區內的分布簡圖Fig.12 Distribution contours of Ma and distribution sketches for propagation of shock wave, expansion waves and weak compression wave in reaction zone at t=0.14 ms

圖13 ΔT=100 K時3個區域內觀測點的Ma變化Fig.13 Variation of Ma in sampling points from three zones at ΔT=100 K

4 結 論

本文對溫度擾動在斜爆震燃燒室的傳播進行了數值模擬，分析了在溫度上升和溫度下降2種擾動作用下的擾動傳播形態，得出以下結論：

1) 在溫度擾動下，ODW能夠順利完成結構的調整。引入溫降擾動后，ODW內部不穩定性被進一步釋放，爆震區呈現出多層次的橫波結構，胞格結構更為清晰，波后物性參數的波動也更為劇烈，而且對比爆震區和爆燃區內Ma在擾動傳播過后的數值波動，發現爆震區內的不穩定性較爆燃區劇烈，這側面印證了Dupré等[24]的結論，即不穩定性對爆震波的穩定傳播起著本質作用。

2) 在溫降擾動下，ODW的結構調整可分為3個階段：弓形CJ ODW及雙三波點結構的形成；雙三波點的匯聚及新三波點的形成；激波、膨脹波和弱壓縮波的組合傳播。擾動在傳播過程中不斷衰減，在誘導區和未反應區，擾動以膨脹波的形式傳播。 在反應區，激波會呈現出弓形激波、馬赫反射、規則反射和近乎垂直于壁面的正激波4種形態。

對ODW結構起主要作用的激波和膨脹波。在爆燃區中激波在上游，膨脹波在下游；在爆震區中則相反，膨脹波在上游，激波在下游。這并不意味著在2個區域形成了不同的激波和膨脹波的組合傳播形式，這種分布位置矛盾的原因是激波下游的膨脹波在爆燃區內強度較大，對氣流狀態的影響明顯，但其傳播到在爆震區內時強度衰減明顯，不足以顯著改變氣體狀態，而激波上游膨脹波的強度分布則相反，其在爆震區中的強度要明顯大于在爆燃區中的強度。

3) 在溫升擾動下，激波、膨脹波和弱壓縮波對ODW結構調整都起到了重要作用。通過分析擾動的傳播形態的分布情況，發現與溫降擾動下，弱壓縮波對氣流狀態的影響很弱不同，弱壓縮波在溫升擾動下強度較大，對ODW的結構調整起到了較大的作用。此外，溫降擾動沿壁面向下游傳播，而溫升擾動則沿ODW面向下游傳播，造成這種形態上巨大差異的原因尚不清楚，還需進一步研究。