張 挺， 譚志新， 張 恒， 范佳銘， 楊志強(qiáng)(. 福州大學(xué) 土木工程學(xué)院， 福州 3506； . 臺(tái)灣海洋大學(xué) 河海工程系， 臺(tái)灣 基隆 04)

管道系統(tǒng)廣泛運(yùn)用于生活中的各個(gè)領(lǐng)域，管內(nèi)流體與管壁結(jié)構(gòu)間的耦合振動(dòng)可能會(huì)產(chǎn)生噪音污染，導(dǎo)致結(jié)構(gòu)失效，甚至爆管等嚴(yán)重的工程問(wèn)題，因此管道振動(dòng)長(zhǎng)久以來(lái)都是一個(gè)研究的熱點(diǎn)。早期學(xué)者們對(duì)輸流管道軸向振動(dòng)響應(yīng)的研究計(jì)算以“經(jīng)典水錘理論”為主，該模型假定管道是固定不動(dòng)的，沒(méi)有考慮管道與流體間的耦合作用，對(duì)管內(nèi)壓力有一定的預(yù)測(cè)作用，因簡(jiǎn)化為較簡(jiǎn)單問(wèn)題，被廣泛應(yīng)用于工程實(shí)際中。然而在實(shí)際問(wèn)題中，水錘作用下管內(nèi)流體的壓力脈動(dòng)會(huì)造成管壁的收縮與膨脹，從而引起管道振動(dòng)，管道振動(dòng)又反作用影響流體壓力，這就是流固耦合(Fluid-Structure Interaction，F(xiàn)SI)現(xiàn)象。隨后，不少學(xué)者從流體與管道結(jié)構(gòu)的動(dòng)量和連續(xù)性的角度出發(fā)，在流體和管道間的流固耦合作用方面進(jìn)行深入研究。Otwell[1]較早提出了輸流直管軸向振動(dòng)流固耦合-四方程模型。Wiggert等[2]接著深入研究流體壓力波的基礎(chǔ)，對(duì)Otwell提出的四方程模型進(jìn)行了改進(jìn)，并對(duì)輸流管道的流固耦合現(xiàn)象作了較為全面的概述與討論。

流固耦合作用是引起輸流管道振動(dòng)的主要原因，管道系統(tǒng)的流固耦合振動(dòng)問(wèn)題被稱為“典型的動(dòng)力學(xué)問(wèn)題”[3]。目前，對(duì)管道系統(tǒng)動(dòng)力學(xué)特性的分析主要有時(shí)域[4-5]和頻域[6-7]兩個(gè)方面，其研究分析的關(guān)鍵技術(shù)就是尋求高精度、高效可靠的數(shù)值計(jì)算方法。其中，主要有特征線法(Method of Characteristics，MOC)[8-10]、有限元法(Finite Element Method，F(xiàn)EM)[11-12]，以及特征線-有限元法(MOC-FEM)[13-14]等。對(duì)于管道系統(tǒng)中的偏微分方程(波方程)，采用特征線法將方程組轉(zhuǎn)化為一組特殊的常微分方程，可直接用于求解計(jì)算，但計(jì)算過(guò)程中存在多特征線插值問(wèn)題；對(duì)于空間復(fù)雜管系的結(jié)構(gòu)模態(tài)分析與響應(yīng)計(jì)算，有限元法具有明顯的優(yōu)越性；特征線-有限元法將兩種方法結(jié)合起來(lái)，可發(fā)揮了兩種方法的優(yōu)勢(shì)。其不足之處在于，在每一步計(jì)算時(shí)刻都需要進(jìn)行液體和管道間的數(shù)據(jù)傳輸，以滿足相容和耦合條件，計(jì)算量相對(duì)較大。

