張 挺， 林震寰， 郭曉梅， 張 恒， 范佳銘

(1. 福州大學(xué) 土木工程學(xué)院， 福州 350116； 2. 臺(tái)灣海洋大學(xué) 河海工程系，臺(tái)灣 基隆 20224)

管道在輸流過(guò)程中，因內(nèi)部流體流動(dòng)會(huì)在管道邊壁上施加力的作用，造成振動(dòng)，并最終導(dǎo)致管道的長(zhǎng)期疲勞失效，如輸流直管的軸向耦合振動(dòng)[1]。而對(duì)于管道偏離中心軸線造成的橫向振動(dòng)，不同支撐條件對(duì)其振動(dòng)特性也有較大的影響。因此，對(duì)不同支撐條件下輸流管道動(dòng)態(tài)響應(yīng)的研究顯得尤為重要。

目前，輸流管道研究模型主要有梁模型和殼模型兩種，在管道內(nèi)徑遠(yuǎn)小于管道長(zhǎng)度時(shí)，采用梁模型來(lái)研究是簡(jiǎn)單而有效的。不同學(xué)者針對(duì)該模型下的輸流管道都有著廣泛研究，Paidoussis[2]基于梁模型，采用牛頓法對(duì)管道單元和流體單元進(jìn)行受力分析，得到較為完整的描述輸流直管橫向振動(dòng)微分方程，并研究了不同支撐條件下輸流直管的穩(wěn)定性問題；隨后，又有眾多學(xué)者對(duì)該方程進(jìn)行不斷改進(jìn)和不同角度的研究。

由于輸流直管橫向振動(dòng)微分方程中具有對(duì)空間和時(shí)間物理量的高階偏導(dǎo)項(xiàng)，因此其精確解的獲取較為困難，目前對(duì)輸流管道振動(dòng)響應(yīng)特性研究采用的數(shù)值分析方法主要有伽遼金法(Galerkin Method)[3-5]、微分變換法(Differential Transformation Method，DTM)[6-7]、微分求積法(Differential Quadrature Method，DQM)[8-9]、精細(xì)積分法(Precise Integration Method，PIM)[10-11]、廣義積分變換法(Generalized Integral Transform Technique，GITT)[12-13]等，這些都可以有效地模擬此類問題。

隨著計(jì)算機(jī)技術(shù)和數(shù)值方法的不斷發(fā)展，無(wú)網(wǎng)格法(Meshless Method)在近幾年有了快速的發(fā)展，其中廣義有限差分法(Generalized Finite Differential Method，GFDM)屬于區(qū)域型的無(wú)網(wǎng)格法，采用泰勒級(jí)數(shù)展開搭配移動(dòng)最小二乘法，將每個(gè)點(diǎn)位上的微分量以鄰近點(diǎn)的物理量線性累加表示。Benito等[14]提出GFDM法的離散公式，針對(duì)其中的一些參數(shù)與特性進(jìn)行一系列的研究，之后將這一方法改良后與其它無(wú)網(wǎng)格方法作比較[15]，并使用該方法求解不同類型的偏微分方程，如對(duì)流擴(kuò)散方程[16]、拋物線型方程[17]。同時(shí)，F(xiàn)an等[18]利用此方法求解Burgers方程式，Li等[19]用于求解淺水波方程，均得到了良好的結(jié)果。Zhang等[20]采用此方法運(yùn)用到工程實(shí)際問題中，例如數(shù)值波浪水槽中非線性波浪的傳播、非線性自由液面的液體晃動(dòng)問題[21]、緩坡方程的求解[22]；Gu等[23]將GFDM方法成功應(yīng)用于三維熱傳導(dǎo)的逆時(shí)變?cè)磫栴}的精確求解。從這些數(shù)值案例模擬的比對(duì)結(jié)果可以看出，GFDM法在求解高階偏微分方程問題上具有很大的潛力。

本文針對(duì)空間含有四階偏導(dǎo)項(xiàng)和時(shí)間含有二階偏導(dǎo)項(xiàng)的兩端支撐輸流直管橫向運(yùn)動(dòng)微分方程，采用GFDM法和Houblot方法分別對(duì)微分方程的空間項(xiàng)和時(shí)間項(xiàng)進(jìn)行離散，建立一種新的高階精度數(shù)值模式，通過(guò)與前人的數(shù)值結(jié)果對(duì)比，驗(yàn)證本研究所提出的數(shù)值模型的準(zhǔn)確性和可行性，在此基礎(chǔ)上，分析不同支撐條件對(duì)輸流直管模型橫向振動(dòng)響應(yīng)特性的影響。

