任意弹性边界支承输流管路系统耦合振动特性分析

马 腾， 杜敬涛， 许得水， 代 路， 张 赟

(1.哈尔滨工程大学动力与能源工程学院, 黑龙江 哈尔滨 150001；2.武汉第二船舶设计研究所热能动力技术重点实验室， 湖北 武汉 430205)

引 言

输流管路广泛存在于船舶与海洋工程、石油化工、生物工程、航空航天、核工业等诸多领域。流固耦合问题是引起管路振动过大以及疲劳断裂的重要因素，全面理解输流管路耦合振动特性是正确设计相关复杂管路系统的前提和基础。为此，国内外众多学者针对输流管路耦合振动特性分析开展了大量的研究[1-3]，针对输流管路动力学特性分析提出了各种建模分析方法，例如：分离变量法、有限元法、传递矩阵法、微分求积法等。

梁波等[4]采用有限元法建立了简支条件下管路系统的动力学特性分析模型，研究了流体流速、质量比等参数对固有频率的影响；金基铎等[5]通过数值模拟的方法探讨了脉动流作用下两端固定支承管路系统的参数共振问题；Panda和Kar[6]对脉动流作用下两端铰支管路系统的内共振参数问题进行了分析。上述研究多是针对悬臂、铰支、固定等理想边界条件，从工程角度，管路系统支承难以做到严格的刚性情况，边界支承方式更加倾向于弹性支承，为此，弹性支承边界管路耦合振动受到更多的关注。Noah和Hopkin[7]以一端固定，一端有线性弹簧和扭转弹簧支承的管路系统为模型，对系统的失稳形式进行了研究。马小强等[8]采用传递矩阵法对置于弹性地基上，中间有弹性支承的管路系统模型进行了稳定性分析。Chellapilla和Simha[9]针对双参数弹性地基输流管路，分析了失稳临界流速问题。包日东等[10]采用微分求积法对弹性支承的输流管路进行了研究，分析了不同系统参数对于管路系统固有频率及稳定性的影响。

对于固定、铰支、自由等经典边界下的管路振动系统四阶线性常微分方程而言，通常可以得到固有振动特性分析的系统超越特征方程，但是，当边界条件改变时，需要对问题进行重新推导，对于弹性边界约束模型，求解则更为困难。最近，杜敬涛等[11-13]提出了一套改进傅里叶级数方法进行任意弹性边界板壳结构面内与弯曲振动分析，数值结果对所提出方法的有效性与可靠性进行了充分验证。本文针对输流管路耦合振动问题，将所提出的改进傅里叶级数方法拓展至管路流固耦合振动特性分析。管路横向振动位移场采用标准的傅里叶级数和边界光滑函数组合展开，进而精确地满足输流管路耦合振动方程与弹性边界条件，通过联立控制方程和边界条件得到系统特征方程，所有的模态参数可以通过求解一个标准的矩阵特征值问题而统一得到。随后，通过给出数值算例对所建立模型的正确性和有效性进行了验证，在此基础上，讨论分析了边界约束刚度对耦合振动特性的影响规律。

1 输流管路耦合振动方程与弹性支承边界

考虑如图1所示的细长输流管路耦合振动问题，通过在两端边界引入4种约束弹簧实现对边界条件的模拟，任意经典边界条件及其组合可以通过设置约束刚度系数得到。忽略非线性项及摩擦力的管路系统横向振动微分方程为

(1)

式中 前4项项分别为系统惯性力、流体的离心力、流体作用的陀螺力和弹性恢复力[14]；W为管路横向振动位移函数；M为单位长度流体的质量；m为单位长度管路的质量；U=const(定常流)为流体流速；E为弹性模量；μ为黏弹性系数；I为管路截面的惯性矩。

图1 复杂边界支承管路系统振动分析模型Fig．1 The vibration analysis model of pipeline system with arbitrary elastic boundary supports

利用弹性边界约束所对应的力平衡和位移协调关系，可以得到边界条件方程为

(2)

为了使运动微分方程不受单位影响而具有广义性，引入如下定义无量纲条件

