傳遞矩陣法在計算輸流管路高頻振動時的穩定性改進

曹銀行, 柳貢民, 胡 志

(1. 哈爾濱工程大學 動力與能源工程學院,哈爾濱 150001; 2. 哈爾濱工程大學 煙臺研究生院,山東 煙臺 264006)

在管路系統的流固耦合振動特性研究中,建立合適的數學模型和配套的計算方法是關鍵[1]。目前為止,常采用的計算模型包括基于一維梁及流體平面波理論的4方程模型[2-4]、8方程模型[5-6]和14方程模型[7-9]等。計算方法則分為時域方法和頻域方法[10]。其中頻域傳遞矩陣法(transfer matrix method, TMM)只需要將每個管路元件作為一個整體單獨處理,不需要根據其結構尺寸進行空間離散,涉及的自由度少[11]。此外,TMM還具有易于編程,計算效率高等優點。

多篇文獻提及,基于14方程模型和TMM計算大跨度管路高頻振動響應時存在數值不穩定的現象,認為這嚴重限制了TMM的使用,有必要加以改進[12-13]。De Jong的研究表明,直管14方程模型中,軸向、橫向和扭轉運動可以相互解耦,而計算失穩現象只存在于橫向振動的計算中。有鑒于此,本文的研究只考慮管路在y-z平面內的橫向運動,不考慮其在x-z平面內的橫向運動以及軸向和扭轉運動。

Riccati傳遞矩陣法(riccati transfer matrix method,RTMM)將傳統TMM中狀態向量的傳遞轉變為Riccati矩陣的傳遞,將管路振動特性的求解從一個兩點邊值問題轉化為一點初值問題,在降低計算矩陣維度的同時有效提高了數值解的穩定性[14]。但是RTMM可能會引入奇異頻率和奇異截面兩個病態問題[15];同時RTMM必須根據邊界條件將狀態向量分為已知元素和數目相等的互補元素以確保Riccati矩陣是一個方陣;最后RTMM在計算過程中需要先求得Riccati矩陣的初始值。RTMM很難處理復雜邊界條件下的振動計算,如彈性支撐和多激勵邊界、柔性邊界等。

全局傳遞矩陣法(global transfer matrix method, GTMM)是另一種常用的改善TMM計算穩定性的方法。為避免過大的“特征長度”,GTMM將大跨度的管路劃分為若干個子單元,而全局矩陣由各子單元的平行處理得到。由于劃分的子單元數目決定了全局矩陣的維度,并進而影響到GTMM的計算速度,De Jong提出了子單元劃分的“條件數”標準。但其建立的傳遞矩陣和條件數推導過程在考慮流動和壓力帶來的離心力和科氏力的影響后會將變得非常困難。

綜上所述,本文將重點討論TMM的改進方法,以使其能夠基于更先進的數學模型,穩定、快速地計算大跨度輸流管路系統的橫向振動。本文提出了基于無量綱化計算結果得到的子單元劃分準則的GTMM、混合能傳遞矩陣法(hybrid energy transfer matrix method,HETMM)和結合變精度算法的傳遞矩陣法(hybrid energy transfer matrix method,VPA-TMM)等三種改進方法并指明了其各自的優缺點。

1 De Jong的子單元劃分標準

De Jong采用的輸流管路面內橫向振動的頻域方程與Timoshenko梁的彎曲方程一致

(1)

(2)

(3)

(4)

符號和下標的含義如表1所示。

表1 符號和下標的含義Tab.1 Meaning of symbols and subscripts

式(1)～式(4)的特征方程為

(5)

在低頻段(ξ?1;χ?1),式(5)的特征值滿足

(6)

式中,kB和kN分別為傳播波和近場波。

直管傳遞矩陣的條件數與exp(kNL)處于同一數量級。在PC中,計算值分辨率為10-16,因此當kNL≥ln(1013)時,舍入誤差變得很重要(大于0.1%)并使計算結果出現不穩定的現象。在計算頻率上限不變的情況下,De Jong提出利用劃分子單元的方法降低單段管路的“特征長度”,對于這些子單元,其條件數都要小于一個預定的值。即

