CN105260806A - 一种管路系统的流固耦合动力学特性预测方法 - Google Patents

Abstract

本发明涉及流固耦合预测领域，具体涉及一种管路系统的流固耦合动力学特性预测方法，其包括，根据管路系统中每个管路元件的属性，将管路系统离散为直管元件、弯管元件和管路附件元件三种类型；根据离散后的管路元件的类型，建立各管路元件相应的场传递流固耦合矩阵或点传递流固耦合矩阵；将离散后的管路系统重新组合为若干传递矩阵单元体，且根据所述管路系统预测的管壁弯曲波最小波长或者预测的上限频率计算传递矩阵单元体的最大特征长度；预测管路系统任意位置的流固耦合动力学特性。本发明能够消除管路系统流固耦合动力学问题数值溢出和预测结果不稳定问题，预测过程简便，非常有利于编程计算，并且预测精度高。

Description

一种管路系统的流固耦合动力学特性预测方法

技术领域

本发明涉及一种管路系统的流固耦合动力学特性预测方法，更具体的，本发明涉及一种用于含有较大场传递管路元件的大型复杂管路系统的流固耦合动力学问题预测。

背景技术

大型复杂管路系统广泛应用于船舶与海洋工程、石油化工、能源与电力工业等诸多领域中，由于管路系统中存在泵、阀门、分支管路等内部激励以及基础或铺板振动等外部激励，因而会诱发管路结构及管内流体的振动，这些管路结构及流体的振动还会相互耦合而产生流固耦合(FSI)振动，因此这些看似安静的管路系统在输送流体的同时，还将振动能量沿管壁和管内流体传播到管路系统的各个位置，造成管路及与之相连的其他元件和精密仪器的破坏，影响管路系统安全和管路动力系统的正常运行，严重时会造成巨大的经济损失，如何实现管路系统流固耦合动力学特性的快速高精度预测一直是国内外研究的热点和难点。

目前管路系统的流固耦合动力学特性的预测方法较多，常见的方法有四种：(1)特征线方法(MOC)，MOC是一种时域数值求解方法，适用于简单管路的时域响应计算，特别是流体压力波的瞬态响应分析，但MOC需要在时间和空间中离散求解，较为复杂，且很难考虑管道的弹性支撑条件，多段管路计算时还存在插值误差和不同特征线的相交等问题；(2)有限元法(FEM)，FEM是通过节点、单元构建管道及流体的运动形式，计算效率取决于模型尺寸的大小或分析频率的高；(3)特征线-有限元法方法(MOC-FEM)，MOC-FEM是通过MOC求解管内流体运动，通过FEM求解管道结构运动，然后根据接触面的平衡关系实现FSI问题的求解，MOC-FEM为输流管路动力学问题的求解提供了一个新的途径，但该方法的求解效率较MOC低得多；(4)传递矩阵方法(TMM)，TMM是通过波动方程直接描述管路系统的运动。与其他方法相比，TMM求解过程简单，易于编程且计算效率高，因而在管路系统结构振动预测，管内流体压力脉动的传播特性分析及流固耦合动力学预测等问题中具有较多应用，但是目前此方法尚难以应用于具有较长管路元件的大型管路系统的流固耦合计算中，原因之一是TMM在计算长管路时会出现计算数值溢出和预测结果不稳定的现象。

针对TMM的计算不稳定、积累数值误差大及数值溢出等“病态问题”，芮筱亭等人在《多体系统传递矩阵法及其应用》一书中集中介绍和讨论了传递矩阵的无量纲化，特征值问题的改进算法以及Riccati传递矩阵方法等几种常见的处理方法(芮筱亭，贠来峰，陆毓琪，何斌，王国平.多体系统传递矩阵法及其应用[M].北京：科学出版社，2008)。胡培民通过一个一维数组和两个变量来存储一个数据的方法解决了超精度数据的存储问题(胡培民.传递矩阵法在高频振动分析中的应用[J].振动与冲击，1996,15(4):50-52.)。Liu和Li(LiuGongmin,LiYanhua.Vibrationanalysisofliquid-filledpipelineswithelasticconstraints[J].JournalofSoundandVibration,2011,330(13):3166-3181.)根据管路系统的流固耦合动力学特性，推导得到了12个一元四阶常微分方程和2个一元二阶常微分方程，进而求解得到了管路系统的频域解析解，该方法在一定程度上能够减少由于矩阵变化带来的计算误差，上述这些的改进措施，在一定程度上均提高了传递矩阵在大型复杂管路系统流固耦合动力学特性预测精度和预测结果的稳定性，但利用这些方法在计算管路系统的流固耦合响应时并无优势，并且难以实现模块化编程预测。