為了提供輸流管道振動(dòng)問(wèn)題簡(jiǎn)單且準(zhǔn)確的數(shù)值仿真模式，本文將分離系數(shù)矩陣差分法(Split-Coefficient Matrix Finite Difference Method，SCM-FDM)[15-17]應(yīng)用于輸流管道流固耦合-四方程模型，并配合隱式歐拉法(Implicit Euler Method，IEM)建立一種簡(jiǎn)單易行的數(shù)值計(jì)算模式，通過(guò)與前人的數(shù)值結(jié)果和經(jīng)典水錘理論進(jìn)行對(duì)比，驗(yàn)證本研究所提出數(shù)值模型的準(zhǔn)確性、穩(wěn)定性，在此基礎(chǔ)上研究泊松耦合和連接耦合對(duì)輸流直管軸向振動(dòng)響應(yīng)特性的影響。

1 數(shù)值模型及初邊值條件

1.1 數(shù)值模型

本文研究的水箱-管道-閥門系統(tǒng)(Reservoir- Pipeline-Valve，RPV)，如圖1所示。管長(zhǎng)L，管道軸向(z)從上游水箱端向下游閥門端為正，假定管道為水平直管，管壁為均勻、各向同性的彈性材料，不考慮橫向慣性力、彎曲變形及摩擦的影響，研究閥門瞬時(shí)關(guān)閉時(shí)，管道在水錘作用下的軸向振動(dòng)響應(yīng)特性。因此若管道較長(zhǎng)，可假想管道按一定間距設(shè)置理想滑動(dòng)支撐，不影響管道軸向振動(dòng)特性。

基于以上假定，輸流直管軸向振動(dòng)流固耦合-四方程模型可表達(dá)為

(1)

(2)

(3)

(4)

圖1 水箱-管道-閥門系統(tǒng)

1.2 邊界條件

對(duì)于水箱-管道-閥門系統(tǒng)，邊界條件需針對(duì)流體和管道在水箱端和閥門端分別給出，并配合控制方程構(gòu)成兩個(gè)端點(diǎn)關(guān)于未知量的方程組。管道上游端與水箱底部固定相連，假定水箱水位恒定，則上游水箱端邊界條件可表示為

(5)

式中：P0為一給定值。

管道下游端與閥門連接，閥門瞬時(shí)關(guān)閉時(shí)，管道與流體間的相對(duì)速度為零。假定閥門為輕質(zhì)閥門，當(dāng)其邊界約束條件不同時(shí)，管道表現(xiàn)出的耦合振動(dòng)特性也將不同。

若閥門為固支，管道在閥門處受到強(qiáng)制理想約束，軸向管道運(yùn)動(dòng)速度為零，管道與流體間的耦合振動(dòng)以泊松耦合為主，此時(shí)下游閥門端邊界條件可表達(dá)為

(6)

當(dāng)閥門為簡(jiǎn)支時(shí)，管道在閥門處能自由滑動(dòng)，軸向不受外力，管道與流體間的耦合振動(dòng)，除受泊松耦合影響外，還受連接耦合影響，此時(shí)下游閥門端邊界條件可表達(dá)為

(7)

式中：Af和As分別為管道流通截面積和管道環(huán)形截面積。

1.3 初始條件

閥門瞬時(shí)關(guān)閉前，管內(nèi)流體流動(dòng)為穩(wěn)態(tài)，流速為常數(shù)，軸向管速及管道應(yīng)力都為零，即：

v=v0

(8)

P=P0

(9)

(10)

σz=0

(11)

2 數(shù)值算法

將控制方程組(1)～(4)表示成如下矩陣矢量形式

(12)

其中，cf和cs分別為流體壓力波速和管道應(yīng)力波速

(13)

(14)

對(duì)于式(12)來(lái)說(shuō)，系數(shù)矩陣A可逆，且4階方陣A-1B存在4個(gè)互異的特征值，所以存在變換矩陣T，可使Φ(z,t)=Τψ(z,t)。于是，可作如下轉(zhuǎn)化

(15)

將上式左右兩邊同時(shí)左乘Τ-1A-1

(16)

