段塞流輸流管道的穩(wěn)定性分析

吳明明，楊翊仁，李 鵬

（西南交通大學力學與航空航天學院，成都 611756）

引 言

輸流管道的動力學行為是一個典型的流固耦合現(xiàn)象，Paidoussis 等[1]認為輸流管道已經(jīng)成為結(jié)構(gòu)動力學和穩(wěn)定性研究的新模型，并歸納總結(jié)了輸流直管、曲管的線性及非線性動力學問題，進一步對管道穩(wěn)定性的研究做了系統(tǒng)化、理論化的闡述。

單相流輸流管道是早期學者們的主要研究對象。Chang 等[2]研究了兩端固支的輸流層合圓柱殼的固有頻率，發(fā)現(xiàn)其隨流速的增大而降低。張瑞平等[3]通過有限差分法，分析了軸向力對兩端支撐管道穩(wěn)定性的影響，認為軸向力的增大會導致失穩(wěn)臨界流速的增大。初飛雪[4]進一步發(fā)現(xiàn)，兩端支撐管道的固有頻率會隨著管內(nèi)壓力增大而減小，并且會隨管長的增大而明顯減小。張挺等[5]采用有限差分法對管道振動微分方程進行離散，研究了輸流直管在不同邊界條件下的振動及響應特性。Huang 等[6]使用伽遼金消元法得到了不同邊界條件下管道固有頻率的計算公式。易浩然等[7]基于哈密頓原理建立了含有集中質(zhì)量的懸臂輸流管道的動力學模型，研究了集中質(zhì)量對管道振動特性以及穩(wěn)定性的影響。周坤等[8]和Yu 等[9]研究了由不同材料組成的周期性管道的穩(wěn)定性和非線性動力學行為，認為當周期數(shù)達到某一臨界值時，管道的失穩(wěn)臨界流速將趨于定值。劉穎等[10]對粘彈性輸流管道的動力學微分方程進行了Laplace變換，研究了系統(tǒng)參數(shù)變化對管道自由端穩(wěn)態(tài)響應的影響。王忠民等[11-12]采用有限差分法研究了彈性地基上輸流管道的流速、地基密度等參數(shù)對管道固有頻率的影響。馬濤等[13]通過數(shù)值方法求解了功能梯度懸臂管道的振動方程，研究了材料的冪律指數(shù)、無量綱尺度參數(shù)及質(zhì)量比對管道穩(wěn)定性的影響。

與單相流近似穩(wěn)態(tài)的流動不同，段塞流管道中氣相、液相交替流動，成為了研究的難點。Wang等[14]研究了左端固支、右端簡支的段塞流管道的響應特性，認為流體速度影響著離心力和科氏力，液塞段出管時對管道的響應有很大影響。Liu 等[15-16]分別用伽遼金方法和有限單元法對懸臂和左端固支、右端簡支邊界下段塞流管道的固有頻率進行了研究，認為固有頻率既受液體表觀速度的影響，又受氣體表觀速度的影響。Khudayarov 等[17]建立了兩相段塞流的粘彈性輸流直管的動力學模型，發(fā)現(xiàn)增大氣塞長度會降低振動響應的幅值和管道的固有頻率，彈性地基密度會影響管道的失穩(wěn)臨界流速。Meng 等[18]研究了在外流作用下的兩相段塞流柔性立管的渦激振動問題，計算了立管的響應,并通過實驗驗證了理論模型的正確性。此外，孔令錢等[19]使用ANSYS Workbench 仿真軟件對流固耦合作用下的輸流管道進行了仿真，計算了管道的頻率及振型。

上述學者相繼研究了段塞流管道的頻率和響應，在穩(wěn)定性方面卻鮮有研究。因此，本文在段塞流管道固有頻率的研究基礎上，進一步對管道穩(wěn)定性的內(nèi)容進行補充。基于氣、液兩相段塞流管道的理論模型，采用伽遼金方法對系統(tǒng)振動方程離散，并進行數(shù)值求解，依據(jù)弗洛凱理論對系統(tǒng)穩(wěn)定性進行判定，討論各參數(shù)變化對管道固有頻率均方根及失穩(wěn)臨界流速的影響，通過理論分析對段塞流管道的穩(wěn)定性、可靠性進行評估，預防事故的發(fā)生。

1 氣液兩相段塞流管道理論模型

氣液兩相段塞流管道和段塞單元的模型如圖1所示，圖中藍色區(qū)域表示液體，白色區(qū)域表示氣體，L為管道長度，x方向為管內(nèi)兩相流體的運動方向，y方向為管道的橫向運動方向，V為兩相流流速，lu為段塞單元長度，lf、lg分別為段塞單元內(nèi)液塞段和氣塞段的長度。

圖1 氣液兩相段塞流管道及段塞單元模型

對管道結(jié)構(gòu)和管內(nèi)流體作如下假設：（1）管道具有較大的長細比，可用歐拉-伯努利梁模型來描述；（2）管道軸線不可伸長；（3）忽略管材及流體的阻尼；（4）管內(nèi)氣體、液體定常不可壓縮，且保持相同的運動速度；（5）忽略重力、軸向壓力的影響。