国内外目前利用传递矩阵方法进行管路系统动力学预测的流程是：首先分析管路系统的构成形式，然后将各管路元件的数理模型带入管路系统中，并根据各管路元件连接处平衡及传递关系，将其组合为整体传递矩阵形式，进而利用定解条件，预测管路系统的动力学特性，但该方法在预测管路系统的流固耦合振动响应时会出现计算数值溢出和预测结果不稳定的现象。

发明内容

针对现有技术中存在的缺陷，本发明的目的是提供一种精确、稳定、高效的管路系统的流固耦合动力学特性预测方法，其包括以下步骤：

步骤一、根据管路系统中每个管路元件的属性，将管路系统离散为直管元件、弯管元件和管路附件元件三种类型；

步骤二、根据离散后的管路元件的类型，建立各管路元件相应的场传递流固耦合矩阵或点传递流固耦合矩阵；

步骤三、将离散后的管路系统重新组合为若干传递矩阵单元体，且根据所述管路系统预测的管壁弯曲波最小波长或者预测的上限频率计算传递矩阵单元体的特征长度范围；

步骤四、结合管路系统边界处流体和结构已知条件组成边界条件矩阵和外部激励列向量，利用传递矩阵方法，预测管路系统任意位置的流固耦合动力学特性。

在上述技术方案的基础上，离散管路系统时，管路系统包括场传递元件和点传递元件，所述场传递元件包括直管元件和弯管元件，所述点传递元件为管路附件元件。

在上述技术方案的基础上，所述管路附件元件包括管路支撑、分支接头和阀门。

在上述技术方案的基础上，根据所述管路系统预测的管壁弯曲波最小波长计算传递矩阵单元体的最大特征长度的计算方法为：

其中，L_max表示重新组合后的单元体的最大特征长度，C表示比例系数，表示与管壁杨氏模量(E)，管路内径(R_i)，管路外径(R_o)，管内流体密度(ρ_f)，管壁密度(ρ_p)和泊松比(υ)相关的计算式，λ_min为预测的管壁的弯曲波最小波长。

在上述技术方案的基础上，根据所述管路系统预测的上限频率计算传递矩阵单元体的最大特征长度的计算方法为：

其中，表示组合后的单元体的最大特征长度，C表示比例系数，表示一个与管壁杨氏模量(E)，管路内径(R_i)，管路外径(R_o)，管内流体密度(ρ_f)和管壁密度(ρ_p)相关的计算式，s_max表示可以预测的上限频率。

在上述技术方案的基础上，对于任意的管路系统,其任意考察位置的状态向量的预测方法为：

${\mathrm{\Φ}}_{arb}=U_{arb}{U}_{tq}{\left[\begin{array}{c}D\\ U_{tr}\end{array}\right]}^{-1}{F}_{exi}$

其中，Φ_arb表示在管路系统中选取的任意考察位置的状态向量，U_arb表示从任意考察位置所在的单元体的始端到考察位置的传递矩阵，U_tq表示提取矩阵,其作用是将考察位置所在单元体的始端的状态矩阵从管路系统的整体传递矩阵中提取出来,D表示由流体已知条件和结构已知条件组成的边界条件矩阵,U_tr表示由各单元体共同组成的整体传递矩阵，[]^-1表示矩阵的逆阵，F_exi表示管路系统的外部激励列向量。

特别地,当管路系统由单个单元体组成时，计算管路系统任意考察位置的状态向量的方法为：

${\mathrm{\Φ}}_{arb}=\left[\begin{array}{c}{0}_{N\×N}\\ I_{N\×N}\\ 0_{N\×N}\end{array}\right]{\left[\begin{array}{ccc}{D}_{N\×3N}& & \\ I_{N\×N}& U{1}_{N\×N}& 0_{N\×N}\\ 0_{N\×N}& I_{N\×N}& U{2}_{N\×N}\end{array}\right]}^{-1}{F}_{3N\×1}$

其中，N＝2，3，4，5，6，7……，式中，Φ_arb表示在管路系统中选取的任意考察位置的状态向量，0_N×N表示N×N维的零矩阵，I_N×N表示N×N维的单位矩阵，D_N×3N表示N×3N维的边界条件矩阵，U1_N×N表示从管路始端到考察位置的N×N维的传递矩阵，U2_N×N表示从考察位置到管路末端的N×3N维传递矩阵，[]^-1表示矩阵的逆阵，F_3N×N表示3N×N维的激励列向量，N表示用于描述管路系统流固耦合振动的变量数，其值为大于等于2的自然数。

