Rijke管熱聲振蕩的穩定性切換行為研究?

采用數值方法模擬了強弱兩種阻尼條件下傳熱遲滯時間對一維Rijke管熱聲系統穩定性的影響,發現Rijke管系統存在穩定性切換現象.在推導了無量綱形式的管內聲波動量方程和能量方程之后,利用Galerkin方法對控制方程進行展開并在時間域內數值求解.分析了強阻尼和弱阻尼條件下,給定熱源的Rijke管熱聲振蕩的穩定性與傳熱遲滯時間的關系.結果顯示:在兩類阻尼條件下,持續增大傳熱與速度的遲滯時間,系統均呈現出穩定性切換現象,即系統在穩定和不穩定兩個狀態間持續轉變;但弱阻尼系統的不穩定區域寬于強阻尼系統的不穩定區域,系統最大振幅相對增大,且系統熱聲振蕩的主模態在不同模態之間發生轉換.最后,通過求解系統各階模態極限環幅值隨傳熱遲滯時間的變化,發現Rijke管熱聲振蕩穩定性切換現象與遲滯時間存在近似周期性關系.

1 引 言

熱聲振蕩是指燃氣輪機和航空發動機等設備由于不穩定熱釋放與壓力脈動耦合作用導致的低頻大振幅自激振蕩現象,通常還伴隨高分貝的低頻噪音.以燃燒領域為例,在燃氣輪機、燃氣發動機、固體火箭發動機等燃燒器中,當火焰面的熱釋放脈動與燃燒室的聲場之間形成正反饋機制時,就可能引發強烈的熱聲振蕩.這種不必要的振蕩常帶來噪音、熄火、工作點偏移和污染物排放等問題,對燃燒設備的安全和高效運行造成影響.另一方面,合理地利用熱聲現象,可制造熱聲發動機、熱聲驅動脈沖管制冷機等熱聲轉換設備,該類設備結構簡單,可靠性高、壽命長且環保性高,其相關研究引起了越來越多的關注.

1878年,Rayleigh首先對熱聲振蕩的產生機理給出描述,向一振蕩的氣團周期性地加入或取出熱量,所產生的效果取決于加熱或者散熱與振蕩的相位關系:當熱量在壓力最高點加入或者最低點取出,則振蕩加強;反之,振蕩減弱.至今,Rayleigh準則仍是被廣泛接受的熱聲振蕩現象產生和維持的合理解釋.1963—1983年,Rott發表了一系列文章,建立了經典線性熱聲理論,在理論上闡明熱聲效應中存在著熱和功的相互轉化,奠定了現代線性熱聲理論的基礎,該理論是目前熱聲研究中最有效、運用最為廣泛的理論.

瑞利準則和線性理論能夠闡述熱聲振蕩的產生和維持的機理,但卻無法描述起振、跳頻、遲滯以及聲壓飽和等非線性現象,這是因為熱聲振蕩是一個非常復雜的非線性問題,而線性理論是對小振幅弱非線性現象的近似,其中考慮了系統的線性,但忽略了熱源函數本身的非線性.于是,發展非線性熱聲理論來解釋此類現象,描述本質上非線性的熱聲自激振蕩的整個過程就顯得十分迫切[1].

Rijke管是研究熱聲振蕩最方便最典型的系統,國內外學者以Rijke管熱聲系統為模型展開了大量研究.關于Rijke管內的非線性熱聲振蕩現象,目前存在著幾種可能的解釋,包括但不限于非線性聲學效應、非線性對流換熱、熱聲非正交性以及非線性管口損失等四類原因.1990年,Heckl[2]對Rijke管中的非線性效應進行了理論和實驗研究,指出Rijke管中的非線性效應主要在于非線性對流換熱和非線性管口損失,前者在速度擾動和主流速度量級相當時作用明顯,導致換熱率下降,是振蕩幅值限制的關鍵原因;后者在壓力振動幅值很高時凸顯出管口損失增大的作用,但作用相對小得多.國內韓飛等[3,4]通過研究Rijke管熱聲相互作用的非線性和管口末端的非線性輻射聲阻,指出非線性效應的作用是限制振幅的增長和激發高階諧波的出現.2003年,Matveev[5]在其博士論文中指出:在大量的系統中,聲音強度導致的非線性聲學損失并不足以成為熱聲振蕩非線性飽和現象的主要原因.Balasubramanian和Subramanian等[6,7]首次提出Rijke管內的非正交熱聲現象,通過數值模擬總結出,即使在沒有阻尼的情況下非正交性也可以導致幅值飽和.

