基于直曲組集算法的復雜液壓管路固有頻率分析

韓 濤， 劉 偉， 張子駿， 趙通來， 劉永壽

(西北工業大學 力學與土木建筑學院， 西安 710129)

航空液壓管道的流致振動是液壓系統“跑冒滴漏”的重要原因之一。但航空管道系統走向、布局復雜。目前航空液壓管路的動力學分析方法存在以下兩個問題：① 主要針對單個管路(包括直管、曲管等)的流固耦合振動分析，難以反映管路局部系統或者整體系統的動態品質，難以滿足航空工程復雜管路分析的需求；② 傳統液壓系統動態特性分析將管路簡化為一維模型，偏重液壓系統功能性質分析，但是難以反映空間管道結構的響應和振動問題。隨著未來戰機減重和提高機動性與生存力的要求，液壓系統向高壓化發展，這樣對航空管道的強度和穩定性要求越來越高。因此，有必要對復雜管路系統流致振動的分析方法展開進一步研究。

輸流單管(包括直管、曲管等)的流固耦合振動問題已有大量的研究，Paisoussis等[1]推導出了輸流直管振動方程并對其進行穩定性分析。Housner[2]研究了兩端支撐直管的動力學特性，他采用Euler-Bernoulli梁模型，忽略重力、結構阻尼、外部拉壓力和流體壓力效應，不考慮流體黏性、可壓縮性，獲得了輸流直管的橫向振動控制方程。輸流曲管的研究過程中一直存在軸線可伸長和不可伸長兩種假定，其主要區別在于是否考慮軸向力作用。Chen[3-4]應用Hamilton原理建立了軸線不可伸長曲管的運動方程。Misra等[5-6]分別基于軸線可伸長和不可伸長假定對曲管平面內和平面外振動進行了分析。在此基礎上，很多學者對單個管路的流固耦合振動特性進行了深入研究，包括線性穩定性研究[7-8]、非線性振動研究[9-10]、模型的數值解法研究[11]以及管道的振動控制研究[12-13]等，為輸流單管的流致振動分析建立了較為完善的理論體系。

然而，在實際工程應用中，輸流管路系統的幾何形狀往往比較復雜。對于這些復雜管路，很難建立其振動分析的解析方法，因此一些數值方法和半解析方法逐漸被研究者采用，如有限元法[14-15]、動力剛度矩陣法[16-21]、傳遞矩陣法[22-23]、子結構法[24]等。其中有限元法通用性最強，但其計算量很大。動力剛度矩陣法的計算精度不依賴于單元數量，對于尺寸很大的單元它依然能夠保證精度，因而可適用于復雜管路的振動分析。Koo等采用Euler-Bernoulli梁理論建立了直管單元的動力剛度矩陣，并采用傳遞矩陣法推導了曲管單元的動力剛度矩陣，完成了復雜管路基于直管Euler梁模型的組集，為復雜管路的動力學分析提供了重要思路。隨后，Koo 等采用同樣方法對液態金屬爐熱管系統的動態特性進行了研究。Shen等[25]研究了由周期復合材料結構組成的復雜管路的減振問題，將管路中的曲管部分劃分為一系列直管單元來近似，基于直管Timoshenko梁模型完成組集。Dai等[26]在輸流直管控制方程中加入軸向組合力項，結合傳遞矩陣法建立管路系統動力剛度矩陣，并分析了流速、軸向組合力對單管和復雜管路頻率響應函數的影響。但是，已有的研究中在構建管路系統動力剛度矩陣時，都是基于單直管模型，即“直管-直管”組集，未曾實現“直管-曲管”組集[27]。另外，算法中對曲管部分常采用有限數目的直管單元來近似，再結合傳遞矩陣法來推導，過程比較繁瑣，還可能導致動力剛度矩陣維度增加。

