一种计算管道结构超声导波频散特性的方法

技术领域

本发明涉及超声导波检测技术，尤其涉及一种计算管道结构超声导波频散特性的方法，属于无损检测技术领域。

背景技术

管道设备作为重要的运输工具，在输送流体和气体中起到重要作用。如今，随着油气使用的增加，管道运输领域也得到了更大的发展。而在长期服役过程中，由于复杂的环境因素，管道可能会出现损伤或者裂纹，从而引起气液泄露的风险。选择有效的检测技术对管道进行安全检测至关重要。作为无损检测技术之一的超声导波检测具有远距离传播、传播速度快以及传播范围广等优点，可适用于圆管类结构的安全检测。

管中的超声导波可以分为轴向导波和周向导波两类模式。轴向导波主要应用于管的长距离检测，多用于检测沿管道周向分布的缺陷和裂纹，对管道中纵向分布的缺陷并不敏感。而周向导波是沿着圆管的圆周方向传播的导波，对管道纵向缺陷较为敏感，应用周向导波检测管道中纵向缺陷可以有很好的效果。在特定的频率下，超声导波可以沿着圆管长度方向或圆周方向传播，且接收到的信号完整，并且在遇到缺陷或结构不连续时，会在回波信号中反映出来，所以可以据此有效判断结构的损伤信息。而选择一个合适的激励频率对管道进行超声导波检测，则需要研究管道的频散特性。

常见的管道超声导波频散特性研究方法中，传统的解析方法适合求解各向同性结构，对于复杂的多层复合结构，解析方法难以获得其解析解，并且会面临丢根和计算不稳定的问题。而有限元法则需要建立管道的三维模型，求解的自由度大，计算成本高。半解析有限元法结合了解析法和有限元法的优点，广泛应用于导波频散特性的研究，但大多数的半解析有限元法通过对管道的二维模型进行求解，虽然较大地降低了求解的自由度，但当求解高频厚积下管道的频散曲线时，仍需要大量的运算时间。若能将管道结构简化为一维模型，求解的速度将会极大地提升。

发明内容

有鉴于此，本发明提供了一种计算管道结构超声导波频散特性的方法，将管道简化为径向的一维模型，并且用统一推导的波动控制方程弱形式求解管道轴向和周向导波的频散特性。

一种计算管道结构超声导波频散特性的方法，包括以下步骤：

1.对声波在圆管中不同方向的质点振动，采用统一的位移解表示沿管道轴向和周向的波传播，推导柱坐标系下各向异性管道轴向和周向导波传播的波动控制方程弱形式的统一表达式：

1.1.一般各向异性材料中的波动控制方程为：

其中，σ

1.2.将表达式(1)波动方程的两边乘以测试函数v

1.3.对表达式(2)应用散度定理，方程的弱形式可以定义如下：

1.4.表达式(3)中，在边界Γ上的积分趋向于0，所以方程的弱形式可以表达为：

1.5.在各向异性弹性介质中，应力与应变的关系

σ

其中，C

1.6.在柱坐标系下以位移场的正交分量表示的应变表达式：

1.7.根据导波在圆管中不同方向的质点振动：纵波的质点运动主要在r和(或)

z轴方向；扭转波质点运动主要在θ轴方向；弯曲波在周向、径向和轴向三个方向均具有质点振动；而周向导波为沿着管道周向传播的导波。采用统一的位移解表示沿管道轴向和周向的波传播：

上式中，m为周向阶次；K

1.8.把质点的位移表达式代入波动控制方程的弱形式中，三维域的求解问题简化为一维域求解，结合表达式(5)、(6)，各向异性材料的管道波动控制方程的弱形式展开式为：

2.利用管道的对称性，构建管道半径方向的一维有限元模型。

3.利用弱形式偏微分方程模块，将表达式(8)中的被积函数转换成有限元软件可以读取的弱形式被积函数表达式。

4.在参数栏中定义弱形式表达式中的参数：求解L(0,n)纵向模态和T(0,n)扭转模态时，令K

5.对模型网格离散化后，添加参数化扫描，对波数或角频率进行扫描以搜索特征值，并设定扫描区间以及所需的特征值数。

6.特征值求解得到波数和角频率，根据相速度公式V

其中，Δf为同一模态下相邻频率点之差；Δk为相邻波数之差。

7.在有限元软件中，添加一维绘图组分别绘制不同方向的波结构图。

本发明具有以下优点：

1)本方法基于半解析有限元理论，只需建立管道的一维有限元模型，运算量极少，求解速度快。

2)本方法推导了统一表示管道轴向和周向导波传播的波动控制方程的弱形式，只需简单设置弱形式方程中的参数，便可分别求解管道中不同导波类型的频散特性，并且利用有限元软件求解，计算出来的管道导波模态完整，无模态缺失。

附图说明

