基于分離系數(shù)矩陣差分法的輸流管道軸向耦合響應(yīng)特性研究

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

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

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

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

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

1.1 數(shù)值模型

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

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

(1)

(2)

(3)

(4)

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

1.2 邊界條件

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

(5)

式中：P0為一給定值。

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

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

(6)

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

(7)

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

1.3 初始條件

閥門瞬時關(guān)閉前，管內(nèi)流體流動為穩(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)

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

(15)

將上式左右兩邊同時左乘Τ-1A-1

(16)

令Τ-1A-1BΤ=Λ，則對角矩陣Λ的對角元素即為矩陣A-1B的特征值，該特征值即為|B-λA|=0對應(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)

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

(19)

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

(20)

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

(21)

(22)

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

3 模型驗證

3.1 案例1

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

表1 Wilkinson 和Curtis實驗參數(shù)

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

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

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

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

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

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

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

(a) z=20 m

(b) z=10 m

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

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

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

(a) z=20 m

(b) z=10 m

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

(a) z=20 m

(b) z=10 m

(a) 管道軸向振動速度

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

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

(a) z=20 m

(b) z=10 m

5 結(jié) 論

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

對水箱-管道-閥門系統(tǒng)而言，不同約束條件管道主要表現(xiàn)出的流固耦合強度不同，閥門固支時以泊松耦合為主，閥門簡支時以泊松耦合和連接耦合的疊加效應(yīng)為主，邊界條件對輸流管道流固耦合作用影響顯著。對比經(jīng)典水錘方程計算結(jié)果發(fā)現(xiàn)，兩種耦合作用對管道軸向振動特性的影響不容忽視，泊松耦合主要影響振動響應(yī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)振動研究進展[J]. 力學(xué)進展, 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] 范士娟, 楊超. 輸水管道流固耦合振動數(shù)值計算[J]. 噪聲與振動控制, 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]. 核動力工程, 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] 姬賀炯,白長青,韓省亮．輸流管耦合動力學(xué)特性分析[J]．噪聲與振動控制,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)的動力學(xué)分析及參數(shù)影響[J]. 北京工業(yè)大學(xué)學(xué)報, 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ù)對輸液管流固耦合振動的影響[J]. 振動與沖擊, 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)振動噪聲的有限元分析方法[J]. 振動工程學(xué)報, 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.