基于CEL 方法的管道内检测器过弯机理分析

2022-10-25吴昱达陈新汉李涛屹

焊管 2022年10期

吴昱达，郑莉，陈新汉，刘洋，柴泽，李涛屹

（1. 北京华航无线电测量研究所，北京 100013;2. 中国石油大学机械与储运工程学院，北京 100100）

0 前言

管道运输作为目前油气运输安全、合理的方式之一，具有连续、快捷、输送量大、成本低等诸多优点。定期对管道进行检测和清管，是管道运营维护必不可少的流程[1-2]。管道内检测器作为长输管道在线无损检测的有效设备，具有安全、高效等优点，其在跟随管道输送介质运动的同时，可采集管道的相关信息，检测完成后，通过分析存储数据来确定管道内外壁的缺陷信息并对缺陷进行定位[3]。目前，管道机器人通过性主要局限在直管且多数的仿真研究都局限于二维模型，无法表征管道内检测器运行过程中真实的运动状态。因此，进行在线油气管道内检测器三维仿真研究对管道内检测器的实用化具有重要意义[4-6]。

压差式管道内检测器，主要采用驱动舱段作为动力装置，在管道内采用皮碗密封前后形成压差，支撑并推动检测器实现长距离行走，从而实现油气管道的在线内检测[7]。为了获得准确的检测结果，需要确保检测器运行速度在一个稳定的范围内。对于输气管道的检测来说，检测器在管道内正常的行走速度范围为2～7 m/s，对于输送液体的管道，正常的运行速度范围在1～5 m/s[8-9]，速度波动对检测结果会有重要的影响。

双节管道内检测器在管道内运行过程，主要是靠液体介质推动动力舱段，在动力舱段牵引下沿着管道运动。双节管道内检测器主要在万向节的作用下实现检测器的过弯动作，而在过弯过程中万向节的尺寸设计直接影响其在发挥作用时的动作状态，进而影响检测器的运动稳定性。为此，本研究建立双节检测器的有限元模型，基于耦合欧拉-拉格朗日（CEL）方法，建立流-固耦合模型，研究检测器在管道内运动时万向节转动角度的变化规律，对确定最佳的万向节结构以及尺寸具有重要的指导意义。

1 基于CEL 方法的流-固耦合模型建立

1.1 耦合欧拉-拉格朗日（CEL）方法

ABAQUS 中的耦合欧拉-拉格朗日（CEL）方法结合了欧拉方法和拉格朗日方法各自的优点。 ABAQUS 可提供一种EC3D8R 体积单元来描述欧拉流体材料模型。欧拉流体可部分或者完全被流体材料填充[10]，在管道内检测器流-固耦合仿真中，流体充满管道，但是检测器占据的欧拉体积为空。欧拉体积分数（EVF）代表了每一个欧拉单元填充的情况，如果欧拉单元完全被材料充满，则欧拉体积分数为1，相反，如果是空的单元，则欧拉体积分数为0。 ABAQUS/CAE提供的体积分数划分工具，方便了用户定义欧拉材料初始边界[11-12]。在进行CEL 分析时，要注重网格的划分精度、大小等因素，以增加模拟结果的准确性。

管道内检测器在运行过程中，会产生与流体的作用，以及与管道内壁的摩擦作用。ABAQUS 提供了一种基于罚函数的一般接触算法，来描述欧拉流体与拉格朗日固体之间的相互作用。在分析过程中，拉格朗日固体网格与欧拉流体网格之间建立了联系，并在计算过程中自动计算和跟踪。欧拉网格必须要进行较为精细的划分，从而保证可以有效地捕捉耦合界面间的接触。

1.2 双节检测器过弯机理

检测器在运行过程中需要配备其他检测设备，多为多节的柔性结构，以便更好地通过不同形状的管道。而在检测器过弯过程中，万向节的长度L 以及允许的转角α6对检测器的运行状态起决定作用。双节舱段参数化分析简化模型如图1 所示。

1.3 有限元模型

管道内检测器的流-固耦合有限元模型如图2 所示。检测器模型与管道模型均采用线性减缩积分单元C3D8R 来描述，考虑到橡胶皮碗的超弹性和非线性行为，皮碗的网格划分密度需大于中心舱体，本研究设置皮碗的网格大小为5 mm，中心舱体网格大小为10 mm。欧拉流体材料采用8 节点线性欧拉实体单元EC3D8R 来描述，这种单元是基于C3D8R 拉格朗日单元扩展而来的，支持欧拉流动状态的分析。

在数值模拟过程中，欧拉材料的初始材料定义采用ABAQUS 中的The Volume Fraction Tool，同时结合材料的预定义场进行控制。

1.4 材料属性

皮碗材料采用应用广泛的聚氨酯[13]，聚氨酯具有强度高、回复性好、耐磨性能优良的特点。本研究采用的聚氨酯材料通过拉压试验获得性能数据，最后通过M-R 模型获得该材料的模型常数。

应变能密度函数可以表示为

式中： U——应变能；

Cij——常数项；

I1， I2——柯西应力的第一、二不变量；

Di——材料可压性，如Di=0，则表示材料不可压，皮碗材料性能参数见表1。

表1 皮碗材料性能参数

在CEL 方法中，管道内的流体材料采用Mie-Grüneisen 和Hugoniot 结合的状态方程（EOS）[13]来表达，该方程定义了流体材料的压力与密度之比，通用的流体Mie-Grüneisen 状态方程表示为

式中： pH——Hugoniot 压力的能量， Pa；

Em——Hugoniot 单位质量的能量， J；

ρ0——参考密度， kg/m3；

Γ0——材料常数；

η——名义体积压缩应变。

通常满足Hugoniot 曲线的压力方程表示为