图1为基于弱形式偏微分方程理论计算管道轴向超声导波频散关系的步骤框图；

图2为无限长空心圆管中超声导波沿波导结构轴向传播的三维模型示意图；

图3为超声导波沿管道结构周向传播的三维模型示意图；

图4为有限元软件中管道的一维模型选取示意图；

图5为利用本方法求解得到的管道纵向模态和扭转模态的相速度频散曲线与理论频散曲线对比示意图；

图6为利用本方法求解得到的管道纵向模态和扭转模态的群速度频散曲线与理论频散曲线对比示意图；

图7为利用本方法求解得到的管道弯曲模态的相速度频散曲线与理论频散曲线对比示意图；

图8为利用本方法求解得到的管道弯曲模态的群速度频散曲线与理论频散曲线对比示意图；

图9为利用本方法求解得到的管道周向导波各模态的相速度频散曲线与理论频散曲线对比示意图；

图10为利用本方法求解得到的管道周向导波各模态的群速度频散曲线与理论频散曲线对比示意图；

图11为本方法求解得到的纵向模态L(0,1)在0.492MHz频率下的径向位移波结构图；

图12为本方法求解得到的纵向模态L(0,1)在0.492MHz频率下的轴向位移波结构图；

图13为利用本方法求解得到的周向类兰姆波CLT0模态在0.514MHz频率下的径向位移波结构图；

图14为利用本方法求解得到的周向SH波CSH0模态在0.710MHz频率下的轴向位移波结构图。

具体实施方式

实施例1

结合本发明方法的内容提供以下管道轴向超声导波频散关系计算方法实例，并与理论计算的结果对比验证本方法的正确性，具体步骤图1所示：

1)对导波在圆管中不同方向的质点振动，采用统一的位移解表示沿管道轴向和周向的波传播，推导得到柱坐标系下各向异性管中轴向导波的波动控制方程的弱形式展开式：

2)利用商用有限元软件COMSOL Multiphysics 6.0中的弱形式偏微分方程模块进行分析，选择特征值研究。构建管道半径方向的一维模型，管道模型的内半径为12mm，外半径为15mm。三维管道模型及波传播示意图如图2所示，一维有限元模型选取示意图如图4所示。

3)在参数栏中定义各向异性材料弹性刚度常数、管道的径向尺寸、波数、角频率、系数K

4)在弱形式偏微分方程模块的其中一个输入栏中输入推导的弱表达式，另外两个设置为0。将角频率定义为特征值，波动控制方程弱形式中的被积函数转换为COMSOL可以读取的弱表达式：

2*pi*r*((u3r+Kz*k*u1*i)*(C25*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C45*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C56*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C15*conj(test(u1r))+C55*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C35*Kz*conj(test(u3))*conj(k)*i)+(u2r-u2/r+

(m*u1*i)/r)*(C26*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C46*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C66*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C16*conj(test(u1r))+C56*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C36*Kz*conj(test(u3))*conj(k)*i)+u1r*(C12*(conj(test(u1))/r-

(m*conj(test(u2))*i)/r)-C14*(Kz*conj(test(u2))*conj(k)*i+

(m*conj(test(u3))*i)/r)-C16*(conj(test(u2))/r-conj(test(u2r))+

(m*conj(test(u1))*i)/r)+C11*conj(test(u1r))+C15*(conj(test(u3r))-

Kz*conj(test(u1))*conj(k)*i)-C13*Kz*conj(test(u3))*conj(k)*i)+

((m*u3*i)/r+Kz*k*u2*i)*(C24*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C44*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C46*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C14*conj(test(u1r))+C45*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C34*Kz*conj(test(u3))*conj(k)*i)+(u1/r+

(m*u2*i)/r)*(C22*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C24*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C26*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C12*conj(test(u1r))+C25*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C23*Kz*conj(test(u3))*conj(k)*i)-lambda^2*rho*conj(test(u1))*u1-

lambda^2*rho*conj(test(u2))*u2-lambda^2*rho*conj(test(u3))*u3+

Kz*k*u3*(C23*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C34*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C36*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C13*conj(test(u1r))+C35*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C33*Kz*conj(test(u3))*conj(k)*i)*i)

其中，conj为复共轭函数；m为周向阶次；lambda为所求的特征值，波数参数化扫描时，特征值即角频率；test()函数表示测试函数。

5)对有限元模型进行网格划分，大小选择“较细化”。添加参数化扫描，对波数k进行扫描，扫描范围为[0，5*pi/a]，其中a为管道厚度。特征值求解中，设置所需特征值数为20。

6)在COMSOL软件的“结果”中，得到波数和角频率，其中角频率ω＝2π/f，利用相速度计算公式V

7)添加有限元软件中的一维绘图组分别绘制不同方向的波结构图。归一化处理后的波结构图如图11、12所示。