因此，本文對曲管采用曲管單元的精確模型，提出復雜管路系統模態分析的直曲組集算法。首先基于Euler-Bernoulli梁理論，在局部坐標系下建立直管單元和曲管單元自由振動的離散模式及動力剛度矩陣；然后在全局坐標系中實現直曲組集，建立復雜管路的系統動力剛度矩陣和特征方程；最后采用所提算法分析管路布局對“Z”形管路固有頻率的影響規律，建立經驗公式，并加以試驗驗證。

1 輸流直管的動力剛度矩陣

輸流直管模型，如圖1所示。管道長度為L，管道各處橫截面相同，彈性模量為E，剪切模量G，單位長度質量為mp。流體相對管壁的流速恒定且為U，流體單位長度質量為mf。定義一個直角坐標系xyz，x軸沿著未變形的管道軸線，y軸，z軸均垂直于未變形的管道軸線。

圖1 輸流直管模型

基于以下假設：① 忽略管道材料阻尼影響；② 忽略重力影響；③ 流體無黏、不可壓縮，那么輸流直管的流固耦合振動方程可表示如下[28-29]

EAp?2wx?x2-mfU?2wx?x?t-(mp+mf)?2wx?t2=0

(1)

EI?4wk?x4+(mfU2+PiAi)?2wk?x2+2mfU?2wk?x?t+

(mp+mf)?2wk?t2=0, (k=y,z)

(2)

GJ?2φx?x2-Ix?2φx?t2=0

(3)

式中:wx為直管的軸向位移；wk為直管的橫向位移；φx為直管沿軸線方向的轉角；t為時間；Pi為內壓；Ai為流體的橫截面積；EAp,EI,GJ分別為管道的拉伸剛度、彎曲剛度和扭轉剛度。

式(1)～式(3)分別具有以下形式的解[30-31]

wx(x,t)=wx(x)exp(iωt)

wk(x,t)=wk(x)exp(iωt)

φx(x,t)=φx(x)exp(iωt)

(4)

式中：ω為圓頻率；將其代入振動方程式(1)～式(3)，消去時間項exp(iωt)，得到輸流直管在頻域內的振動方程，其解的形式可以表達為

wx(x)=Aexp(kax)

wk(x)=Bkexp(kbx)

φx(x)=Cexp(kcx)

(5)

式中：ka,kb,kc分別為軸向波數，橫向波數和扭轉波數，由下列頻散方程決定

(6)

(mp+mf)ω2=0

(7)

(8)

式(6)和式(8)有兩個根，式(7)有四個根，根據Euler梁理論，建立直管單元的離散模式，見表1。其中φy(x),φz(x)分別為沿著y軸，z軸方向的轉角；Nx(x)為軸向力；Ny(x),Nz(x)分別為剪力沿著y軸，z軸方向的分量；Mx(x)為沿著x軸方向的扭矩；My(x),Mz(x)為彎矩沿著y軸，z軸方向的分量。

定義直管單元兩端的位移狀態矢量Ws和力狀態矢量Fs如下

(9)

(10)

Cs=(A1A2By1By2By3By4

Bz1Bz2Bz3Bz4C1C2)T

(11)

式中：上標L為管道單元的左端；R為管道單元的右端；下標s為直管單元，結合表1和式(9)～式(11)，則有

Ws=W1sCs

(12)

Fs=D2sCs

(13)

由式(12)和式(13)得到輸流直管的動力學關系

Fs=DsWs

(14)

表1 輸流直管自由振動的離散模式

2 輸流曲管的動力剛度矩陣

輸流曲管模型，如圖2所示。管道的曲率半徑為R，管道各處橫截面相同流體單位長度質量為mf。定義一個曲線坐標系xyz，x軸與未變形的曲管軸線相切，y軸垂直于未變形的曲管中線且處于曲管平面內，z軸垂直于未變形曲管所在的平面。

圖2 輸流曲管模型

根據Paidoussis的研究，對于輸流曲管，如果不考慮管道材料阻尼的影響，管內流體無黏、不可壓縮，忽略重力的影響，那么其三維振動方程可表示如下