劃分后的子單元應利用GTMM平行處理。

2 無量綱模型得到的子單元劃分標準和GTMM

本章采用較簡單的拉普拉斯變換法獲取傳遞矩陣,并從計算結果的角度建立新的子單元劃分準則,而De Jong則從數值不穩定的原因出發,且實踐證明是更保守的。最后本節還將給出GTMM的具體實現措施。式(7)～式(10)如下

(7)

(8)

(9)

(10)

2.1 基于拉普拉斯變換的傳遞矩陣推導

包含式(7)～式(10)中四個獨立變量的時域狀態向量可表示為

之后式(7)～式(10)可以寫為矩陣方程的形式

(11)

式中,B,C為4×4的系數矩陣。令y(z,0)=04×1,式(11)可以通過拉普拉斯變換ζ轉換到頻域進行求解

(12)

用U和V分別表示A-1B的特征值矩陣和特征向量矩陣,并令v(z,s)=V-1Φ(z,s),式(12)可以簡化為

(13)

式(13)是一個一階微分方程,其解可以由式(14)直接給出。

v(z,s)=E(z,s)v0(s)

(14)

其中,

將v(z,s)=V-1Φ(z,s)代入式(14)可得

Φ(z,s)=VE(z,s)v0(s)

(15)

令z=0,則E(0,s)=I4×4,再代入式(15)可以得到v0(s)=V-1Φ(0,s)。進而式(14)可以轉化為

Φ(0,s)=VE-1(z,s)V-1Φ(z,s)

(16)

式中,VE-1(z,s)V-1=T為橫向振動的場傳遞矩陣。

基于拉普拉斯變換的傳遞矩陣推導避免了建立特征方程時的“消元”過程,因此更加簡單,不易出錯。同時它只需要處理一階微分方程式(13)。

2.2 TMM的高頻失穩現象

輸流管路橫向振動的TMM計算公式為

(17)

式中:Ds,De為具有相似形式的始末端邊界條件矩陣;Fs,Fe分別為作用于管道兩端的頻域激勵。

式(17)可以簡化為

DoveΦove=Fove

(18)

式中,Dove,Φove,Fove分別為整體傳遞矩陣、整體狀態向量和整體激勵向量。通過求解式(18)可以得到Φove,其第一部分是管路始端的狀態向量;任意一點的狀態向量都進而可以通過式(16)得到。

如圖1所示的Dundee管道為例,其參數如表2所示。 管路兩端自由,內部流體被無質量的堵頭封閉,所以流速等于0,初始壓力也為0。

圖1 Dundee管Fig.1 Dundee pipe

表2 Dundee管的材料與流體特性參數Tab.2 Material and fluid properties of dundee pipe

在B點施加一個單位力或力矩的脈沖激勵,計算得到的B點橫向振動和轉動響應如圖2和圖3所示。TMM計算結果在高頻段出現了明顯的計算失穩現象,這是因為傳遞矩陣VE-1(z,s)V-1=T中也有類似于exp(kNL)的元素,即exp[z/λ(s)],同時也包含那些數值很小的元素,所以計算結果數值不穩定的原因和通過子單元劃分構建全局矩陣的解決方法可以和De Jong的研究相同。

圖2 B點的橫向振動速度響應計算結果Fig.2 Calculated results of the transverse vibration velocity response at point B

圖3 B點的橫向轉動速度響應計算結果Fig.3 Calculated lateral torsional velocity response at point B

2.3 根據無量綱化計算結果得到的子單元劃分標準

其中,μ為管路材料泊松比。用上述無量綱參數替換式(7)～式(10)中的相應項可得到無量綱的橫向振動4方程模型

(19)

(20)

(21)

(22)

基于無量綱4方程模型的傳遞矩陣推導和輸流管路橫向振動計算與2.1節和2.2節相同,這里不做贅述。

