三維方腔內(nèi)導(dǎo)電流體第1次Hopf分叉

張敬奎，常家鵬，崔苗，李琦芬，任洪波，楊涌文

1.上海電力大學(xué) 能源與機械工程學(xué)院，上海 200090

2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析優(yōu)化與CAE軟件全國重點實驗室，大連 116024

近年來，在天體物理、熱核聚變、磁流體推進、電磁冶金、磁流體發(fā)電等領(lǐng)域中，磁流體動力學(xué)（Magnetohydrodynamics， MHD）理論獲得發(fā)展和應(yīng)用，其中流動轉(zhuǎn)捩問題受到廣泛關(guān)注和研究。自1883年雷諾［1］在著名的圓管流動實驗中發(fā)現(xiàn)流體存在截然不同的流動狀態(tài)開始，流體流動轉(zhuǎn)捩問題［2-3］便成為了流體力學(xué)研究的焦點問題。外加磁場作用下導(dǎo)電磁流體轉(zhuǎn)捩引起的流動不穩(wěn)定性，由于受多物理場耦合影響，更成為流體力學(xué)研究的難點。本文選取三維腔體內(nèi)磁流體Hopf分叉引起的不穩(wěn)定轉(zhuǎn)變開展研究。

在工程和自然領(lǐng)域，普遍存在非線性系統(tǒng)參數(shù)演化產(chǎn)生的自激振蕩，振蕩平衡點的失穩(wěn)則可能引起Hopf分叉。Hopf分叉理論在工程中具有重要應(yīng)用，例如，在航空航天領(lǐng)域，Hopf分叉臨界參數(shù)被用于預(yù)測超聲速飛行器壁板發(fā)生振顫的臨界速度和振顫頻率；在流體機械領(lǐng)域，Hopf分叉理論被應(yīng)用于壓氣機旋轉(zhuǎn)失速平衡點的判斷。在雷諾數(shù)Re較低的情況下，流體保持穩(wěn)態(tài)層流的流動狀態(tài)。增加Re至某一臨界值，穩(wěn)態(tài)的流動轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)周期性振蕩流動，稱為第1次Hopf分叉（the first Hopf bifurcation）。繼續(xù)增加Re至另一臨界值，使得非穩(wěn)態(tài)周期性振蕩流動轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)準(zhǔn)周期性振蕩流動，稱為第2次Hopf分叉（the second Hopf bifurcation）。再繼續(xù)增加至某一更高的臨界值，非穩(wěn)態(tài)準(zhǔn)周期性振蕩流動將轉(zhuǎn)變?yōu)橥牧?，稱為第3次Hopf分叉（the third Hopf bifurcation）。本文選取三維腔體內(nèi)導(dǎo)電流體由穩(wěn)態(tài)流動轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)周期性振蕩流動的不穩(wěn)定轉(zhuǎn)變開展研究。這種不穩(wěn)定轉(zhuǎn)變已在多項數(shù)值研究和實驗研究工作中（例如文獻［4-5］）觀察到，是由次臨界狀態(tài)下的第1次Hopf分叉引起，具有唯一主導(dǎo)的振蕩頻率。

近年來，眾多研究者對第1次Hopf分叉引起的不穩(wěn)定轉(zhuǎn)變開展了大量研究工作。大部分研究工作針對二維腔體內(nèi)流體流動的不穩(wěn)定性轉(zhuǎn)變，例如文獻［6-17］。這些研究工作主要采用了線性穩(wěn)定性分析方法，獲得了二維方腔內(nèi)第1次Hopf分叉失穩(wěn)對應(yīng)的臨界雷諾數(shù)Recr（范圍為7 402～10 500）和臨界無量綱角頻率（范圍為0.439 9～0.61）。盡管從第1篇預(yù)測二維腔體內(nèi)第1次Hopf分叉臨界參數(shù)的研究工作［18］至今，已持續(xù)30年，但如上述數(shù)據(jù)，臨界雷諾數(shù)和臨界振蕩頻率的預(yù)測結(jié)果仍具有很大分歧。近十多年，也有研究者開始針對三維腔體內(nèi)流體第1次Hopf分叉的不穩(wěn)定轉(zhuǎn)變開展研究。相對二維流動Hopf分叉不穩(wěn)定性的研究，三維流場下采用廣泛使用的線性穩(wěn)定性分析方法難于求解特征值方程，而采用直接數(shù)值求解則需要大量的網(wǎng)格節(jié)點和計算時長。因此，三維腔體內(nèi)Hopf分叉流動的研究結(jié)果也相對較少。