EI?4wy?s4+1R?3wx?s3+??s(AiPi-Nx)?wy?s+wxR+

1R(AiPi-Nx)+mfU2?2wy?s2+1R?wx?s+1R+

2mfU?2wy?t?s+1R?wx?t+(mp+mf)?2wy?t2=0

(15)

EI?4wz?s4-1R?2φx?s2-GJR?2φx?s2+1R?2wz?s2+

??s(AiPi-Nx)?wz?s+mfU2?2wz?s2+

2mfU?2wz?t?s+(mp+mf)?2wz?t2=0

(16)

EIR?3wy?s3+1R?2wx?s2-??s(AiPi-Nx)+1R×

(AiPi-Nx)?wy?s+wxR+mfU2R?wy?s+wxR-

mfU?2wx?t?s-1R?wy?t-(mp+mf)?2wx?t2=0

(17)

-GJ?2φx?s2+1R?2wz?s2+EIRφxR-?2wz?s2+Ix?2φx?t2=0

(18)

式中:wx,wy,wz分別為曲管的切向位移、徑向位移以及垂直于曲管平面的位移；φx為曲管軸線方向的轉角；s為自然坐標；t為時間；Ai為流體的橫截面積；Pi為流體內壓；管道的拉伸剛度、彎曲剛度、扭轉剛度分別用EAp，EI，GJ表示，則管道的軸向力Nx表達為

Nx=EApε

(19)

在軸線不可伸長的情況下，ε=0；在軸線可伸長的情況下，ε=?wx?s-wyR。不難看出：式(15)和式(17)為曲管的面內振動，式(16)和式(18)為曲管的面外振動。

式(15)～式(18)解的形式設為

wx(s,t)=wx(s)exp(iωt)

wy(s,t)=wy(s)exp(iωt)

wz(s,t)=wz(s)exp(iωt)

φx(s,t)=φx(s)exp(iωt)

(20)

將式(20)代回式(15)～式(18)，消去時間項exp(iωt)，可得到輸流曲管在頻域內的振動方程，略去非線性項，在軸線可伸長情況下得到輸流曲管的自由振動方程如下

EI?4wy(s)?s4+1R?3wx(s)?s3+AiPi?2wy(s)?s2+1R?wx(s)?s-

本文根據試驗結果，通過對不同條件下干燥時間的分析并結合文獻[14]得出辣椒干燥時間的評價方法，如表2所示。

EAPR?wx(s)?s-wy(s)R+mfU2?2wy(s)?s2+1R?wx(s)?s+

2iωmfU?wy(s)?s+wx(s)R-(mp+mf)ω2wy(s)=0

(21)

EI?4wz(s)?s4-1R?2φx(s)?s2-GJR?2φx(s)?s2+1R?2wz(s)?s2+

AiPi?2wz(s)?s2+mfU2?2wz(s)?s2+2iωmfU?wz(s)?s-

(mp+mf)ω2wz(s)=0

(22)

EIR?3wy(s)?s3+1R?2wx(s)?s2+AiPiR+mfUR×

?wy(s)?s+wx(s)R+EAp?wx(s)?s-wy(s)R-

iωmfU?wx(s)?s-wy(s)R+(mp+mf)×

(23)

-GJ?2φx(s)?s2+1R?2wz(s)?s2+EIR×

φx(s)R-?2wz(s)?s2-Ixω2φx(s)=0

(24)

式(21)～式(24)具有如下形式的解：

wx(s)=Aexp(kas)

wy(s)=A′exp(kas)

wz(s)=Bexp(kbs)

φx(s)=B′exp(kbs)

(25)

式中：ka為面內振動的波數，既包含軸向波數，也包含面內彎曲波數，由式(26)決定

deta11a12

a21a22=0

(26)

式中:

(AiPi+mfU2)/R2

kb為面外振動的波數，既包含扭轉波數，也包含面外彎曲波數，由式(27)決定

(27)

式中：

(mp+mf)ω2

