基于分離系數矩陣差分法的輸流管道軸向耦合響應特性研究

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

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

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

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

1 數值模型及初邊值條件

1.1 數值模型

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

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

(1)

(2)

(3)

(4)

圖1 水箱-管道-閥門系統

1.2 邊界條件

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

(5)

式中：P0為一給定值。

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

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

(6)

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

(7)

式中：Af和As分別為管道流通截面積和管道環形截面積。

1.3 初始條件

閥門瞬時關閉前，管內流體流動為穩態，流速為常數，軸向管速及管道應力都為零，即：

v=v0

(8)

P=P0

(9)

(10)

σz=0

(11)

2 數值算法

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

(12)

其中，cf和cs分別為流體壓力波速和管道應力波速

(13)

(14)

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

(15)

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

(16)

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

(17)

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

Λ++Λ-

(18)

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

(19)

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

(20)

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

(21)

(22)

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

3 模型驗證

3.1 案例1

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

表1 Wilkinson 和Curtis實驗參數

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

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

表2 Delft水力學基準問題A模型參數

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

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

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

(a) z=20 m

(b) z=10 m

4 兩種不同耦合作用對響應的影響

4.1 流體速度響應

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

(a) z=20 m

(b) z=10 m

4.2 管道軸向振動速度響應

(a) z=20 m

(b) z=10 m

(a) 管道軸向振動速度

(b) 管道平均應力

4.3 管壁應力響應

(a) z=20 m

(b) z=10 m

5 結 論

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

對水箱-管道-閥門系統而言，不同約束條件管道主要表現出的流固耦合強度不同，閥門固支時以泊松耦合為主，閥門簡支時以泊松耦合和連接耦合的疊加效應為主，邊界條件對輸流管道流固耦合作用影響顯著。對比經典水錘方程計算結果發現，兩種耦合作用對管道軸向振動特性的影響不容忽視，泊松耦合主要影響振動響應幅值，而連接耦合不僅會影響振動幅值，同時也會影響振動頻率。

[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] 任建亭, 姜節勝. 輸流管道系統振動研究進展[J]. 力學進展, 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] 范士娟, 楊超. 輸水管道流固耦合振動數值計算[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] 席志德, 馬建中, 孫磊. 考慮流-固耦合效應的空間管道水錘方法研究[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] 姬賀炯,白長青,韓省亮．輸流管耦合動力學特性分析[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] 葉紅玲, 邵沛澤, 陳寧, 等. 流固耦合輸流管系統的動力學分析及參數影響[J]. 北京工業大學學報, 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] 楊超, 范士娟. 管材參數對輸液管流固耦合振動的影響[J]. 振動與沖擊, 2011, 30(7): 210-213.

YANG Chao，FAN 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] 孫玉東, 王鎖泉, 劉忠族,等. 液-管耦合空間管路系統振動噪聲的有限元分析方法[J]. 振動工程學報, 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.