(3)

管路系统横向振动无量纲挠度响应表达式为

w(x,t)=w(x)ejwt

(4)

式中w(x)为管路系统横向振动的无量纲挠度分布函数，ejωt为简谐时间因子，ω为无量纲固有频率。为了后续计算方便，取Ω=jω。由于黏弹性系数μ的数值非常小，对于结果的影响微乎其微，后文推导计算中暂时忽略该项。

利用上述无量纲定义，对方程(1)进行简化，得到无量纲化后的振动微分方程

(5)

相应的弹性边界条件式(2)变换为

(6)

2 改进傅里叶级数展开与求解

对于经典边界条件下管路横向振动无量纲挠度通常可以采用傅里叶级数法展开

(7)

式中λn=nπ。可以发现，该经典形式的1阶导数和3阶导数在管路两端将恒为零，然而，对于弹性边界约束，由式(6)可知，力平衡关系中该项并不应该恒为零。为了克服经典傅里叶级数展开微分在管路两端的不连续，在该式基础上进一步引入边界光滑函数，从而使得横向位移函数在整个求解域内[0, 1]足够光滑[12, 15]，此处，在原挠度方程的基础上增加4个辅助函数，即

(8)

式中ζ1(x),ζ2(x),ζ3(x),ζ4(x)为4个辅助函数。通过辅助函数的构造，位移函数在端点处的一阶导数，3阶导数将不再会恒等于零，从而具备表示复杂约束下边界情况的能力，进一步改善解的收敛性和精确性。此处，辅助函数构造为:

以x=0和x=1两端支承边界条件上的辅助函数为例，容易证明

‴4(1) = 1

而在边界上所有的其他1阶和3阶导数均等于零，从而成功地克服了振动位移在求解域内展开为傅里叶级数过程中在x=0和x=1边界上的不连续性。这里需要指出，边界光滑辅助函数的构造形式并不唯一，只要能够解决边界条件的求导失效问题即可。本文选取如上辅助函数简单形式，方便后续推导和计算。

将所构造的位移方程(8)代入方程(5)，并对左右两端分别乘以2cosλmx后在0到1上积分，化简后可写成如下形式

(10)

其中：

(11a)

(11b)

(11c)

(11d)

(11e)

其他方程中系数的表达式同理可得下面以第二项为例给出计算结果：

(12a)

(12b)

(12c)

(12d)

(12e)

将系数矩阵写成矩阵形式，同样以第二项为例：

Q2=2βu·

(13)

(14)

同理可以得到Q1，P1，Q3，P3，Q4，P4各项矩阵的形式。

管路系统的振动方程可写成如下形式

Ω2(Q1A+P1B)+Ω(Q2A+P2B)+

(Q3+Q4)A+(P3+P4)B=0

(15)

将位移方程(8)带入边界条件(6)中，可得：

(16a)

(16b)

(16c)

(16d)

进一步写为矩阵形式

HA=FB

(17)

式中

(18)

(19)

(20)

(21)

从而得到所构建改进傅里叶级数中标准级数系数与边界光滑函数系数之间的约束关系矩阵表达式为

F-1HA=B

(22)

将式(22)代入式(15)中得到

(K+ΩC+Ω2M)A=0

(23)

其中:

K=(Q3+Q4)+(P3+P4)F-1H

C=Q2+P2F-1H

M=Q1+P1F-1H

通过采用状态空间法可以对式(23)进行特征值问题求解，由于流体作用的原因，导致ω为复数，ω=Ω/j，其中ω的实部代表管路耦合振动系统的固有频率，虚部代表指数衰减。

3 数值结果与分析

本节中，将采用MATLAB对上述理论模型进行编程仿真，由于流体作用导致阻尼项的存在，应用状态空间法来求解管路系统的固有频率及特征模态方程。将首先验证本文方法在复杂边界条件管路系统振动特性方面的准确性和可靠性，在此基础上，讨论分析边界约束刚度变化对管路耦合振动固有特性的影响规律。