令Τ-1A-1BΤ=Λ，則對(duì)角矩陣Λ的對(duì)角元素即為矩陣A-1B的特征值，該特征值即為|B-λA|=0對(duì)應(yīng)的一系列特征根λj(j=1、2、3、4)，則T為λj相應(yīng)的特征矢量矩陣。

(17)

定義，λ+1=max(0,λ1)，λ-1=min(0,λ1)，同理定義λ+2、λ-2，λ+3、λ-3，λ+4、λ-4。將Λ作如下拆分

Λ++Λ-

(18)

式中:Λ+表示往正方向(向上游)傳遞的波，而Λ-表示往負(fù)方向(向下游)傳遞的波。并將ψ(z,t)=Τ-1Φ(z,t)，代回式(16)，可得：

(19)

經(jīng)由特征值與特征矢量的推導(dǎo)，式(12)可以轉(zhuǎn)換為式(19)，式(19)中成功地將波傳遞的信息包含在系數(shù)矩陣中，因此可以根據(jù)式(19)中不同的波傳遞方向，而采用不同的差分公式，以發(fā)展準(zhǔn)確且可行的數(shù)值仿真模式。采用有限差分法對(duì)上式進(jìn)行離散，并采用隱式歐拉法進(jìn)行時(shí)間項(xiàng)的離散：為保證特征值λj(j=1、2、3、4)特征方向與差分方向一致，?Φ(z,t)/?z與正特征值矩陣Λ+匹配時(shí)，采用后項(xiàng)差分；與負(fù)特征值矩陣Λ-匹配時(shí)，采用前項(xiàng)差分，整理可得：

(20)

式中:上標(biāo)n為計(jì)算時(shí)間層，初始時(shí)刻n=1；下標(biāo)i為管道從上游到下游均布的計(jì)算節(jié)點(diǎn)，總節(jié)點(diǎn)數(shù)為N；Δt、Δz分別為時(shí)間步長(zhǎng)和空間步長(zhǎng)，其中Δz=L/(N-1)。

(21)

(22)

式(20)～式(22)左邊為n+1時(shí)刻待求未知量的關(guān)系式，右邊為n時(shí)刻已知量的關(guān)系式?？梢?jiàn)，本研究所提分離系數(shù)矩陣差分法避開了特征線法復(fù)雜的時(shí)間或空間插值，只根據(jù)波的傳播方向選擇差分公式進(jìn)行計(jì)算，配合上隱式歐拉法，使得計(jì)算簡(jiǎn)單易行。

3 模型驗(yàn)證

3.1 案例1

本研究第一個(gè)驗(yàn)證案例為Wilkinson 和 Curtis所做的薄壁直鋼管水錘實(shí)驗(yàn)，管道末端為簡(jiǎn)支，應(yīng)用上述數(shù)值模型進(jìn)行模擬，本案例模型參數(shù)如表1。初始流速v0=5 m/s，其余初始條件都為0。

表1 Wilkinson 和Curtis實(shí)驗(yàn)參數(shù)

圖2為計(jì)算得到的z=6.10 m和z=1.20 m兩處管內(nèi)壓強(qiáng)響應(yīng)曲線，并將計(jì)算結(jié)果與Tijsseling[18]采用特征線法得到的數(shù)值結(jié)果進(jìn)行對(duì)比，吻合良好，準(zhǔn)確地捕捉到管道內(nèi)水錘壓力的非線性變化過(guò)程。為了進(jìn)一步驗(yàn)證本文數(shù)值模式的穩(wěn)定性和有效性，圖中分別給出了三組不同總點(diǎn)數(shù)N(圖2(a))和三組不同時(shí)間步長(zhǎng)Δt(圖2(b))的數(shù)值結(jié)果對(duì)比，可見(jiàn)隨著總布點(diǎn)數(shù)的增加或時(shí)間步長(zhǎng)的縮短，分離系數(shù)矩陣差分法與特征線法得到的數(shù)值結(jié)果越一致，表明本文所提出的數(shù)值計(jì)算模式是可行且正確的。

