基于改進伽遼金法解釋簡支輸流管顫振誤判

丁 明，范祖相，孟 帥

（1.上海交通大學海洋工程國家重點實驗室，上海 200240；2.上海齊耀動力技術(shù)有限公司，上海 201203）

引言

輸流管系統(tǒng)已被廣泛應用于航空航天、核工程、石油化工及海洋工程等領(lǐng)域。輸流管流固耦合問題可歸結(jié)為典型的無窮維連續(xù)陀螺系統(tǒng)動力學模型。Pa?doussis 在其著 作中詳 細闡述 了內(nèi)流效應［1］。內(nèi)流效應分為兩類：一類是保守內(nèi)流效應，輸流管在振動過程中不會從內(nèi)流獲取或損失能量。另一類是非保守內(nèi)流效應，輸流管可從內(nèi)流中損失能量（阻尼作用）或者吸取能量導致系統(tǒng)失去穩(wěn)定性。

輸流管內(nèi)流效應主要取決于內(nèi)流的方向、內(nèi)流是單相流還是多相流，輸流管形狀、邊界條件以及周圍環(huán)境等［1］。當前學者們的主要精力是將內(nèi)流效應理論應用于解決實際工程問題［2-12］。孫玉東等［4］以實際管路系統(tǒng)為研究對象，考慮內(nèi)流和管路結(jié)構(gòu)的相互作用，建立了液-管耦合振動噪聲動力學分析一體化有限元模型，分別對消聲壓力筒和管路撓性元件的噪聲、振動衰減作用進行研究。郭海燕等［5］基于內(nèi)流效應對海洋立管響應特性的影響，對現(xiàn)用的兩部經(jīng)典海洋立管設計標準（API RP2RD 和DNV-OS-F201）提出改進建議。陳正翔等［6］開展了實際工程常見的受多個彈性支座支承的輸液管系統(tǒng)穩(wěn)定性分析。張智勇等［7］進行了充液直管固-液耦合振動響應研究。齊歡歡等［8］采用了時滯主動控制方法對輸液管顫振失穩(wěn)進行控制以提高其臨界流速。張挺等［9］基于有限積分法研究瞬時關(guān)閥時輸流直管軸向耦合振動響應特性。段金龍等［10］分析了剪切流作用下頂張式海洋立管在不同內(nèi)流速度和密度下橫向渦激振動響應特性。鮑健等［11］分析了內(nèi)外流對細長海洋彈性管振動特性的作用機制。針對水動力段塞流誘導的柔性立管振動響應問題，高岳等［12］開展了段塞流誘導的懸鏈線型柔性立管模型振動響應測試。

內(nèi)流效應尚有眾多基礎科學問題未解決［1，13-27］。例如，對于自由端含有點質(zhì)量的懸臂輸流管，利用伽遼金法計算系統(tǒng)固有頻率時通常采用兩種方法：（1）將集中質(zhì)量作為邊界條件；（2）利用狄拉克函數(shù)將點質(zhì)量嵌入輸流管振動方程。當內(nèi)流速度為零時兩種方法計算結(jié)果一致。但是當內(nèi)流速度不為零時兩種方法計算結(jié)果顯著不同，此時第二種方法正確，該問題發(fā)生原因未可知［13］。對于空氣環(huán)境懸臂噴管，當內(nèi)流速度達到臨界流速時，輸流管會發(fā)生顫振而失去穩(wěn)定性。研究者通過現(xiàn)有內(nèi)流效應理論可很好地預測臨界內(nèi)流速度，但是將同樣方法應用在浸沒在靜水環(huán)境中的懸臂吸管時，理論上系統(tǒng)會在較小內(nèi)流速度下失穩(wěn)，但實驗中臨界內(nèi)流速度非常大［14-16］。學者們通過修正懸臂吸管的數(shù)學模型研究該問題。例如，Pa?doussis 等［15］提出在懸臂管的低端吸口處需考慮外流環(huán)境負增壓效應。Kuiper 等［16-17］研究發(fā)現(xiàn)外部水曳力對吸管穩(wěn)定性有著重要影響。Kuiper等［18］通過模型實驗發(fā)現(xiàn)，當內(nèi)流速度超過臨界值后，懸臂吸管將會呈現(xiàn)復雜的動態(tài)響應特征，且整個實驗中沒有觀察到平面振動。Giacobbi 等［19］采用解析法、數(shù)值計算和模型實驗對上述現(xiàn)象進行了深入的探討。盡管該問題尚未解決，Adiputra 等［20-21］基于現(xiàn)有內(nèi)流效應理論分析了海洋溫差發(fā)電裝置懸臂吸管系統(tǒng)的穩(wěn)定性。在新型吸管工程應用背景下，Butt等［22-23］對懸臂吸管系統(tǒng)在內(nèi)流和軸向外流聯(lián)合下的穩(wěn)定性進行理論及實驗分析。