3.1 不同边界支承下管路系统固有频率的计算

表1 不同流速下两端铰支时前4阶固有频率的结果

Tab.1 Results of the first four natural frequencies of pinned-pinned pipe with different flow velocities

流速ω1ω2ω3ω4u=09.869639.478588.8270157.9160精确解[16]9.869639.478488.8264157.9137微分求积法[17]9.871639.486388.8442157.9454u=19.309739.022588.3984157.5009u=27.453637.632287.1050156.2520u=32.783135.229984.9233154.1579u=3.535(失稳临界流速)033.473083.3769152.6838u=4031.632381.8084151.1982u=5026.358477.6825147.3410u=6017.137672.4090142.5384u=710.108310.108365.7240136.7179u=810.676910.676857.0092129.7664u=910.020410.020443.9360121.4954u=10027.795127.7951111.5560

图2 两端铰支下前4阶固有频率随流速的变化关系Fig.2 Relationship between the first four orders of natural frequencies and flow velocities of pinned-pinned pipe

前4阶固有频率随流速的变化关系如图2所示。在稳定区域内，随着流速增加，固有频率随之降低，当流速u=3.535时，第1阶固有频率ω1=0，此时的流速uc称为失稳临界流速；当流速u=6.465时，第1阶固有频率与第2阶固有频率发生耦合；当u=9.616时，第2阶固有频率与第3阶固有频率发生耦合；当u=12.73时，第3阶固有频率与第4阶固有频率出发生耦合。

当边界条件为固定支承时，各边界支承的刚度系数分别为k0= 1010，K0=1010，k1=1010，K1=1010，得到不同流速下前4阶固有频率，结果如表2所示。

表2 不同流速下两端固定时前4阶固有频率的结果

Tab.2 Results of the first four natural frequencies of clamped-clamped pipe with different flow velocities

流速ω1ω2ω3ω4u=022.373361.6728120.9034199.8594精确解[16]22.373361.6728120.9034199.8594微分求积法[17]22.377861.6852120.9276199.8995u=122.020761.3318120.5508199.5020u=220.945060.3047119.4906198.4282u=319.085158.5784117.7157196.6333u=416.307856.1271115.2131194.1093u=512.295252.9050111.9623190.8437u=65.658648.8263107.9316186.8187u=6.465(失稳临界流速)046.5948105.7788184.6817u=7043.7139103.0717182.0091u=8037.107197.3017176.3790u=92.163926.773790.4760169.8761u=1018.960918.960982.2932162.4205

图3 两端固定下前4阶固有频率随流速的变化关系Fig.3 Relationship between the first four natural frequencies and flow velocities of clamped-clamped pipe

前4阶固有频率随流速的变化关系如图3所示。在稳定区域内，随着流速增加，固有频率随之降低，当流速u=6.465时，第1阶固有频率ω1=0，此时的流速uc称为失稳临界流速；当流速u=9.495时，第1阶固有频率与第2阶固有频率发生耦合；当u=12.3时，第2阶固有频率和第3阶固有频率发生耦合。

当边界条件为悬臂支承时，各边界支承的刚度系数分别为k0=1010，K0=1010，k1=0，K1=0，得到不同流速下前4阶固有频率，结果如表3所示。

表3 不同流速下悬臂支承时前4阶固有频率的结果

Tab.3 Results of the first four natural frequencies of cantilevered pipe with different flow velocities

流速ω1ω2ω3ω4u=03.516022.034561.6976120.9035精确解[16]3.516022.034561.6972120.9109微分求积法[17]3.516722.038961.7096120.9260u=13.323721.765961.3597120.5553u=22.691220.956160.3421119.5085u=31.227319.592558.6320117.7557u=3.273(失稳临界流速)019.121258.0419117.1529u=4017.650156.2043115.2836u=5015.067853.0126112.0708u=63.0263111.838848.9649108.0833u=75.781411.018943.8451103.2667u=87.429013.283236.940197.5285u=99.001520.953625.533090.6957u=104.756221.032328.553782.3931

