一种含柱状缺陷的涡流无损检测磁场的半解析计算方法

摘要

本发明公开了一种含柱状缺陷的涡流无损检测磁场的半解析计算方法，利用数学物理和电磁场原理，得出拉氏域下磁矢量势的精确解析表达式，并基于MATLAB对磁矢量势级数形式的复杂表达式进行了数值法拉普拉斯逆变换，得到含柱状缺陷导体中涡流无损检测信号的时域解。本发明可准确高效地计算含柱状缺陷导体涡流无损检测磁场，为涡流无损检测技术在缺陷的定量识别领域提供了理论基础。

说明书

技术领域

本发明属于金属试件缺陷检测技术领域，具体涉及一种含柱状缺陷的涡流无损检测磁场的半解析计算方法的设计。

背景技术

在金属试件的缺陷检测领域，涡流无损检测方法相对于其他检测方法具有检测精度高、造价低廉、无需耦合剂、对人体和环境安全绿色等诸多优势。而探究电磁场信号同激励源信号和缺陷信息之间关系，即涡流无损检测的正向问题的研究是涡流缺陷识别的基础，对实际工程中检测仪器开发具有理论指导意义。

目前涡流缺陷检测正向问题研究有三种方法：解析法、数值仿真法和实验法。在缺陷识别研究中，目前大多采用数值仿真法和实验法，其中实验法对环境和仪器要求高且只能得出空间某一点的电磁场信号；数值仿真法对计算机配置要求高且运算时间长。而解析或半解析的方法具有结果精确、运算耗时少等特点。

解析法根据麦克斯韦方程组，建立以磁矢量势为求解目标的偏微分方程，转化至拉氏域后添加磁场边界条件，求解出磁矢量势的函数解析表达式，再进行拉普拉斯逆变换可以得出时域下的电磁场信号。然而，由于缺陷的存在，求解模型的边界条件变得复杂，导致偏微分方程求解困难；而且拉氏域的磁矢量势函数表达式形式复杂，进行拉普拉斯逆变换比较困难。因此如何求解偏微分方程和如何对磁矢量势表达式进行拉普拉斯逆变换是需要解决的核心问题。

针对被测体中含缺陷的涡流检测问题，现有方法采用数值仿真技术，根据检测激励信号和模型材料及其几何尺寸等求电磁场磁矢量势。其技术方案是采用有限元法将原模型求解域划分成为有限个单元，单元之间以节点为连接点，假设单元间的相互作用只能通过节点传递。这样就将原来具有无限自由度的连续变量微分方程和边界条件转化为只包含有限个节点变量的代数方程组，以利用ANSYS、COMSOL等商业数值计算软件求解出数值解。

现有的利用商业数值计算软件对涡流缺陷检测中电磁场物理量计算主要存在以下几个缺点：

(1)有限元仿真获得的数值解无法得知模型材料、缺陷尺寸同电磁场信号的显式联系；

(2)有限元仿真过程中对模型划分网格时，对模型做了离散化处理，导致求解结果存在计算误差；

(3)利用有限元仿真软件进行仿真时，由于计算量很大，对计算机配置要求很高，而且计算时间很长，计算效率低。

发明内容

本发明的目的是为了解决现有技术中利用商业数值计算软件对涡流缺陷检测中电磁场物理量计算存在的上述问题，提出了一种含柱状缺陷的涡流无损检测磁场的半解析计算方法。

本发明的技术方案为：一种含柱状缺陷的涡流无损检测磁场的半解析计算方法，包括以下步骤：

S1、按照材料以及缺陷的分布，对缺陷检测模型划分求解区域；

S2、根据电磁学原理对求解区域建立偏微分方程，并运用数学物理方法求解拉氏域偏微分方程通解的形式；

S3、对模型求解区域添加内边界条件和外边界条件，利用贝塞尔函数正交性建立线性方程组，从而求出偏微分方程通解的系数，得出拉氏域模型磁矢量势级数形式的解析解；

S4、编写拉普拉斯逆变换的运算程序，并将拉氏域模型磁矢量势的解析解录入运算程序，封装为MATLAB的m文件；

S5、给出m文件的输入值，并在MATLAB中运行运算主函数，得出给定输入值对应的输出结果。

本发明的有益效果是：

(1)本发明可以获得拉氏域含柱状缺陷磁矢量势的级数形式解析解。

(2)本发明大幅提高了含柱状缺陷模型电磁场计算的求解效率，据统计相同计算条件下，计算效率为数值计算效率的5倍。

(3)本发明可为缺陷涡流无损检测提供理论、技术及数据等支持。

附图说明

