基于廣義有限差分法的輸流直管振動響應(yīng)特性研究

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

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

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

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

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

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

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

1 控制方程及邊界條件

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

(1)

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

引入無量綱變量

代入式(1)，整理可得兩端支撐輸流直管的無量綱橫向運動微分方程

(2)

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

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

(3)

(b)一端固支一端簡支(如圖1(b))

(4)

(c)兩端簡支(見圖1(c))

(5)

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

2 數(shù)值方法

2.1 廣義有限差分法

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

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

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

(6)

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

(7)

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

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

A·Dη=b

(8)

其中，

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

b=BQ

(9)

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

(10)

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

(11)

(12)

(13)

2.2 Houbolt法

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

(14)

(15)

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

(16a)

(16b)

2.3 方程離散

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

(17)

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

(1)兩端固支

(18)

(2)一端固支一端簡支

(19)

(3)兩端簡支

(20)

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

(21)

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

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

η=φeλτ

(22)

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

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

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

(23)

其特征方程可表述為

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

(24)

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

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

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

(25)

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

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

(26)

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

3 數(shù)值模型驗證

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

3.1 兩端支撐梁模型

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

當管道運動速度u=0.5時，給定初始條件為

(27)

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

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

圖3 兩端固支梁受迫振動中點處振幅比較(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)受迫振動，本節(jié)數(shù)值模型參數(shù)分別取N=604，ns=20，Δt=0.001 0，同時給定初始條件為

(28)

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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