基于上述假設及文獻[14]，可以得到兩相流管道的振動方程如下：

式中：E為管道彈性模量；I為管道的截面二次矩；S為管道橫截面面積；w(x,t)為管道橫向振動位移；x為管道沿軸向到入口處的距離；t為時間；ρp為管道密度；ρ(x,t)為管道內(nèi)單位長度上兩相流體的密度，其值根據(jù)時間t和空間x的不同，要么等于液體密度ρf，要么等于氣體密度ρg。

假設管內(nèi)存在n個完整的段塞單元，以管道右端出口處的段塞單元為研究對象，考慮該單元從管道內(nèi)完全流出的過程，依次取圖2 所示t1～t5的5 個不同時刻。t1時刻，該段塞單元中的氣塞段準備出管；t3時刻，段塞單元中的氣塞段完全出管，液塞段準備出管；t5時刻，整個段塞單元完全出管。t2是介于t1和t3的任意時刻，此時段塞單元中的氣塞段部分出管；t4是介于t3和t5的任意時刻，此時段塞單元中的液塞段部分出管。t5時刻后管道中兩相段塞流的運動狀態(tài)又將回到t1時刻，即上述過程構(gòu)成段塞流系統(tǒng)的一個運動周期，換言之，該系統(tǒng)是一個周期時變的線性系統(tǒng)。

圖2 不同時刻的管內(nèi)流動狀態(tài)

圖2 中，x0表示在一個運動周期內(nèi)管道入口與第一個完整段塞單元之間的距離，其表達式為：

假設t1為初始時刻，即要求t1= 0，根據(jù)已知的兩相流流速V、段塞單元長度lu和氣塞長度lg，可以得到t3、t5的表達式分別為：

選擇簡支梁的振型函數(shù)φ為基函數(shù)，并采用伽遼金方法對式（1）進行離散，通過對振動方程在管道長度L上進行積分，便可得到積分形式的氣液兩相段塞流管道的橫向振動方程。

在一個運動周期內(nèi)，當t1≤t

將上述方程改寫成矩陣形式的表達式為：

其中：M(t)、C(t)、K(t)分別為系統(tǒng)中與時間相關的質(zhì)量矩陣、阻尼矩陣和剛度矩陣，q、q?、q?分別為系統(tǒng)的廣義位移、廣義速度和廣義加速度。

定義如下矩陣和向量：

將其代入式（7）中，即可得到狀態(tài)方程的標準形式：

令式（10）的解為：

其中：ω=ωr+iωi，ωi表示管道的固有頻率，ωr表示與阻尼相關的衰減項。將解代入式（10）后，得到系統(tǒng)的特征方程為：

通過上述分析，便可將動力學方程轉(zhuǎn)化為廣義復特征值的問題進行求解。

2 理論模型的驗證

在管道內(nèi)徑ri=5 cm，外徑ro=6 cm，L=15 m,E=210 GPa，ρp=7850 kg/m3，ρf=1000 kg/m3，ρg=1.293 kg/m3，lu=5 m，lf=3 m，V=10 m/s 的條件下，通過對管道振動方程的求解，可以得到管道的一階固有頻率f1和二階固有頻率f2隨時間的變化曲線如圖3 所示，圖中可見，管道固有頻率是隨時間周期變化的。

圖3 管道前兩階固有頻率隨時間的變化曲線

當段塞單元中的氣塞長度為零時，即段塞單元長度等于液塞長度，這表示氣液兩相段塞流模型完全退化為單相流管道模型，周期時變系統(tǒng)也變?yōu)槎ǔＯ到y(tǒng)，此時管道的固有頻率為恒定的常數(shù)。計算本模型在該條件下的固有頻率，并通過與Huang等[6]關于單相流管道系統(tǒng)固有頻率的計算結(jié)果進行對比，即可完成本文中氣液兩相段塞流理論模型的驗證。

在系統(tǒng)參數(shù)保持不變的條件下，控制兩相流流速V從2 m/s 增加至34 m/s，計算得到本文模型與Huang模型[6]的一階固有頻率f1隨流速V的變化曲線如圖4 所示，圖中可見，二者的計算結(jié)果高度吻合，驗證了本文理論模型的正確性。

圖4 兩種模型的一階頻率隨流速的變化曲線對比

3 臨界流速預測與頻率參數(shù)研究

對于兩端簡支的單相流輸流管道，通常會基于線性定常系統(tǒng)的穩(wěn)定性準則對管道進行穩(wěn)定性分析。考慮到圖3 中的一階頻率的周期波動幅值在0.0001以下，因此在對頻率的精度要求不高時，完全可以將一階頻率視為常數(shù)，故可以對本文的周期時變系統(tǒng)進行擬定常化的分析，缺點在于會產(chǎn)生微小的誤差。為了便于分析，考慮用單周期內(nèi)的頻率均方根值來表征系統(tǒng)的固有頻率。