图1为本发明提供的一种含柱状缺陷的涡流无损检测磁场的半解析计算方法流程图。

图2为本发明实施例的缺陷检测模型求解区域划分示意图。

图3为本发明实施例的拉普拉斯逆变换子函数计算流程图。

图4为本发明实施例的主程序函数计算流程图。

具体实施方式

下面结合附图对本发明的实施例作进一步的说明。

本发明提供了一种含柱状缺陷的涡流无损检测磁场的半解析计算方法，如图1所示，包括以下步骤：

S1、按照材料以及缺陷的分布，对缺陷检测模型划分求解区域。

对于电磁场问题，求得磁矢量势后理论上可以推导出电磁场所有物理量，因此电磁场计算问题往往以磁矢量势作为最初求解对象。本发明实施例针对图2(a)的缺陷检测模型对其检测物理量进行计算，首先以磁矢量势作为求解目标量。为求解需要，先考虑图2(b)只有单匝检测线圈的情况，并根据材料以及缺陷的分布，将缺陷检测模型划分为5个求解区域，其中无缺陷区域包括1区、2区、4区和5区，1区、2区和5区为空气域，4区为导体域；缺陷区域为3区，3区既包含导体域又包含空气域，如图2所示。

S2、根据电磁学原理对求解区域建立偏微分方程，并运用数学物理方法求解拉氏域偏微分方程通解的形式。

本发明实施例中，由于缺陷检测模型具有轴对称性，该模型的磁矢量势只含有周向分量，根据电磁场理论和麦克斯韦方程组，各求解区域的磁矢量势A满足式(1)所示的偏微分方程：

式中ρ表示径向坐标，z表示轴向坐标，p表示一个参数变量，对于空气域(1区、2区和5区)，p²＝0；对于导体域(3区和4区)，p²＝-jωμ₀μ_rσ；为求解方便将式(1)转为拉氏域下，此时导体域中p²＝-sμ₀μ_rσ。其中j为复数标记，ω表示角频率，μ₀表示真空磁导率，μ_r表示导体相对磁导率，σ表示导体电导率，s表示拉氏域下表达式的频率参数。

对无缺陷区域(1区、2区、4区和5区)的磁矢量势进行有限汉克尔变换得到：

式中表示经过有限汉克尔变换后的磁矢量势表达式，表示频域下的磁矢量势表达式，J₁表示第一类一阶贝塞尔函数，λ_i表示满足边界条件的正常数，b表示求解区半径。

将式(2)代入式(1)后，可将偏微分方程转化为常微分方程求解问题，即求解式(3)所示的常微分方程：

求解此常微分方程，并令缺陷检测模型上无穷远和下无穷远处的边界函数值为0，即当z＝+∞和z＝-∞时磁矢势为0，可得式(4)-(7)：

式(4)-(7)中表示第i区经过有限汉克尔变换后的磁矢量势表达式，i＝1,2,4,5；e表示自然常数，γ_i满足等式C₁₁、C₂₁、C₂₂、C₄₁、C₄₂、C₅₁均为待定系数。

对于含缺陷的3区，采用分离变量法求解偏微分方程，得出：

式中表示第3区频域下的磁矢量势表达式；q_i表示满足边界条件的正常数；p_i满足等式B₃、C₃、E、F为待定系数，J₁(q_iρ)表示第一类一阶贝塞尔函数，Y₁(q_iρ)表示第二类一阶贝塞尔函数，b表示求解区半径，c表示缺陷半径。

S3、对模型求解区域添加内边界条件和外边界条件，利用贝塞尔函数正交性建立线性方程组，从而求出偏微分方程通解的系数，得出拉氏域模型磁矢量势级数形式的解析解。

模型边界条件分为两类，模型的外边界条件和各求解区域间的内边界条件。其中外边界条件只包含狄里克莱边界条件，内边界条件既包含狄里克莱边界条件又包含诺依曼边界条件。为便于求解依次添加各个求解区域的边界条件，因此该步骤具体包括以下分步骤：

S31、对3区中ρ＝c处添加内边界条件。

3区中ρ＝c处几何形状为一个圆柱的侧面，首先设置磁感应强度ρ分量B_ρ在该处连续，即满足：

EJ₁(q_iρ)|_ρ＝c＝J₁(q_iρ)+FY₁(q_iρ)|_ρ＝c(9)

其次令磁场强度z分量H_z在该处连续，即满足：

通过式(9)可以将式(8)变换形式为：

将式(11)带入式(10)，可解得式(11)中：

式中J₀表示第一类零阶贝塞尔函数，Y₀表示第二类零阶贝塞尔函数。

S32、对1区和2区的交界处以及4区和5区的交界处添加内边界条件。