對于兩端簡支輸流管系統(tǒng)，其內(nèi)流效應保守［1］。利用傳統(tǒng)伽遼金法線性分析時，在超臨界區(qū)預測到模態(tài)耦合顫振現(xiàn)象，但沒有能量輸入以維持該振動。Ch'ng 等［24］基于輸流管非線性模型發(fā)現(xiàn)，兩端固定特別是兩端簡支的輸流管道在有限（伽遼金）近似的數(shù)值積分中得出穩(wěn)定的極限環(huán)顫振運動是可能的。Holmes［25］基于有限維非線性模型得到的研究結(jié)果與線性模型一致，同樣預測到了顫振現(xiàn)象。其后，Holmes［26］基于有限維非線性模型進行無限維度分析，通過李雅普諾夫第二法評估輸流管道在超臨界區(qū)的穩(wěn)定性，研究證明顫振不能發(fā)生。Pa?doussis［1］和Sadeghi 等［27］基于非線性模型分析發(fā)現(xiàn)，系統(tǒng)在超臨界區(qū)發(fā)生屈曲失穩(wěn)，且隨著內(nèi)流速度進一步增加，輸流管屈曲響應幅值增加。Pa?doussis 對該研究歷程進行了總結(jié)［1］。馬騰等［28］通過在輸流管的兩端引入邊界支承彈簧，采用一種改進傅里葉級數(shù)的方法建立了輸流管耦合振動分析模型，依然預測到模態(tài)耦合顫振現(xiàn)象。于是產(chǎn)生了第一個疑問：既然對于非線性模型，采用不同求解方法會影響超臨界區(qū)響應預測，那么對于線性模型可否通過修正求解方法避免顫振誤判呢？考慮到伽遼金法是加權(quán)殘值法的一種，本研究嘗試從加權(quán)殘值法求解角度對線性模型產(chǎn)生顫振誤判進行探討。第二個疑問：既然非線性模型可以避免顫振誤判，探究線性模型發(fā)生顫振誤判的機理有工程意義嗎？作者認為該研究具有重要的工程應用背景。例如，隨著海洋開發(fā)步入深海區(qū)，海洋立管長徑比的大幅增加致使立管的柔性增強，柔性立管開始得到應用（尤其是在深海采礦領(lǐng)域），且新型柔性立管系統(tǒng)涌現(xiàn)（例如，為克服海床多變地形環(huán)境，有學者提出一種懸浮柔性輸流管系統(tǒng)［29］），柔性立管的內(nèi)流效應突顯，且可能進入超臨界區(qū)［30-32］。在海洋立管響應特性仿真計算中，因基于線性模型的伽遼金法求解簡便實用且可對激發(fā)的模態(tài)進行分析等優(yōu)勢被一直廣泛應用，尤其在試驗數(shù)據(jù)處理中很難采用非線性模型［33-36］。例如，在深海立管（內(nèi)部傳輸高速流體）渦激振動試驗數(shù)據(jù)處理中，一般先建立海洋立管的結(jié)構(gòu)動力線性模型，然后采用伽遼金法，利用實驗數(shù)據(jù)，分析海洋流體力以及激發(fā)的立管模態(tài)［33］。

1 顫振悖論

1.1 數(shù)學模型

簡支輸流管系統(tǒng)如圖1 所示，且采用二維坐標系oxy。原點o為輸流管左端點，y為橫向，x為未發(fā)生形變時的軸向，L為輸流管的長度。忽略重力效應，輸流管在恒預張力作用下的橫向振動控制方程為［1］：

圖1 簡支輸流管系統(tǒng)Fig.1 A simply-supported fluid-conveying pipe

式中M為單位長度內(nèi)流質(zhì)量；m為輸流管道單位長度的質(zhì)量；w為橫向形變；EI為輸流管道橫向彎曲剛度；C0為輸流管道阻尼系數(shù)；U為內(nèi)流單相流速；T為預張力；上角標表示對x的導數(shù)；上標表示對時間t的導數(shù)。