(a) 不同空間步長(zhǎng)(Δt=10-7 s)

(b) 不同時(shí)間步長(zhǎng)(N=40 000)

Fig.2 Response curves of pressure by using different numbers of total nodes (Δt=10-7s) and different time increments (N=40 000)

3.2 案例2

Tijsseling[18]針對(duì)Delft水力學(xué)基準(zhǔn)問(wèn)題A采用特征線法模擬水箱-管道-閥門系統(tǒng)的水錘壓力特性，本研究將以所提出的分離系數(shù)矩陣差分法模擬同一案例，并將計(jì)算結(jié)果與Tijsseling數(shù)值模擬結(jié)果進(jìn)行對(duì)比。模型參數(shù)如表2，初始流速v0=1 m/s，其余初始條件都為0，閥門考慮固支與簡(jiǎn)支兩種工況。計(jì)算時(shí)間步長(zhǎng)Δt=10-7s，總節(jié)點(diǎn)數(shù)N=40 000。

表2 Delft水力學(xué)基準(zhǔn)問(wèn)題A模型參數(shù)

圖3分別給出了管道閥門處(z=20 m)和中點(diǎn)處(z=10 m)管內(nèi)壓強(qiáng)的響應(yīng)曲線。將計(jì)算結(jié)果與特征線法得到的數(shù)值結(jié)果進(jìn)行對(duì)比，吻合良好，表明本文的數(shù)值模式具有良好的準(zhǔn)確性和穩(wěn)定性。同時(shí)，將模擬結(jié)果與經(jīng)典水錘理論計(jì)算結(jié)果對(duì)比發(fā)現(xiàn)，耦合作用不僅影響壓力幅值，還影響振動(dòng)頻率，其對(duì)管道振動(dòng)的影響不容忽視。

當(dāng)閥門固支(即只考慮泊松耦合)時(shí)，隨著閥門的瞬時(shí)關(guān)閉，液體流動(dòng)受阻，產(chǎn)生壓力波，其以cf在管道內(nèi)傳播，壓力脈動(dòng)所到之處又引起管道的徑向脹縮，從而在管壁內(nèi)產(chǎn)生了軸向應(yīng)力波，其以cs在管壁內(nèi)傳播，由于壓力波與應(yīng)力波不同速，且相互耦合，導(dǎo)致管內(nèi)脈動(dòng)壓力幅值較不考慮耦合作用時(shí)增大，且脈動(dòng)曲線出現(xiàn)較多局部突變，高頻振動(dòng)成分非常明顯，表現(xiàn)出壓力波主導(dǎo)，管壁應(yīng)力波疊加的強(qiáng)烈振動(dòng)。此外，考慮泊松耦合的計(jì)算結(jié)果相位與經(jīng)典水錘理論的計(jì)算結(jié)果相位幾乎同步。

當(dāng)閥門簡(jiǎn)支(即同時(shí)考慮泊松耦合和連接耦合)時(shí)，管內(nèi)脈動(dòng)壓力幅值進(jìn)一步增大，振動(dòng)曲線更加不規(guī)則，且出現(xiàn)明顯的相位延遲。與泊松耦合相比，邊界約束條件的減弱，振動(dòng)周期變長(zhǎng)，高頻振動(dòng)成分不再明顯。

(a) z=20 m

(b) z=10 m

4 兩種不同耦合作用對(duì)響應(yīng)的影響

4.1 流體速度響應(yīng)

