采用改進的CE/SE方法模擬方管中氫氧爆轟波的穩定傳播結構

沈洋，劉凱欣，陳璞，張德良

1.北京大學 工學院，北京 100871 2. 中國科學院 力學研究所 高溫氣體動力學國家重點實驗室，北京 100190

爆轟是一個包含了復雜化學反應的流體動力學過程，由于在建筑工程、能源生產、國防建設等領域的重要應用，爆轟波的起爆和傳播機制，特別是胞格結構形成機理，仍然是當下熱門的研究方向。目前關于一維爆轟波的理論模型已構建完善，對于平面氣相爆轟波傳播的研究也已經相當成熟，學者們把目光越來越多地聚焦到爆轟波三維結構的研究上面。Hanana等[1]在考察方管中爆轟波傳播壓力結構時，使用煙熏法得到了管壁上胞格結構的清晰圖樣。經過細致的分析，他們認為方管中爆轟波傳播具有兩種穩定的結構：對角模式和直角模式。除此之外，周凱等[2]利用爆轟驅動技術與膨脹管的結合，設計并初步研究了高速氣流爆轟驅動膨脹管。Zheng和Wang[3]通過實驗考察了低溫等離子體對爆燃轉爆轟過程的加速作用。這些都是實驗研究三維爆轟波形成與傳播的重要嘗試。

由于實驗技術的限制，想要捕捉三維管道內部爆轟波復雜結構的細節信息目前來講仍有些困難。相比實驗手段，數值方法的最大優勢是可以獲得任意物理量在任意時刻的全局分布，但由于三維模擬需要的網格數以千萬計，計算量也是一個不容忽視的問題。在模擬爆轟波傳播過程中的一個關鍵問題是選擇合理的化學反應模型。Oran 等[4]提出了基元反應模型并用來考察平面氫氧爆轟的胞格結構。Sichel等[5]在Taki和Fujiwara[6]的二步模型基礎上引入了考慮氣體組分變化的參數，并在氫氧爆轟系統的算例中給出了具體的公式，提高了二步模型的計算精度。

Tsuboi等[7-8]采用基元反應和non-MUSCL(non-Monotonic Upwind Scheme for Conservation Laws)格式對Hanana的方管實驗進行了數值模擬，發現能夠得到有效數值結果的管道尺寸比爆轟實際的胞格尺寸要小，另外他們還對在圓管中的旋轉爆轟進行了數值模擬[9-10]。竇華書[11-12]和王成[13-14]等用一步模型和五階WENO(Weighted Essentially Non-Oscillatory)格式對三維爆轟波在方管中的傳播進行了進一步的數值分析，發現爆轟波在較小尺寸的方管中也存在類似于圓管中的旋轉傳播模式。另外，翁春生和Gore[15]、申華[16]、Ivanov[17]、蔡曉東[18]以及黃玥[19-20]等同樣對方管中爆轟波的傳播做了詳盡的數值分析并且得到了類似的結論。然而，絕大部分的數值模擬工作需要高階精度的數值格式和細致的化學反應模型來獲得較好的胞格圖案，否則得到的結果將不夠清晰，而本文采用的時空守恒元/解元(CE/SE)算法是兼顧計算效率與質量的一種高精度數值格式。

時空守恒元/解元算法是近年來興起的一種全新的求解雙曲守恒型方程的計算方法，將時間和空間統一起來同等對待，巧妙地定義時空間的守恒元和解元使得局部和全局都嚴格保證守恒率。NASA Lewis研究中心的Chang[21]首先提出這種差分格式，并將其推廣到二維情況[22]。Wang[23]針對含有激波的氣動聲學問題，對CE/SE算法的計算精度進行了分析。Zhang等[24]基于四邊形網格劃分方案改進了守恒元和解元的結構，并由此順利推廣到了三維情況。

本文針對三維氫氧爆轟問題，使用基于四邊形網格劃分方案的改進時空守恒元/解元格式和Sichel新型二步化學反應模型，以較小的計算代價得到了方管中3種穩定傳播模式的爆轟波結構，并著重于討論截面尺寸在兩種傳播模式下對管道中爆轟波結構穩定性的影響。