引入無量綱參數(shù)：

可得輸流管系統(tǒng)橫向振動無量綱振動控制方程為：

1.2 加權(quán)殘值法

基于加權(quán)殘值法［37］，可假設：

式中?j(ξ)為采用的試函數(shù)，試函數(shù)需要滿足正交以及邊界條件；N為基函數(shù)的數(shù)目；qj(τ)為對應的廣義時間坐標，得：

由于η(ξ，τ)只是近似解，設定η(ξ，τ)與精確解之間的殘差為R[η，ξ，τ]，即：

加權(quán)殘值法需要選定一組權(quán)函數(shù)ψi(i=1，…，N)，使得殘差的加權(quán)積分為0，即：

采用不同權(quán)函數(shù)ψi，加權(quán)殘值法有不同實現(xiàn)版本。若采用的權(quán)函數(shù)為基函數(shù)本身，則稱為伽遼金法。

1.3 伽遼金法

可采用簡支梁模態(tài)函數(shù)為試函數(shù)，即：

若采用伽遼金法，權(quán)函數(shù)和試函數(shù)相同即ψj=?j。利用簡支梁模態(tài)函數(shù)正交性條件：

最后得到輸流管系統(tǒng)結(jié)構(gòu)動力方程為：

式中M為質(zhì)量矩陣，且M=IN×N，IN×N為N階單位陣；C為阻尼矩陣，且C=C1+C2，C1為系統(tǒng)結(jié)構(gòu)阻尼矩陣，C1=csIN×N，cs為結(jié)構(gòu)阻尼系數(shù)，C2為內(nèi)流科里奧利力矩陣為剛度矩陣，Kij=(iπ)4+(iπ)2(Γ-u2)；q=(q1q2…qN)T。

矩陣C1為反對稱矩陣，故C為反對稱矩陣，可令即轉(zhuǎn)化為標準特征值問題：

式中λ為系統(tǒng)特征值當u≠0 時λ為復數(shù)。Ω實部Re(Ω)為固有頻率，虛部Im(Ω)與阻尼有關(guān)［1］。

利用MATLAB 自編程，已完成對輸流管道固有頻率計算程序的驗證（詳見參考文獻［30］）。參考文獻［1］中的某柔性簡支輸流管系統(tǒng)，設β=0.1，Γ=0，忽略結(jié)構(gòu)阻尼c=0，取N=4，系統(tǒng)固有頻率隨著內(nèi)流速度的變化曲線如圖2 所示，與文獻［1］中的前三階固有頻率數(shù)據(jù)吻合較好。隨著內(nèi)流速度的增加，系統(tǒng)的第一階固有頻率逐步減小。當u=π，Re(Ω1)=0，第一階模態(tài)發(fā)生屈曲失穩(wěn)。當u=2π，Re(Ω2)=0，第二階模態(tài)發(fā)生屈曲失穩(wěn)。隨著u的持續(xù)增加，第一階和第二階模態(tài)耦合顫振，于是產(chǎn)生一悖論：內(nèi)流效應保守，系統(tǒng)不會從內(nèi)流中吸取能量來維持顫振，為什么會發(fā)生顫振誤判？

圖2 輸流管系統(tǒng)固有頻率隨u 變化曲線Fig.2 The varying curves of the natural frequencies of the pipe with increase of u

2 悖論解釋

2.1 伽遼金法的局限

伽遼金法中，若采用N個試函數(shù)，則結(jié)構(gòu)動力方程（13）是由N個加權(quán)殘值方程構(gòu)成的，其中第i個加權(quán)殘值方程（上標N為采用的基函數(shù)數(shù)目）為：

當基函數(shù)增加到N+1 時，第i個加權(quán)殘值為：

對于任 意i(i=1，…，N)，根據(jù)式(16) 和(17)得：

假設基于伽遼金法求出的解收斂，則當N→∞時，對于任意的i(i=1，…，N) 有：

因此，采用伽遼金法隨基函數(shù)數(shù)目N的增加不存在收斂解，這顯然不滿足加權(quán)殘值法的收斂條件［37］。同時依據(jù)式的差異與內(nèi)流速度u成正比。因此，若u足夠大，采用伽遼金法求得的輸流管固有頻率是不可信的，甚至是錯誤的。