1区和2区的交界(z＝d)处以及4区和5区的交界(z＝-h)处几何形状为水平平面，其内边界条件表示为如式(12)-(13)的形式：

式(12)-(13)中下标i和j代表不同求解区域，和可写为分别代表1、2、4、5区频域下经过有限汉克尔变换后的磁矢量势表达式；μ_i和μ_j可写为μ₁、μ₂、μ₄、μ₅，分别代表1、2、4、5区的磁导率；可写为对于1区和2区的交界处有表示频域下的激励电流表达式，r表示激励线圈半径；对于4区和5区的交界处有通过改变本发明可适用于多种激励方式。

S33、对2区和3区的交界处以及3区和4区的交界处添加内边界条件。

2区和3区的交界(z＝0)处以及3区和4区的交界(z＝-l)处的内边界需满足磁矢势的连续和磁感应强度ρ分量连续，即：

式(14)-(15)中下标i和j代表不同求解区域，和可写为分别代表第2、3、4区频域下的磁矢量势表达式。

对于2区和4区，带入式(14)和式(15)之前需要进行汉克尔逆变换：

S34、添加模型求解区域的外边界条件。

外边界条件要求模型求解区域外部边界磁矢量势为0，模型上无穷远处和下无穷远处在求解偏微分方程时已满足该外边界条件。

为了获取级数形式的解，将模型侧面外边界处的磁矢量势设置为0，即：

带入各求解区域，即满足式(18)-式(19)的所有正根：

J₁(bλ_i)＝0(18)

J₁(bq_i)[J₁(q_ic)Y₀(q_ic)-μ_rY₁(q_ic)J₀(q_ic)]+(μ_r-1)J₁(q_ic)J₀(q_ic)Y₁(bq_i)＝0(19)

S35、求解待定系数。

将式(4)-(7)及式(11)带入式(12)-(16)的内边界条件，利用贝塞尔函数的正交性和贝塞尔-傅里叶变换，消去含ρ变量，可得一个含有8个方程的线性方程组，将其整理为矩阵形式，可得：

Kx＝w(20)

式中K为一个已知矩阵，x为待定系数组成的未知列向量，w为已知列向量；

由矩阵的形式可以得出，矩阵为一个满秩矩阵，因此式(20)有且仅有一组解，式(20)中：

x＝[C₁₁>21>22>3>3>41>42>51]^T

其中d、l、h表示求解分区的轴向边界坐标。

求解式(20)的线性方程组后，可得出各求解域拉氏域的磁矢量势级数形式的解析表达式，需要注意的是，此结果为图2(b)所示单匝线圈激励模型的结果，为获得图2(a)所示圆柱形线圈的结果，还需进行式(21)的积分运算：

式中表示圆柱形激励线圈激励下的磁矢量势频域表达式，s表示圆柱形激励线圈的横截面，N_d表示线圈绕线密度。

S4、编写拉普拉斯逆变换的运算程序，并将拉氏域模型磁矢量势的解析解录入运算程序，封装为MATLAB的m文件。

为获得时域下的结果，还需要进行拉普拉斯逆变换。由于得出的拉氏域的磁矢量势表达形式复杂，解析方法无法对其进行拉普拉斯逆变换，因此本发明实施例中采用基于MATLAB进行拉普拉斯逆变换的数值算法。

为了便于计算，本发明将一个复杂函数式进行拉普拉斯逆变换数值计算的算法程序封装为一个子函数，嵌套在主函数内进行计算，该子函数首先将连续的时间段转化为离散的时间轴向量，再采用循环结构，每次循环将原复杂表达式转化为含符号s的离散向量，由于离散向量可利用数值积分的方法得出数值解，计算流程如图3所示。

S5、给出m文件的输入值，并在MATLAB中运行运算主函数，得出给定输入值对应的输出结果。

由于解的形式是级数形式，为了保证计算精度，使结果累加m次，主函数计算流程图如图4所示。

图4中输入数据为拉普拉斯逆变换后时域信号起始时间、终止时间、计算点数及检测位置信息；gen_xu_lie子函数可以获取式(18)的正跟序列并使结果输入cha_xu_lie子函数可获得计算所需的系数向量进行计算前数据初始化；式(19)的正跟序列通过图解法获取进行计算前数据初始化；G矩阵为m×n的矩阵，n为用户输入的计算点数，m为累加的项数，计算时将级数解每一项进行拉普拉斯逆变换后的数值结果存储在G矩阵的每一行，最终需对G矩阵每一行乘对应系数相加，本发明使G矩阵左乘数据初始化阶段得到的系数列向量完成此运算，最终得到向量形式的数值解。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。