1 基本方程組和化學反應模型

爆轟波傳播是爆轟間斷以超聲速在流體中傳播的過程，在數值模擬時，一般情況下可以不考慮熱傳遞和黏性效應，因此流體力學基本方程組實際上就是歐拉方程和描述化學反應的方程組成的方程組。

在三維笛卡兒坐標系中流體力學基本方程組和二步反應模型可以表示為

(1)

式中：流矢量U=[ρρuρvρweραρβ]T，ρ為密度，u、v、w分別為速度矢量在笛卡兒坐標系中沿x、y、z3個方向的分量，α和β分別為誘導和放熱反應的進行度(初始為0，反應完全為1)，e為單位體積的總能；E=[ρuρu2+pρuvρuw(e+p)uραuρβu]T，p為壓強；F=[ρvρuvρv2+pρvw(e+p)vραvρβv]T；G=[ρwρuwρvwρw2+p(e+p)wραwρβw]T；剛性源項S=[0 0 0 0 0ωαωβ]T，ωα和ωβ分別為誘導和放熱反應的速率。流矢量E、F、G均為U的函數，e的表達式為

(2)

式中：Q為單位質量的氫氧燃燒熱；γ為氫氣的比熱容比。

Sichel等在氫氧混合系統下構建的改進二步模型[5]分為兩個步驟。第1步描述了一個以自由基形成為特征，無顯著能量釋放的誘導期，經過詳細比較，在諸多現有的模型中選擇了Burks和Oran提出的反應模型，得到了擬合式(3)。第2步描述了產生燃燒產物與釋放大量熱量的放熱反應期，在這里，通用Arrhenius公式式(4)被認為是對氫氧混合氣體放熱反應的最準確描述。通過再現基元反應模型所描述的化學動力學模型的性質，可擬合得到Arrhenius公式具體的輸入參數。

(3)

(4)

2 數值方法和差分格式

爆轟波是含有化學反應的強間斷流動，對于數值模擬的精度，特別是強間斷面附近的模擬精度要求很高。CE/SE方法用于三維爆轟波的數值模擬有著得天獨厚的優勢：首先，相比傳統差分方法，它從控制方程的時空積分形式出發，在時間和空間上都能夠很好地保證物理量的守恒性；其次，通過巧妙的守恒元解元的設置以及權函數的引入，其在捕捉爆轟波的強間斷方面具有良好的效果；最后，通過改進的守恒元解元劃分，格式很容易推廣到三維情況。總的來說，CE/SE格式是一種保持高精度和低計算量的具有極高性價比的數值格式。下面在三維一階Taylor展開CE/SE格式[16]的基礎上推導二階Taylor展開下三維CE/SE遞推格式。

由于二步模型的化學反應弛豫時間比CFL(Courant, Friedrichs and Lewy)時間步長小1 ～ 2個數量級，因此采用化學反應與爆轟推進解耦的方法，在推導格式時，不必考慮源項的作用。根據散度定理,守恒方程可以寫成積分形式：

(5)

式中：張量H={U,E,F,G}為時空間上流矢量的組合；S(V)為任意封閉的時空體域V的邊界；ds=dσ·n，其中dσ和n分別為邊界S(V)的面積和外法向。利用解元[16]在網格基點的二階Taylor展開，并讓展開式在守恒元[16]上進行積分，代入式(5)，最終可以得到U的半步遞推格式：

(6)

式中：算符±a、 ±b、 ±c是相互獨立取正負的；而算符±a和?a意味著其中一個取正號時，另一個則取負號；對a、b、c的求和表示對算符±a、 ±b、 ±c窮舉所有正負；符號“∧”代表一個重定義的函數。流矢量N重定義函數的第m項分量在點A的展開為

(7)

(8)

如此，就能利用遞推格式得到下半步空間流矢量U的數值解。進一步利用解元物理量在交界點處的連續性，可得到下半步U的x、y、z方向導數的迎風和逆風表達，再利用權函數得到方向導數的一個加權結果，就能將此遞推格式循環下去。