与现有技术相比，本发明以重组优化后的传递矩阵单元体为计算对象，可快速地建立管路系统的流固耦合动力学问题的预测模型，通过划分单元体分段计算的技术，能够消除管路系统流固耦合动力学问题数值溢出和预测结果不稳定问题，预测过程简便，非常有利于编程计算，并且预测精度高。

本发明所涉及的预测技术通过管路系统预测的管壁弯曲波最小波长或者预测的上限频率计算传递矩阵单元体的特征长度范围，其在保证了预测精度的同时，降低计算量，使得本发明的管路系统的流固耦合动力学特性预测方法更为简单、可以实现程式化建模、具有预测精度高和速度快等显著特点。

附图说明

图1为大型复杂管路系统流固耦合动力学预测步骤框图；

图2为管路模型示意图；

图3为传统传递矩阵方法计算的管路系统横向流固耦合振动响应；

图4为本发明预测的管路系统横向流固耦合振动响应；

图5为复杂管路系统示意图；

图6为管路系统A点x方向振动加速度响应曲线；

图7为管路系统B点x方向振动加速度响应曲线；

图8为管路系统C点x方向振动加速度响应曲线。

具体实施方式

以下结合附图及实施例对发明作进一步详细说明。

参见图1所示，本发明管路系统的流固耦合动力学特性预测方法具体包括以下步骤：

步骤一、根据管路系统中每个管路元件的属性，将管路系统离散为直管元件、弯管元件和管路附件元件三种类型；离散管路系统时，管路系统包括场传递元件和点传递元件，所述场传递元件包括直管元件和弯管元件，所述点传递元件为管路附件元件。

步骤二、根据离散后的管路元件的类型，建立各管路元件相应的场传递流固耦合矩阵或点传递流固耦合矩阵；其具体为，对场传递元件建立场传递流固耦合矩阵，对点传递元件建立点传递流固耦合矩阵。

步骤三、将离散后的管路系统重新组合为若干传递矩阵单元体，且根据所述管路系统预测的管壁弯曲波最小波长或者预测的上限频率计算传递矩阵单元体的特征长度范围；

当已知预测的管壁弯曲波最小波长时，则采用最小波长方法计算传递矩阵单元体的特征长度范围；根据所述管路系统预测的管壁弯曲波最小波长计算传递矩阵单元体的最大特征长度的计算方法为：

其中，L_max表示重新组合后的单元体的最大特征长度，C表示比例系数，表示与管壁杨氏模量(E)，管路内径(R_i)，管路外径(R_o)，管内流体密度(ρ_f)，管壁密度(ρ_p)和泊松比(υ)相关的计算式，λ_min为预测的管壁的弯曲波最小波长。为了取得更好的计算效果，传递矩阵单元体的特征长度为：时为最佳。

当已知预测的上限频率时，则采用上限频率方法，根据所述管路系统预测的上限频率计算传递矩阵单元体的最大特征长度的计算方法为：

其中，L_max表示组合后的单元体的最大特征长度，C表示比例系数，表示一个与管壁杨氏模量(E)，管路内径(R_i)，管路外径(R_o)，管内流体密度(ρ_f)和管壁密度(ρ_p)相关的计算式，s_max表示可以预测的上限频率。为了取得更好的计算效果，传递矩阵单元体的特征长度为：时为最佳。

以重组优化后的传递矩阵单元体为计算对象，可快速地建立管路系统的流固耦合动力学问题的预测模型，通过划分单元体分段计算的技术，能够消除管路系统流固耦合动力学问题数值溢出和预测结果不稳定问题，预测过程简便，非常有利于编程计算，并且预测精度高。

步骤四、结合管路系统边界处流体已知条件和结构已知条件组成边界条件矩阵和外部激励列向量，利用传递矩阵方法，预测管路系统任意位置的流固耦合动力学特性。其具体方法为：

对于任意的管路系统,其任意考察位置的状态向量的预测方法为：

${\mathrm{\Φ}}_{arb}=U_{arb}{U}_{tq}{\left[\begin{array}{c}D\\ U_{tr}\end{array}\right]}^{-1}{F}_{exi}$

其中，Φ_arb表示在管路系统中选取的任意考察位置的状态向量，U_arb表示从任意考察位置所在的单元体的始端到考察位置的传递矩阵，U_tq表示提取矩阵,其作用是将考察位置所在单元体的始端的状态矩阵从管路系统的整体传递矩阵中提取出来,D表示由流体已知条件和结构已知条件组成的边界条件矩阵,U_tr表示由各单元体共同组成的整体传递矩阵，[]^-1表示矩阵的逆阵，F_exi表示管路系统的外部激励列向量。