在非線性熱聲振蕩的研究領域內,馬大猷[8]對熱聲振蕩問題進行過系統研究,并根據瑞利準則推導出了詳細的Rijke管方程的嚴格解,給出了相應的管內非線性行波和駐波解.Yoon等[9]以Rijke型火箭發動機為對象,描述了非線性速度敏感型熱聲不穩定系統演化過程中的自舉(bootstrapping)現象,即系統第一階模態的能量先傳遞給第二階模態,激發第二階模態幅值增長,而第一階模態幅值下降;隨后第二階模態再將能量傳遞給第一階模態,導致第一階模態幅值隨之增長的現象.李國能等[10]在對一端開口一端封閉的Rijke型預混燃燒器的研究中發現熱聲不穩定的起振過程存在著頻率跳躍:系統首先激發低階的熱聲振蕩,然后低階熱聲振蕩逐步消退,激發起更高階的熱聲不穩定,依次類推,直到激發起適合當前燃燒器結構的穩定持續的壓力振蕩.黃鑫等[1]總結了Rijke型熱聲自激振蕩的研究進展,指出目前還沒有能完全解釋熱聲振蕩機理的理論,已有的理論只適用于弱非線性,應該在管內自激振蕩非線性現象以及建立并完善非線性模型兩個方面做進一步研究.2014年,Sayadi等[11]在其研究中指出,在熱釋放率較小時,系統主要振蕩頻率為線性化后不穩定模態的頻率,但當熱釋放率增大到一定程度時,系統將激發出其他高階頻率,而在線性分析中,這些頻率卻可能是穩定的.Kashinath[12]對燃燒G方程和聲波方程耦合求解,獲得了Rijke熱聲振蕩通向混沌的兩種途徑:倍周期分岔及Ruelle-Takens-Newhouse途徑,發現火焰皺褶和夾斷是產生周期性聲波的原因.2017年,Li等[13]研究了時間遲滯、聲學損失以及燃燒-流動耦合對Rijke管穩定性的影響,熱源選取1990年Fleifil等[14]提出的n-τ模型,根據加熱功率和系統阻尼將穩定性區間分為遲滯無關區域和穩定性切換區域.在遲滯無關區間,改變遲滯時間的大小,對系統穩定性沒有影響;但在穩定性切換區域,系統穩定性隨著遲滯時間的增大在穩定和失穩之間轉變.但該文以線性模型為基礎,而且僅選取了第一階模態,因此所得結果實際上應該是線性結果.

本文以一維水平電熱絲網加熱的Rijke管為研究對象,采用非線性熱源模型,從Navier-Stokes(N-S)方程出發,引入非線性傳熱模型和阻尼模型,推導出管內聲波擾動量控制方程組.采用Galerkin方法逼近控制方程并數值求解,取系統前十階模態進行計算,研究了在強阻尼和弱阻尼條件下,傳熱遲滯參數對系統穩定性的影響.結果發現,除公認的兩個重要因素熱源位置和熱源大小之外,阻尼系數和傳熱遲滯時間也可能對系統穩定性造成影響,不僅系統穩定性隨傳熱遲滯時間在穩定和不穩定狀態間切換,而且系統熱聲振蕩的主模態也和遲滯時間相關.最后,通過求解系統各階模態極限環幅值隨傳熱遲滯時間的變化,發現這種穩定性切換現象與傳熱遲滯時間存在周期性關系.

2 數學物理模型及求解

2.1 Rijke管模型及控制方程

采用如圖1所示的一維水平電熱絲熱源的Rijke管熱聲模型,兩端均為壓力開口邊界.其中,Rijke管長度為L,熱源位置,來流速度.