3 數值模擬通用條件

實驗表明爆轟波具有非常復雜的三維結構，主要包括前導激波、馬赫桿(MS)、橫波(TW)等，它們相互作用形成了爆轟波陣面。爆轟波陣面具有不穩定性同時又具有時空周期性，爆轟波陣面的周期性變化形成了所謂的胞格結構，體現了爆轟的很多特征參數，但是其機理仍未完全探明，因此爆轟波傳播形成的胞格結構是爆轟波研究的一個熱點。

本文針對方管中爆轟波傳播的物理問題進行數值模擬。影響方管內爆轟波穩定傳播結構的因素有很多，本文重點研究截面尺寸在兩種傳播模式下對爆轟波穩定傳播以及形成的胞格結構的影響，對于其余控制變量采用統一的設置。以下是本文程序所使用的通用模擬條件。

1) 初始條件：在管道左側5%的區域設置p=50p0,ρ=ρ0的高壓區，其余空間則充滿了p0、T0、ρ0的氫氧混合氣。其中p0=105Pa,T0=298 K,ρ0=p0/(RT0)，在本文中氫氣與氧氣以2:1的摩爾比相混合，因此混合氣體常數R=689 J/(kg·K)。整個空間處于靜止的狀態，在計算開始時隔板消失。

2) 邊界條件：與管道方向平行的4個壁面以及左邊界都使用反射邊界條件，右邊界使用自由邊界條件。實際上右邊界條件的設置對計算結果沒有任何影響。

3) 網格設置：為了找到合適的計算網格，首先對網格尺寸d=1/10,1/25,1/50 mm 3種網格密度的方管算例進行了試算，然后對數值結果進行了比較。不同網格密度下側壁胞格結構如圖1所示，所謂的胞格結構在數值上即為壓力歷史極值在壁面上的分布。結果顯示，當截面尺寸(D)在2 mm及以上時，d=1/25 mm的網格密度已經足夠滿足計算需求；但當截面尺寸為1 mm時，d=1/25 mm得到的胞格結構相比d=1/50 mm 的結果稍顯模糊，因此將網格密度加倍。

圖1 不同截面和網格尺寸下歷史壓力極值在方管上表面的分布(對角模式)Fig.1 Distribution of maximum pressure histories on upper surface of square duct with different cross-sectional sizes and grid sizes (diagonal mode)

4) 方管長度設置：當截面尺寸較小時(≤2 mm)，爆轟波形成穩定傳播模式可能會出現較長的過渡區，此時方管長度用160 mm為宜，一般情況下80 mm即能出現穩定的爆轟波結構。

5) 初始擾動設置：為了形成橫波，在初始30個 時間步的化學反應區內(0<β<0.99)對內能e給予小幅度的正弦擾動。擾動分為沿著邊界的直角模式和沿著對角線的對角模式，其無量綱化的擾動幅度e′/e在方管yz截面的分布函數如圖2所示(δ為擾動幅度極值)。

圖2 方管中兩種擾動模式在yz截面上的振幅分布Fig.2 Amplitude distribution of two kinds of perturbation modes on yz section of square duct

圖3 不同擾動幅度下歷史壓力極值在xy側壁的分布(對角擾動，D=2 mm，d=1/25 mm)Fig.3 Distribution of maximum pressure histories on xy sides with different amplitudes of perturbation (diagonal mode, D=2 mm, d=1/25 mm)

為了確定適合本文數值模擬的擾動振幅δ，對圖1(a)中的算例(網格尺寸為1/25 mm)分別給予不同的擾動幅度(δ=0.01, 0.05, 0.50)進行試算，結果如圖3所示。可以看到當擾動幅度δ= 0.50時，爆轟波已經失穩，模擬產生了非物理的結果；而當擾動幅度δ= 0.01時，形成的爆轟胞格圖案相對比較模糊。最終選擇δ= 0.05作為本文方管爆轟數值模擬的振幅。

6) 并行設置：對方管軸向進行分割，相鄰的計算塊保留界面附近的兩層數據以作通訊，利用MPICH2接口的通訊函數形成多線程并行計算。