1 控制方程及邊界條件

如圖1所示，考慮兩端支撐輸流直管，管道長(zhǎng)度為L(zhǎng)，管道軸線為x軸，管道橫向?yàn)閥軸。忽略重力、內(nèi)部阻尼和流體壓力的影響，則其橫向運(yùn)動(dòng)微分方程可表述為

(1)

式中：Y為管道的橫向位移；E為管道的彈性模量；I為管道的橫截面慣性矩；mf和mp分別為單位長(zhǎng)度管內(nèi)流體的質(zhì)量和管道的質(zhì)量；U為流體流速。

引入無(wú)量綱變量

代入式(1)，整理可得兩端支撐輸流直管的無(wú)量綱橫向運(yùn)動(dòng)微分方程

(2)

針對(duì)不同的管道支撐條件，其無(wú)量綱化的邊界條件可表述為:

(a)兩端固支(見圖1(a))

(3)

(b)一端固支一端簡(jiǎn)支(如圖1(b))

(4)

(c)兩端簡(jiǎn)支(見圖1(c))

(5)

圖1 兩端支撐輸流直管模型Fig.1 Schematic diagram of fluid-conveying pipe

2 數(shù)值方法

2.1 廣義有限差分法

輸流直管橫向運(yùn)動(dòng)微分方程式(2)，在空間坐標(biāo)上最高具有四階偏導(dǎo)項(xiàng)，本文采用廣義有限差分法進(jìn)行離散，其方法是基于移動(dòng)最小二乘法與泰勒級(jí)數(shù)四階展開。首先在整個(gè)計(jì)算區(qū)域內(nèi)布N個(gè)點(diǎn)，再將每個(gè)點(diǎn)位上的偏微分項(xiàng)轉(zhuǎn)換成由子區(qū)域內(nèi)各點(diǎn)物理量與權(quán)重系數(shù)乘積的線性累加。對(duì)于區(qū)域內(nèi)的第i點(diǎn)而言，選擇ns個(gè)最鄰近點(diǎn)，形成一個(gè)子區(qū)域，如圖2所示。

圖2 子區(qū)域中選擇臨近點(diǎn)示意圖Fig.2 Schematic diagram of nodes in local region

當(dāng)?shù)趇點(diǎn)的子區(qū)域形成后，在該子區(qū)域內(nèi)以第i點(diǎn)為中心進(jìn)行泰勒級(jí)數(shù)展開，因式(2)對(duì)空間項(xiàng)的微分最高階數(shù)為四階，從而略去四階以上各項(xiàng)，并定義一個(gè)函數(shù)B(η)

(6)

式中：j為子區(qū)域內(nèi)的節(jié)點(diǎn)編號(hào)；δij=ξi-ξi,j為沿著布點(diǎn)方向上第i點(diǎn)與第j點(diǎn)的距離；ξi和ξi,j分別為第i點(diǎn)和第j點(diǎn)的坐標(biāo)值；ηi,j為第i個(gè)子區(qū)域中的第j個(gè)點(diǎn)的物理量；w(δij)為權(quán)重函數(shù)

(7)

式中：dmi為第i點(diǎn)與子區(qū)域內(nèi)最遠(yuǎn)點(diǎn)的距離。

根據(jù)移動(dòng)最小二乘法的思想，將函數(shù)B(η)分別對(duì)?η/?ξ，?2η/?ξ2，?3η/?ξ3和?4η/?ξ4求極小值，可得

A·Dη=b

(8)

其中，

從式(8)可知，系數(shù)矩陣A是一個(gè)對(duì)稱矩陣，是由第i點(diǎn)與其子區(qū)域內(nèi)ns個(gè)點(diǎn)的物理量計(jì)算得到，而矩陣b是由子區(qū)域內(nèi)各節(jié)點(diǎn)物理量和空間坐標(biāo)構(gòu)成，因此可將矩陣b分解為

b=BQ

(9)

式中：Q=[ηi,ηi,1,ηi,2,ηi,3,…,ηi,ns]T為子區(qū)域內(nèi)第i點(diǎn)和與其相鄰ns點(diǎn)的物理量。從而，每一點(diǎn)位上的前四階偏微分項(xiàng)Dη可表示為

(10)

