郗雪辰,楊鵬飛,滕宏輝
(1.北京理工大學宇航學院,北京 100081;2.山西警察學院治安系,山西太原 030401;3.中國科學院力學研究所高溫氣體動力學國家重點實驗室,北京 100190;4.中國科學院大學工程科學學院,北京 100049)
與傳統(tǒng)的爆燃燃燒相比,爆轟燃燒依靠強激波壓縮進行點火并實現(xiàn)自持傳播,具有能量轉(zhuǎn)化速率快、波傳播速度高的特點?;诒Z燃燒的航空航天推進系統(tǒng)可望實現(xiàn)較高的熱循環(huán)效率,且燃燒室尺寸小、結(jié)構(gòu)質(zhì)量小[1-2],在高超聲速推進領(lǐng)域具有潛在的工程應(yīng)用價值。氣相爆轟波涉及流動、燃燒及其耦合作用等復雜的物理-化學過程,爆轟波面結(jié)構(gòu)復雜多變,波面失穩(wěn)后形成的非定常波系對燃燒組織和發(fā)動機穩(wěn)定工作構(gòu)成挑戰(zhàn)[3]。波面失穩(wěn)源于化學反應(yīng)和前導激波不斷變化的耦合關(guān)系,宏觀上表現(xiàn)為多道橫波與入射波共同構(gòu)成的復雜波面。弱耦合的入射激波與強耦合的Mach桿,通過橫波聯(lián)結(jié),三者的交點被稱為三波點。研究者利用煙跡技術(shù)記錄了三波點的運動軌跡[4],通常稱為胞格結(jié)構(gòu),并用來表征爆轟波的動力學過程。過去幾十年,研究者圍繞爆轟的胞格結(jié)構(gòu)開展了大量的實驗研究,獲取了氫氣、甲烷、乙炔、乙烯等氣體燃料在不同狀態(tài)下的胞格特性[5]。隨著計算機技術(shù)的發(fā)展,利用數(shù)值模擬圍繞胞格結(jié)構(gòu)開展了大量研究,深化了對氣相爆轟的認識[6-9]。
雖然真實爆轟波均為三維,但是一維爆轟波在很多情況下也是一種很好的宏觀近似,反映了激波和放熱在傳播方向上的耦合關(guān)系,可作為理解多維爆轟波不穩(wěn)定的基礎(chǔ),因此得到了廣泛的研究。Fickett等[10]針對一維爆轟波開展了開創(chuàng)性的研究,Erpenbeck[11]對一維和二維爆轟波進行了穩(wěn)定性分析,利用漸進理論得到了爆轟波的穩(wěn)定性邊界。He等[12]研究了活化能對一維爆轟波的影響,發(fā)現(xiàn)隨著活化能的增加,爆轟波從穩(wěn)定傳播發(fā)展為間歇爆轟波,最后由于振蕩幅度過大,爆轟波不再能穩(wěn)定傳播。以上研究假設(shè)波前氣體狀態(tài)固定不變。事實上,在爆轟發(fā)動機研究的過程中發(fā)現(xiàn),高速來流中爆轟波的波系結(jié)構(gòu)對來流參數(shù)的變化非常敏感[13],且來流狀態(tài)的變化還可能引起波前組分不均勻,而這種不均勻性也被認為是造成發(fā)動機預測性能與實測性能之間差異的重要原因之一[14-16]。因此研究爆轟波在非定常來流中的傳播就有了非常重要的實際意義。
非定常來流問題之前多圍繞來流狀態(tài)突變和瞬時離散擾動等進行[17-20]。近來,Kasimov等[21]采用Burgers模型研究了一維爆轟波對連續(xù)擾動的響應(yīng)過程。結(jié)果顯示,連續(xù)小擾動不會影響一維穩(wěn)定爆轟波的動力學過程,但會激發(fā)脈沖爆轟波的內(nèi)在不穩(wěn)定性,使其呈現(xiàn)出更為復雜的動力學特征。Burgers模型獲得的結(jié)果與單步反應(yīng)模型的結(jié)果定性一致。Kim等[22]采用單步反應(yīng)模型對一維脈沖爆轟波在連續(xù)密度擾動中的傳播進行了數(shù)值模擬,研究了擾動波長對爆轟波振蕩行為的影響。結(jié)果發(fā)現(xiàn)對于自身不穩(wěn)定的爆轟波,波長在一定范圍內(nèi)的擾動,反而會抑制一維爆轟的不穩(wěn)定性。這一結(jié)果不僅驗證了非線性理論在爆轟動力學研究領(lǐng)域應(yīng)用的可行性,也為研究利用微擾動技術(shù)控制動態(tài)爆轟參數(shù)等問題提供了新的思路,具有非常重要的啟發(fā)意義。
此前的研究采用簡化的單步反應(yīng)模型,該模型在討論爆轟波氣體動力學和化學動力學之間的相互作用時,存在一定的局限性。兩步反應(yīng)模型可更好地模擬真實爆轟中的誘導-放熱反應(yīng)過程,且具有比單步反應(yīng)更多的自由參數(shù),更利于開展數(shù)值研究。本文采用兩步誘導-放熱反應(yīng)模型,模擬一維爆轟波在連續(xù)正弦密度擾動中的傳播,主要關(guān)注擾動波長對穩(wěn)定爆轟波動力學特性的影響。模擬獲得了一些新的現(xiàn)象,以此為基礎(chǔ)對擾動爆轟波的傳播機理進行了分析,深化了對擾動來流中一維爆轟波傳播機理的認識。
控制方程采用的是忽略黏性的一維Euler方程
假設(shè)氣體為具有固定比熱比的理想氣體,總能量和狀態(tài)方程為
式中,ρ,u,p,e,q,T分別代表流體的密度、速度、壓力、總能量、化學反應(yīng)放熱量和氣體溫度?;瘜W反應(yīng)放熱量q的表達式為
q=ηQ
其中η表示放熱反應(yīng)進程變量。上述所有變量均使用波前未反應(yīng)氣體狀態(tài)參數(shù)進行無量綱化
其中,下角標0表示波前的氣體參數(shù)。上述方程與化學動力學模型耦合,用來描述爆轟波的結(jié)構(gòu)。采用Ng等[23]在Short等[24]基礎(chǔ)上修正的兩步鏈分支反應(yīng)模型,用兩個化學速率控制方程模擬化學反應(yīng)過程。第1步是熱中性誘導區(qū)域點火過程,反應(yīng)速率是對溫度敏感的Arrhenius形式
其中,ξ表示誘導區(qū)反應(yīng)過程的變量,其變化范圍是1 → 0;誘導區(qū)定義為爆轟波前導激波與最大釋熱速率位置之間的距離,對于兩步誘導-放熱總包反應(yīng)模型而言,化學反應(yīng)進程變量ξ=0對應(yīng)著最大釋熱速率的位置;kI表示誘導區(qū)的化學反應(yīng)速率常數(shù),令kI=-Uvn,Uvn是以一維CJ爆轟波為坐標系前導激波后氣流的流動速度;EI表示誘導區(qū)活化能;Ts表示氣體經(jīng)過一維CJ正爆轟波前導激波之后的溫度;H(1-ξ)是階躍函數(shù)
第2步是緩慢放熱之后的快速熱量釋放過程,其反應(yīng)速率方程如下
其中,η是放熱反應(yīng)進程變量,其變化范圍是0 → 1;kR表示放熱區(qū)反應(yīng)速率常數(shù);ER表示放熱區(qū)活化能。上述公式中的EI和ER都通過RT0無量綱化,Ts通過T0無量綱化。
本文采用基于Mach數(shù)分裂的AUSMPW+格式[25]求解通量,空間差分采用的是具有3階精度的MUSCL格式,時間推進采用3階Runge-Kutta方法。計算域網(wǎng)格尺寸為無量綱的 0.05,誘導區(qū)寬度內(nèi)的平均網(wǎng)格密度為20。將穩(wěn)態(tài)Zel′dovich-von Neumann-D?ring(ZND)爆轟參數(shù)作為初始點火源設(shè)置在計算域最左側(cè),計算開始后爆轟波自左向右傳播。模擬的總體長度不同算例有差別,最少為100倍擾動波長距離,以保證壓力峰值統(tǒng)計結(jié)果的可靠性。為提高計算效率,采用了移動網(wǎng)格技術(shù),計算域隨著爆轟波面不斷前移。計算域的長度取決于擾動波長和誘導區(qū)長度,保證不小于擾動波長且至少達到200倍的誘導區(qū)長度。根據(jù)計算經(jīng)驗,此長度能夠有效避免動網(wǎng)格后邊界截斷對爆轟傳播的影響。為了避免初始數(shù)值振蕩的影響,爆轟波首先在定常流動中穩(wěn)定傳播一段時間之后,再施加連續(xù)的密度擾動,其擾動形式為
其中A和λ分別表示擾動的振幅和波長,振幅被設(shè)置為波前靜止氣體壓力的10%,波長是利用誘導區(qū)無量綱化后的長度。計算采用的無量綱參數(shù)為:放熱量Q=50,比熱比γ=1.2,誘導區(qū)活化能EI=6.0Ts,放熱區(qū)活化能ER=1.0Ts。放熱反應(yīng)速率常數(shù)kR作為可調(diào)參數(shù),用來控制爆轟波的內(nèi)在不穩(wěn)定性。
可燃氣體經(jīng)過爆轟波前導激波面的絕熱壓縮,壓力和溫度迅速升高,并引發(fā)快速的化學反應(yīng)。反應(yīng)過程包含誘導階段和緊隨其后的放熱階段。誘導階段的壓力和溫度基本保持不變,放熱階段釋放能量的同時驅(qū)使反應(yīng)產(chǎn)物膨脹,表現(xiàn)為溫度上升且壓力下降。采用ZND模型可以更直觀地表達上述過程。圖1顯示了利用兩步誘導-放熱反應(yīng)模型獲得的穩(wěn)態(tài)ZND結(jié)構(gòu)。爆轟波設(shè)置在流場最左側(cè),并自左向右傳播,跨過激波面的壓力和溫度瞬間上升至von Neumann狀態(tài),隨后進入誘導反應(yīng)階段,誘導結(jié)束后開始放熱反應(yīng),壓力逐漸下降,溫度繼續(xù)上升至CJ狀態(tài)。兩步反應(yīng)模型中的放熱反應(yīng)速率由參數(shù)kR控制,從圖1中可以看出,放熱階段的溫度和壓力變化會隨kR的增大加劇,然而并不影響其最終狀態(tài)。對于一維爆轟波而言,kR只控制放熱進程中溫度和壓力等狀態(tài)參數(shù)變化的快慢,對一維爆轟波的靜態(tài)ZND結(jié)構(gòu)影響較小。
圖1 穩(wěn)態(tài)ZND爆轟結(jié)構(gòu)Fig.1 Steady ZND detonation structure
雖然一維爆轟的ZND結(jié)構(gòu)受kR的影響很小,但其動態(tài)傳播特性會隨kR的調(diào)整發(fā)生顯著變化。傳播特性可以通過記錄前導激波后的壓力變化來描述,具體表現(xiàn)為壓力在傳播方向上的振蕩。研究者根據(jù)不同的振蕩模態(tài),將一維爆轟波分為穩(wěn)定爆轟波和不穩(wěn)定爆轟波,后者也被稱為脈沖爆轟波。穩(wěn)定爆轟波的激波后壓力在傳播過程中不發(fā)生變化,而在不穩(wěn)定爆轟波中,壓力則表現(xiàn)為一定范圍內(nèi)的振蕩,且隨著kR的變化呈現(xiàn)出不同的振蕩模態(tài)。圖2顯示了kR對前導激波后壓力變化過程的影響。波后壓力在kR= 1.4時始終保持von Neumann狀態(tài)不變,此時為穩(wěn)定爆轟波。將參數(shù)調(diào)整到kR=1.8后,轉(zhuǎn)變?yōu)槊}沖爆轟波,其特征是波后壓力具有單一振幅的周期性振蕩。繼續(xù)增大kR,則演化成更加復雜的振蕩模態(tài),單個周期內(nèi)有兩個不同的壓力極大值,即所謂的雙周期振蕩。進一步增加kR,壓力振蕩的振幅明顯增加且周期性消失,具有高度不規(guī)則的特征,這意味著爆轟波在傳播過程中的振蕩行為變得不可預測。
(a) kR=1.4
利用相空間方法,可以描述壓力振蕩狀態(tài)的演變過程。圖3展示了根據(jù)圖2計算結(jié)果繪制的脈沖爆轟波的平面相軌跡,它可以反映傳播過程中前導激波后壓力的變化趨勢。kR=1.8時,壓力的變化具有固定的振幅和循環(huán)周期。kR=1.9時,壓力變化分裂成兩個不同振幅的循環(huán)周期,且兩個周期的長度不一致。kR=2.1時壓力的變化軌跡變得復雜,無顯著的周期性。
(a) kR=1.8,1.9
上述計算結(jié)果表明,在給定其他兩步反應(yīng)參數(shù)的情況下,隨著kR逐漸增大,前導激波后壓力的變化會從定值過渡為周期性振蕩,最終演變成高度不規(guī)則振蕩。為了更詳細地反映演化的過程,進一步計算了不同kR參數(shù)的一維爆轟波,并對各參數(shù)下的壓力峰值進行了統(tǒng)計,繪制出局部壓力峰值隨kR變化的分岔圖譜,如圖4所示。其中,局部峰值壓力用Pmax表示。該圖譜可以明確展示各振蕩模態(tài)相應(yīng)kR的取值區(qū)間。圖中的點表示取不同kR時統(tǒng)計的壓力峰值。一組kR參數(shù)內(nèi),出現(xiàn)兩個點則表示雙周期振蕩,多個點即多周期或不規(guī)則振蕩。
圖4 壓力峰值隨kR的變化Fig.4 Evolutions of local maximum pressure with kR
在放熱反應(yīng)速率較低的區(qū)間內(nèi)(kR≤ 1.45),前導激波后壓力在傳播過程中不發(fā)生變化,始終等于von Neumann狀態(tài)的壓力(Pvn=42.1),屬于穩(wěn)定模態(tài)。放熱速率加快(1.45 根據(jù)圖4的統(tǒng)計結(jié)果,認為kR<1.45都屬于穩(wěn)定爆轟波。為了研究穩(wěn)定爆轟在擾動影響下的振蕩特性,首先模擬了kR=1.0時一維穩(wěn)定爆轟波在連續(xù)變化密度擾動中的傳播,主要關(guān)注擾動波長對爆轟傳播動力學過程的影響,這種影響通過傳播過程中壓力振蕩模態(tài)的變化來表達。為了反映振蕩模態(tài)與擾動波長之間的關(guān)系,采用了與圖4相同的方法,對壓力峰值隨擾動波長的變化進行了統(tǒng)計,結(jié)果如圖5所示。擾動波長λ=0表示未施加擾動,此時穩(wěn)定爆轟的波面壓力不發(fā)生變化。根據(jù)振蕩模態(tài)的演化情況,將擾動對爆轟動力學的影響過程分為以下幾個階段:無影響階段(λ<90),表現(xiàn)為正弦擾動引起的單周期振蕩。調(diào)整波長只會略微改變振幅,并不會引起振蕩模態(tài)的變化。此過程與圖4中kR在1.45~1.86區(qū)間內(nèi)的變化趨勢相似;相互作用階段(90≤λ<1 000),原始的穩(wěn)定爆轟波開始響應(yīng)擾動波,因此表現(xiàn)出相對豐富的現(xiàn)象。擾動與爆轟波的相互作用過程,具有以下兩個特征:(1)隨著λ增大,振蕩的周期數(shù)先經(jīng)過倍周期分岔,后逐漸回歸單周期振蕩模態(tài)。此外,在λ=800附近,壓力振蕩的波峰處出現(xiàn)短暫的小幅高頻振蕩,因此在統(tǒng)計圖譜中出現(xiàn)了兩個峰值。小幅振蕩并不影響爆轟波的振蕩模態(tài)。(2)伴隨擾動波長的增加,峰值壓力的最大值經(jīng)歷了線性上升和下降;擾動支配階段(λ≥1 000),此階段振蕩模態(tài)被鎖定且振幅不再隨λ變化,擾動對爆轟傳播動力學的影響失效。計算獲得了與單步反應(yīng)不同的結(jié)果,發(fā)現(xiàn)對于穩(wěn)定的爆轟波,在特定波長范圍內(nèi),擾動也可以導致波后壓力振蕩。 圖5 擾動波長對局部壓力峰值的影響Fig.5 Influence of perturbation wavelength on local maximum pressure 為進一步觀察穩(wěn)定爆轟波的參數(shù)區(qū)間內(nèi),擾動對爆轟波動力學過程的影響趨勢是否一致,統(tǒng)計了kR=1.4時壓力峰值隨擾動波長的變化,如圖6所示。擾動波長λ<100時,激波后壓力只存在一個峰值,此時為單周期振蕩;λ增大,發(fā)現(xiàn)壓力峰值的數(shù)量顯著增多,說明波后壓力轉(zhuǎn)變?yōu)椴灰?guī)則振蕩。隨著λ進一步增大,壓力振幅越來越大,峰值壓力的最大值線性升高。當λ>550之后,振幅開始下降,最終演化為擾動支配的單周期振蕩。調(diào)整擾動波長,壓力振蕩的整體演化趨勢與kR=1.0時基本一致,但很顯然kR=1.4時擾動對爆轟動力學過程的影響更顯著。結(jié)果表明,擾動對一維穩(wěn)定爆轟波動力學過程的影響,與放熱反應(yīng)速率有關(guān)。反應(yīng)速率越快,爆轟波對擾動的響應(yīng)程度越大。 圖6 擾動波長對局部壓力峰值的影響Fig.6 Influence of perturbation wavelength on local maximum pressure 利用功率譜密度(power spectral density,PSD)方法,可以更直觀地反映整體演化趨勢的一致性。圖7展示了kR分別在1.0和1.4條件下,擾動波長λ=40,400,800對應(yīng)的功率譜密度。其中,λ分別為40和800時,有一個明顯的主頻信號,說明能量都集中在這個頻率上,因此導致振蕩具有很強的規(guī)則性,即只有一個峰值的單周期振蕩。需要注意的是,在λ=800的情況下,除主頻信號外,還存在一些倍數(shù)于基礎(chǔ)頻率的諧波,但這并不表示振蕩有多個周期。kR=1.0,λ=400時,主頻信號數(shù)量增加,此時為多周期振蕩。相同擾動波長條件下將kR調(diào)整到1.4,沒有明顯的主頻功率,這意味著能量在頻譜范圍內(nèi)的分布相對更廣。PSD結(jié)果也說明擾動對穩(wěn)定爆轟波動力學過程的影響只在適當波長范圍內(nèi)生效。 圖7 前導激波壓力的功率譜密度Fig.7 Power spectral density of the leading shock pressure 與kR=1.0不同的是,在kR=1.4時峰值壓力的最大值并未線性下降,而是在λ>600之后出現(xiàn)振蕩幅度的突然下降,這說明當擾動波長大于該值后,波前擾動對穩(wěn)定爆轟波的影響會被大幅削弱。圖8展示了擾動波長分別為600和650時,波后壓力的振蕩歷程??梢钥吹?,λ=600時波面壓力在振蕩的過程中會間歇性下降至Pvn的一半后重新上升,說明爆轟波在傳播的過程中,激波和化學反應(yīng)面會在擾動的干預下經(jīng)歷短暫的解耦,且重新耦合后會隨機出現(xiàn)壓力的突躍。擾動波長增大至λ=650后,這種隨機的壓力突躍消失,因此在圖6中表現(xiàn)為最大壓力的突然下降。對圖8的結(jié)果進行功率譜密度分析,如圖9所示??梢钥吹溅?600時能量在頻譜范圍的分布相對更廣。λ=650時頻率在0.01 附近存在一個非常明顯的主頻信號,說明此時振蕩的規(guī)則性更強。 (a) λ=600 (a) λ=600 前文在給定擾動幅值A(chǔ)=10%的條件下,研究了擾動波長的影響,發(fā)現(xiàn)對于不同的放熱反應(yīng)速率,穩(wěn)定爆轟波隨擾動波長的演化過程呈現(xiàn)出不同的特征。本節(jié)初步研究了擾動幅值的影響,計算結(jié)果如圖10所示。 (a) kR=1.0,A=20% kR=1.0,A=20%的擾動圖譜顯示,擾動對爆轟波的影響與圖5相比更加顯著,局部壓力峰值隨λ的演化過程與圖6相似。另外,注意到擾動波長從200增大到260時,不穩(wěn)定性并沒有增加,反而呈現(xiàn)出更穩(wěn)定的狀態(tài)。這與單步反應(yīng)中獲得的重新穩(wěn)定區(qū)相似,但這種現(xiàn)象在兩步反應(yīng)模型中是否具有普遍性,尚待進一步驗證。kR=1.0,A=5%時繪制的擾動圖譜,相較于圖6的結(jié)果,擾動對爆轟波動力學過程的影響被顯著削弱,表現(xiàn)為響應(yīng)區(qū)間縮小和振蕩模態(tài)的周期性變化。上述結(jié)果說明,擾動幅值越大,對爆轟波在擾動來流中傳播的影響越大。 現(xiàn)有的研究結(jié)果發(fā)現(xiàn),穩(wěn)定爆轟波和不穩(wěn)定爆轟波對擾動的響應(yīng)過程有明顯的區(qū)別[21]。前文已經(jīng)討論了一維穩(wěn)定爆轟的擾動振蕩特性,觀察到區(qū)別于單步反應(yīng)模型的現(xiàn)象。接下來討論擾動對一維不穩(wěn)定爆轟波的影響。圖11展示了給定擾動幅值A(chǔ)=10%的情況下,對不穩(wěn)定爆轟波施加不同波長擾動后統(tǒng)計的壓力變化結(jié)果。由圖2可知,kR=1.8無擾動情況下,前導激波后壓力的變化呈單一振幅的周期性振蕩。施加擾動后,在爆轟波的內(nèi)在不穩(wěn)定性和擾動波的共同作用下,前導激波后壓力波動增加,局部峰值壓力分布在更寬的范圍。擾動波長增大,壓力的整體變化趨勢與穩(wěn)定爆轟的結(jié)果相似。區(qū)別在于持續(xù)增大擾動波長λ,壓力振蕩行為也不會轉(zhuǎn)變?yōu)楹唵蔚闹芷谡袷帯R=2.1時,壓力隨λ的變化趨勢與低kR的結(jié)果不同。隨著λ的增大,逐漸演化為高頻振蕩的過程中,隨機出現(xiàn)一些振幅很高的低頻振蕩。 (a) kR=1.8 本文采用兩步誘導-放熱反應(yīng)模型分析了連續(xù)密度/溫度擾動對一維爆轟波傳播特性的影響,重點關(guān)注了擾動波長、擾動幅值對不同化學反應(yīng)活性的穩(wěn)態(tài)爆轟振蕩特征的影響規(guī)律,并與連續(xù)擾動作用下單周期振蕩和多周期振蕩模態(tài)的不穩(wěn)定爆轟波進行了對比。 對于穩(wěn)定爆轟波,擾動波長會在一定范圍內(nèi)觸發(fā)爆轟波的內(nèi)在不穩(wěn)定性,使前導激波后壓力呈現(xiàn)出多模態(tài)振蕩的特征,且振蕩模態(tài)隨著反應(yīng)物化學反應(yīng)活性的增加而趨于復雜。超過一定的波長(本文中波長臨界值接近800~1 000),爆轟波恢復為單峰振蕩的模態(tài)??傮w上,爆轟波的振蕩峰值壓力隨著擾動波長的增加呈現(xiàn)先增加后減小的趨勢。同時,觀察到在某些擾動波長下,爆轟波的周期性解耦-再起爆過程被削弱,導致振蕩峰值壓力出現(xiàn)突降。 通過對擾動幅值影響的分析,發(fā)現(xiàn)擾動幅值越大越容易觸發(fā)爆轟波的內(nèi)在不穩(wěn)定性。進一步對比擾動作用下單周期和多周期振蕩模態(tài)的結(jié)果可知,其振蕩模態(tài)更加復雜多變,外界施加的周期性擾動難以讓其恢復到初始的振蕩模態(tài)。因此,連續(xù)擾動作用下的爆轟波傳播特性主要是由兩方面決定:放熱區(qū)/激波面的耦合關(guān)系,以及擾動幅值和波長。前者決定了爆轟波的基礎(chǔ)振蕩模態(tài)和穩(wěn)定特性,后者通過對波前氣體干擾影響燃燒-流動的耦合關(guān)系,實現(xiàn)對爆轟波傳播模態(tài)的調(diào)整。 致謝本文得到了國家自然科學基金(11822202)資助。2.2 一維穩(wěn)定爆轟的擾動振蕩特性
2.3 擾動幅值對擾動振蕩的影響
2.4 一維不穩(wěn)定爆轟的擾動振蕩特性
3 結(jié)論