基于計算發(fā)現(xiàn)，隨著N增加不會消除耦合顫振誤判，同時隨著模態(tài)提高和內(nèi)流速度的增加，預測結(jié)果差別增大。舉例說明，設輸流管參數(shù)β=0.1，Γ=0，c=0，基于伽遼金法在N=4，5 時計算出的輸流管固有頻率隨內(nèi)流速度的變化曲線如圖3所示。

圖3 基于伽遼金法采用不同N 預測輸流管系統(tǒng)固有頻率隨u 的變化Fig.3 The predicted natural frequencies of the pipe vary with the increase of u at different N via Galerkin method

2.2 改進伽遼金法

可否找到合適的權(quán)函數(shù)ψi確保滿足加權(quán)殘值法收斂性條件呢？經(jīng)過驗算，提出權(quán)函數(shù)為：

式中?N+i(ξ)依然為簡支梁的模態(tài)函數(shù)；N仍然為式（3）中采用的基函數(shù)數(shù)目。利用正交性條件：

M為系統(tǒng)質(zhì)量矩陣，且M=IN×N；C為阻尼矩陣，且C=C1+C2，系統(tǒng)結(jié)構(gòu)阻尼矩陣C1=csIN×N。需要強調(diào)說明的是，利用式（22）內(nèi)流引入科里奧利力矩陣C2的元素為：

基于權(quán)函數(shù)和試函數(shù)正交性，科里奧利力項C2消失了；K為系統(tǒng)的剛度矩陣，且K=K1+(Γu2)K2，其中：

若Φ矩陣可逆，則式（26）可轉(zhuǎn)化為：

可以發(fā)現(xiàn)，式（31）由N個獨立加權(quán)殘值方程構(gòu)成。M，C，K矩陣的非對角元素皆為0，從而實現(xiàn)了系統(tǒng)結(jié)構(gòu)動力方程全解耦。其中，第i個加權(quán)殘值方程為：

觀察式(32)可知，加權(quán)殘值與基函數(shù)的數(shù)量N無關(guān)，這意味著當基函數(shù)的數(shù)量增加到N*(N*＞N) 時：

對任意i(i=1，…，N)皆成立。換言之，隨著N增加，利用式（31）求解能保證收斂。由于此方法屬于加權(quán)殘值法的一種，而且與伽遼金法采用相同的基函數(shù)，這里定義為改進伽遼金法。

忽略結(jié)構(gòu)阻尼，則可以直接求出第i階固有頻率為：

采用簡支輸流管的參數(shù)β=0.1，Γ=0，首先計算發(fā)現(xiàn)Φ可逆，然后取N=4 計算輸流管系統(tǒng)固有頻率隨流速的變化曲線如圖4 所示。隨著內(nèi)流速度增加，系統(tǒng)的第一階固有頻率逐步減少。u=π 時Re(Ω1)=0，發(fā)生第一階模態(tài)屈曲失穩(wěn)。u=2π 時Re(Ω2)=0，發(fā)生第二階模態(tài)屈曲失穩(wěn)。很明顯，在亞臨界區(qū)與文獻［1，26-27］數(shù)據(jù)較好吻合。但是在內(nèi)流超臨界區(qū)未發(fā)生模態(tài)耦合顫振誤判。同時發(fā)現(xiàn)，當增加試函數(shù)數(shù)目N時，計算結(jié)果具有一致收斂性且始終未發(fā)生模態(tài)耦合顫振誤判。舉例說明，當N=5 時，計算結(jié)果如圖4 所示，與N=4 時的仿真曲線較好吻合。

圖4 基于改進伽遼金法采用不同試函數(shù)數(shù)量（N=4 和N=5）預測輸流管系統(tǒng)固有頻率隨u 的變化Fig.4 The natural frequencies of the pipe vary with the increase of u at N=4 and N=5 via the modified Galerkin method

2.3 比較分析

內(nèi)流速度u=0 時，內(nèi)流產(chǎn)生的科氏力項消失?；谫み|金法和改進伽遼金法得到的結(jié)構(gòu)動力方程即式(13)和(26)都可以完全解耦，任一加權(quán)殘值方程與采用基函數(shù)數(shù)量N及輸流管系統(tǒng)參數(shù)無關(guān)。此時，兩種方法都能滿足加權(quán)殘值法的收斂性條件，采用兩種方法計算輸流管系統(tǒng)固有頻率的結(jié)果一致。