從式(26)和式(27)解得面內振動包含六個波數，面外振動也包含六個波數，那么結合Euler梁理論，可以建立曲管單元的離散模式，如表2所示。表中系數αj,βj可結合式(26)和式(27)得到

(28)

(29)

結合表2，同樣可以得到

Wc=D1cCc

(30)

Fc=D2cCc

(31)

式中:Wc為位移狀態矢量，與式(9)相同，Fc為力狀態矢量，與式(10)相同，下標c為曲管單元，系數列陣為

Cc=(A1A2A3A4A5A6

B1B2B3B4B5B6)T

(32)

聯立式(30)～式(32)得到

Fc=DcWc

(33)

3 基于動力剛度矩陣法的直曲組集

對于復雜管路系統，其直曲組集計算過程，如圖3所示。按照“先分解再組集”的思路，將復雜管路系統分解為單管，依據Euler梁理論，建立單管的局部動力剛度矩陣，然后運用轉換矩陣建立全局坐標系下單管的動力學關系，再按照節點進行組集，建立起系統的動力剛度矩陣，結合邊界條件構建系統特征方程，最后進行求解。

以管路模型為例進行詳細說明，如圖4所示。首先將管路分解為單直管和單曲管，單管動力剛度矩陣的建立過程如前文所述，第i個管道單元的動力剛度矩陣記為Di，根據式(14)和式(33)，單管在局部坐標系下的動力學關系為

表2 輸流曲管自由振動的離散模式

Fi=DiWi

(34)

圖3 算法流程圖

建立全局坐標系x0y0z0和局部坐標系x1y1z1，x2y2z2,…,xiyizi,…,xmymzm，m為管路系統中直管單元的數目。那么依據局部坐標系xiyizi相對于全局坐標系x0y0z0的方向余弦矩陣ti，可以得到管道單元i的轉換矩陣Ti，將管道單元i的力狀態矢量和位移狀態矢量轉換到全局坐標系下

Fig=TiFi

(35)

Wig=TiWi

(36)

圖4 管路組集示意圖

將式(35)和式(36)代入式(34)，得到

(37)

式中：下標g為在全局坐標系下，由于轉換矩陣Ti為正交矩陣，故管道單元i在全局坐標下的動力剛度矩陣為

(38)

那么，對于單元i，其在全局坐標系下的動力學關系為

Fig=DigWig

(39)

將其分解為

(40)

同理，單元i+1的動力學關系可表示為

(41)

管道單元i與單元i+1的交點記為節點i，其受力如圖5所示。在工程應用中，直管單元與曲管單元在節點處相切或近似相切的情形比較常見，但也有不相切的情形，這兩種情形在全局坐標系下的受力分析可以得到統一，也就是說，無論節點i在局部坐標系下受力情況如何，都可以通過轉換矩陣將其所受力轉換為沿著全局坐標軸x0,y0,z0方向的分量，節點i處的力矩分析也是如此。

圖5 節點i受力分析圖

在全局坐標系下，根據該節點處的力平衡條件和位移連續性條件，可以得到

(42)

(43)

(44)

至此，便完成了“直管-曲管”的組集過程，依次類推，依據節點建立整個管路系統的動力學關系如式(45)所示

(45)

FK=-KW(ω)

(46)

K為彈性系數，節點i處受力平衡，得到

(47)

依據上式將式(46)疊加到式(45)中，即可得到含彈性支撐的系統動力剛度矩陣。

4 算 例

4.1 計算驗證

直曲組集算法可用于含有“直管-曲管”典型管道元件的復雜管路計算，但是目前關于這類管路的理論成果較少，缺乏對比性，而關于單管的計算結果較多。如果所提算法是正確的，那么必然也可適用于單管的計算。所以先采用該法對單曲管進行計算，并與文獻[18]進行對比，對所提算法進行初步驗證。