圖4分別給出了管道閥門處(z=20 m)和管道中點(diǎn)處(z=10 m)在考慮泊松耦合和同時(shí)考慮兩種耦合時(shí)的流體流速響應(yīng)曲線?？梢?jiàn)，當(dāng)閥門理想固支，此時(shí)只考慮泊松耦合，閥門瞬時(shí)關(guān)閉將導(dǎo)致閥門處(z=20 m)水流流速由初始流速瞬間變?yōu)榱悖也辉俨▌?dòng)；而當(dāng)閥門簡(jiǎn)支時(shí)，即同時(shí)考慮泊松耦合和連接耦合，此時(shí)閥門處水流流速伴隨管道的自由伸縮作周期性振動(dòng)，見(jiàn)圖4(a)。在管道中點(diǎn)處(z=10 m)，泊松耦合和經(jīng)典水錘理論的水流流速響應(yīng)曲線相位幾乎同步，且振動(dòng)規(guī)律極為相似。而同時(shí)考慮泊松耦合和連接耦合的流速響應(yīng)曲線發(fā)生相位延遲，且響應(yīng)曲線較為不規(guī)則，見(jiàn)圖4(b)。

(a) z=20 m

(b) z=10 m

4.2 管道軸向振動(dòng)速度響應(yīng)

(a) z=20 m

(b) z=10 m

(a) 管道軸向振動(dòng)速度

(b) 管道平均應(yīng)力

4.3 管壁應(yīng)力響應(yīng)

(a) z=20 m

(b) z=10 m

5 結(jié) 論

本研究將分離系數(shù)矩陣差分法配合隱式歐拉法應(yīng)用于輸流管道軸向振動(dòng)流固耦合問(wèn)題，即提出一種準(zhǔn)確可行的數(shù)值仿真模式，以應(yīng)用于四方程模型的數(shù)值計(jì)算中，研究水錘激勵(lì)下輸流直管耦合軸向振動(dòng)響應(yīng)特性。該數(shù)值模式有效結(jié)合特征線法的概念與有限差分法，避開了特征線法復(fù)雜的插值計(jì)算，只根據(jù)波的傳播方向進(jìn)行差分，計(jì)算簡(jiǎn)單易行。計(jì)算結(jié)果與前人的數(shù)值結(jié)果對(duì)比，吻合良好，表明其具有較高的適應(yīng)性、穩(wěn)定性和準(zhǔn)確性。

對(duì)水箱-管道-閥門系統(tǒng)而言，不同約束條件管道主要表現(xiàn)出的流固耦合強(qiáng)度不同，閥門固支時(shí)以泊松耦合為主，閥門簡(jiǎn)支時(shí)以泊松耦合和連接耦合的疊加效應(yīng)為主，邊界條件對(duì)輸流管道流固耦合作用影響顯著。對(duì)比經(jīng)典水錘方程計(jì)算結(jié)果發(fā)現(xiàn)，兩種耦合作用對(duì)管道軸向振動(dòng)特性的影響不容忽視，泊松耦合主要影響振動(dòng)響應(yīng)幅值，而連接耦合不僅會(huì)影響振動(dòng)幅值，同時(shí)也會(huì)影響振動(dòng)頻率。

[1] OTWELL R S. The effect of elbow restraint on pressure transients[D]. East Lansing, MI: Michigan State University, 1984.

[2] WIGGERT D C, TIJSSELING A S. Fluid transients and fluid-structure interaction in flexible liquid-filled piping[J]. Applied Mechanics Reviews, 2001, 54(5): 455-481.

[3] 任建亭, 姜節(jié)勝. 輸流管道系統(tǒng)振動(dòng)研究進(jìn)展[J]. 力學(xué)進(jìn)展, 2003, 33(3): 313-324.

REN Jianting, JIANG Jiesheng. Advances and trends on vibration of pipes conveying fluid[J]. Advances in Mechanics, 2003, 33(3): 313-324.

[4] 范士娟, 楊超. 輸水管道流固耦合振動(dòng)數(shù)值計(jì)算[J]. 噪聲與振動(dòng)控制, 2010, 30(6): 43-46.

FAN Shijuan，YANG Chao. Numerical calculation of fluid-structure coupling vibration of a water pipeline[J]. Noise and Vibration Control, 2010, 30(6): 43-46.

[5] 席志德, 馬建中, 孫磊. 考慮流-固耦合效應(yīng)的空間管道水錘方法研究[J]. 核動(dòng)力工程, 2013, 34(2): 1-4.