前4阶固有频率随流速的变化关系图如图4所示。在稳定区域内，随着流速增加，固有频率随之降低，当流速u=3.273时，第1阶固有频率ω1=0，此时的流速uc称为失稳临界流速。

图4 悬臂支承下前4阶固有频率随流速的变化关系Fig.4 Relationship between the first four natural frequencies and flow velocities of cantilevered pipe

文献[16-17]中给出了两端铰支，两端固定，悬臂支承条件下，流速为0时前4阶固有频率的精确解和微分求积法的结果，同时给出不同边界条件下，前4阶固有频率随流速的变化情况。通过上述结果对比，可以验证本文方法的正确性；同时对比数据可以发现，相比于微分求积法，改进傅里叶级数法得到的结果更加精确。后续工作中，进一步结合模态测试与识别技术，可以将本文所建立模型应用于实际管路系统特性分析。

3.2 复杂约束对管路系统振动特性影响

3.2.1 自由状态到两端铰支

为了分析管路系统从自由状态到两端铰支时的振动特性，通过改变不同的弹簧刚度系数来模拟不同的边界支承，控制2个扭转刚度系数K0=0，K1=0，2个支承刚度系数k0，k1从0变到1010。选取前2阶固有频率进行分析，图5，6分别为第1阶固有频率，第2阶固有频率随刚度系数k0，k1的变化关系。

图5 第1阶固有频率随支承刚度系数的变化关系Fig.5 Relationship between the first natural frequency and bearing stiffness coefficient

图6 第2阶固有频率随支承刚度系数的变化关系Fig.6 Relationship between the second-order natural frequency and bearing stiffness coefficient

可以看出，当2个扭转刚度系数K0=0，K1=0时，管路系统的固有频率随支承刚度系数k0，k1的增加而增大，弹性支承刚度系数在大约0～1000的范围内时，固有频率随刚度系数的变化比较明显，即可以定义为影响固有频率变化的敏感区域；当弹性刚度系数大于1000以上时，固有频率的变化则相对微小。同时左右两端的弹性支承刚度系数k0，k1对系统的影响是近似相等的。因此给出在刚度变化敏感范围内固有频率随刚度系数、流速的变化情况，以第1阶固有频率为例，如图7所示。可以发现随着弹性支承刚度系数的增大，第1阶固有频率随之增大，随着流速增大，第1阶固有频率逐渐减小，直到达到失稳临界流速时，第1阶固有频率ω1=0；但失稳临界流速不随刚度系数的变化而发生改变。

图7 第1阶固有频率随支承刚度系数和流速变化情况Fig.7 Relationship among the first natural frequency, bearing stiffness coefficient and flow velocities

3.2.2 两端铰支到两端固定

为了分析管路系统从两端铰支到两端固定时的振动特性，令2个弹性支承刚度系数k0=1010，k1=1010。2个扭转刚度系数K0，K1从0变到1010。以第1阶固有频率为例进行分析，图8为第1阶固有频率随扭转刚度系数K0，K1的变化情况。

图8 第1阶固有频率随扭转刚度系数的变化关系Fig.8 Relationship between the first natural frequency and torsional stiffness coefficient

图9 第1阶固有频率随扭转刚度和流速变化情况Fig.9 Relationship among the first natural frequency, torsional stiffness coefficient and flow velocities

可以看出，当2个支承刚度系数k0=1010，k1=1010时，管路系统的固有频率随扭转刚度系数K0，K1的增加而增大，当扭转刚度系数在大约0-1000的范围内时，随着扭转刚度系数的变化，管路系统固有频率的变化很明显；当扭转刚度系数增大到1000以上时，固有频率的变化很微小。同时可以看出K0，K1对固有频率的影响是近似相等的。给出敏感区域内固有频率随扭转刚度系数和流速的变化情况，以第1阶固有频率为例，如图9所示。可以发现随着扭转刚度系数的增大，第1阶固有频率随之增大；随着流速增大，第1阶固有频率随之减小；同时可以看出，随着流速增大到一定程度，当扭转刚度系数较小时，会出现失稳的情况；随着扭转刚度系数的增大，失稳临界流速也会随之增大。因此，在一定范围内，增大扭转刚度系数会有利于管路系统的稳定性。