表1總結(jié)了三維方腔內(nèi)流體流動第1次Hopf分叉的臨界參數(shù)預(yù)測結(jié)果。從表中可以看到，現(xiàn)有研究獲得的三維流動第1次Hopf分叉失穩(wěn)的臨界雷諾數(shù)范圍為1 700～2 250，對應(yīng)的臨界無量綱角頻率范圍為0.472～0.651。文獻［4-5，19-21］給出了具體的臨界參數(shù)預(yù)測結(jié)果，結(jié)果差異不大。而文獻［22-24］則給出了較大差異的臨界參數(shù)預(yù)測范圍。二維和三維流動第1次Hopf分叉失穩(wěn)對應(yīng)的臨界雷諾數(shù)存在更大的差異。這是由于展向壁面邊界條件的影響，并且腔體內(nèi)第1次Hopf分叉失穩(wěn)對應(yīng)的臨界雷諾數(shù)與展向尺寸的長度呈負(fù)相關(guān)［24-25］。這也意味著，展向尺寸和展向壁面邊界強化了三維腔體內(nèi)流體流動的穩(wěn)定性。

表1 三維方腔內(nèi)由第1次Hopf分叉的臨界參數(shù)預(yù)測結(jié)果Table 1 Prediction of critical parameters by the first Hopf bifurcation in a three-dimensional square cavity

從當(dāng)前國內(nèi)外的研究工作中可以看到，研究者對二維和三維腔體內(nèi)流體第1次Hopf分叉失穩(wěn)的研究工作還不完善。由于采用的方法不同，不同研究者所預(yù)測的流動狀態(tài)發(fā)生轉(zhuǎn)變的臨界參數(shù)還沒有取得一致的結(jié)果。而針對磁場作用下導(dǎo)電流體第1次Hopf分叉失穩(wěn)引起的流動不穩(wěn)定性，目前還未有研究工作涉及。磁場作用下腔體內(nèi)導(dǎo)電流體由次臨界流動狀態(tài)下的準(zhǔn)穩(wěn)態(tài)流動轉(zhuǎn)變?yōu)橹芷谛苑欠€(wěn)態(tài)振蕩流動狀態(tài)的速度、振幅和頻率演化特性有待揭示。本文采用自己開發(fā)的配置點譜方法（Spectral Collocation Method，SCM）和人工壓縮法（Artificial Compressibility Method，ACM）相結(jié)合的高精度數(shù)值方法SCM-ACM［26-27］，直接求解三維流動控制方程，精確捕捉次臨界流動狀態(tài)下的速度演化過程，研究一定磁場強度對振幅衰減和振蕩頻率的影響規(guī)律，預(yù)測一定哈特曼數(shù)Ha條件下第1次Hopf分叉失穩(wěn)的臨界雷諾數(shù)，采用Fourier分析法獲得臨界無量綱角頻率，運用Richardson拓展法提升臨界參數(shù)預(yù)測的準(zhǔn)確度。

1 數(shù)學(xué)方程

1.1 人工壓縮格式的控制方程

