李云飞，陈振茂

（1.中国工程物理研究院总体工程研究所，四川 绵阳 621900；2.西安交通大学航天航空学院，陕西 西安 710049）

Ar瞬态方程非线性涡流检测信号数值模拟研究

李云飞1，陈振茂2

（1.中国工程物理研究院总体工程研究所，四川 绵阳 621900；2.西安交通大学航天航空学院，陕西 西安 710049）

针对一个简化的三维铁磁性材料数值计算模型，对其非线性涡流检测信号进行数值模拟研究。基于Ar瞬态方程的基本原理改良并开发用于计算非线性涡流检测信号的数值模拟程序，得出非线性涡流检测信号与磁性材料因塑性变形而产生的磁导率变化之间的相关性。对比非线性涡流数值模拟结果和实验检测结果表明：塑性变形导致的磁导率变化是非线性涡流检测信号变化的直接诱因。

非线性涡流；磁性材料；数值模拟；塑性变形

0 引 言

碳素钢等铁磁材料因其良好的机械性能和耐腐蚀性被用作核电站中主要结构材料。因地震等强载荷产生的塑性变形等损伤可能导致结构材料断裂和失效，进而引发核电事故。因此，结构材料塑性变形等机械损伤的定量无损检测和评价成为保证核电安全的关键技术。

国内对于塑性变形损伤的传统评价方法主要有基于金相分析、硬度压痕测试和X射线衍射法等，属于破坏性测试或者无法在现场使用。国外基于非线性涡流检测的研究：日本M.SHIWA[1]研究发现，非线性涡流信号频谱分析后的三次谐波分量幅值随铬钼钢焊后热处理温度发生显著变化；日本东北大学T.TAKAGI和T.UCHIMOTO[2-4]研究过非线性涡流信号与石墨铸铁中白口含量以及铸铁内石墨形态特征的相关性。

作者前期采用非线性涡流检测对碳素钢材料的塑性损伤评价做过一定研究[5]，而国内外关于非线性涡流检测的数值模拟鲜见报道。本文作为非线性涡流数值模拟的一个尝试，基于Ar瞬态方程和Crank-Nicholson直接积分法首次开发了非线性涡流数值模拟程序，并且对比数值模拟结果与实验检测结果。

1 Ar瞬态方程

磁准静态场（MQS）的麦克斯韦方程忽略位移电流密度项，其微分形式的基本方程组为

式中：H——磁场强度，A/m；

J——电流面密度，A/m2；

E——电场强度，V/m；

B——磁感应强度，T；

D——电通密度，C/m2；

ρ——电荷体密度，C/m3。

由于磁场B是无源场，所以有▽·B=0，引入一个矢量函数A（称为磁矢量）使得B=▽×A，推导可得：

对全体分析区域进行划分，并定义各个划分区域（见图1）。

图1 Ar法计算区域划分

全体区域：衰减区域+混合区域+正常区域。

衰减区域：仅由变形矢量磁位Ar记述的区域Vr。

混合区域：Ar和仅依赖电流源的As相混合的一层有限元区域Vtr。

正常区域：由通常矢量磁位A=（Ar+As）记述的区域Vt。

导体区域：包含在正常区域中的导体区域。

其中A=Ar+As，Ar称为变形磁矢位，As是仅仅由电流源所产生的磁矢位。由于在Ar法中仅依赖电流源的场As和仅依赖涡电流源的场Ar是独立考虑的，所以不划分电流源的网格。

利用六面体棱边元方法对Ar方法的控制方程式（9）～式（12）进行有限元离散可得伽辽金有限元方程：

式中：[N]——矢量形状函数；

n——单位法线矢量。

可简写为

2 数值模拟方法与程序开发

本文将非线性涡流数值模拟作为瞬态问题处理，基于Ar方法的低频涡流场控制方程式（14）可写成

式中：[K]，[C]，[M]——全局系数矩阵；

I（t）——依赖于时间函数的激励电流，本文采用正弦激励电流；