特别地,当管路系统由单个单元体组成时，计算管路系统任意位置的状态向量的方法为：

${\mathrm{\Φ}}_{arb}=\left[\begin{array}{c}{0}_{N\×N}\\ I_{N\×N}\\ 0_{N\×N}\end{array}\right]{\left[\begin{array}{ccc}{D}_{N\×3N}& & \\ I_{N\×N}& U{1}_{N\×N}& 0_{N\×N}\\ 0_{N\×N}& I_{N\×N}& U{2}_{N\×N}\end{array}\right]}^{-1}{F}_{3N\×1}---\left(3\right)$

其中，N＝2，3，4，5，6，7……，式中，Φ_arb表示在管路系统中选取的任意考察位置的状态向量，0_N×N表示N×N维的零矩阵，I_N×N表示N×N维的单位矩阵，D_N×3N表示N×3N维的边界条件矩阵，U1_N×N表示从管路始端到考察位置的N×N维的传递矩阵，U2_N×N表示从考察位置到管路末端的N×3N维传递矩阵，[]^-1表示矩阵的逆阵，F_3N×N表示3N×N维的激励列向量，N表示用于描述管路系统流固耦合振动的变量数，其值为大于等于2的自然数。

最后，利用数字信号处理技术和MATLAB程序，对计算结果进行后处理并绘制管路系统流固耦合动力学响应曲线，实现大型复杂管路系统的流固耦合动力学特性的高精度可视化预测。

本发明所涉及的预测技术通过管路系统预测的管壁弯曲波最小波长或者预测的上限频率计算传递矩阵单元体的特征长度范围，其在保证了预测精度的同时，降低计算量，使得本发明的管路系统的流固耦合动力学特性预测方法更为简单、可以实现程式化建模、具有预测精度高和速度快等显著特点。

下面通过具体的实施例来说明本发明管路系统的流固耦合动力学特性预测方法与现有技术的技术优势。

实施例1：

如图2所示，一根两端自由的充水长管路，管路系统的特性参数如表1所示，在管道B端施加沿y方向的单位激励。

表1管路结构及内部流体的特性参数

如果将管路作为一根直管整体计算(即一个单元体)，计算得到的B端的流固耦合响应曲线如图3所示，在255Hz以后，管道横向振动响应结果开始变得不稳定。

图3出现不稳定现象是由于计算单元体的特征长度超出公式(1)或公式(2)所规定的最大特征长度，根据图3，如果将分析频率范围设定为0-1000Hz，则对应的分析的上限频率为1000Hz，根据公式(2)，需要将图2中的管路分成AC和BC两个单元体进行分段计算，利用本发明计算得到的B端的流固耦合响应曲线如图4所示，与图3相比，通过本发明计算得到的图4在高频段的计算精度明显提高，说明依据上限频率方法中公式(2)确定的最大特征长度，将管路通过分段计算的方式消除流固耦合振动不稳定的问题是有效可行的。

实施例2：

本例以图5所示的复杂空间管路系统为实施对象，该管路结构材料的基本参数为：弹性模量为200GPa，泊淞比为0.3，密度7800kg/m³。各边界端固支约束，各弹性支撑垂向(x方向)悬吊，刚度均为6.78×10³N/m。图中1为DN100的闸阀，2为DN100的截止阀，3为DN65的闸阀，4为DN65的截止阀。N₁、N₂、N₃三处均连接有六自由度的弹簧，各向刚度为1×10⁹N/m(近似固支边界条件)，N₄处连接有三个自由度的线性弹簧，各向弹簧刚度为8×10⁶N/m,管路系统预测的管壁弯曲波最小波长为17m。

在分支点处施加激励，任意选取管路上的A、B、C三点作为的振动响应考察点，首先采用本发明最小波长方法中的公式(1)计算传递矩阵单元体的特征长度，然后根据计算的特征长度划分单元体并预测考察点的振动响应，预测如图6、图7和图8所示，新发明的预测方法求解复杂管路系统的振动响应与有限元方法(FEM)仿真计算结果吻合良好，说明了利用最小波长方法确定的最大特征长度，将管路通过分段计算的方式能够避免流固耦合振动预测不稳定的问题，同时也进一步验证了本发明在复杂空间管路系统流固耦合动力学预测的适用性。

通过本发明在以上例子的实施结果可以看出，本发明可用于大型复杂管道流固耦合动力学问题频域响应预测，并且具有较高的预测精度。