3.2.3 弹性支承刚度和扭转刚度同时变化

通过上文分析可知，无论是支承刚度还是扭转刚度，管路左右两侧的支承条件对系统特性的影响是近似相等的，为了研究边界支承更普遍的情况，取支承刚度系数k0=k1=k，扭转刚度系数K0=K1=K。图10为管路系统第1阶固有频率随刚度系数的变化情况。

图10 第1阶固有频率随支承和扭转刚度系数的变化情况Fig.10 Relationship among the first natural frequency, bearing stiffness coefficient and torsional stiffness coefficient

可以看出固有频率随刚度系数的增加而增大，当刚度系数增大到一定程度后，固有频率变化很小。同时可以看出，对于管路系统的固有频率，在频率较低时，弹性支承刚度系数k要比扭转刚度系数K的影响大。

4 结 论

本文采用一种改进傅里叶级数方法获得了任意弹性边界支承输流管路耦合振动系统模态参数预报的精确解。通过在管路系统两端引入边界约束弹簧，任意边界条件及其组合可以设置相应的刚度系数而统一得到。管路横向振动位移采用标准傅里叶级数附件边界光滑函数进行构建，从而使得弹性边界约束方程所需要的1阶和3阶位移导数在整个求解域内足够光滑，联立输流管路耦合振动方程和边界条件方程，获得输流管路模态特性分析的系统矩阵方程。

在数值算例中，应用MATLAB语言编程仿真，计算了不同弹性边界支承下管路系统的固有频率。与精确解和解析解的对比充分验证了本文方法的可靠性和有效性。相比文献中其他方法，本文模型统一考虑边界约束刚度设置，当边界约束刚度发生改变时，无需对理论模型或程序进行重新修改，进行相应设置边界刚度系数即可。在此基础上，计算分析了不同流速、边界支承刚度等重要参数对管路系统固有频率的影响规律。结果表明，管路系统由自由状态变为两端铰支、再由两端铰支变为两端固定的过程中，当流速达到失稳临界流速之前，即稳定区域内，管路系统的固有频率随流速的增大而减小；当第1阶固有频率ω1变为0时，此时的流速uc称之为失稳临界流速；当流速u增大到一定值时，第1阶固有频率ω1开始增大；当流速u达到某一定值后，第1阶固有频率ω1和第2阶固有频率ω2发生耦合。同时发现，管路系统的边界支承刚度系数对固有频率的影响存在一个敏感区域，刚度系数大约为0～1000左右。当刚度系数越过此区域时，将不再具有显著影响；在敏感区域内，固有频率会随着弹簧刚度系数的增大而增大；在固有频率较低时，相比于扭转刚度系数的影响，支承刚度系数对管路系统固有频率的影响更大。同时当管路系统由自由状态变到两端铰支的过程中，弹簧刚度系数的变化将不会影响到管路系统的失稳临界流速；而当边界条件由两端铰支到两端固定的过程中，失稳临界流速将会随着扭转刚度系数的增大而增大。因此，在一定范围内，增大扭转刚度系数会有利于管路系统的稳定性。

[1] PaïDoussis M P, Li G X. Pipes conveying fluid: a model dynamical problem[J]. Journal of Fluids & Structures, 1993, 7(2):137—204.

[2] 黄玉盈, 邹时智. 输液管的非线性振动,分叉与混沌——现状与展望[J]. 力学进展, 1998, 28(1):30—42.

Huang Yuying, Zou Shizhi. Advances and trends of nonlinear dynamics if pipes conveying fluid[J]. Advances in Mechanics, 1988, 28(1):30—42.

[3] 张立翔, 黄文虎, TIJSSELING A S. 输流管道流固耦合振动研究进展[J]. 水动力学研究与进展, 2000, 15(3):366—379.