式中： c0——流体声速， m/s；

s——待定常数；

Us——流体冲击速度， m/s；

Up——流体质点速度， m/s。

将公式（7）代入（6），可以表示为

其中， ρ0c20＞0 代表弹性体最小公称应变的弹性模量。

在模拟水的情况下， ABAQUS 规定s=0，Γ0=0， η=1-ρ0/ρ[13]。，可以得到一个简单的EOS 模型公式，即

综上分析，设定流体材料的各项参数分别为：密度ρ0=1 000 kg/m3，黏度μ=0.001 Pa·s，c0=1 483 m/s， s=0， Γ0=0。

1.5 边界条件

管道内检测器的有限元模拟边界条件主要包括管道的入口及出口边界[14-17]，设置入口边界为速度入口，同时设置出口为自由边界。限制管道空间所有自由度，并定义其为刚体。初始欧拉材料边界为流体充满管道，对于检测器占有的空间，欧拉流体材料为空，即EVF=0。同时设置管道内壁与皮碗的摩擦系数为0.4，设置全局摩擦系数为0。

2 管道内检测器过弯特性分析

基于非线性有限元软件ABAQUS 建立有限元模型，分析双节管道内检测器在管道运行过程中万向节的运动状态变化，进而分析不同的万向节长度对机器人通过性的影响及规律。

本研究采用的管道为152.4 mm （6 in）弯管，双直角拐弯，管道的曲率半径分别为1.5D、 3D （D 为管道的内径），在保证检测器单舱段各项几何尺寸不变的前提下，调整万向节的长度分别为135 mm、 142.5 mm、 150 mm、157.5 mm 和165 mm，通过有限元计算模拟分析，对比检测器在完整通过管道过程中万向节的前铰链、后铰链的角度变化以及前后舱段的轴线角度（锐角）变化情况。

双舱段仿真机段中以20 帧为单位提取舱段运行过程中万向节前、后铰接处弯曲角度，取前后舱段中轴线所成夹角绘制角度变化曲线。

万向节轴向长度为135 mm 时的角度变化曲线如图3 所示，由图3 可知，通过3D 大弯管时万向节前铰接处最大弯曲角度为13.2°，后铰接处最大弯曲角度为14.1°；通过1.5D 小弯管时万向节前铰接处最大弯曲角度为19.5°，后铰接处最大弯曲角度为28.6°。

万向节轴向长度为142.5 mm 时的角度变化曲线如图4 所示，由图4 可知，通过3D 大弯管时万向节前铰接处最大弯曲角度为11.2°，后铰接处最大弯曲角度为16.5°；通过1.5D 小弯管时万向节前铰接处最大弯曲角度为24°，后铰接处最大弯曲角度为31.8°。

万向节轴向长度为150 mm 时的角度变化曲线如图5 所示，由图5 可知，通过3D 大弯管时万向节前铰接处最大弯曲角度为10.4°，后铰接处最大弯曲角度为19.2°；通过1.5D 小弯管时万向节前铰接处最大弯曲角度为24°，后铰接处最大弯曲角度为31.3°。

万向节轴向长度为165 mm 时的角度变化曲线如图6 所示，由图6 可知，通过3D 大弯管时万向节前铰接处最大弯曲角度为13.2°，后铰接处最大弯曲角度为18.9°；通过1.5D 小弯管时万向节前铰接处最大弯曲角度为25.6°，后铰接处最大弯曲角度为33.8°。

由以上计算可知，随着万向节长度的增加，前铰链的弯曲角度变化峰值先减小后增大，而后铰链的角度变化峰值先增大后减小。

双节管道检测器在通过直角弯管时，由于万向节的转向作用，双节检测器可以顺利的通过拐弯位置。图7 所示为不同长度万向节条件下，前后舱段的轴线夹角变化曲线，由图7 可知，万向节夹角出现两次波峰和波谷，且在检测器通过1.5D 位置时的角度变化峰值明显大于通过3D 拐弯位置时的角度，且在小曲率半径位置时，万向节的角度变化更为剧烈，因此，检测器在通过小曲率管道位置时，由于角度剧烈的变化，会导致该位置产生剧烈的振动，因而该位置的检测精度会受到很大的影响。在通过3D 弯管位置时，随着万向节长度的增加，在万向节从L=135 mm 到L=165 mm 变化时，前后舱段的轴线夹角峰值分别为27°、 24°、 20°和30°，说明随着万向节长度的增加，舱段轴线的夹角峰值先减小后增大；而当通过1.5D 弯管时，当万向节从L=135 mm 到L=165 mm 变化时，前后舱段的轴线夹角峰值分别为50°、 58°、 56°、55°。由此可知，随着万向节长度的增加，舱段的轴线夹角呈现先增大后减小的趋势。在检测器经过管道的3D 曲率位置时，当L=142.5 时，角度波动出现二次峰值，随着万向节长度的增加；当L=165 mm 时，仅出现一次波动峰值，这可能是由于随着万向节的长度增加，前后舱段的皮碗均能稳定地贴合管道内壁，导致该位置的万向节角度变化平稳，波动减小。

为分析万向节在管道内运行过程中的应力变化规律，将万向节进行刚体化处理，进而获取在检测器过弯过程中的应力云图，以长度为135 mm的万向节为例提取检测器的应力云图如图8 所示。由图8 可以看出，在检测器运行时，万向节受到应力作用，实际中，万向节会受到拉力、压力以及扭转力的综合作用。当检测器通过曲率半径1.5D 弯道位置时，万向节所受应力最大，达到2 870 MPa，同时最大应力位置发生在铰接的配合位置，因此，在运行过程中，该位置会发生较大的应力集中，在设计万向节的结构时应重点考虑该位置的加强设计。