本发明不局限于上述实施方式，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围之内。本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (7)

1.一种管路系统的流固耦合动力学特性预测方法，其特征在于：其包括以下步骤，

步骤一、根据管路系统中每个管路元件的属性，将管路系统离散为直管元件、弯管元件和管路附件元件三种类型；

步骤二、根据离散后的管路元件的类型，建立各管路元件相应的场传递流固耦合矩阵或点传递流固耦合矩阵；

步骤三、将离散后的管路系统重新组合为若干传递矩阵单元体，且根据所述管路系统预测的管壁弯曲波最小波长或者预测的上限频率计算传递矩阵单元体的特征长度范围；

步骤四、结合管路系统边界处流体已知条件和结构已知条件组成边界条件矩阵和外部激励列向量，利用传递矩阵方法，预测管路系统任意位置的流固耦合动力学特性。

2.如权利要求1所述的一种管路系统的流固耦合动力学特性预测方法，其特征在于：离散管路系统时，管路系统包括场传递元件和点传递元件，所述场传递元件包括直管元件和弯管元件，所述点传递元件为管路附件元件。

3.如权利要求2所述的一种管路系统的流固耦合动力学特性预测方法，其特征在于：

所述管路附件元件包括管路支撑、分支接头和阀门。

4.如权利要求1所述的一种管路系统的流固耦合动力学特性预测方法，其特征在于，

根据所述管路系统预测的管壁弯曲波最小波长计算传递矩阵单元体的最大特征长度的计算方法为：

其中，L_max表示重新组合后的单元体的最大特征长度，C表示比例系数，表示与管壁杨氏模量(E)，管路内径(R_i)，管路外径(R_o)，管内流体密度(ρ_f)，管壁密度(ρ_p)和泊松比(υ)相关的计算式，λ_min为预测的管壁的弯曲波最小波长。

5.如权利要求1所述的一种管路系统的流固耦合动力学特性预测方法，其特征在于：根据所述管路系统预测的上限频率计算传递矩阵单元体的最大特征长度的计算方法为：

其中，L_max表示组合后的单元体的最大特征长度，C表示比例系数，表示一个与管壁杨氏模量(E)，管路内径(R_i)，管路外径(R_o)，管内流体密度(ρ_f)和管壁密度(ρ_p)相关的计算式，s_max表示可以预测的上限频率。

6.如权利要求1所述的一种管路系统的流固耦合动力学特性预测方法，其特征在于：对于任意的管路系统,其任意考察位置的状态向量的预测方法为：

${\mathrm{\Φ}}_{arb}=U_{arb}{U}_{tq}{\left[\begin{array}{c}D\\ U_{tr}\end{array}\right]}^{-1}{F}_{exi}$

其中，Φ_arb表示在管路系统中选取的任意考察位置的状态向量，U_arb表示从任意考察位置所在的单元体的始端到考察位置的传递矩阵，U_tq表示提取矩阵,其作用是将考察位置所在单元体的始端的状态矩阵从管路系统的整体传递矩阵中提取出来,D表示由流体已知条件和结构已知条件组成的边界条件矩阵,U_tr表示由各单元体共同组成的整体传递矩阵，[]^-1表示矩阵的逆阵，F_exi表示管路系统的外部激励列向量。

7.如权利要求6所述的一种管路系统的流固耦合动力学特性预测方法，其特征在于：当管路系统由单个单元体组成时，通过下述的方法计算管路系统任意考察位置的状态向量：

${\mathrm{\Φ}}_{arb}=\left[\begin{array}{c}{0}_{N\×N}\\ I_{N\×N}\\ 0_{N\×N}\end{array}\right]{\left[\begin{array}{ccc}{D}_{N\×3N}& & \\ I_{N\×N}& U{1}_{N\×N}& 0_{N\×N}\\ 0_{N\×N}& I_{N\×N}& U{2}_{N\×N}\end{array}\right]}^{-1}{F}_{3N\×1}$

其中，N表示用于描述管路系统流固耦合振动的变量数，其值为大于等于2的自然数，Φ_arb表示在管路系统中选取的任意考察位置的状态向量，0_N×N表示N×N维的零矩阵，I_N×N表示N×N维的单位矩阵，D_N×3N表示N×3N维的边界条件矩阵，U1_N×N表示从管路始端到考察位置的N×N维的传递矩阵，U2_N×N表示从考察位置到管路末端的N×3N维传递矩阵，[]^-1表示矩阵的逆阵，F_3N×N表示3N×N维的激励列向量。