因輸流直管橫向運(yùn)動(dòng)微分方程式(2)中只含有對(duì)空間物理量的一階、二階和四階偏導(dǎo)數(shù)，故提取式(10)中每一個(gè)點(diǎn)位i上未知物理量(位移)的一階、二階和四階偏微分量的表達(dá)式，即

(11)

(12)

(13)

2.2 Houbolt法

由于兩端支撐輸流直管橫向運(yùn)動(dòng)微分方程式(2)，在時(shí)間坐標(biāo)上最高具有二階偏導(dǎo)項(xiàng)，本文采用Houbolt法對(duì)時(shí)間項(xiàng)進(jìn)行離散。該法屬于四點(diǎn)格式的隱式時(shí)間積分法，具有二階精度且無(wú)條件穩(wěn)定，即通過(guò)對(duì)n-2，n-1，n和n+1四個(gè)時(shí)刻的位移η進(jìn)行三次插值來(lái)近似表示其一階時(shí)間導(dǎo)數(shù)?η/?τ和二階時(shí)間導(dǎo)數(shù)?2η/?τ2，其表達(dá)式為[24]

(14)

(15)

因Houbolt法在求解未知時(shí)間層物理量ηn+1時(shí)，需已知前兩個(gè)時(shí)間層的物理量ηn-1和ηn-2，從而需要起步條件，本文采用Euler法進(jìn)行起步，即

(16a)

(16b)

2.3 方程離散

首先，使所有內(nèi)部點(diǎn)滿足控制方程式，采用GFDM對(duì)控制方程式中的空間變量偏微分進(jìn)行離散，即將式(11)～式(13)代入式(2)中，可得

(17)

其次，使所有邊界點(diǎn)滿足對(duì)應(yīng)的邊界條件，采用GFDM法對(duì)邊界條件進(jìn)行離散可得：

(1)兩端固支

(18)

(2)一端固支一端簡(jiǎn)支

(19)

(3)兩端簡(jiǎn)支

(20)

式(17)結(jié)合邊界條件式(18)～式(20)中的一種，可定義一種支撐條件下輸流直管的動(dòng)力學(xué)方程組，即：

(21)

式中：M，C，K分別為離散系統(tǒng)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣；η為待求物理量未知矩陣。

對(duì)式(21)可進(jìn)行模態(tài)分析，設(shè)方程式的特解為

η=φeλτ

(22)

式中：φ為N階位移幅值列陣。

將式(22)代入式(21)中，可得兩端支撐輸流管道的廣義特征值問題

(λ2M+λC+K)φ=0

(23)

其特征方程可表述為

|λ2M+λC+K|=0

(24)

上式是關(guān)于λ的2N次實(shí)系數(shù)代數(shù)方程，設(shè)無(wú)重根，通過(guò)求解可得2N個(gè)共軛對(duì)形式的互異特征值λ，其值通常為復(fù)數(shù)，虛部即表示輸流管道振動(dòng)頻率。

同時(shí)，采用Houbolt法對(duì)式(21)中的時(shí)間項(xiàng)一階和二階微分進(jìn)行離散，可得每一內(nèi)部點(diǎn)位i上的兩端支撐輸流直管橫向振動(dòng)微分方程的離散形式，即

(i=3,4,…,N-2)

(25)

式(25)結(jié)合邊界條件式(18)～式(20)中的一種，可定義一種支撐條件下輸流直管的線性代數(shù)方程組為

[C]N×N·{η}N×1={f}N×1

(26)

式中：C為稀疏矩陣；f為控制方程式與邊界條件離散后的非齊次項(xiàng)；η為待求未知矩陣。結(jié)合不同案例的初始條件，通過(guò)式(26)即可求得不同支撐條件下輸流直管在不同時(shí)刻，每一點(diǎn)位上的物理量。

3 數(shù)值模型驗(yàn)證

為了驗(yàn)證本文所提出的數(shù)值模式的準(zhǔn)確性和魯棒性，本節(jié)對(duì)兩端支撐軸向運(yùn)動(dòng)梁模型和兩端固支輸流直管模型進(jìn)行數(shù)值模擬，并與前人所做研究結(jié)果進(jìn)行對(duì)比，其中當(dāng)流體流速u用管道橫向運(yùn)動(dòng)速度代替，式(2)即轉(zhuǎn)化為兩端支撐的梁模型。

3.1 兩端支撐梁模型