為了方便對比，管道材料參數及邊界條件設置與文獻[18]相同：楊氏模量E=210 GPa，截面慣性矩Ip=9.41×10-8m4，管截面面積Ap=4.05×10-4m2，流體截面面積Af=1.26×10-3m2，管道單位長度質量mp=3.18 kg/m，流體單位長度質量mf=1.26 kg/m，管道內壓P=10 MPa，流體速度V=15 m/s，曲管的約束條件設為兩端固支。計算結果，如表3所示。

從表3可知：直曲組集方法計算的一階、三階、五階頻率與文獻[18]給出的前三階頻率相近，相對誤差2%，但卻多出了兩階頻率，是參考文獻中所沒有的。這是因為文獻[18]中計算的是曲管的平面振動，而本文采用的是曲管的三維振動模型，不僅包含面內振動，也包含面外振動，為了進一步說明，求得曲管的前五階振型，如表4所示。可以更為清楚地看出，一階、三階、五階頻率是曲管面內振動頻率，二階、四階頻率是曲管面外振動頻率。采用所提算法求得的面內振動前三階頻率與參考文獻中平面振動前三階固有頻率相吻合，初步證明所提算法是正確的。

表3 曲管固有頻率對比結果

表4 曲管的前五階振型

為了進一步證明直曲組集算法可用于復雜管路，分別采用該法與有限元法對“Z”形管路進行計算。管路參數設置如下：楊氏模量E=200 GPa，泊松比υ=0.3，密度ρp=7.93×103kg/m3，管道外徑do=9.9 mm，壁厚t=1.2 mm，中間直管單元長520 mm，兩端直管單元的長度均為250 mm，曲管曲率半徑均為40 mm。管內流體密度ρf=898 kg/m3，內壓P=10 MPa，流速V=0.5 m/s，約束條件設為兩端固支。

直曲組集計算結果，如圖6所示，圖6中下尖點對應的橫坐標即為管路的固有頻率，分別用W1,W2,W3,W4,W5來表示。然后采用有限元法計算了不同單元數目下“Z”形管路的前三階固有頻率。對比結果，如圖7所示。

圖6 “Z”形管路的前五階頻率圖

圖7 兩種方法計算結果對比

從圖7可知，隨著劃分單元數目的增多，有限元法計算結果逐漸趨于穩定，并與直曲組集計算結果逐漸吻合，說明直曲組集算法是正確的。而且，直曲組集算法計算時僅用了5個單元，而有限元法在計算單元數目達到30個時才趨于穩定，說明直曲組集算法在保證計算精度的同時，減少了計算單元數目，實現了大尺寸單元組集計算。

4.2 管路布局對管路固有頻率的影響

管路布局對管路的振動特性有著重要影響，因為管路布局發生變化時，管路的固有頻率也會隨之發生變化，而固有頻率是判斷管路穩定性的一個重要參數。研究管路布局對管路固有頻率的影響規律，有助于管路設計時，提高管路的穩定性，避開激振頻率，預防發生共振等。所以，接下來以航空管路系統中常見的“Z”形管路為例，探討管路布局對其固有頻率的影響規律。

“Z”形管A端、B端位置固定，如圖8所示。定義兩個幾何參數：曲管至A端的距離l以及曲管的曲率半徑r。那么管路布局的變化可通過改變l和r的值來實現。為了方便尋找規律，采用控制變量法，即先給定曲管曲率半徑r，通過改變l來觀察其對“Z”形管路固有頻率的影響；然后再給定l，改變r值，觀察“Z”形管路固有頻率的變化規律。

圖8 “Z”形管路布局圖

管道材料參數設置同上，AB兩端點之間的橫向距離為600 mm，豎向距離為580 mm，管內流體密度為ρf=1 000 kg/m3，內壓Pi=21 MPa，流速U=12.48 m/s。根據航空管路設計標準，曲管的曲率半徑不小于管道外徑的4倍，即r≥4do，故取r=40 mm，分別計算了l=80 mm,150 mm,220 mm,290 mm四種工況下管路的前三階固有頻率，如圖9所示。

圖9 不同l值下管路的前三階頻率