本文同樣選取表1中研究者所采用的三維方腔頂蓋驅(qū)動流為對象，分析磁場作用下三維腔體內(nèi)導(dǎo)電流體流動的第1次Hopf分叉失穩(wěn)。如圖1所示，三維方腔的邊長為L，針對次臨界流動狀態(tài)下導(dǎo)電流體的三維流動，在X軸方向施加一定的磁場。引入人工壓縮法后的流動控制方程如式（1）和式（2）所示。人工壓縮法由Chorin［28］首次提出，適用于穩(wěn)態(tài)和非穩(wěn)態(tài)的不可壓縮或弱可壓縮流體流動的求解，具有形式簡單、易于實施和收斂的特點。

式中：p為壓力；u為速度矢量；ρ為流體密度；v為運動黏度；c為人工壓縮因子；t為偽時間變量；J為電流密度，J=ε(u×B)，其中ε為電導(dǎo)率；B為磁場。

采用式（3）的變換，對式（1）和式（2）進行無量綱化處理，獲得如下形式的無量綱化控制方程：

式中：Re為雷諾數(shù)；Ha為哈特曼數(shù)；U0為頂蓋移動速度；μ為動力黏度；B0為X軸方向施加的磁場；物 理 區(qū) 域 為Ω={E：(X，Y，Z)∈[0，1]×[0，1]×[0，1]}，這里E表示區(qū)域。

1.2 邊界條件和初始條件

三維方腔各壁面均設(shè)置無滑移邊界條件，頂部以恒定速度沿X軸方向移動，無量綱速度為1。為了獲得次臨界流動狀態(tài)下的準(zhǔn)穩(wěn)態(tài)結(jié)果，合適的初始條件對快速獲得準(zhǔn)穩(wěn)態(tài)結(jié)果至關(guān)重要。在低雷諾數(shù)條件下，例如當(dāng)Re=100時，除邊界條件外的速度和壓力初始值均給定0，即可以快速獲得三維穩(wěn)態(tài)流動的結(jié)果。然而，當(dāng)Re增至次臨界流動狀態(tài)時，這種初始條件很難在短時間內(nèi)計算獲得準(zhǔn)穩(wěn)態(tài)的結(jié)果，相關(guān)研究和分析參見文獻［21］。本研究通過逐次增加Re的方法，以前次Re條件下的速度和壓力結(jié)果作為增大Re后的初始條件，從而可以快速獲得次臨流動狀態(tài)的數(shù)值結(jié)果。

2 數(shù)值方法

2.1 空間離散方法

本文采用配置點譜方法SCM對人工壓縮格式的控制方程進行離散化處理。這種方法于1944被提出，從20世紀(jì)90年代開始逐漸進入主流的科學(xué)計算方法序列，目前已廣泛應(yīng)用于流體力學(xué)、傳熱學(xué)、電磁學(xué)等領(lǐng)域的數(shù)值求解，并被證明具有高精度和指數(shù)收斂的特性。作者已采用SCM與ACM相結(jié)合，開發(fā)了直角坐標(biāo)和柱坐標(biāo)系統(tǒng)下求解不可壓縮流動問題的數(shù)值方法SCMACM［26-27］，并檢驗了其求解二維和三維流動問題的精確性和收斂特性。采用SCM離散控制方程的主要過程如下：

選用切比雪夫-高斯-洛巴托配置點（Chebyshev-Gauss-Lobatto， CGL）［29］，計算所有配置點N中第i配置點ri的位置，以X方向為例，配置點位置如式（9）所示。配置點在標(biāo)準(zhǔn)譜空間取值，因此 需 要 將 實 際 物 理 區(qū) 域{E：(X)∈[Xmin，Xmax]}轉(zhuǎn) 化 為 譜 空 間 區(qū) 域{E：(r)∈[-1，1]}，如式（10），則空間偏微分項可轉(zhuǎn)化為式（11）所示的形式。

通用變量φ(r)可以通過所有配置點的級數(shù)進行近似，即

式中：h(r)為拉格朗日插值多項式，即