圖1 簡化的水平Rijke管模型Fig.1.Simplified Schematic of the horizontal Rijke tube.

模型的控制方程從N-S方程出發,假設氣體為理想流動,忽略氣體的熱傳導和黏性損失.由于Rijke管內空氣溫度范圍和壓強范圍介于240 K

采用聲流分析法,將變量分為穩態流動量和擾動聲學量部分,令p0,ρ0,u0及Q0分別為速度和熱源熱釋放量的平均量;分別為對應的擾動量,即

由于速度和速度擾動項相對音速都很小,因此,在推導過程中,含有速度和速度擾動項都忽略不計.得到關于各擾動項的方程組

對上述方程進行無量綱化,令:

其中x,t,u′,p′以及xf分別表示無量綱后的距離、時間、擾動速度、擾動壓力和熱源位置;a0為流動平均音速;M為平均流馬赫數.將(12)式代入方程(10),(11)中,得到無量綱形式的管內聲場動量和能量方程:

2.2 阻尼模型

Rijke管熱聲系統的阻尼損失主要分為兩部分,分別是邊界層損失和管口末端的聲能輻射損失,Howe[15]曾提出如下形式的阻尼模型:

式中ξ是管內總阻尼系數,j代表系統的第j階聲波模態,ξj是第j階聲模態的阻尼系數;A和P分別為Rijke管截面積和截面周長;ν和κ流體的動力黏性系數和熱擴散率,ωj=jπ為系統的第j階聲模態的無量綱頻率.令

則有

將阻尼模型引入方程(13),得到含阻尼的管內擾動方程

2003年,Matveev[5]在其博士論文中曾采用Howe模型對其實驗系統的阻尼進行計算并用于數值分析,根據Matveev的實驗參數可得其系統阻尼參數為c1=0.028,c2=0.0001,Subramanian等[7]曾使用該阻尼系數對Matveev的實驗模型進行數值模擬,給出了較好的對比結果.此外,Subram anian等[16],Juniper[17]以及Sayadi等[11]對熱聲不穩定的研究中也使用到Howe阻尼模型,且均取c1=0.1,c2=0.06.下文將采用這兩組阻尼參數對Rijke管熱聲系統的穩定性進行研究,并分別稱之為“弱阻尼”和“強阻尼”.

2.3 熱源模型

1914年,King提出了電加熱絲非定常熱釋放率的速度時滯模型,該模型雖然為非線性模型,但只能預測速度擾動幅值大于主流速度的非線性現象,即僅對熱源處發生回流的情況有效.1990年,Heckl[2]指出,K ing模型不符合擾動速度為主流速度的1/3時就出現非線性的實驗觀測結果,他以K ing模型為基礎提出了改進的非線性模型,稱為Heckl模型.Heckl模型的時均值及擾動速度小于1/3主流速度時的預測結果與King模型一致,當擾動速度大于1/3主流速度,更符合實驗的非線性結果.本文熱源函數采用Heckl模型,取熱源擾動項如下:

其中Lw,dw及Tw為電熱絲的長度、直徑和溫度.由于熱慣性的存在,傳熱和流場速度之間存在傳熱遲滯時間τ.

傳熱遲滯時間可利用Lighthill[18]提出的經驗公式進行計算,

從中可知對電熱絲網加熱的Rijke管,影響傳熱遲滯時間的因素為電熱絲直徑和管內氣體平均流動速度.

無量綱化的傳熱遲滯時間為

2.4 擾動方程的Galerk in逼近

擾動方程組為時間和空間的偏微分方程組,這里采用Galerkin方法將控制方程表達成頻域和時域函數的疊加,基函數的選取并不惟一,原則是必須滿足邊界條件,本文選取線性化的系統自共軛部分的特征函數作為基函數.對兩端開口的Rijke管,在x=0和x=1處,忽略管口損失,有p′(x,t)=0,?u′(x,t)/?x=0,可將聲波速度場和壓力場分別表示如下:

考慮到計算可行性,只能采用有限數量的模態數目.2014年,Selimefendigil和?ztopb[19]以熱機為研究對象發現系統前兩階和前十一階模態分別占據了97%和99.9%的流體動能.2017年,Sui等[20]在對Rijke型熱聲不穩定的實驗研究中指出,標準化的第一階特征值占總空間平均壓力脈動的92%,第二階模態是極限環振動的關鍵因素,其中第二、第三階的脈動能量占5%和2.3%.馮建暢等[21]在對Rijke管熱聲不穩定的分岔分析中發現:當取前九階模態和前十階模態時,得到的數值解差異可以忽略不計,說明了十階聲學模態的收斂性.因此在后續計算中,將取前十階Galerkin模態進行計算和分析.

將熱源方程(19)代入方程(18),利用(22)式和(23)式展開并投影到基函數系,得到時域內常微分方程組

其中j=1,2,3…N,N=10,且

3 計算結果與討論

從(21)及(22)式中可以得出:系統的可變參數包括熱源位置xf、熱源強度K、阻尼系數c1和c2以及傳熱遲滯時間τ.通常在實驗條件下熱源位置xf,熱源強度K可以精確調整和測量,獲得研究較多也較為深入;而系統阻尼和傳熱遲滯時間τ則難以精確測量和改變,因此前人研究較少.本文采用MATLAB軟件中求解延遲微分方程的dde23函數對方程進行數值求解,研究阻尼和遲滯參數對系統穩定性的影響.

3.1 傳熱遲滯對強阻尼系統穩定性的影響:穩定性切換

首先計算在強阻尼條件下傳熱時滯參數對Rijke管熱聲振蕩的影響,選取計算參數為:xf=0.25,K=0.8,c1=0.1,c2=0.06,初始條件為P1=1,Pj=0,?j=1及Uj=0,?j=1,…N.在此條件下,持續增大傳熱遲滯時間τ,觀察系統各階模態振蕩幅值隨時間演化情況,所得結果如圖2所示.由于系統壓力振蕩的高階幅值非常小,因此圖2僅顯示前三階壓力分量的變化情況.

圖2 強阻尼系統不同遲滯下的振蕩波形圖Fig.2.Time evolution of the heavily damped system with different time delay.

當τ非常小,如τ=0.05時,系統快速地收斂到穩定狀態,如圖2(a)所示.系統第一個臨界點為τ=0.119,越接近該值,系統衰減得越慢,圖2(b)在τ=0.118的條件下,當t>140后,系統才呈較明顯的衰減趨勢;當τ大于0.119后,系統穩定性發生改變,振蕩幅值不再隨時間衰減,而是維持在特定幅值振蕩,形成圖2(c)所示的極限環振蕩.此后系統的極限環振蕩幅值先隨τ增大,當τ=0.49時達到極大值,如圖2(d)所示;之后系統振蕩幅值隨τ的τ增大而減小,直到如圖2(e)中所示,系統都將處于極限環振蕩的狀態.而當τ增大到臨界點τ=0.94后,系統將隨時間逐漸趨于穩定,如圖2(f)示.繼續增大遲滯時間τ,系統達到平衡點的時間越來越短,圖2(g)中當τ=1.5時,系統在t=20范圍內就達到平衡.第三個臨界點為τ=2.01,在0.936<τ<2.01時,系統最終都將趨于穩定狀態,而τ>2.01后,系統進入新的不穩定區域,初始擾動會最終發展成周期性的極限環振蕩,如圖2(h)所示.并且在此之后,系統的穩定性依舊隨著遲滯參數τ的增大不斷地在穩定和不穩定之間切換.

這種系統穩定性隨著傳熱遲滯的增加不斷切換的現象稱為穩定性切換現象,該現象廣泛存在于各類時滯系統中.可用Rayleigh準則來解釋Rijke管內的穩定性切換現象:改變傳熱遲滯時間,即改變系統熱釋放與波動速度之間的相位差,亦改變了熱釋放與波動壓力的相位差.對系統任一階振蕩模態而言,在一個振蕩周期內,由熱源、阻尼以及其他模態作用共同導致的能量變化為正則振蕩加強;為負值則振蕩減弱;當達到平衡時,則該模態或處于穩定狀態,或處于特定幅值下的極限環振蕩狀態.