4 數值結果分析

4.1 不同截面尺寸的算例比較

不同的傳播模式下，截面大小對胞格結構尺寸的影響不盡相同。利用如圖2所示的初始擾動形成直角和對角傳播模式。在控制其他變量不變的情況下方管的截面尺寸在1～8 mm范圍內改變。截面尺寸為2～8 mm時，網格尺寸均為1/25 mm；截面尺寸為1 mm時，為了保持最終結果的精度，網格尺寸加密為1/50 mm，所有算例的計算參數如表1所示。將表1中所有算例的xy側壁截面胞格結構繪制在一起，如圖4所示。

觀察圖4(a)中對角模式算例1-1～算例1-6的模擬結果可以發現，當截面尺寸在3～8 mm浮動時，穩定的爆轟波結構在壁面均投影出一個胞格寬度，胞格數目沒有改變，胞格尺寸與截面尺寸相對應。當橫截面尺寸減小為1 mm 時，爆轟波結構將會從對角模式逐漸轉變為螺旋模式。

表1 不同算例的模擬參數Table 1 Simulation parameters in different cases

如果給予初始直角擾動，通過截面尺寸的變化得到算例2-1～算例2-5，其模擬結果和對角模式結果有一些差別。當橫截面尺寸增加至8 mm時，較大的截面尺寸拉伸爆轟波傳播結構，使得橫波分裂并逐漸累積相位差，最終胞格結構分裂成為兩個較小的胞元。另外，當橫截面尺寸減小為2 mm乃至1 mm時，爆轟波結構將會從直角模式轉變為螺旋模式。

以上數值結果指出，方管尺寸變化雖然可以使胞格尺寸隨之產生適應性的變化，但存在一定的極限。通常來講，胞格尺寸越小，橫向胞格個數越少，胞格結構越穩定。但保持小胞格尺寸和維持橫向胞格數是矛盾的，若橫向胞格數目增加一倍，則胞格尺寸必然減小一半，而隨著截面尺寸的變化，這兩個因素的主導地位也會隨之產生變化。

當方管截面尺寸只增加一點時，維持橫向胞格數目成為維持胞格結構穩定性的主導因素，所以胞格尺寸會相應的增大以保證橫向胞格數目不發生變化。

圖4 不同截面尺寸算例典型區域的xy截面胞格結構Fig.4 Cellular patterns in typical areas on xy sides with different cross-sectional sizes

當方管截面尺寸增加過大，保持小的胞格尺寸則成為了維持胞格結構穩定性的主導因素，此時胞格圖樣分裂成更多的橫向胞格以保證每個胞元的小尺寸。

另外，幾何性質決定胞格結構總是在壁面邊界擁有整數或者半整數個胞格，也即魚鱗狀或者半魚鱗狀的圖案；而直角模式中，由于橫波平行于壁面進行掃掠，因而永遠不會出現半整數的胞格結構。

最后，當截面尺寸足夠小時(1 mm)，則不管初始擾動條件如何，穩定的爆轟波傳播結構都會轉變為螺旋模式，旋轉方向是順時針還是逆時針取決于初始擾動形式。臨界情況下(2 mm)，對角模式會在計算區域后半段衰減為半胞結構，而直角模式則會在區域中部形成四頭螺旋的中間傳播態，并最終轉化為螺旋模式。關于螺旋模式旋轉方向以及臨界尺寸下爆轟波結構變化將在下文進行詳細的討論。

在考察截面尺寸對方管氫氧爆轟波傳播結構的影響時，擾動形式是一個不可忽略的因素。乍一看這似乎只是一個數學游戲，沒有過多的物理意義，然而數值模擬結果表明，無論初始給予什么擾動，隨機擾動，非對稱擾動，正弦擾動還是矩形擾動，最終的傳播模式總是有規律可循，或者是對角模式，或者是直角模式，或者是螺旋模式。究其原因，可能是初始擾動橫波在數學上可以表達為不同頻率和振幅的傅里葉級數，但是在擾動傳播過程中，與方管尺寸不相容的頻率項都受到強阻尼而衰減了，最終留下特定頻率的擾動與爆轟波相互正作用形成主要項，這類擾動項最終形成了所觀察到的爆轟波傳播結構。