當(dāng)管道運(yùn)動(dòng)速度u=0時(shí)，兩端支撐梁模型退化為簡(jiǎn)單的直梁模型，其前N階固有頻率具有精確解，采用總點(diǎn)數(shù)N=604，時(shí)間步長(zhǎng)Δt=0.001，表1給出了GFDM法計(jì)算得到的不同支撐條件下直梁模型前四階固有頻率，與DTM、DQM以及Thomson的精確解[25]研究成果進(jìn)行對(duì)比，吻合良好，說(shuō)明該數(shù)值模式具有相當(dāng)高的精度。為了進(jìn)一步說(shuō)明其計(jì)算效率，本文也將三種不同數(shù)值方法的計(jì)算時(shí)間列于表1，雖然GFDM的總布點(diǎn)數(shù)N達(dá)604，但所用計(jì)算時(shí)間均在0.85 s左右，低于DTM(N=17)，部分高于DQM(N=60)，這是由于應(yīng)用GFDM法對(duì)控制方程進(jìn)行離散后，所得到的代數(shù)方程組的系數(shù)矩陣為稀疏矩陣，因此可提高計(jì)算效率。

當(dāng)管道運(yùn)動(dòng)速度u=0.5時(shí)，給定初始條件為

(27)

以管道兩端固支為例，計(jì)算得到管道中點(diǎn)處振幅η的時(shí)間歷程，如圖3所示?？梢?，應(yīng)用GFDM法得到的結(jié)果與An等研究中采用GITT法得到的結(jié)果吻合良好。圖3中分別采用了不同總點(diǎn)數(shù)N(見圖3(a))、不同時(shí)間步長(zhǎng)Δt(見圖3(b))和不同選點(diǎn)數(shù)ns(見圖3(c))進(jìn)行數(shù)值模擬，結(jié)果表明，隨著總點(diǎn)數(shù)N和選點(diǎn)數(shù)ns的增加或時(shí)間步長(zhǎng)Δt的減小，GFDM的計(jì)算結(jié)果與An等研究中的數(shù)值結(jié)果越接近，說(shuō)明本文提出的數(shù)值模式具有良好的穩(wěn)定性。

表1 兩端支撐梁固有頻率(u=0)

圖3 兩端固支梁受迫振動(dòng)中點(diǎn)處振幅比較(u=0.5)Fig.3 Comparison of amplitudes at midpoint of forced vibration of clamped beam at u=0.5

3.2 兩端固支輸流直管模型

同樣以兩端固支輸流直管模型為例，在流體流速u的作用下，輸流直管將出現(xiàn)受迫振動(dòng)，本節(jié)數(shù)值模型參數(shù)分別取N=604，ns=20，Δt=0.001 0，同時(shí)給定初始條件為

(28)

圖4為不同流體流速u和不同質(zhì)量比β輸流直管中點(diǎn)處位移的時(shí)間歷程。將數(shù)值結(jié)果與Gu等的研究成果進(jìn)行對(duì)比，結(jié)果也是非常吻合,進(jìn)一步說(shuō)明了本文所提出的數(shù)值模型具有較高的精確度。從圖4可知，在相同流速u情況下，隨著質(zhì)量比β的增加，兩端固支輸流直管中點(diǎn)處的振動(dòng)速率加快，而振幅無(wú)明顯變化，見圖4(a)和圖4(b)所示。在相同質(zhì)量比β情況下，隨著流體流速u的減小，兩端固支輸流直管中點(diǎn)處振動(dòng)幅值減小，但振動(dòng)速率加快，見圖4(a)和圖4(c)所示。

圖4 不同流體流速u和質(zhì)量比β兩端固支輸流直管中點(diǎn)振幅比較Fig.4 Comparison of amplitudes at midpoint of differential fluid velocity and mass ratio of clamped-clamped pipe conveying fluid

4 不同支撐條件下輸流直管振動(dòng)響應(yīng)

當(dāng)管道在輸流過(guò)程中，由于端部約束限制條件不同，其振動(dòng)響應(yīng)特性也將不同。為了進(jìn)一步研究支撐條件對(duì)輸流直管橫向振動(dòng)響應(yīng)特性的影響，本節(jié)針對(duì)三種不同支撐條件(兩端固支、兩端簡(jiǎn)支和一端固支一端簡(jiǎn)支)下輸流直管進(jìn)行模擬。數(shù)值模型參數(shù)仍取總點(diǎn)數(shù)N=604，時(shí)間步長(zhǎng)Δt=0.001，子區(qū)域選點(diǎn)數(shù)ns=20，流體流速u=1.5，質(zhì)量比β=0.5。