3.2 傳熱遲滯對弱阻尼系統穩定性的影響:穩定性切換及主模態改變

選取計算參數為:xf=0.25,K=0.8,c1=0.028,c2=0.0001,P1=1,Pj=0,?j=1及Uj=0,?j=1,…N. 在此條件下,持續增大傳熱遲滯時間τ,觀察系統各階壓力振蕩幅值隨時間的演化情況,所得結果列于圖3中.與上文一致,圖3只顯示了系統前三階壓力幅值的變化情況.

系統穩定性改變的第一個臨界點在τ=0.02.如圖3(a)所示,當τ<0.02時,系統隨時間演化最終收斂到平衡點;而當τ增大到0.02時,系統最終將進入圖3(b)所示的極限環振蕩狀態;繼續增大τ,從圖3(c)及圖3(d)可以發現,隨著第一階壓力幅值的增長,當τ=0.2時,第二階壓力幅值也不斷地增長甚至略微大于第一階幅值,系統振蕩的主模態從第一階變為第二階;繼續增大τ,系統的極限環幅值開始減小,且第二階模態幅值減小的速度要大于第一階幅值,圖3(e)中當τ=1.0時,第二階幅值減小為第一階幅值的1/4左右.

圖3 弱阻尼系統不同遲滯下的振蕩波形圖Fig.3.Time evolution of the weakly damped system with different time delay.

圖3 弱阻尼系統不同遲滯下的振蕩波形圖(續)Fig.3.Time evolution of the weakly damped system with different time delay.

當τ=1.2時,如圖3(f)所示,系統幅值并不呈單調遞增的趨勢,而是第一階模態先迅速衰減,然后第三階模態幅值迅速增大到約0.6,與圖3(d)中第二階幅值占主導的情況不同的是,圖3(d)中系統各階模態都是呈單調遞增的趨勢上升的極限環狀態,而圖3(f)經歷了一個第一階模態衰減而后第三階模態幅值上升到極限環振蕩的過程,可以推測當τ=1.2時,第一階和第二階壓力模態的相位與熱釋放相位不同,只有第三階模態相位與熱釋放相位相同.

在0.02<τ<1.27范圍內,系統振蕩的主模態在第一、第二及第三階之間切換,但一直處于不穩定的極限環狀態,而當達到第二個臨界點τ=1.27時,系統再次進入穩定區間,從圖3(g)和圖3(h)可以觀察到該轉變.

從圖3(i)和圖3(j)可看到,第三個臨界值為τ=1.73,在τ>1.73后,系統再次進入不穩定區間,此后系統一直處于極限環振蕩狀態,且振蕩主模態在第一階、第二階及第三階之間切換.系統第4個臨界點為τ=3.3,當τ>3.3后,系統再次進入穩定區間,圖3(k)和圖3(l)中給出了系統從不穩定到穩定的轉變.

3.3 強、弱阻尼系統穩定性對比

將兩種不同阻尼條件下系統穩定性區間進行對比,結果列于表1.

表1 兩種不同阻尼條件下系統的穩定性區間Table 1.Comparison of stability region of differently damped systems.

從表1可以發現,當阻尼從c1=0.1,c2=0.06減小到c1=0.028,c2=0.0001后,系統的第一個不穩定區間范圍從0.119<τ<0.93擴大為0.02<τ<1.27,第二個不穩定區間從2.01<τ<3.07擴大到1.73<τ<3.3,而穩定區間范圍則分別從τ<0.119,0.94<τ<2.01及3.07<τ<3.90縮小到τ<0.02,1.27<τ<1.73及3.3<τ<3.7,說明減小系統阻尼,系統的每個不穩定區域的范圍增大,穩定區域范圍縮小.這是因為阻尼越大,能量耗散越快,系統就越趨于穩定.

其次,對比圖2和圖3還可以發現:強阻尼條件下,不穩定區域內始終是第一階模態的幅值遠高于其他高階成分;而弱阻尼條件下,系統的主模態并不固定在第一階模態,而是在第一階、第二階和第三階之間轉換,這是由于所選取的阻尼模型對高階成分的衰減作用更加強烈導致的.