表2給出了不同截面尺寸在兩種傳播模式下對爆轟波胞格結構的綜合性影響。

表2 截面尺寸在兩種傳播模式下對爆轟波胞格結構的影響

Table 2 Effect of cross-sectional size with two kinds of propagation modes on cellular patterns of detonation wave

D/mm胞格結構對角模式直角模式1(很小)順時針螺旋逆時針螺旋2(臨界情況)對角→半胞結構直角→四頭螺旋→單頭螺旋>2(適中)對角直角

4.2 方管中的3種穩定結構

為了驗證CE/SE算法搭配Sichel改進二步模型的可靠性，從表1的所有算例中選取了4個典型算例(算例1-1，算例1-4，算例2-1，算例2-4)，用來考察方管中的3種穩定傳播模式：對角模式，直角模式和螺旋模式，并用本文的模擬結果與前人的結果進行比較。圖5給出了方管中3種穩定傳播模式在一個周期內的壓力等值面演化過程。

在圖5(a)和圖5 (b)中，計算區域為80 mm×4 mm×4 mm，初始分別給予對角擾動和直角擾動。從等值面演化過程中可以看出，爆轟波在傳播過程中，前導激波與橫波(TW)交匯形成一些三波線(TL)，其軌跡在壁面形成魚鱗狀的圖案，即胞格結構；前導激波波陣面的各區域被三波線交替分割為凸離爆轟方向的馬赫桿(MS)和相對凹陷的入射波(IS)。相鄰三波線對撞、分離使得新的馬赫桿不斷產生，舊的馬赫桿衰減為入射波。在對角模式圖5(a)中，兩組共8條三波線形成兩個封閉的長方形，彼此相互垂直且與通道壁面呈45°。類似的，在直角模式圖5(b)中同樣存在兩組共4條三波線，兩組平行線彼此相互垂直，且各自平行于壁面。在每個胞格周期內，當一組三波線掃掠至平行的壁面附近時，會在另一組三波線形成的壁面胞格結構中顯示出橫向的亮線，稱為拍波(Slapping wave)，以此作為直角模式和對角模式胞格圖案最為顯著的不同之處。在圖5(c)和圖5(d)中，分別給予和圖5(a)、圖5(b)相同的初始擾動，但是把截面尺寸減少至1 mm。計算結果顯示存在兩種類似的螺旋結構，一種沿順時針螺旋，另一種沿逆時針螺旋。此時只有一對相互垂直且正交于壁面的三波線存在。和圖5(b)的直角模式不同之處在于這一對三波線存在大約π/4的相差，失去了軸對稱性質，如圖6所示。如果觀察此時的三維爆轟波結構，會看到一條明亮的壓力極值帶沿著管壁進行螺旋運動，在管道內部則沒有明顯的壓力集中區域，如圖7所示(圖中pmax為歷史壓力極值)。

圖5 約一個周期內的壓力等值面演化過程Fig.5 Evolutionary process of pressure iso-surfaces during about one period

對于螺旋模式的旋轉方向，有順時針的模擬結果，也有逆時針的結果。竇華書等[11]使用隨機擾動和一步反應模型，得到的模擬結果是順時針螺旋；王成等[13]使用沿著對角線的正弦擾動和一步反應模型，得到的模擬結果是逆時針螺旋，而Tsuboi等[7]使用基元反應模型并設置了不對稱擾動，得到了逆時針旋轉的結果。

圖6 兩種螺旋模式下三波線運動方式示意圖Fig.6 Sketch of motion of triple point lines in two kinds of spinning mode

圖7 兩種螺旋模式下歷史壓力極值等值面Fig.7 Iso-surfaces of maximum pressure histories in two kinds of spinning modes

本文的模擬結果顯示，只要設置了圖2(a)所示直角模式的擾動，不管截面尺寸和網格密度多少，模擬結果均為逆時針螺旋；而設置圖2(b)所示的對角擾動則一定會產生順時針螺旋。