Zhang Lixiang, Huang Wenhu, TIJSSELING A S. Review of FSI analysis of fluid-conveying pipes[J]. Journal of Hydrodynamics, 2000, 15(3):366—379.

[4] 梁 波, 唐家祥. 输液管道动力特性与动力稳定性的有限元分析[J]. 固体力学学报, 1993,14(2):167—170.

Liang Bo, Tang Jiaxiang. Analysis of the dynamic characteristics and stability of pipes conveying fluid by finite element method[J]. Acta Mechanica Solida Sinica, 1993,14(2):167—170.

[5] 金基铎, 宋志勇, 杨晓东. 两端固定输流管道的稳定性和参数共振[J]. 振动工程学报, 2004, 17(2):190—195.

Jin Jiduo, Song Zhiyong, Yang Xiaodong. Stability and parametric resonances of a clamped-clamped pipe conveying fluid[J]. Journal of Vibration Engineering, 2004, 17(2):190—195.

[6] Panda L N, Kar R C. Nonlinear dynamics of a pipe conveying pulsating fluid with parametric and internal resonances[J]. Nonlinear Dynamics, 2007, 49(1):9—30.

[7] Noah S T, Hopkins G R. Dynamic stability of elastically supported pipes conveying pulsating fluid[J]. Journal of Sound & Vibration, 1980, 71(1):103—116.

[8] 马小强, 向 宇, 黄玉盈. 求解弹性地基上任意支承输液直管稳定性问题的传递矩阵法[J]. 工程力学, 2004, 21(4):194—198.

Ma Xiaoqiang, Xiang Yu, Huang Yuying. A transfer matrix method for solving stability of pipes conveying fluid on elastic foundation with various end supports[J].Engineering Mechanics, 2004, 21(4):194—198.

[9] Chellapilla K R, Simha H S. Critical velocity of fluid-conveying pipes resting on two-parameter foundation[J]. Journal of Sound & Vibration, 2007, 302(1-2):387—397.

[10] 包日东, 闻邦椿. 微分求积法分析弹性支承输流管道的稳定性[J]. 东北大学学报(自然科学版), 2007, 28(7):1017—1020.

Bao Ridong, Wen Bangchun. Differential quadrature method to analyze stability of elastically-supported fluid conveying pipelines[J]. Journal of Northeastern University (Natural Science), 2007, 28(7):1017—1020.

[11] Du J, Li W L, Jin G, et al. An analytical method for the in-plane vibration analysis of rectangular plates with elastically restrained edges[J]. Journal of Sound & Vibration, 2007, 306(3-5):908—927.

[12] Li W L, Zhang X, Du J, et al. An exact series solution for the transverse vibration of rectangular plates with general elastic boundary supports[J]. Journal of Sound & Vibration, 2009, 321(1-2):254—269.

[13] 李文达,杜敬涛,杨铁军,等. 基于改进傅里叶级数方法的旋转功能梯度圆柱壳振动特性分析[J]. 哈尔滨工程大学学报,2016,03:388—393.

Li Wenda, Du Jingtao, Yang Tiejun, et al. Vibration characteristics analysis of the rotating functionally graded cylindrical shell structure using an improved Fourier series method[J]. Journal of Harbin Engineering University, 2016,03:388—393.

[14] 倪 樵, 黄玉盈, 陈贻平. 微分求积法分析具有弹性支承输液管的临界流速[J]. 计算力学学报, 2001, 18(2):146—149.

Ni Qiao, Huang Yuying, Chen Yiping. Differential quadrature method for analyzing critical flow velocity of the pipe conveying fluid with spring support[J]. Chinese Journal of Computational Mechanics, 2001, 18(2):146—149.

[15] W L， L I. Free vibrations of beams with general boundary conditions[J]. Journal of Sound & Vibration, 2000, 237(4):709—725.

[16] Thomson W. Theory of Vibration with Applications[M]. CRC Press, 1996.

[17] Ni Q, Zhang Z L, Wang L. Application of the differential transformation method to vibration analysis of pipes conveying fluid[J]. Applied Mathematics & Computation, 2011, 217(16):7028—7038.