從圖9可知，在曲管曲率半徑一定的情形下，隨著曲管位置距離A端越遠，即l值越大，“Z”形管路的前兩階頻率均隨之減小，第三階頻率先減小后增大。

一階固有頻率是判斷管路靜態失穩的重要參數，此處重點研究“Z”形管路的一階頻率隨l值的變化規律。定義一個比例系數α，α=l/yAB，yAB代表圖8中A端，B端的豎向距離，逐步改變l值，得到“Z”形管路的一階固有頻率隨α的變化，如圖10所示。

從圖10可知，“Z”形管一階頻率W1隨著比例系數α的增大先減小后增大，在α=0.5時達到最小值，這與管路的幾何形狀有關，管路A端、B端約束條件相同，當曲管距離B端的距離達到與原來距離A端的距離相等(例如l=430 mm與l=150 mm，圖10中α=0.741與α=0.259)時，管路幾何形狀是對稱的，故一階頻率是相同的，振型也是對稱的。觀察圖形可知：W1與α近似呈余弦函數關系，不妨假設

W1=Acosωα+θ+B

(48)

為了進一步確定式中待定系數，定義：“Z”形管路幾何形狀呈中心對稱(即α=0.5)時，管路的一階固有頻率為W0.5；α0=r0/yAB，r0表示給定的曲率半徑；當α=α0時，管路的一階固有頻率記為Wα0。那么結合圖10，可以將式(48)改寫為

(49)

圖10 一階固有頻率隨比例系數α的變化規律

至此，只需將一組計算數據代入式(49)，即可確定ω值，這樣就得到了“Z”形管路一階固有頻率關于比例系數α的表達式。為了驗證式(49)是否正確，接下來對其進行試驗驗證，如圖11所示。

對比結果，如圖12所示。式(49)與計算結果基本吻合，與試驗結果雖有誤差，但變化趨勢基本相符。而且通過數值分析發現，試驗結果與式(49)計算結果的相對誤差<5%，說明在誤差許可的范圍內，式(49)可用于“Z”形管路的設計計算。

圖12 不同比例系數下的結果對比

接下來研究曲管曲率半徑r發生變化時對“Z”形管路固有頻率的影響。給定l=290 mm，管道材料參數同上，分別計算了r=40 mm,80 mm,120 mm,160 mm四種工況下管路的前三階固有頻率，如圖13所示。

圖13 不同r值下管路的前三階固有頻率

由圖13可知，隨著曲管曲率半徑r的增大，“Z”形管路的前三階固有頻率均會隨之增大，一階頻率增大得較為緩慢，二階、三階增大地較為明顯。同樣，這里重點研究曲管曲率半徑對“Z”形管路一階固有頻率的影響規律。

選取不同r值進一步計算，“Z”形管路一階固有頻率的變化趨勢，如圖14所示。隨著曲管曲率半徑的增大，“Z”形管路的一階頻率增大速率越來越快。采用最小二乘法對計算結果分別進行二次多項式擬合和三次多項式擬合，發現擬合曲線近乎重合，故可近似認為“Z”形管路的一階固有頻率是曲率半徑r的二次函數，即

W1=a0+a1r+a2r2

(50)

對式(50)進行試驗驗證，分別選取了r1=40 mm,r2=80 mm,r3=120 mm,r4=160 mm,r5=200 mm五種尺寸進行試驗，結果如圖15所示。二次擬合結果與計算結果基本吻合，試驗結果與擬合曲線的變化趨勢基本相同，但相對擬合結果偏小，這是因為管道上貼有傳感器等附件，給管道增添了附加質量，所以導致試驗結果偏小。數值分析顯示，擬合結果與試驗結果的相對誤差<5%，如果考慮附加質量的影響，相對誤差會更小，說明擬合結果是正確的，在誤差許可的范圍內，式(50)可用于“Z”形管路的設計計算。

圖14 一階固有頻率隨曲率半徑的變化規律

圖15 不同曲率半徑下的結果對比