圖4 B點的橫向轉動速度響應計算結果Fig.4 Calculated lateral torsional velocity response at point B

(23)

2.4 新的子單元劃分標準的驗證和GTMM的實現

圖5 7.2 m管路分析模型Fig.5 7.2 m pipe analysis model

圖6 E點的橫向轉動速度響應計算結果Fig.6 Calculated lateral torsional velocity response at point E

根據式(23)可知,只有最大的管道子單元長度取小于等于3.6 m的值時才能獲得1～1 000 Hz內穩定的計算結果;7.2 m長的管道應先劃分為兩個子單元CE=DE=3.6 m再進行GTMM計算,全局矩陣的維度為12×12,即

直接利用TMM計算式(17)與二單元GTMM計算結果的對比見圖6。如圖6所示,基于子單元劃分的GTMM明顯保證了TMM在高頻下的計算穩定性,同時也說明依據式(23)來劃分大跨度的管路是有效可行的。

圖5的算例中,流體速度和壓力均等于0,式(1)～式(4)和式(7)～式(10)可視為等價的,但根據De Jong的研究,1～1 000 Hz內TMM計算穩定時管路子單元的最大劃分長度為3.3 m,所以圖5中的計算模型應劃分為3個子單元并采用16×16的全局矩陣,這顯然是更為保守和復雜的,計算所需時間也會因子單元數目的增加而增加。

3 管路振動計算的混合能傳遞矩陣法

高強等人建立了一種求解層狀介質中波傳播問題的穩定方法,并將其命名為混合能矩陣法[17-18]。本節首次在流體管路系統橫向振動計算的TMM中引入混合能矩陣法的思想,給出了程式化的計算過程,解決了TMM高頻響應計算失穩的問題。首先仍然利用拉普拉斯變換的方法來獲取傳遞矩陣,但是這里取Φ′(z,s)=T′Φ′(0,s)。即

(24)

引入兩個混合能向量[qi,pi-1]T,[qi-1,pi]T,它們分別結合了管道子單元一端的速度變量及其另一端的力(力矩)變量,并且根據式(24)可得

(25)

式中,G,Q可分別看作動柔度矩陣和動剛度矩陣,根據式(25)還可以得到

(26)

結合式(25)和式(26)則可以得到

(27)

其中,

Mi-1,i+1=Mi,i+1(I2×2+Gi-1,iQi,i+1)-1Mi-1,i

(28a)

(28b)

(28c)

Ni-1,i+1=Ni-1,i(I2×2+Qi,i+1Gi-1,i)-1Ni,i+1

(28d)

重復使用式(25)～式(28),可以得到當大跨度的管路被劃分為n個子單元時,始末端混合能向量之間的關系滿足

(29)

對于圖1和圖5中的計算模型,A(D)點和B(E)點的管路自由邊界條件滿足

將邊界條件與式(29)結合可以得到

(30)

將Dundee管劃分為2～4個子單元時HETMM的橫向振動速度響應計算結果見圖7。同時如圖8所示,在最短的計算時間內,HETMM利用少量的子單元就可以得到在800 Hz之前與GTMM,RTMM和TMM一致,且在1～1 000 Hz內穩定的計算結果。

圖7 B點的橫向振動速度響應的HETMM計算結果Fig.7 HETMM calculation of the transverse vibration velocity response at point B

圖8 B點的橫向振動速度響應計算結果Fig.8 Calculation results of the transverse vibration velocity response at point B

4 結合高精度算法的傳遞矩陣法

De Jong的研究指出了TMM在高頻段計算不穩定的根本原因是由于計算機使用大約16位小數來表示一個實數而引起的舍入誤差問題。因此,TMM計算結果出現不穩定的現象是一個數學問題,而不是過大的“特征長度”這一物理問題。但是,無論De Jong本人還是后來的研究者都對長跨度的管道進行了一定程度的子單元劃分,都可以看作是從物理角度而非從根本上來解決TMM計算數值不穩定的問題。