此外,弱阻尼條件下,當τ=1.2時,系統幅值并不呈單調遞增的趨勢,而是先經歷了一個低頻模態衰減,然后高頻模態快速增長達到飽和狀態的過程.而強阻尼條件下卻不出現該現象,這說明阻尼也是決定系統振蕩模態的一個重要因素.

以上討論表明,系統阻尼和傳熱時滯參數不僅影響到系統是否穩定,還影響到系統極限振蕩時的主模態.在之前的研究中,學者們一般認為系統主模態只和熱源位置有關,而該結果則說明,不僅僅是熱源的位置,系統的阻尼以及傳熱遲滯時間都是決定系統主振蕩模態的重要參數.

圖4和圖5分別給出了兩類阻尼下系統前三階模態對傳熱遲滯τ的分岔圖譜.其中不同顏色的點分別代表不同壓力模態在系統達到極限環振蕩或者穩定狀態時的峰值,黑色表示一階模態,藍色表示二階模態,紅色表示三階模態.從中可以發現系統穩定性與時滯參數存在近似周期關系,且變化周期約為2,這符合系統第一階模態的周期,即系統各模態中的最大周期.

圖4 強阻尼系統對時間遲滯τ的分岔圖Fig.4.Bifurcation plot of heavily damped system for variation of time lag(τ).

圖5 弱阻尼系統對時間遲滯τ的分岔圖Fig.5.Bifu rcation plot of weakly damped system for variation of time lag(τ).

圖4 和圖5中,當τ=0.49時,系統第一階模態的幅值達到極大值,根據瑞利準則可以推測,當τ=0.49時,第一階模態壓力振蕩和放熱脈動相位差最小.

對比圖4和圖5中的最大振蕩幅值可發現,在強阻尼條件下,系統振蕩最大幅值約為1,而弱阻尼系統的最大振蕩幅值為2.83.這是因為強阻尼系統不僅能量耗散的更多,且高階模態成分被抑制程度較高,因此,輸入系統的熱能轉化為聲能的部分更少,相對集中在一階模態.而弱阻尼系統能量耗散較少,且阻尼對高階模態的抑制作用也減小,因此,輸入系統的熱能更多的轉化為聲能,不僅第一階模態幅值變大,第二階和第三階模態得到的能量亦有所增多,幅值也相對增大.

圖4中系統第一階模態的幅值占據主要地位,遠大于第二階、第三階模態幅值,與圖2符合得很好.而圖5中的不穩定區域則分為兩類,以0—2區間內為例,當0.02<τ<1.05時,系統第一階、第二階和第三階模態都處于不穩定狀態,且第一階和第二階模態幅值量級相當,第三階模態幅值較小.而1.08<τ<1.27時,系統只有第三階模態處于不穩定狀態.這說明在該范圍內,只有第三階模態的壓力振蕩相位和熱釋放脈動相位相同.

4 結 論

以一維Rijke管系統熱聲振蕩為研究對象,通過分離擾動項的方法得到了管內聲場壓力和速度控制方程,采用Galerkin方法對控制方程進行了數值求解,分析了強阻尼和弱阻尼系統動力學特性與傳熱遲滯時間τ的關系,得出如下結論:

1)在給定的熱源位置和熱源強度下,對強阻尼系統和弱阻尼系統,增大熱源相對速度的遲滯時間,系統的穩定性都將在穩定和不穩定兩個狀態間轉變,即系統存在穩定性切換現象;

2)在給定的熱源位置和熱源強度下,弱阻尼系統不穩定區域大于強阻尼系統的不穩定區域,且由于更多的熱能被轉化為聲能,系統極限環振蕩最大幅值亦從1增大到2.83;

3)在給定的熱源位置和熱源強度下,強阻尼系統熱聲振蕩的主模態始終是第一階模態,而弱阻尼條件下,系統熱聲不穩定的主模態在第一階、第二階和第三階模態間轉換,這意味著當阻尼較弱時,即使在相同的加熱位置,系統也可能發生聲波頻率改變或者復頻率聲波的現象;

4)系統穩定性與時滯參數存在近似周期性切換關系,且變化周期約為2,與系統一階模態周期相等.