基于中心线测量的山地管道力学分析软件开发

兰才富，赵 飞，赵 荣，吴承睿，郭家瑞，邓宝信，王 琳

（1.西南管道南宁输油气分公司，广西南宁 530200；2.西南石油大学，四川成都 610500）

0 引言

我国西南地区油气管道沿线地质条件复杂、地灾频发，极易引起管道变形失效，因此需要对管道变形、失稳等状态进行安全监测。然而，通过外部地质环境监测推测管道本体安全状态的方法，沿线工作量大、难度高，且不能直接体现管道本体的真实变形；传统的应力应变监测系统只能针对局部管段进行监测，全线应用的代价较大，而且多为事后监测，隐患监控滞后。因此，亟需开发对管道全线进行管道应力应变状态计算、风险发展趋势及安全评估的手段和工具[1-3]。

在此背景下，基于中心线测量结果，开发可对管道全线应变在空间上连续检测的管道本体安全状态评价方法：①利用单次测量结果，计算全线应力应变分布状态；②利用多次测量结果，分析管道应力应变发展趋势；③评价管道全线在不同工况下的可靠性，寻找出管道力学状态的危险环节[4-6]。研究基于管道中心线检测数据的管道应力应变算法，构建反映长输管道全线应力应变状态及其风险发展趋势的管道安全评估方法，开发相应工具。研究结果可实现油气管道的中心线、位移、应变应力的自主检测和风险等级评估，为管道运行和维修维护决策提供适用、可靠、高效的技术工具，保障管道安全。

1 数学模型及求解方法

通过db4 小波的6 层软阈值分解对角速度进行降噪处理，软件的核心算法是管道力学分析的应力应变算法[7]。实现有限元分析标准化和规范化的载体为单元，其搭建出的复杂结构即为有限元分析的对象。有限元分析最主要的内容就是研究单元，首先给出单元的节点位移和节点力，然后基于单元节点位移与节点力的相互关系，获得相应的刚度系数，进而得到单元的刚度方程，再将单元组装为整体刚度方程，即可得到整体结构的基于节点位移的整体平衡方程。

基本分析流程为：节点编号和单元划分→计算各单元刚度方程→组装各单元刚度方程→处理边界条件并求解→求支反力→求各单元其他力学量。

由于六面体单元的几何形状与管道结构相协调，因此与四面体相比，六面体构成环形单元的节点数目更少、网格质量更高，更便于在管道力学计算中实现对网格规模的控制。采用有限元法分析三维空间问题时，可将研究对象划分为如图1 所示的八节点六面体单元[8]。六面体单元模型由8 个节点组成，每个节点有3 个自由度，单元共计24 个自由度。

图1 八节点六面体单元模型

其节点位移列阵qe与节点力列阵Pe分别表示为：

八节点单元的每个方向的位移场可以设定8 个待定系数，根据确定位移模式的基本原则（从低阶到高阶、唯一确定性），选取该单元的位移模式：

由节点条件确定待定系数（ai，bi，ci），i=0，1，2，……7。代入上式，可计算整理得到该单元的形状函数矩阵：

得到该单元的形状函数矩阵后，按照有限元分析的标准过程推导相应的几何矩阵、刚度矩阵、节点等效载荷矩阵以及刚度方程。

由弹性力学平面问题的几何方程，可得到单元应变场：

由弹性力学中平面问题的物理方程，可得到单元的应力表达，然后计算单元的势能，就能得到单元的刚度矩阵与等效节点载荷矩阵：

将单元势能对节点位移取一阶极值，即可得到单元的刚度方程：

由六面体环形单元构造出管段的刚度矩阵K、弹性矩阵D、几何矩阵B 后，进一步得管道总体的刚度方程：

管段总体位移向量后，由弹性力学空间问题的物理方程可以得到应力和应变场：

2 软件设计理念及特色

为了保证软件的实用性，详细调研了国内外相关软件（Abaqus、ANSYS）功能，明确软件的设计目标为：吸收权威软件的优势功能，围绕管道安全检/监测管理的业务需求扩充（消弱）其部分功能，保证计算精度可靠的情况下，增强软件的适用性；保证软件操作简便、直观。为实现该目标，突破了常规力学分析软件思路，引入新的软件设计理念：采用模块化编程，有机地吸取有限元软件与信号分析软件各自的优势，实现计算软件的精确性与大数据处理能力的有机结合。依据该理念，完成了西南管道基于中心线测量的山地管道力学分析软件设计。