以上結果表明，不同的擾動形式的確對旋轉方向產生了影響，且同一擾動形式在不同截面尺寸或者其他初始條件變化情況下只會產生一種旋轉方向，不具有隨機性。

為了驗證改進CE/SE算法和Sichel二步模型的計算效果，把本文計算的胞格結構和前人的工作進行了對比，如圖8、圖9所示。可以看到對于直角和對角模式來說，模擬結果能夠與實驗很好的吻合；對于螺旋模式來說，數值模擬得到的爆轟波結構與其他高精度數值格式相差不大。

圖8 不同數值模擬方法計算方管中對角模式和直角模式胞格結構與實驗圖樣的比較Fig.8 Comparison of cellular patterns in diagonal or rectangular mode of square ducts using different numerical methods with experimental results

圖9 不同數值模擬方法計算方管中螺旋模式胞格結構的比較Fig.9 Comparison of cellular patterns in spinning mode of square ducts using different numerical methods

一般認為，化學反應模型決定了胞格結構的形狀和特征尺度，數值模型決定了胞格圖案的精細程度和對比度。比較的結果表明，使用Sichel二步模型計算得到的胞格可以和基元反應模型相媲美，略優于一步模型；使用改進CE/SE格式得到的胞格結構細節清晰，不遜色于高階精度的WENO格式。而使用改進CE/SE格式和Sichel二步模型所花費的計算資源卻相當少，一般1 000萬三維網格配置和5 000個時間步的方管算例用普通的4核CPU計算，在一周內即可得到完整結果。

4.3 臨界情況結構變化

在臨界情況算例1-2中，截面尺寸為2 mm，在計算區域前半部分，爆轟波仍然具有一個橫向胞格，然而此時爆轟波結構處于不穩定的狀態，胞格形狀會出現一定的變形。經過一段較短的轉化區之后，爆轟波穩定傳播，此時橫波掃掠前導激波形成的紡錘體沿xy,xz平面被剖分成1/4，投影在壁面上顯示為半胞結構。

與典型的對角模式胞格結構相比，此時前導激波波陣面上三波線的數量減少為一組4條，形成一個封閉的長方形，如圖10所示。可以想見，過小的截面尺寸壓迫爆轟波橫向結構，使得兩組橫波的相位差逐漸縮小，直至合并成一組橫波，形成在這個截面尺寸下更為穩定的傳播結構。

和算例1-2類似，算例2-2也接近轉化為螺旋模式的臨界情況，此時直角模式下并不存在如圖10所示的半胞結構，而需要經過一個更長的轉化區。在此種情況下，一開始爆轟胞格圖樣按照初始擾動形式產生一個橫向胞格。接下來兩組三波線會逐漸在計算區域中段累積π/8弧度的相差。如圖11所示，三波線的運動模式與算例1-4中看到的情形類似，但是兩組三波線交匯形成的封閉圖樣是長方形而不是正方形，這是π/8相差產生的直觀結果。另外考察圖11中yz截面歷史壓力極值的分布，可以清楚地看到4個壓力集中的點(三波點)螺旋式地前進，這種波陣面運動模式就是圖6中4個單頭螺旋以不同相位的組合，因此可以稱這種胞格傳播模式為四頭螺旋模式。進一步觀察圖4(b)中截面尺寸為2 mm的xy壁面胞格結構(即歷史壓力極值分布)可以看到,與一般直角模式相比，四頭螺旋的拍波由于相位差不同,并不會出現在魚鱗狀胞格結構的頂點連線附近,而是會出現在頂點連線和中心匯聚點的中間位置。

圖10 算例1-2中一個周期內壓力等值面和三波線的運動Fig.10 Pressure iso-surfaces and motion of triple point lines in one period in Case 1-2

圖11 過渡區間壓力等值面演化和三波線運動以及yz截面歷史壓力極值變化(D=2 mm，初始直角擾動)Fig.11 Pressure iso-surfaces evolution and motion of triple point lines as well as variation of maximum pressure histories on yz sections in transition zone (D=2 mm, initial rectangular perturbation)

