曹銀行, 柳貢民, 胡 志

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

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

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

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

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

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

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

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

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

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

(7)

(8)

(9)

(10)

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

包含式(7)～式(10)中四個(gè)獨(dú)立變量的時(shí)域狀態(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)是一個(gè)一階微分方程,其解可以由式(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)避免了建立特征方程時(shí)的“消元”過程,因此更加簡單,不易出錯。同時(shí)它只需要處理一階微分方程式(13)。

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

輸流管路橫向振動的TMM計(jì)算公式為

(17)

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

式(17)可以簡化為

DoveΦove=Fove

(18)

式中,Dove,Φove,Fove分別為整體傳遞矩陣、整體狀態(tài)向量和整體激勵向量。通過求解式(18)可以得到Φove,其第一部分是管路始端的狀態(tài)向量;任意一點(diǎn)的狀態(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點(diǎn)施加一個(gè)單位力或力矩的脈沖激勵,計(jì)算得到的B點(diǎn)橫向振動和轉(zhuǎn)動響應(yīng)如圖2和圖3所示。TMM計(jì)算結(jié)果在高頻段出現(xiàn)了明顯的計(jì)算失穩(wěn)現(xiàn)象,這是因?yàn)閭鬟f矩陣VE-1(z,s)V-1=T中也有類似于exp(kNL)的元素,即exp[z/λ(s)],同時(shí)也包含那些數(shù)值很小的元素,所以計(jì)算結(jié)果數(shù)值不穩(wěn)定的原因和通過子單元劃分構(gòu)建全局矩陣的解決方法可以和De Jong的研究相同。

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

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

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

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

(19)

(20)

(21)

(22)

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

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

(23)

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

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

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

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

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

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

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

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

(24)

引入兩個(gè)混合能向量[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個(gè)子單元時(shí),始末端混合能向量之間的關(guān)系滿足

(29)

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

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

(30)

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

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

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

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

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

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

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

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

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

5 GTMM,HETMM與VPA-TMM的比較

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

(1) 全局傳遞矩陣法GTMM

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

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

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

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

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

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

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

6 結(jié) 論

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

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