（1）具有较好的通用性，其数据预处理、中心线分析、特征识别、力学分析、安全评估模块的功能可以单独使用，满足不同用户的多方面需求。

（2）应用领域及用户包括专业设计人员、专业技术人员、完整性管理人员及现场工人等用户。

（3）使用逻辑与用户需求相契合，满足不同层次用户需求，既可以通过导入数据进行快速计算得到结果，也可以用于对管道力学性能的详细计算和分析。

（4）具有良好的人机界面，在自主开发的界面进行数据集成与处理，采用大量图形和表格引导输入数据和显示计算结果，多页面、积木式窗口设计，操作简便，与各类常规软件（Word、Excel、Access 等）之间可建立数据接口。

（5）采用高效的计算方法，计算效率要优于通用有限元软件，并且操作界面简洁、结果展示详尽，具备一键报告生成功能。

（6）软件界面有操作提示，并且对用户的每一步操作都进行操作日志记录。

3 软件应用效果分析

软件能否得以广泛应用并发挥其效能，关键是中心线数据预处理、管道特征识别、应力应变算法的准确性。

3.1 数据预处理效果

小波变化作为一种有效的时频分析方法，能够在保留时间信息的前提下对不同频率的信号进行分析[9-10]。“db”系小波在惯性导航的降噪上被广泛的应用，通过“db4”小波6 层低频重构得到的角速度信息与原始角速度信息对比如图2 所示。

图2 降噪前后角速度对比

3.2 特征识别方法使用效果

连续敷设管道采用弯头来改变管道的走向，以满足实际的运输需求，在地质环境复杂的山区，弯头的使用更为频繁。为避免弯头引起的姿态信息改变进而影响管道的应力应变计算，开发基于检测数据的弯头识别方法[11]。经过弯头后，检测器的姿态信息发生显著变化，在角速度的信号中体检为局部的明显凸起，通过选取适当的阈值能够有效的对直管和弯头进行区分、通过曲率可以对弯头的类型进行区分。

表1 为50 km 实际管道本软件方法与第三方漏磁检测的弯头检测结果对比，在50 km 管道，本方法共识别出722 处冷弯及342 处热弯，共计1064 处弯头，同管段第三方漏磁检测结果为720 处冷弯、341 处热弯；本方法的冷弯合计722 处、热弯合计342 处，第三方检测的冷弯合计720 处、热弯合计341 处，误差仅为0.65%。

表1 软件识别弯头及第三方检测弯头对比

3.3 应力应变算法验证

采用本课题方法与Abaqus 计算结果进行对比验证。管材弹性模量210 GPa，密度为7850 kg/m3，泊松比0.3，管径为1016 mm，壁厚为16 mm，管长10 m。模拟集中载荷作用于水平管道中间顶部，表2 为不同工况下的集中载荷。

表2 不同工况下的集中载荷和最大应力

采用六面体环形单元计算得到的管道应力分布，与Abaqus 计算得到的管道应力分布结果吻合。

3.4 前50 km 安全评估结果

软件的部分界面如图3 所示，根据应力应变计算结果筛选得到的关注点共6 处，根据管道的应变均值及弯曲管道的长度，将弯曲管段分为局部弯曲和整段弯曲，长度划分界限为2/3 选取管道长度。除了应力应变计算结果外，通过匹配内检测数据还可以得到关注点的具体坐标信息，方便后续的实地考察及维护工作。

图3 软件界面

4 结论

基于MATLAB，开发了基于管道中心线检测数据的管道应力应变评价软件，能够实现以下功能：

（1）对原始姿态信息进行降噪处理，以此提高应力应变计算结果的准确性。

（2）有效识别中心线数据中的弯头数据，通过中心线数据能够识别弯头的类型，曲率半径及弯头的走向。

（3）采用六面体法沿线计算管道的应力应变。

（4）通过应力应变计算结果、弯曲变形管段长度对发生变形的危险管道进行识别及分类。

通过与原始姿态信号、第三方检测结果、Abaqus仿真计算的对比，分别验证了该软件数据预处理模块、特征识别模块、应力应变计算模块的正确性。该软件可通过中心线检测结果，快速给出全线的管道应力应变水平，及时发现管道中有应力应变集中的管段，为长输管道的安全评估提供有效的方法。