當內(nèi)流速度u＞0 時，基于伽遼金法得到的結(jié)構(gòu)動力方程因科氏力項為反對稱矩陣導致無法完全解耦。在本文提出的改進伽遼金法中，利用選取的權(quán)函數(shù)和試函數(shù)的正交性讓科氏力項為零，進而使結(jié)構(gòu)動力方程完全解耦。計算發(fā)現(xiàn)，基于改進伽遼金法可消除顫振誤判，同時發(fā)現(xiàn)，采用伽遼金法可能高估了系統(tǒng)的高階固有頻率，且隨模態(tài)階數(shù)增加和內(nèi)流速度加快更加顯著。

舉例說明，選取輸流管參數(shù)β=0.1，Γ=0，c=0，基于伽遼金法和改進伽遼金法計算出系統(tǒng)固有頻率隨內(nèi)流速度的變化曲線如圖5 所示。在u=0 時，兩種方法計算的結(jié)果一致，即頻率曲線起點相同，但隨著內(nèi)流速度的增加，第二階、第三階和第四階固有頻率開始分離，且第四階分離最為顯著。

圖5 基于伽遼金法與改進伽遼金法采用N=4 預測某個輸流管系統(tǒng)的固有頻率隨u 的變化Fig.5 The predicted natural frequencies vary with the increase of u based on Galerkin method and modified Galerkin method at N=4

3 結(jié)論

兩端簡支輸流管系統(tǒng)的內(nèi)流效應保守，即系統(tǒng)不會從內(nèi)流中獲取或損失能量。利用伽遼金法通過線性方法計算得到兩端簡支輸流管的固有頻率，在亞臨界區(qū)與試驗數(shù)據(jù)擬合較好，但是當內(nèi)流速度超過臨界流速后，預測到一階模態(tài)和二階模態(tài)耦合顫振現(xiàn)象。盡管學者們認定顫振現(xiàn)象不可能發(fā)生，但其發(fā)生的機理尚不清楚。考慮到伽遼金法屬于加權(quán)殘值法的一種，本研究基于加權(quán)殘值法對該悖論進行探討。

從線性角度看，內(nèi)流效應引入離心力、科氏力和慣性力三項。采用伽遼金法求解時，因科氏力項為反對稱矩陣，由加權(quán)殘值方程構(gòu)成的輸流管結(jié)構(gòu)動力方程不能完全解耦，這導致采用伽遼金法求解不符合加權(quán)殘值法收斂條件。提出一組新的權(quán)函數(shù)，可利用權(quán)函數(shù)和試函數(shù)的正交性使科氏力項消失，進而使輸流管系統(tǒng)結(jié)構(gòu)動力方程完全解耦，并滿足加權(quán)殘值法收斂條件。該方法屬于加權(quán)殘值法的一種，定義為改進伽遼金法?；诟倪M伽遼金法計算兩端簡支輸流管固有頻率，在亞臨界區(qū)與文獻［1，26-27］吻合較好，在超臨界區(qū)未發(fā)生顫振誤判，進而完成對模態(tài)耦合顫振悖論的解釋。同時發(fā)現(xiàn)，基于傳統(tǒng)伽遼金法可能高估輸流管系統(tǒng)高階的固有頻率，且隨著模態(tài)階數(shù)增加和內(nèi)流速度加快更加嚴重。由于固有頻率對于輸流管動態(tài)特性分析與穩(wěn)定性預判有著重要影響，本研究具有工程應用價值。例如，海洋立管正常作業(yè)時，仿真計算通常設定立管為兩端簡支。需要準確預測立管系統(tǒng)固有頻率以防止/抑制渦激振動［29］、躲避參激不穩(wěn)定區(qū)［30］等。

需要指出的是，式（26）轉(zhuǎn)化到式（31）的前提是Φ可逆。盡管通過MATLAB 計算發(fā)現(xiàn)，當N=1～100 時，Φ皆可逆。未來研究將對Φ的可逆性進行嚴格證明。同時，采用的權(quán)函數(shù)是基于簡支梁邊界條件提出的，并不適用于其他邊界條件輸流管系統(tǒng)（比如懸臂輸流管等）。但是本研究為實現(xiàn)各種邊界條件輸流管系統(tǒng)結(jié)構(gòu)動力方程的全解耦，尤其為處理內(nèi)流科氏力項提供了一種思路。同時，關(guān)于是否有其他因素影響顫振誤判需進一步研究。