XI Zhide, MA Jianzhong, SUN Lei. Method to study water hammer with fluid-structure interaction in spatial pipe[J]. Nuclear Power Engineering, 2013, 34(2): 1-4.

[6] 姬賀炯,白長(zhǎng)青,韓省亮．輸流管耦合動(dòng)力學(xué)特性分析[J]．噪聲與振動(dòng)控制,2013,33(5): 10-14．

JI Hejiong, BAI Changqing, HAN Shengliang． Analysis of dynamic characteristics of fluid-structure interaction in fluid-filled pipes[J]． Noise and Vibration Control，2013，33(5): 10-14．

[7] 葉紅玲, 邵沛澤, 陳寧, 等. 流固耦合輸流管系統(tǒng)的動(dòng)力學(xué)分析及參數(shù)影響[J]. 北京工業(yè)大學(xué)學(xué)報(bào), 2015, 41(2): 167-173.

YE Hongling, SHAO Peize, CHEN Ning, et al. Dynamic analysis and parameters’ influences on fluid-structure interaction in a fluid-filled pipes system[J]. Journal of Beijing University of Technology, 2015, 41(2): 167-173.

[8] TIJSSELING A S, LAVOOIJ C S W. Waterhammer with fluid-structure interaction[J]. Applied Scientific Research, 1990, 47(3): 273-285.

[9] LAVOOIJ C S W, TIJSSELING A S. Fluid-structure interaction in liquid-filled piping systems[J]. Journal of Fluids and Structures, 1991, 5(5): 573-595.

[10] 楊超, 范士娟. 管材參數(shù)對(duì)輸液管流固耦合振動(dòng)的影響[J]. 振動(dòng)與沖擊, 2011, 30(7): 210-213.

YANG Chao，F(xiàn)AN Shijuan. Influence of pipe parameters on fluid-structure coupled vibration of a fluid-conveying pipe[J]. Journal of Vibration and Shock, 2011, 30(7): 210-213.

[11] 孫玉東, 王鎖泉, 劉忠族,等. 液-管耦合空間管路系統(tǒng)振動(dòng)噪聲的有限元分析方法[J]. 振動(dòng)工程學(xué)報(bào), 2005, 18(2):149-154.

SUN Yudong, WANG Suoquan, LIU Zhongzu,et al.Unified finite element method for analyzing vibration and noisein 3D piping system with liquid-pipe coupling[J]. Journal of Vibration Engineering, 2005, 18(2):149-154.

[12] SREEJITH B, JAYARAJ K, GANESAN N, et al. Finite element analysis of fluid-structure interaction in pipeline systems[J]. Nuclear Engineering and Design, 2004, 227(3): 313-322.

[13] WANG Z M, TAN S K. Vibration and pressure fluctuation in a flexible hydraulic power systemon an aircraft[J]. Computers and Fluids, 1998, 27(1): 1-9.

[14] AHMADI A, KERAMAT A. Investigation of fluid-structure interaction with various types of junction coupling[J]. Journal of Fluids and Structures, 2010, 26(7): 1123-1141.

[15] CHEN Y G, PRICE W G. Numerical simulation of liquid sloshing in a partially filled container with inclusion of compressibility effects[J]. Physics of Fluids, 2009, 21(11): 112-105.

[16] OUYANG L B, AZIZ K. Transient gas-liquid two-phase flow in pipes with radial influx or efflux[J]. Journal of Petroleum Science and Engineering, 2001, 30(3): 167-179.

[17] LU D M, SIMPSON H C, GILCHRIST A. The application of split-coefficient matrix method to transient two phase flows[J]. International Journal of Numerical Methods for Heat and Fluid Flow, 1996, 6(3): 63-76.

[18] TIJSSELING A S. Exact solution of linear hyperbolic four-equation system in axial liquid-pipe vibration[J]. Journal of Fluids and Structures, 2003, 18(2): 179-196.