假設系統(tǒng)參數(shù)lu、lf、L、ρg、ρf、ρp的參考值分別取1 m、0.5 m、10 m、1.293 kg/m3、1000 kg/m3、7850 kg/m3，式（12）的特征值實、虛部的均方根隨兩相流流速的變化曲線如圖5 所示，圖中可見，系統(tǒng)的一、二階頻率均隨流速的增大而降低，而且一階頻率會率先降低至零。依據(jù)線性定常系統(tǒng)的穩(wěn)定性準則，可以推斷當流速約為83.3 m/s 時，特征值的一階虛部均方根等于零、一階實部均方根大于零，表明管道系統(tǒng)發(fā)生失穩(wěn)。

圖5 系統(tǒng)特征值虛、實部均方根隨兩相流流速的變化曲線

通過改變某一參數(shù)，其他參數(shù)分別取對應的參考值時，得到管道的一、二階頻率f1和f2隨參數(shù)的變化曲線如圖6所示。由圖6（a）、圖6（d）可以看出，段塞單元長度和氣體密度的變化對管道一、二階固有頻率的影響較小；從圖6（b）～6（f）中可以發(fā)現(xiàn)，管道的一、二階固有頻率均隨液塞及管道長度、液體及管道密度的增大而明顯降低，且當管道長度增大到某一臨界值時，頻率均方根會逐漸降低至零。

圖6 管道固有頻率隨系統(tǒng)參數(shù)的變化曲線

4 穩(wěn)定性及其參數(shù)研究

4.1 穩(wěn)定性判定

周期性時變系統(tǒng)[20]的特征值是隨時間周期變化的，分析時需要將其轉(zhuǎn)化為定常系統(tǒng)，再根據(jù)定常系統(tǒng)的系數(shù)矩陣特征值模的大小來判斷原系統(tǒng)的穩(wěn)定性。

式（10）中系數(shù)矩陣E(t)的周期為T，且滿足E(t) =E(t+T)，假設該系統(tǒng)的標準基解矩陣為Z(T)，即要求Z(T)是如下微分方程組的解：其中I為單位矩陣。

引入李雅普諾夫變換：

其中L(t)為李雅普諾夫變換矩陣。將式（14）代入式（10），即可將周期時變系統(tǒng)轉(zhuǎn)化為定常系統(tǒng)：

式中H為該定常系統(tǒng)的系數(shù)矩陣。通過計算系數(shù)矩陣特征值的模，并將計算結(jié)果與數(shù)值1進行比較，即可完成原系統(tǒng)穩(wěn)定性的判斷。

當系統(tǒng)參數(shù)取參考值時，分別計算流速為83.2 m/s、83.3 m/s 和83.4 m/s 時矩陣H的特征值γ1和γ2的模，計算結(jié)果見表1。

表1 不同流速下矩陣H特征值的模對應的系統(tǒng)狀態(tài)

由表1 可見，當流速大于83.3 m/s 時，系統(tǒng)特征值的模有一個大于1，表明管道已經(jīng)發(fā)生了失穩(wěn)，這與第3 節(jié)中管道失穩(wěn)臨界流速的預測結(jié)果一致，驗證了通過固有頻率均方根來預測失穩(wěn)臨界流速的可行性。

4.2 穩(wěn)定性的參數(shù)研究

為了研究參數(shù)lu、lf、L、α、β、γ變化對管道穩(wěn)定性的影響，通過改變其中某一參數(shù)，其他參數(shù)分別取對應的參考值，得到管道失穩(wěn)臨界流速隨參數(shù)的變化曲線如圖7 所示。從圖7（a）與圖7（d）中可以發(fā)現(xiàn)，段塞單元長度和氣體密度變化對管道臨界流速的影響較小；從圖7（b）～7（e）中可以發(fā)現(xiàn)，臨界流速會隨著液塞及管道長度、液體密度的增大而呈現(xiàn)顯著下降的趨勢。從圖7（f）中可以看出，臨界流速不隨管道密度的改變而發(fā)生變化。

圖7 失穩(wěn)臨界流速隨系統(tǒng)參數(shù)的變化曲線

5 結(jié)論

基于氣、液兩相段塞流管道的理論模型，采用伽遼金方法對系統(tǒng)振動方程離散并進行數(shù)值求解，討論了參數(shù)變化對管道振動特性的影響，并運用弗洛凱理論確定了管道的失穩(wěn)臨界流速，分析了參數(shù)變化對段塞流輸流管道穩(wěn)定性的影響。主要結(jié)論如下：

（1）氣液兩相段塞流管道的固有頻率是隨時間周期性變化的，且波動周期與管內(nèi)流速和段塞單元長度有關。

（2）管道固有頻率的均方根隨流速、液塞及管道長度、液體及管道密度增大而減小，受段塞單元長度和氣體密度的影響較小。

（3）管道的失穩(wěn)臨界流速隨液體密度、液塞及管道長度的增大而減小，受段塞單元長度和氣體密度的影響較小，且不隨管道密度的變化而發(fā)生改變。