在Hanana等[1]的實驗中可以得到這種非同相(Partially out of phase)的直角模式。如圖12所示，利用煙熏法得到的非同相管壁胞格圖樣中，垂直管壁的亮條紋即數值模擬中的拍波，此時亮條紋位于相鄰兩個橫波交匯點的中間區域，與算例2-2中四頭螺旋模式下的胞格結構十分類似，進一步驗證了計算結果的合理性。

圖12 Hanana[1]實驗中不同相的實驗結果Fig.12 One of typical soot records with partially out of phase in Hanana’s experiments[1]

當然，四頭螺旋結構也只是計算區域中段的過渡結構，最終處于平行狀態的兩組橫波會逐漸縮小相位差乃至合并成一組，形成如圖6所示的單頭螺旋結構。

為了進一步比較直角模式和對角模式的穩定性，分析臨界尺寸時波陣面附近壓力脈沖極值的變化與爆轟波結構變化的關系，本文提取了臨界情況算例歷史壓力極值在中心線(y= 1 mm,z= 1 mm)處沿著x方向的分布圖，如圖13所示。在起爆階段，兩種傳播模式受初始擾動影響，都暫時形成一個橫向胞格的對角模式和直角模式。在轉化區，兩種模式的爆轟波都會出現中間轉換結構，如前所述，對角模式下爆轟波由一個橫向胞格逐漸轉變為半胞結構。雖然胞格結構的轉化區很短，但是按照壓力脈沖幅度的變化來看，在很長一段區域內每個周期的壓力振幅和峰值并不穩定。在對角模式下的過渡結構中壓力脈沖峰值的振蕩出現不規則的漲落，而在直角模式下則呈現振蕩幅度不斷下降的趨勢，直到兩種模式達到新的穩定結構。因此從這個角度來講，對角模式在此種尺寸下會轉變為半胞結構是合理的。另一方面，在直角模式下，可以發現壓力歷史的包絡線呈現出一個減小的趨勢，直到喪失穩定性。如前所述，此區域內直角模式胞格圖樣會逐漸轉變為四頭螺旋結構。最終穩定區域兩種模式的壓力脈沖極值都形成等幅振蕩，這意味著爆轟波結構形成自持穩定。對角模式的穩定結構是半胞結構，而直角模式在穩定區形成單頭螺旋爆轟，此時三波線匯聚成點沿著壁面做螺旋運動，中心線附近的歷史壓力極值是很小的。

圖13 兩種傳播模式歷史壓力極值沿著x方向中心線的比較(D=2 mm)Fig.13 Comparison of maximum pressure histories along central lines in x-axis direction between two kinds of propagation modes (D=2 mm)

從結果上來說，對角模式胞格結構比直角模式更為穩定，前者三波線累計相位差后仍然可以保持結構的穩定性，后者累計相位差后容易產生中間態的螺旋結構，并最終轉變為單頭螺旋結構。

5 結 論

本文使用改進的三維CE/SE格式和Sichel二步反應模型對方管中氫氧混合爆轟進行了模擬，主要考察了不同截面尺寸對兩種傳播模式下爆轟波結構與胞格結構形成的影響，初步討論了臨界尺寸下爆轟波結構的變化和穩定性,得出如下結論：

1) 當截面尺寸在合理范圍內變化時，爆轟波結構隨之產生適應性的變形；若變化過大，則會出現胞格的合并與分裂、三波線的消失與新生。而當截面尺寸足夠小時，不管是直角模式還是對角模式最終都會轉化為螺旋模式。

2) 臨界尺寸下爆轟波傳播模式的轉化伴隨著三波線運動相位差的逐漸累積，最終形成固定相差，直到出現三波線合并消失后，胞格結構達到新的穩定。

3) 臨界尺寸下直角傳播模式可得到四頭螺旋中間結果，與實驗結果相吻合。

4) 臨界尺寸下爆轟波傳播模式的轉化還伴隨著壓力脈沖峰值的非等幅振蕩，新的穩定傳播結構形成后振蕩變為等幅振蕩。