实施例2

结合本发明方法的内容提供以下管道结构周向导波频散关系计算方法实例，并与理论计算的结果对比验证本方法的正确性，具体步骤如图1所示：1)对声波在圆管中不同方向的质点振动，采用统一的位移解表示沿管道轴向和周向的波传播，推导得到柱坐标系下各向异性管中周向导波的波动控制方程的弱形式展开式：

2)选取有限元分析软件COMSOL Multiphysics 6.0中的弱形式偏微分方程模块进行分析，选择特征值研究。构建管道半径方向的一维模型，管道模型内半径为12mm，外半径为15mm。三维管道模型及波传播示意图如图3所示，一维有限元模型选取示意图如图4所示。

3)在参数栏中定义各向异性材料弹性刚度常数、管道的径向尺寸、波数、角频率、系数K

杨氏模量为2.09GPa，泊松比为0.29，密度为7824kg/m

阵为：

4)在弱形式偏微分方程模块的其中一个输入栏中输入推导的弱表达式，另外两个设置为0。将角频率定义为特征值，将管道波动控制方程的弱形

式中被积函数转换为COMSOL可以读取的弱表达式：

2*pi*r*((u3r+Kz*k*u1*i)*(C25*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C45*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C56*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C15*conj(test(u1r))+C55*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C35*Kz*conj(test(u3))*conj(k)*i)+(u2r-u2/r+

(m*u1*i)/r)*(C26*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C46*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C66*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C16*conj(test(u1r))+C56*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C36*Kz*conj(test(u3))*conj(k)*i)+u1r*(C12*(conj(test(u1))/r-

(m*conj(test(u2))*i)/r)-C14*(Kz*conj(test(u2))*conj(k)*i+

(m*conj(test(u3))*i)/r)-C16*(conj(test(u2))/r-conj(test(u2r))+

(m*conj(test(u1))*i)/r)+C11*conj(test(u1r))+C15*(conj(test(u3r))-

Kz*conj(test(u1))*conj(k)*i)-C13*Kz*conj(test(u3))*conj(k)*i)+

((m*u3*i)/r+Kz*k*u2*i)*(C24*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C44*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C46*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C14*conj(test(u1r))+C45*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C34*Kz*conj(test(u3))*conj(k)*i)+(u1/r+

(m*u2*i)/r)*(C22*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C24*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C26*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C12*conj(test(u1r))+C25*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C23*Kz*conj(test(u3))*conj(k)*i)-lambda^2*rho*conj(test(u1))*u1-

lambda^2*rho*conj(test(u2))*u2-lambda^2*rho*conj(test(u3))*u3+

Kz*k*u3*(C23*(conj(test(u1))/r-(m*conj(test(u2))*i)/r)-

C34*(Kz*conj(test(u2))*conj(k)*i+(m*conj(test(u3))*i)/r)-

C36*(conj(test(u2))/r-conj(test(u2r))+(m*conj(test(u1))*i)/r)+

C13*conj(test(u1r))+C35*(conj(test(u3r))-Kz*conj(test(u1))*conj(k)*i)-

C33*Kz*conj(test(u3))*conj(k)*i)*i)

其中，conj为复共轭函数；lambda为所求的特征值，波数参数化扫描时，特征值即角频率；m为角波数，m＝kR1；k为环向波数；test()函数表示测试函数。

5)对有限元模型进行网格划分，大小选择“较细化”。添加参数化扫描，对波数k进行扫描，扫描范围为[0，5*pi/a]，其中a为管道厚度。特征值求解中，设置所需特征值数为20。

6)在COMSOL软件的“结果”中添加全局计算，得到波数和角频率，其中ω＝2π/f，利用相速度计算公式V

7)添加有限元软件中的一维绘图组分别绘制不同方向的波结构图。归一化处理后的波结构图如图13、14所示。

本发明公开了一种计算管道结构超声导波频散特性的方法，包括：对导波在圆管中不同方向的质点振动，采用统一的位移解表示沿管道轴向和周向的波传播，推导柱坐标系下各向异性管中轴向和周向导波的波动控制方程的弱形式；利用管道的对称性，建立管道半径方向的一维有限元模型；定义有限元软件中的导波波动方程弱形式表达式，通过对弱形式表达式中的参数进行设置，可分别计算沿圆管不同方向传播的导波频散特性；利用有限元软件的特征值求解，得到波数和频率，分别计算管道轴向导波和周向导波各模态的相速度和群速度，绘制频散曲线和波结构图,本发明只需建立管道的一维有限元模型，并通过简单的参数设置即可快速地求解管道中不同方向传播的超声导波频散关系，无需繁琐编程。

最后应说明的是：本发明不限于上述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整和替代而不会脱离本发明的保护范围。因此，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但是不仅仅限于以上实施例，在不脱离本发明构思的情况下，还可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。