式中：T′N(r)為切比雪夫多項式TN(r)的一階偏導(dǎo)數(shù)。TN(r)=cos(icos-1(r))，i=0，1， …，N。

因此，可以獲得變量φ(r)在r=r0處的一階和二階偏導(dǎo)數(shù)，即

式中：D(1)和D(2)為CGL配置點的一階和二階系數(shù)矩陣，其表達式為

在獲得空間偏導(dǎo)數(shù)的系數(shù)矩陣后，即可整理出無量綱控制方程（4）～（7）的空間離散式，即

式中：A、B和C為X、Y和Z方向配置點Nx、Ny和Nz所對應(yīng)的一階系數(shù)矩陣D(1)；E、F和M為X、Y和Z方向配置點所對應(yīng)的二階系數(shù)矩陣D(2)。

2.2 時間離散方法

在時間離散處理上，本文選用四階龍格庫塔方法，這種處理方法對不穩(wěn)定的流動求解具有很好的適用性，由于論文篇幅限制，相關(guān)離散處理可參見文獻［26-27］。

3 結(jié)果與分析

3.1 有效性驗證

為了檢驗數(shù)值結(jié)果的有效性，本文分別選取了方腔頂蓋驅(qū)動流的數(shù)值模擬結(jié)果和實驗結(jié)果進行驗證，如圖2所示。在Re=1 000和Re=1 480的層流狀態(tài)下，本文計算結(jié)果與數(shù)值模擬結(jié)果都吻合良好。圖2（a）中的數(shù)值研究工作［30-33］均采用了高精度的數(shù)值方法，可以作為基準(zhǔn)解進行驗證。圖2（b）為本文計算結(jié)果與Liberzon等［22］開展的數(shù)值模擬和實驗研究結(jié)果的對比驗證。Liberzon等開展了極少被開展的方腔頂蓋驅(qū)動流實驗研究，他們分別以水和甘油為介質(zhì)，采用粒子圖像測速儀（PIV）測量了流體速度分布。本文結(jié)果與Liberzon等的數(shù)值模擬結(jié)果吻合良好，而與他們的實驗結(jié)果相比卻略有差異，這是由于實驗邊界條件的設(shè)置、測量誤差、振動和噪聲等因素所致。值得說明的是，在2種Re條件下本文僅用了363的節(jié)點數(shù)目，遠少于其他數(shù)值模擬所采用的網(wǎng)格數(shù)目。

3.2 磁場對速度演化的影響

在次臨界流動狀態(tài)下，依次計算了不同雷諾數(shù)和磁場條件下的速度演化。三維流場中，不同節(jié)點具有不同的速度和振蕩幅度，但具有相同的振蕩頻率和振幅衰減速率。為了準(zhǔn)確捕捉速度振蕩的演化過程和精確計算速度振幅的衰減，本文選取速度振幅最大區(qū)域的節(jié)點反映速度場振蕩信號的變化。這種選取一個節(jié)點的速度來反映速度場振蕩變化的方法也在其他研究工作中被采用，例如文獻［4-5，22］等。速度U的振幅最大，因此，所有計算工況均選擇U振幅最大區(qū)域的節(jié)點（0.33，0.328，0.5）作為監(jiān)測點。Ha=0時選取點速度U在Re分別1 900、1 910、1 915和1 920的時間演化如圖3所示。在次臨界流動狀態(tài)下，速度隨著時間的增加周期性振蕩，且振幅逐漸衰減。當(dāng)計算時間無限長時，速度將趨于穩(wěn)態(tài)。同時可以看到，隨著Re的增加衰減速率逐漸減小，可以預(yù)見，當(dāng)Re增加至某一臨界值時，速度隨著時間周期性振蕩的衰減速率將減小為0，此時腔體內(nèi)流動狀態(tài)由三維穩(wěn)態(tài)流動轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)周期性振蕩流動，即發(fā)生第1次Hopf分叉。然而，通過直接模擬獲得第1次Hopf分叉的臨界雷諾數(shù)是難以實現(xiàn)的，目前國內(nèi)外的研究多以振蕩衰減速率的變化預(yù)測發(fā)生第1次Hopf分叉的臨界雷諾數(shù)。

