采用改進的CE/SE方法模擬方管中氫氧爆轟波的穩(wěn)定傳播結(jié)構(gòu)

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

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

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

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

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

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

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

1 基本方程組和化學(xué)反應(yīng)模型

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

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

(1)

式中：流矢量U=[ρρuρvρweραρβ]T，ρ為密度，u、v、w分別為速度矢量在笛卡兒坐標系中沿x、y、z3個方向的分量，α和β分別為誘導(dǎo)和放熱反應(yīng)的進行度(初始為0，反應(yīng)完全為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，ωα和ωβ分別為誘導(dǎo)和放熱反應(yīng)的速率。流矢量E、F、G均為U的函數(shù)，e的表達式為

(2)

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

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

(3)

(4)

2 數(shù)值方法和差分格式

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

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

(5)

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

(6)

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

(7)

(8)

如此，就能利用遞推格式得到下半步空間流矢量U的數(shù)值解。進一步利用解元物理量在交界點處的連續(xù)性，可得到下半步U的x、y、z方向?qū)?shù)的迎風(fēng)和逆風(fēng)表達，再利用權(quán)函數(shù)得到方向?qū)?shù)的一個加權(quán)結(jié)果，就能將此遞推格式循環(huán)下去。

3 數(shù)值模擬通用條件

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

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

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

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

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

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

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

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

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

圖3 不同擾動幅度下歷史壓力極值在xy側(cè)壁的分布(對角擾動，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)

為了確定適合本文數(shù)值模擬的擾動振幅δ，對圖1(a)中的算例(網(wǎng)格尺寸為1/25 mm)分別給予不同的擾動幅度(δ=0.01, 0.05, 0.50)進行試算，結(jié)果如圖3所示?？梢钥吹疆敂_動幅度δ= 0.50時，爆轟波已經(jīng)失穩(wěn)，模擬產(chǎn)生了非物理的結(jié)果；而當擾動幅度δ= 0.01時，形成的爆轟胞格圖案相對比較模糊。最終選擇δ= 0.05作為本文方管爆轟數(shù)值模擬的振幅。

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

4 數(shù)值結(jié)果分析

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

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

觀察圖4(a)中對角模式算例1-1～算例1-6的模擬結(jié)果可以發(fā)現(xiàn)，當截面尺寸在3～8 mm浮動時，穩(wěn)定的爆轟波結(jié)構(gòu)在壁面均投影出一個胞格寬度，胞格數(shù)目沒有改變，胞格尺寸與截面尺寸相對應(yīng)。當橫截面尺寸減小為1 mm 時，爆轟波結(jié)構(gòu)將會從對角模式逐漸轉(zhuǎn)變?yōu)槁菪Ｊ健?/p>

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

如果給予初始直角擾動，通過截面尺寸的變化得到算例2-1～算例2-5，其模擬結(jié)果和對角模式結(jié)果有一些差別。當橫截面尺寸增加至8 mm時，較大的截面尺寸拉伸爆轟波傳播結(jié)構(gòu)，使得橫波分裂并逐漸累積相位差，最終胞格結(jié)構(gòu)分裂成為兩個較小的胞元。另外，當橫截面尺寸減小為2 mm乃至1 mm時，爆轟波結(jié)構(gòu)將會從直角模式轉(zhuǎn)變?yōu)槁菪Ｊ健?/p>

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

當方管截面尺寸只增加一點時，維持橫向胞格數(shù)目成為維持胞格結(jié)構(gòu)穩(wěn)定性的主導(dǎo)因素，所以胞格尺寸會相應(yīng)的增大以保證橫向胞格數(shù)目不發(fā)生變化。

圖4 不同截面尺寸算例典型區(qū)域的xy截面胞格結(jié)構(gòu)Fig.4 Cellular patterns in typical areas on xy sides with different cross-sectional sizes

當方管截面尺寸增加過大，保持小的胞格尺寸則成為了維持胞格結(jié)構(gòu)穩(wěn)定性的主導(dǎo)因素，此時胞格圖樣分裂成更多的橫向胞格以保證每個胞元的小尺寸。

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

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

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

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

表2 截面尺寸在兩種傳播模式下對爆轟波胞格結(jié)構(gòu)的影響

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

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

4.2 方管中的3種穩(wěn)定結(jié)構(gòu)

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

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

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

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

圖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

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

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

為了驗證改進CE/SE算法和Sichel二步模型的計算效果，把本文計算的胞格結(jié)構(gòu)和前人的工作進行了對比，如圖8、圖9所示?？梢钥吹綄τ谥苯呛蛯悄Ｊ絹碚f，模擬結(jié)果能夠與實驗很好的吻合；對于螺旋模式來說，數(shù)值模擬得到的爆轟波結(jié)構(gòu)與其他高精度數(shù)值格式相差不大。

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

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

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

4.3 臨界情況結(jié)構(gòu)變化

在臨界情況算例1-2中，截面尺寸為2 mm，在計算區(qū)域前半部分，爆轟波仍然具有一個橫向胞格，然而此時爆轟波結(jié)構(gòu)處于不穩(wěn)定的狀態(tài)，胞格形狀會出現(xiàn)一定的變形。經(jīng)過一段較短的轉(zhuǎn)化區(qū)之后，爆轟波穩(wěn)定傳播，此時橫波掃掠前導(dǎo)激波形成的紡錘體沿xy,xz平面被剖分成1/4，投影在壁面上顯示為半胞結(jié)構(gòu)。

與典型的對角模式胞格結(jié)構(gòu)相比，此時前導(dǎo)激波波陣面上三波線的數(shù)量減少為一組4條，形成一個封閉的長方形，如圖10所示?？梢韵胍?，過小的截面尺寸壓迫爆轟波橫向結(jié)構(gòu)，使得兩組橫波的相位差逐漸縮小，直至合并成一組橫波，形成在這個截面尺寸下更為穩(wěn)定的傳播結(jié)構(gòu)。

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

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

圖11 過渡區(qū)間壓力等值面演化和三波線運動以及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所示，利用煙熏法得到的非同相管壁胞格圖樣中，垂直管壁的亮條紋即數(shù)值模擬中的拍波，此時亮條紋位于相鄰兩個橫波交匯點的中間區(qū)域，與算例2-2中四頭螺旋模式下的胞格結(jié)構(gòu)十分類似，進一步驗證了計算結(jié)果的合理性。

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

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

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

圖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)

從結(jié)果上來說，對角模式胞格結(jié)構(gòu)比直角模式更為穩(wěn)定，前者三波線累計相位差后仍然可以保持結(jié)構(gòu)的穩(wěn)定性，后者累計相位差后容易產(chǎn)生中間態(tài)的螺旋結(jié)構(gòu)，并最終轉(zhuǎn)變?yōu)閱晤^螺旋結(jié)構(gòu)。

5 結(jié) 論

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

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

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

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

4) 臨界尺寸下爆轟波傳播模式的轉(zhuǎn)化還伴隨著壓力脈沖峰值的非等幅振蕩，新的穩(wěn)定傳播結(jié)構(gòu)形成后振蕩變?yōu)榈确袷帯?/p>