5 結 論

利用管單元的動力剛度矩陣，建立“直管-曲管”組合管道系統固有模態組集解法，通過驗證與計算，得到以下結論：

(1) 所提算法對曲管的計算結果與參考文獻吻合，相對誤差<2%。

(2) 所提算法對“直管-曲管”組合管道系統是適用的，實現了大尺寸單元組集計算，較有限元法減少了計算單元數目。

(3) 試驗顯示，“Z”形管基頻關于布局參數的經驗公式誤差<5%，證明所擬公式正確，在誤差許可范圍內，可用于“Z”形管的設計計算。

[1] PAISOUSSIS M P, LI G X. Pipes conveying fluid: a model dynamical problem[J]. Journal of Fluids & Structures, 1993, 7(2):137-204.

[2] HOUSNER G W. Bending vibration of a pipe line containing flowing fluid[J]. Journal of Applied Mechanics, 1952, 19(3):205-208.

[3] CHEN S S. Vibration and stability of a uniformly curved tube conveying fluid[J]. Journal of Acoustic Society of America, 1972, 51(Sup 1):223-232.

[4] CHEN S S. Flow induced in-plane instabilities of curved pipes[J]. Nuclear Engineering and Design, 1972, 23(1): 29-38.

[5] MISRA A K, PAIDOUSSIS M P, VAN K S. On the dynamics of curved pipes transporting fluid. Part I: Inextensible theory[J]. Journal of Fluids and Structures, 1988, 2(3): 221-244.

[6] MISRA A K, PAIDOUSSIS M P, VAN K S. On the dynamics of curved pipes transporting fluid. PartII: Extensible theory[J]. Journal of Fluids and Structures, 1988, 2(3):245-261.

[7] YANG Xiaodong, YANG Tianzhi, JIN Jiduo. Dynamic stability of a beam-model viscoelastic pipe for conveying pulsative fluid[J]. Acta Mechanica Solida Sinica, 2007, 20(4):350-356.

[8] QIAN Qin, WANG Lin, NI Qiao. Vibration and stability of vertical upward-fluid-conveying pipe immersed in rigid cylindrical channel[J]. Acta Mechanica Solida Sinica, 2008, 21(5):331-340.

[9] MENG Dan, GUO Haiyan, XU Sipeng. Non-linear dynamic model of a fluid-conveying pipe undergoing overall motions[J]. Applied Mathematical Modelling, 2011, 35(2):781-796.

[10] WANG L, NI Q. A reappraisal of the computational modelling of carbon nanotubes conveying viscous fluid[J]. Mechanics Research Communications, 2009, 36(36):833-837.

[11] LIANG Feng, WEN Bangchun. Forced vibrations with internal resonance of a pipe conveying fluid under external periodic excitation[J]. Acta Mechanica Solida Sinica, 2011, 24(6):477-483.

[12] YU Dianlong, WEN Jihong, ZHAO Honggang, et al. Vibration reduction by using the idea of phononic crystals in a pipe-conveying fluid[J]. Journal of Sound and Vibration, 2008, 318(1/2):193-205.

[13] 金基鐸, 楊曉東, 張宇飛. 固定約束松動對輸流管道穩定性和臨界流速的影響[J]. 振動與沖擊, 2009, 28(6): 95-99.

JIN Jiduo, YANG Xiaodong, ZHANG Yufei. Analysis of critical flow velocities of pipe conveying fluid under relaxation of boundary conditions[J]. Journal of Sound and Vibration, 2009, 28(6): 95-99.

[14] JUNG D, CHUNG J. In-plane and out-of-plane motions of an extensible semi-circular pipe conveying fluid[J]. Journal of Sound & Vibration, 2008, 311(1/2):408-420.

[15] OLSON L G, JAMISON D. Application of a general purpose finite element method to elastic pipes conveying fluid[J]. Journal of Fluids and Structures, 1997, 11(2): 207-222.