圖3 Ha=0選取點速度U的時間演化Fig.3 Time evolution of velocity U for selected point with Ha=0

選取Re=1 920的次臨界流動，在一定范圍內(nèi)增加Ha（0～5），研究磁場對次臨界振蕩流動特性的影響。選取點速度U在Ha分別為0、1、2和5的時間演化如圖4所示。在次臨界流動狀態(tài)下，速度U振蕩流動的衰減速度對Ha小范圍內(nèi)的變化響應(yīng)迅速。隨著Ha的增加，速度U振蕩的衰減加快，振蕩流動趨于更加穩(wěn)定，從而使得發(fā)生第1次Hopf分叉的臨界雷諾數(shù)增加。此外，每個Re條件下，趨于穩(wěn)態(tài)的平均速度U（絕對值）隨Ha的增加而顯著增加。這是由于磁場強度的增加抑制了主旋渦的流動強度，并使主旋渦中心發(fā)生位移，引起該點X方向速度U與速度矢量夾角減小所致。

圖4 Re=1 920時不同Ha條件下選取點速度U的時間演化Fig.4 Time evolution of velocity U for selected point for different Ha with Re=1 920

3.3 磁場對振幅衰減速率的影響

提取圖4中所有速度振蕩的振幅，在時間軸上的分布如圖5所示。極小的振幅隨時間呈指數(shù)形式衰減。次臨界流動狀態(tài)下的震蕩信號可以通過振幅隨時間變化的擬合函數(shù)描述，如式（22）所示。實部σ表示振幅的增長速率，可以通過式（23）計算獲得。σ＜0時對應(yīng)次臨界流動狀態(tài)，σ=0時為第1次Hopf分叉發(fā)生的臨界狀態(tài)，而當(dāng)σ＞0時為超臨界流動。虛部ω表示振蕩頻率，本文通過Fourier分析將時域振蕩信號轉(zhuǎn)化為頻域振蕩信號。

圖5 Re=1 920時不同Ha條件下選取點速度U的振幅演化Fig.5 Amplitude evolution of velocity U for selected point for different Ha with Re=1 920

圖5中可以看到，Re=1 920時不同Ha條件下，最大速度振幅（Ua）的衰減速率均為負(fù)值，為次臨界流動狀態(tài)下的結(jié)果。隨著Ha的增加，速度振幅的衰減速率急劇增加。Ha從0增加至5時，振幅衰減速率的絕對值從3.16×10-4快速增加至1.66×10-2。相同Re下，速度振幅衰減速率的絕對值隨Ha在一定范圍的增加呈嚴(yán)格拋物線形式增加，如圖6所示。磁場強度顯著抑制了次臨界流動狀態(tài)的速度振蕩，使振幅急速減小，流動更快地趨向穩(wěn)態(tài)流動。從而引起發(fā)生第1次Hopf分叉的臨界雷諾數(shù)隨Ha的增加而顯著增加。

圖6 Re=1 920時不同Ha條件下速度振幅的衰減速率Fig.6 Decay rate of velocity amplitude for different Ha with Re=1 920

3.4 磁場對無量綱角頻率和臨界雷諾數(shù)的影響

采用Fourier分析方法獲得的速度振幅頻域分布如圖7所示。在所有Ha參數(shù)下，速度振蕩均僅有唯一主導(dǎo)的無量綱角頻率，驗證了流動由穩(wěn)態(tài)轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)周期性振蕩流動是由Hopf分叉所引起。盡管磁場對速度振蕩具有顯著的抑制作用，然后一定范圍的磁場變化（Ha=0～5）對主導(dǎo)無量綱角頻率卻無影響，僅使主導(dǎo)角頻率所對應(yīng)的速度振幅有所減小。最大振幅對應(yīng)的無量綱角頻率均為0.575 2，與無磁場作用下的角頻率完全相同。