圖5(a)給出了三種不同支撐條件的輸流直管中點(diǎn)振動(dòng)幅值時(shí)程，從圖5(a)可知，輸流直管均作周期性有規(guī)律的振動(dòng)，當(dāng)兩端簡(jiǎn)支時(shí)輸流直管中點(diǎn)處的振幅最大，一端固支一端簡(jiǎn)支時(shí)次之，而兩端固支時(shí)輸流直管中點(diǎn)處的振幅最小，為兩端簡(jiǎn)支時(shí)的一半；為了進(jìn)一步了解三種不同支撐情況下管道中點(diǎn)處的頻域響應(yīng)，通過(guò)傅里葉變換得到了三種不同支撐條件下輸流直管中點(diǎn)處的頻譜圖，如圖5(b)所示。不同支撐條件下輸流直管均只出現(xiàn)一階主頻，兩端簡(jiǎn)支(f=1.34)時(shí)振動(dòng)頻率最小，兩端固支(f=3.42)時(shí)振動(dòng)頻率最大，而一端固支一端簡(jiǎn)支(f=2.32)時(shí)介于二者之間。

圖5 不同支撐輸流直管中點(diǎn)處振幅比較Fig.5 Comparison of amplitudes at midpoint of pipe conveying fluid for different boundary conditions

圖6為輸流直管在τ=13時(shí)刻全管振動(dòng)幅值曲線?？梢?，在忽略重力、內(nèi)部阻尼和流體壓力影響的條件下，兩端簡(jiǎn)支和兩端固支其振動(dòng)響應(yīng)最大值出現(xiàn)在輸流直管中點(diǎn)處，對(duì)于一端固支一端簡(jiǎn)支支撐條件，因兩端支撐條件不對(duì)稱且簡(jiǎn)支端約束限制較低，致使管道振幅最大值出現(xiàn)位置向右偏移。

圖6 τ=13時(shí)不同支撐輸流直管振幅比較Fig.6 Comparison of amplitudes of pipe conveying fluid for different boundary conditions at τ=13

5 結(jié) 論

本文針對(duì)空間上最高具有四階偏導(dǎo)項(xiàng)和時(shí)間上最高具有二階偏導(dǎo)項(xiàng)的兩端支撐輸流管道橫向運(yùn)動(dòng)微分方程，采用區(qū)域型無(wú)網(wǎng)格法分析了兩端支撐軸向運(yùn)動(dòng)梁模型和兩端固支輸流直管模型的橫向振動(dòng)問題，對(duì)于該問題在空間和時(shí)間上分別采用GFDM法和Houblot法進(jìn)行離散，建立高階精度的無(wú)網(wǎng)格法數(shù)值模式。

(1)通過(guò)與前人研究成果進(jìn)行比較以及方法本身影響因素(總點(diǎn)數(shù)N、時(shí)間步長(zhǎng)△t和子區(qū)域選點(diǎn)數(shù)ns)測(cè)試對(duì)比，結(jié)果吻合良好，表明所提出的數(shù)值模型在求解輸流直管振動(dòng)響應(yīng)問題上具有良好的準(zhǔn)確性和魯棒性。

(2)對(duì)比分析了三種不同支撐(兩端固支、兩端簡(jiǎn)支和一端固支一端簡(jiǎn)支)下輸流直管振動(dòng)響應(yīng)特性，結(jié)果表明，輸流直管均以一階主頻作周期性有規(guī)律的振動(dòng)，在對(duì)稱支撐(兩端簡(jiǎn)支和兩端固支)情況下，振幅與振動(dòng)頻率成反比，即振幅越大，頻率越小；在非對(duì)稱支撐(一端固支一端簡(jiǎn)支)時(shí)，振動(dòng)幅值和頻率介于兩種對(duì)稱支撐之間，且振幅最大值出現(xiàn)位置向右偏移。

本文僅針對(duì)在忽略重力、內(nèi)部阻尼和流體壓力影響的條件下兩端支撐的水平輸流直管進(jìn)行研究，從本文研究測(cè)試的結(jié)果可以看出，由于GFDM法可將每個(gè)點(diǎn)位上高階微分項(xiàng)進(jìn)行快速轉(zhuǎn)換為線性代數(shù)方程組，因此具有較大的潛力，下一步可將其運(yùn)用到輸流管道非線性振動(dòng)和彎曲輸流管道的研究中。