A——矢量磁位，Wb/m。

对于瞬态问题，矢量磁位A关于时间的导数∂A/∂t可由关于时间步的差分（Ak-Ak-1）/Δt近似计算，其中Ak=A（t0+kt），k表示第k次时间步，t表示时间步长，t0表示初始时间。为了提高计算稳定性，矢量磁位A可由Crank-Nicholson直接积分法替换。

其中0≤θ≤1是控制积分稳定性的参数。

将式（16）代入式（15），则当前步的矢量磁位A可通过下式计算：

对于常规涡流法的数值模拟通常将材料的磁导率视为定值来处理，这与铁磁材料实际磁化曲线不符。因控制方程的刚度矩阵与单元的磁导率相关，本文将所有导体单元的磁导率定义为一个数组，以每个单元中心的磁感应强度作为判定参数，满足一定条件则对该单元的磁导率进行重新赋值，然后进入下一个积分时间步。直至所有积分时间步计算完成，输出所有时间步检测线圈的电压检测值。

利用定长时间步的Crank-Nicholson直接积分法改良Ar法的数值模拟程序，在逐步积分法的每个时间步中导入碳素钢磁导率μ和磁感应强度B的非线性磁化曲线，通过非线性磁化曲线对磁导率进行重新赋值，而不同于常规涡流法中将磁导率视为一个定值来处理。图2为碳素钢μ-B近似磁化关系曲线[6]。

图2 碳素钢μ-B关系曲线

将μ-B关系曲线近似分为若干个区间，通过式（17）计算出矢量磁位A，从而可以计算出磁感应强度B。如果某一积分时间步计算所得的单元中心磁感应强度B满足Bi＜B＜Bi+1，则通过下式迭代方法对该单元磁导率重新赋值：

重新赋值后的磁导率进行单元刚度阵组装并进行下一次的积分时间步的信号计算，最后输出所有时间步的检测信号。由此，利用基于Ar方法并改良后的程序可计算得出非线性涡流模拟信号。

3 模型算例与结果分析

利用定长时间步的Crank-Nicholson直接积分法改良Ar法数值模拟程序，在逐步积分法的时间步中导入碳素钢的非线性磁化曲线，可获取非线性涡流的数值模拟信号。图3为非线性涡流法的数值计算模型，为了提高计算速度，采用有限元和棱边元方法对整个分析区域进行离散，将导体外电流源的作用等效在导体的边界上，因此线圈是不需要网格划分的，只需要对导体区域和空气区域进行网格划分。算例中激励线圈和检出线圈直径为15mm，高为3mm，匝数为500匝。

根据导入的非线性磁化曲线，采用改良后的Ar法数值模拟程序计算得到在任意激励频率下的非线性涡流检测信号。在激励线圈中施加一正弦电流，图4为检出线圈中得到的检测信号与正弦激励信号对比图，从图中可知检测信号不完全是一正弦信号，也出现了倍频信号。

图3 非线性涡流数值计算模型

图4 检测线圈的非线性涡流模拟信号

因碳素钢塑性变形后会产生反磁致伸缩效应，并且导致材料的磁各向异性，使磁畴无法自由移动，从而影响材料的电磁特性[7-10]。根据该原理，本文在非线性涡流检测信号的数值模拟中通过改变材料不同的初始磁导率，等效模拟材料不同程度塑性变形的情况。对所得到的检测信号进行傅里叶变换的频谱分析，得到信号的功率频谱图如图5所示，图中可见基波和高次谐波的波峰。

在数值模拟中任意选取10kHz作为激励频率，计算得到不同初始磁导率情况下材料的非线性涡流信号，提取检测信号频谱图中基波、三次谐波对应的幅值参数，可得到不同初始磁导率（即不同塑性变形）下的关系结果，如图6所示。

图5 非线性涡流模拟信号频谱图

图6 幅值与塑性变性关系