Symbolic Math Toolbox提供了sym和vpa兩種方案,通過使用符號化或可變精度的數字來解決計算機的舍入誤差問題。兩種計算方案與計算機默認的雙精度(16位小數,即前文所述PC的計算值分辨率10-16)方案對于sin(π)的計算結果與誤差分析如表3所示。

表3 不同方案下對 sin(π) 的計算Tab.3 The calculation of sin(π) under different schemes

通過設置合理的字節數(這里選取為32字節),并以vpa格式輸入TMM計算程序中的相關參數,不需要對2.1節和2.2節中傳遞矩陣和TMM整體方程的推導和計算過程進行任何修改就可以得到穩定的結果,如圖9所示。

與GTMM和HETMM相比,最高可設置1 024位字節數的VPA-TMM可視為一種能夠一勞永逸地解決TMM長跨度高頻計算數值不穩定問題的改進方法,但其所需內存和計算時間也會大大增加。對于E點的振動響應計算,GTMM需要0.286 379 s,HETMM需要0.181 202 s。當利用vpa函數設置字節數為16～27時,VPA-TMM無法得到1～1 000 Hz內穩定的振動響應計算結果,而當設置28～35位字節時的計算時間如表4所示。

5 GTMM,HETMM與VPA-TMM的比較

GTMM,HETMM與VPA-TMM這三種方法都以傳統TMM為基礎,故而都能以與傳統TMM相同或相似的方式方便統一地處理各種邊界支撐和激勵條件。關于它們的其它優缺點的討論如下。

(1) 全局傳遞矩陣法GTMM

傳統的TMM經過多年的發展已經形成了較為完善的計算和理論體系,各類常見的管路系統元件,如直管、彎管、管路支撐、分支管等的傳遞矩陣模型和各類常見的邊界條件模型及其建立過程都可見于各種文獻,并進而可以輕而易舉地應用于GTMM。但隨著管路系統跨度或元件數量的增加,全局矩陣的形式要不斷更新。即使可以通過更好的子單元劃分標準,矩陣維度和計算時長的增加也是不可避免的。

(2) 混合能量傳遞矩陣法HTEMM

HETMM計算矩陣的維度小,計算時間最短,管路系統子單元數量增加時,計算矩陣的形式也能保持不變,使得計算程序也幾乎無變化,可以無成本地進行試算以確定合適的子單元劃分數目。但作為本文新引進的一種管路系統動力學特性分析方法,HETMM的發展還不如傳統的TMM成熟。HETMM需要對現有的傳遞矩陣和邊界條件模型進行進一步的改造,這是它的缺點,同時也似乎是它的優點。

(3) 變精度算法與傳遞矩陣法相結合的方法VPA-TMM

VPA-TMM可以從根本上解決TMM高頻段計算結果不穩定的問題,而且不需要改變傳遞矩陣的推導過程、形式和整體計算矩陣,但計算時間和成本明顯增加。雖然可以通過設置合理的字節數在一定程度上減少計算時間,但幅度也非常有限。VPA-TMM不適合大型管路系統動力學特性計算或管路系統優化設計這種需要大量樣本計算的情況。

關于3種方法的比較如表5所示。

表5 3種TMM改進方法的比較Tab.5 Comparison of 3 TMM improvement methods

6 結 論

本文提出了三種方法來解決TMM高頻計算結果不穩定的問題,即基于無量綱模型計算結果得到的子單元劃分標準的GTMM,HETMM和VPA-TMM。與RTMM相比,這三種方法都能程式化地處理各種邊界支撐和激勵條件,并且都是基于更全面地考慮了流體速度和壓力影響的4方程模型和更簡單的拉普拉斯變換方法來獲取的傳遞矩陣。

在這三種方法中,HETMM系首次從層狀介質中的波傳播計算擴展到管路系統的振動分析領域,可以直接以較為成熟的TMM為基礎。HETMM計算矩陣的維度和形式不隨子單元數目的增加而改變,計算時間又是最短的。可以認為HETMM是這三種方法中最好的一種。