基于分离系数矩阵差分法的输流管道轴向耦合响应特性研究

张 挺， 谭志新， 张 恒， 范佳铭， 杨志强(. 福州大学 土木工程学院， 福州 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.