[16] TANG Bin. Combined dynamic stiffness matrix and precise time integration method for transient forced vibration response analysis of beams[J]. Journal of Sound and Vibration, 2008, 309(3):868-876.

[17] VIOLA E, RICCI P, ALIABADI M H. Free vibration analysis of axially loaded cracked Timoshenko beam structures using the dynamic stiffness method[J]. Journal of Sound & Vibration, 2007, 304(1/2):124-153.

[18] 李寶輝, 高行山, 劉永壽，等. 輸液曲管平面內振動的波動方法研究[J]. 固體力學學報, 2012, 33(3): 1-7.

LI Baohui, GAO Hangshan, LIU Yongshou, et al. In-plane vibration analysis of curved pipe conveying fluid with wave propagation method[J]. Chinese Journal of Solid Mechanics, 2012, 33(3): 1-7.

[19] KOO G H, PARK Y S. Vibration analysis of a 3-dimensional piping system conveying fluid by wave approach[J]. International Journal of Pressure Vessels & Piping, 1996, 67(3):249-256.

[20] KOO G H, YOO B. Dynamic characteristics of KALIMER IHTS hot leg piping system conveying hot liquid sodium[J]. International Journal of Pressure Vessels and Piping, 2000, 77(11):679-689.

[21] KOO G H, PARK Y S. Vibration reduction by using periodic supports in a piping system[J]. Journal of Sound & Vibration,1998,210(1):53-68.

[22] DONG M L, CHOI M J, OH T Y. Transfer matrix modelling for the 3-dimensional vibration analysis of piping system containing fluid flow[J]. Journal of Mechanical Science & Technology, 1996, 10(2):180-189.

[23] JUNG D, CHUNG J, MAZZOLENI A, et al. Dynamic stability of a semi-circular pipe conveying harmonically oscillating fluid[J]. Journal of Sound & Vibration, 2008, 315(1):100-117.

[24] ALDRAIHEM O J. Analysis of the dynamic stability of collar-stiffened pipes conveying fluid[J]. Journal of Sound & Vibration, 2007, 300(3/4/5):453-465.

[25] SHEN H, WEN J, YU D, et al. The vibrational properties of a periodic composite pipe in 3D space[J]. Journal of Sound and Vibration, 2009, 328(1):57-70.

[26] DAI H L, WANG L, QIAN Q, et al. Vibration analysis of three-dimensional pipes conveying fluid with consideration of steady combined force by transfer matrix method[J]. Applied Mathematics & Computation, 2012, 219(5):2453-2464.

[27] IBRAHIM R A. Overview of mechanics of pipes conveying fluids—Part I: Fundamental studies[J]. ASME Journal of Pressure Vessel Technology, 2010, 132(3):1-32.

[28] PAIDOUSSIS M P. Fluid-structure interactions: slender structures and axial flow[M]. Salt Lake City, UT:Academic Press, 2004.

[29] 包日東, 梁峰. 兩端一般支承裂紋管道的動力學特性[J]. 振動與沖擊, 2016, 35(7):220-224.

BAO Ridong, LIANG Feng. Dynamic characteristics of a cracked pipe conveying fluid with both elastically supported ends[J]. Journal of Vibration and Shock, 2016, 35(7):220-224.

[30] 梁峰，包日東.微尺度輸流管道考慮熱效應的流固耦合振動分析[J].振動與沖擊,2015,34(5):141-144.

LIANG Feng, BAO Ridong. Fluid-structure interaction of microtubes conveying fluid considering thermal effect[J]. Journal of Vibration and Shock, 2015,34(5):141-144.

[31] 王忠民, 鄒德志, 姜全友. 彈性地基上輸流管道主參數共振的主動振動控制[J]. 振動與沖擊, 2016, 35(4):182-187.

WANG Zhongmin, ZOU Dezhi, JIANG Quanyou. Active vibration control for principal parametric resonance of pipes conveying fluid resting on an elastic foundation[J]. Journal of Vibration and Shock, 2016, 35(4):182-187.