可发现非线性涡流模拟信号频谱图中基波幅值和三次谐波随材料的初始磁导率即塑性变形程度的增大出现下降的趋势，这一趋势与此前实验研究中碳素钢的实验检测结果是一致的[5]。

4 结束语

本文基于Ar方法，采用Crank-Nicholson直接积分法改良开发了非线性涡流数值模拟程序。对简化的磁性材料三维模型，以不同初始磁导率等效材料不同塑性损伤，计算得到不同损伤程度下非线性涡流信号。提取信号中基波幅值和三次谐波幅值，得出与塑性损伤之间关系曲线。通过对比，该关系曲线与此前实验研究中碳素钢的实验检测结果是一致的，证明了塑性变形导致的磁导率变化是非线性涡流检测信号变化的直接诱因，也验证了非线性涡流检测对碳素钢等磁性材料塑性变形无损定量评价的可行性。

[1]Shiwa M，Chen G Z，Horii Y，et al.Evaluation of PWHT temperature for Cr-Mo steel welded joint by AC magnetic signal analysis[J].Electromagnetic Nondestructive Evaluation，2006（7）：295-304.

[2]Uchimoto T，Matsukawa J，Abe T，et al.Evaluation of chill contents in flake graphite cast irons using AC magnetization method[J].Electromagnetic Nondestructive Evaluation，2008（11）：62-69.

[3]Uchimoto T，Takagi T.Characterization of matrices and graphite forms of cast irons by electromagnetic nondestructive evaluation[J].Electromagnetic Nondestructive E-valuation，2009（12）：207-214.

[4]Stupakov O，Uchimoto T，Takagi T，et al.Evaluation of ductile cast iron microstructure by magnetic hysteresis and barkhausen noise methods[J].Electromagnetic Nondestructive Evaluation，2009（12）：232-239.

[5]李云飞，陈振茂.基于非线性涡流检测的核电结构材料塑性变形损伤无损评价研究[C]∥第22届全国结构工程学术会议论文集，2013：26-29.

[6]何松，张辉，倪泽钧.常用钢材磁特性曲线速查手册[M].北京：机械工业出版社，2003：44-50.

[7]于庆波，李子林，韦玄，等.低碳合金钢变形奥氏体的相变研究[J].材料热处理学报，2007，28（4）：83-87.

[8]李维娟，李胜利，王方.奥氏体变形条件下低碳微合金钢中的贝氏体相变[J].金属热处理，2005，30（6）：70-73.

[9]李明，任爱，薛飞.铁磁性材料热交换管的远场涡流检测探讨[J].无损检测，2005，27（1）：6-9.

[10]周理平.涡流无损评价的理论与实验方法[J].无损检测，1996，18（7）：186-188.

Numerical simulation method of nonlinear eddy current testing signals based on transient Ar formulation

LI Yun-fei1，CHEN Zhen-mao2
（1.Institute of Systems Engineering，CAEP，Mianyang 621900，China；2.School of Aerospace，Xi’an Jiaotong University，Xi’an 710049，China）

For a simplified three-dimensional numerical model of magnetic materials，numerical simulation research of nonlineareddy currentdetection signal of the model wascarried out. Numerical simulation program used to calculate nonlinear eddy current detection signal was improved and developed based on Ar method principle.Furthermore，relationship between nonlinear eddy currentdetection signaland permeability variation caused by plastic deformation ofmagnetic materials was obtained.It is proved that the permeability variation caused by plastic deformation be the direct incentive of nonlinear eddy current detection signal’s change through comparing the results of numerical simulation with experimental testing results.

nonlinear eddy current；magnetic material；numerical simulation；plastic deformation

TG115.28；O482.54；TN911.7；TM930.114

：A

：1674-5124（2014)06-0006-04

10.11857/j.issn.1674-5124.2014.06.002

2014-01-22；

：2014-03-16

国家自然科学基金项目（50977070，11372295，11302206）

李云飞（1986-），男，广西桂林市人，工程师，硕士，研究方向为电磁无损检测与评价。