傳遞矩陣法在計算輸流管路高頻振動時的穩(wěn)定性改進(jìn)

曹銀行, 柳貢民, 胡 志

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

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

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

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

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

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

1 De Jong的子單元劃分標(biāo)準(zhǔn)

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

(1)

(2)

(3)

(4)

符號和下標(biāo)的含義如表1所示。

表1 符號和下標(biāo)的含義Tab.1 Meaning of symbols and subscripts

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

(5)

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

(6)

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

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

劃分后的子單元應(yīng)利用GTMM平行處理。

2 無量綱模型得到的子單元劃分標(biāo)準(zhǔn)和GTMM

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

(7)

(8)

(9)

(10)

2.1 基于拉普拉斯變換的傳遞矩陣推導(dǎo)

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

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

(11)

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

(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)。進(jìn)而式(14)可以轉(zhuǎn)化為

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

(16)

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

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

2.2 TMM的高頻失穩(wěn)現(xiàn)象

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

(17)

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

式(17)可以簡化為

DoveΦove=Fove

(18)

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

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

圖1 Dundee管Fig.1 Dundee pipe

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

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

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

圖3 B點的橫向轉(zhuǎn)動速度響應(yīng)計算結(jié)果Fig.3 Calculated lateral torsional velocity response at point B

2.3 根據(jù)無量綱化計算結(jié)果得到的子單元劃分標(biāo)準(zhǔn)

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

(19)

(20)

(21)

(22)

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

圖4 B點的橫向轉(zhuǎn)動速度響應(yīng)計算結(jié)果Fig.4 Calculated lateral torsional velocity response at point B

(23)

2.4 新的子單元劃分標(biāo)準(zhǔn)的驗證和GTMM的實現(xiàn)

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

圖6 E點的橫向轉(zhuǎn)動速度響應(yīng)計算結(jié)果Fig.6 Calculated lateral torsional velocity response at point E

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

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

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

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

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

(24)

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

(25)

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

(26)

結(jié)合式(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)

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

(29)

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

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

(30)

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

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

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

4 結(jié)合高精度算法的傳遞矩陣法

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

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

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

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

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

5 GTMM,HETMM與VPA-TMM的比較

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

(1) 全局傳遞矩陣法GTMM

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

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

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

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

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

關(guān)于3種方法的比較如表5所示。

表5 3種TMM改進(jìn)方法的比較Tab.5 Comparison of 3 TMM improvement methods

6 結(jié) 論

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

在這三種方法中,HETMM系首次從層狀介質(zhì)中的波傳播計算擴(kuò)展到管路系統(tǒng)的振動分析領(lǐng)域,可以直接以較為成熟的TMM為基礎(chǔ)。HETMM計算矩陣的維度和形式不隨子單元數(shù)目的增加而改變,計算時間又是最短的?？梢哉J(rèn)為HETMM是這三種方法中最好的一種。