圖7 Re=1 920時不同Ha條件下速度振幅的頻域分布Fig.7 Frequency distribution of velocity amplitude for different Ha with Re=1 920

不同Ha參數(shù)條件，分別逐次增加Re，對次臨界流動狀態(tài)下的流動工況進行計算，獲得不同Re下速度振蕩的衰減速率（Ha=1），如圖8所示。衰減速率的絕對值隨Re的增加呈嚴(yán)格線性減小趨勢。繪制不同Re下衰減速率的曲線，即可獲得衰減速率為0時發(fā)生第1次Hopf分叉的臨界雷諾數(shù)Recr。采用Richardson拓展方法對不同Ha參數(shù)下的Recr進行改進，以獲得更為精準(zhǔn)的Recr，實施過程參見文獻［4，21］。圖9給出了一定Ha范圍內(nèi)（Ha=0～5）發(fā)生第1次Hopf分叉的Richardson拓展方法改進后的Recr。無磁場作用下，發(fā)生第1次Hopf分叉的Recr已在表1中獲得驗證，其中文獻［21］中的結(jié)果即為本文無磁場作用下的Recr。隨著Ha的增加，Recr顯著增加，且呈拋物線增加趨勢。這源于速度振蕩衰減速率的絕對值隨Ha的增加呈拋物線形式增加。隨Ha從0增加至5，Recr從1 916.6以拋物線形式增加至2 040.1。

圖8 Ha=1時不同Re條件下速度振蕩的衰減速率σ及臨界雷諾數(shù)Recr預(yù)測Fig.8 Decay rate of velocity oscillation for different Re with Ha=1， and prediction of critical Reynolds number Recr

圖9 不同Ha條件下的臨界雷諾數(shù)RecrFig.9 Critical Reynolds number Recr for different Ha

4 結(jié) 論

本文采用課題組開發(fā)的高精度數(shù)值方法SCM-ACM，通過逐次增加Re并更新初始條件的方法，直接求解了次臨界流動狀態(tài)下導(dǎo)電流體的控制方程，精確預(yù)測了Hopf分叉失穩(wěn)的臨界參數(shù)。本文預(yù)測三維Hopf分叉失穩(wěn)的方法，避免了求解系統(tǒng)矩陣特征方程，克服了三維系統(tǒng)線性分析方法特征根求解困難的問題，相關(guān)方法和結(jié)果能夠應(yīng)用于超聲速飛行器、風(fēng)機、電磁冶金等領(lǐng)域的工程設(shè)計和設(shè)備運行控制。通過分析次臨界流動狀態(tài)下的流動特性，以及磁場對振蕩速度、振幅衰減、頻譜分布和Recr的影響規(guī)律，獲得如下結(jié)論：

1） 速度振蕩幅度的衰減速率對較小范圍內(nèi)的Ha變化響應(yīng)異常迅速，磁場對速度振蕩具有顯著的抑制作用。隨著Ha在一定范圍（Ha=0～5）內(nèi)的增加，速度振幅的衰減速率呈拋物線形式增加，振蕩流動趨于更加穩(wěn)定。

2） 磁場變化對次臨界流動狀態(tài)下速度振蕩的無量綱角頻率無影響。不同Ha條件下，速度振蕩具有唯一主導(dǎo)的無量綱角頻率，且均為無磁場作用下的無量綱角頻率0.5752，反映磁場作用下流動狀態(tài)由穩(wěn)態(tài)流動轉(zhuǎn)變?yōu)橹芷谛哉袷幜鲃右餐瑯邮怯蒆opf分叉所引起。

3） 隨著Ha的增加，發(fā)生第1次Hopf分叉的Recr顯著增加，且呈拋物線增加趨勢。