ANSYS在瞬变电磁三维正演中的应用

智庆全，武军杰，王兴春，孙怀凤， 杨毅，张杰，邓晓红，陈晓东,杜利明

(1.中国地质科学院 地球物理地球化学勘查研究所，河北 廊坊 065000; 2.自然资源部地球物理电磁法探测技术重点实验室，河北 廊坊 065000; 3.山东大学 岩土与结构工程研究中心，山东 济南 250061; 4.中国冶金地质总局 山东正元地质勘查院，山东 济南 250101)

0 引言

瞬变电磁法是应用最广泛的电磁探测方法之一，近年来在金属矿勘查、煤田地质勘探、隧道超前预报等领域[1-5]发挥了重要作用。瞬变电磁法一般采用不接地回线或接地导线作为发射源，供以稳恒电流，观测电流关断后地下地质体感生的二次电磁场，通过分析二次场的特征来推断地下地质体的电性、位置、深度等信息。利用瞬变电磁法进行勘探工作的关键技术在于根据瞬变电磁响应规律推断地质体信息，而瞬变电磁正演是研究瞬变电磁响应规律的有效途径。随着计算机技术的应用，瞬变电磁正演方法有了快速发展[6-9]，其中，一维层状模型的瞬变电磁正演通常基于频谱法，采用线性数字滤波法进行快速计算，二维、2.5维的瞬变电磁正演主要采用有限单元法和有限差分法，三维瞬变电磁正演方法主要为有限单元法、积分方程法、时域有限差分法等。目前，一维瞬变电磁正演方法经过多年发展已渐趋成熟，基于一维正演的一维反演技术是当前进行瞬变电磁反演的主要方法，而二维、2.5维瞬变电磁方法由于其二维近似的条件常常难以满足，在实际应用中存在诸多限制；三维瞬变电磁正演方法相对一维、二维方法能更精确地模拟实际地质情况，基于三维正演的解释方法能取得更高的解释精度，但一般需要花费很高的时间成本。

瞬变电磁的高精度、快速三维正演技术已成为当前瞬变电磁正演研究的热点，国内外很多学者都针对这一问题展开了研究。目前，瞬变电磁的三维正演技术主要可分为频时转换法和直接时域法两大类。频时转换法首先计算三维频率域电磁响应，然后采用余弦变换、G-S变换等频时技术将频率域响应转换到时间域，进而获得瞬变电磁响应。时域法一般通过时间离散来直接求取时域响应，主要包括时域有限差分法、有限体积法、有限单元法等。时间域有限差分法基于交错网格和显式时间差分格式[10-12]，该方法具有数学概念直观、表达简单的特点，但是需要满足稳定性条件，其网格尺寸和时间步长均受到严格限制。时间域有限体积法和时间域有限元法通过采用后退欧拉时间离散格式对控制方程进行无条件稳定的隐式时间离散[13-17]，极大地放宽了对时间步长的限制，而时间域有限元法可以采用非规则四面体网格进行空间离散，便于处理起伏地形和复杂地电结构的情况。

随着数值计算技术的发展，瞬变电磁三维正演的速度和精度均有了很大程度的提高，但仍然缺乏模型建立和剖分简便、交互便利、便于推广使用的成熟瞬变电磁三维正演软件。在目前已有的商业软件中，ANSYS是比较流行的大型有限元分析软件之一[18]，其Multiphysics模块可以分析解决多学科问题，如：结构、力学、热学、流体、电磁场以及任意两种或两种以上物理场之间的耦合问题，其低频电磁学模块所采用的控制方程与瞬变电磁法正演模拟的控制方程相同(均为扩散方程)[19-22]，其正演算法为时域有限元法，具备网格剖分操作简单、便于处理大型复杂构造问题、计算速度较快的特点，适合用于瞬变电磁法的三维正演模拟。鉴于此，本文尝试利用ANSYS软件进行瞬变电磁正演，并通过算例来分析ANSYS软件的计算精度。

1 ANSYS电磁仿真理论基础

电磁感应现象可由麦克斯韦方程组进行描述。麦克斯韦方程组实际上由4个电磁学定律组成，即安培环路定律、法拉第电磁感应定律、高斯电通定律和高斯磁通定律，其微分形式为：

(1)

式中：E、H分别为电磁场强度矢量，D、B分别为电位移矢量和磁感应强度，ρ为体电荷密度。为了计算简便，常常定义磁通量势和电标量势：

(2)

使得耦合方程(1)转化为关于电场和磁场的两个偏微分方程

(3)

式中：μ、ε分别为磁导率和介电常数。利用有限元方法求解上述偏微分方程，即可获得磁矢势和电势的值，然后根据式(2)即可导出电磁场强度、涡流密度等物理量。

实际上，求解偏微分方程(3)时，必须要先给定初始条件和边界条件，才能形成一个完整的初边问题。一般地，利用ANSYS进行瞬变电磁问题模拟时，初始条件往往要通过一个静态分析获得；边界条件一般使用Dirichlet边界条件(即直接给定磁矢势和电势在边界上的值)或Neumann边界条件(即通量平行或垂直条件)。

2 ANSYS电磁仿真基本流程

ANSYS作为一种有限元仿真工具，求解电磁问题的一般思路就是依照前节给出的理论基础，将所处理区域划分为有限个单元，然后将磁矢势和电势所满足的变分问题极小化以获得节点上的磁矢势和电势值，继而导出计算区域的其他待求量。其基本流程为：

1) 定义物理环境，主要包括单元类型选取、材料定义、坐标系选用等；

2) 进行几何建模，并分配几何体材料特性、单元类型等属性，对求解区域按照一定规则进行剖分；

3) 根据问题的特征，施加边界条件和载荷；

4) 求解和后处理，利用ANSYS求解器求解磁矢势和电势，然后利用其后处理器获得与问题相关的其他量。

3 基于ANSYS的瞬变电磁三维正演算例

为分析和验证利用ANSYS软件进行瞬变电磁三维正演的正确性和精度，分别设计了均匀半空间模型、层状地层模型和单个三维异常体模型进行了仿真模拟。

3.1 均匀半空间模型

3.1.1 模型描述

如图1所示，仿真模型为回线源激发的均匀半空间，回线半径为100 m，均匀半空间电阻率值为100 Ω·m。剖分区域分别由线圈区、空气区和地层区三部分组成， 其外边界为一个半径为10 000 m的球面。在t<0时，回线中通以20 A的电流，t=0时刻，撤去发射电流，观测此后的电磁场变化规律。

图1 瞬变电磁模型示意Fig.1 The diagram of TEM model

3.1.2 几何模型

几何模型按照物理模型及各分析对象的特性建立，使用直角坐标系，其x-y平面与地层和空气的分界面重合，z轴正方向定义为地层深部方向。线圈位于空气层中，与地层—空气层分界面相接。为方便电流密度载荷的施加和剖分，使用较粗的导线围成发射线圈，其横截面积为1 m2，作为真实线圈的近似。应该指出的是，在距离线圈较近时(<5 m)，这种近似的仿真场和细导线产生的电磁场具有较大的偏差，而在较远的区域，这种近似引起的误差可以忽略。同时，对于这种具有明显柱对称性质的问题，将其简化为二维问题可大大减小仿真的计算量，此处考虑到进行三维仿真的普适性，仍将其作为三维问题进行分析。

3.1.3 有限元模型

1)单元类型。ANSYS软件为三维电磁场分析提供了solid97、solid236、solid237等单元类型，其中solid236、solid237为矢量单元类型，基于矢量单元的有限元分析将由ANSYS软件自动采用矢量有限元法进行分析。考虑到基于标量有限单元法求解三维电磁场模拟问题时的“伪解”现象，本次分析使用矢量单元solid236，按照单元选项不同分为2类(表1)。

表1 单元类型

2)材料属性。本次分析地层、线圈、空气3种介质。对应这三种介质，选取3种材料属性，如表2所示。

表2 材料属性

3)网格控制。网格剖分的质量直接影响着计算结果的精度。Solid236单元为一阶单元，其四面体结构属于常应变单元，六面体结构属于梯度单元，当剖分粒度相同时，六面体网格相比于四面体网格具有更高的计算精度[23]，因此应尽量将计算区域剖分为六面体。本例中，线圈使用扫掠方式[24]生成规则网格，地层和空气区域采用自由网格划分，自动细化级别设为4。剖分完成后，共有549 036个单元，92 165个节点。

4)加载与求解。在线圈中加载20 A/m2的电流密度以实现发射电流的加载，所有外边界节点上加载简单狄利克雷边界，即令Az=0，V=0。分析类型选择Transient，求解方法选择FULL。考虑到发射线圈中在t<0时有稳定电流，应先将时间积分效应开关关闭，进行一次静态磁场分析以作为瞬态分析的初始条件，然后打开时间积分效应开关，进行瞬态分析。自动分析采用自动时间步长，最大时间步为1 ms，最小为1 μs，初始时间步设为10 μs。

5) 后处理及分析

ANSYS直接求解的是磁矢势和电势，在求解完成后，可以在后处理模块提取电磁场值、涡流等参数。为验证ANSYS的模拟精度，首先提取中心点时间域磁场响应，并与中心点磁场垂直分量的解析解进行对比。其解析解表达式为

(4)

对比结果如图2所示。容易看出，ANSYS解和解析结果吻合较好，最大相对误差约为3.5%，说明ANSYS模拟结果具有较高的计算精度。本次正演模拟在CPU核心数为4核，主频为2.20 GHz的PC上运行，耗时约为30 min。

图2 均匀半空间ANSYS计算磁场响应和解析结果对比Fig.2 The comparison diagram of ANSYS and analytical solution of homogeneous half-space

3.2 层状介质模型

为分析ANSYS对于层状大地瞬变电磁响应的模拟精度，在均匀半空间中加入一个低阻层，使其成为H型大地模型，然后利用ANSYS进行模拟计算，将计算结果与数字滤波解进行对比。本算例的设计参数为：采用方形回线进行发射，回线边长为100 m，发射电流为1A；H型地层模型的背景电阻率值取为100 Ω·m，低阻层电阻率值为10 Ω·m，深度范围为200～300 m；观测时间为10-5～10-1s，按照对数间隔设置21个时间道。图3中给出了中心点的衰减电压计算结果，并和数字滤波解进行了对比。计算结果表明，对于H型地层，ANSYS的模拟结果和数字滤波解吻合较好，最大相对误差约为4.5%。

图3 H型地层ANSYS计算衰减电压和数字滤波结果对比Fig.3 The comparison diagram of ANSYS and analytical solution of H-type model

3.3 单个异常体模型

在均匀半空间中放置一个异常体，然后计算主剖面上的瞬变电磁响应，以分析单个异常体的瞬变电磁响应规律。如图4所示，模型的设计参数为：采用方形回线进行发射，回线边长为100 m，发射电流为1 A；大地模型的背景电阻率值取为100 Ω·m；异常体电阻率值为1 Ω·m，中心位置为(0 m, 0 m, 250 m)，规模为200 m×200 m×100 m；观测时间为10-5～10-1s，按照对数间隔设置21个时间道；主剖面测线范围为-200～200 m，测点间距50 m。

图4 单个异常体模型示意Fig.4 The diagram of single anomalous body model

图5为主剖面上的多测道曲线。从多测道图中可以看出，当测点通过低阻异常体正上方时，第10～15道的衰减电压出现了“上凸”的特征。这是符合瞬变电磁场理论的：由于在瞬变电磁场传播过程中，遇到低阻体将在其中感生出较强的涡流，使得瞬变电磁场的衰减速率减缓，在多测道曲线上即表现为“上凸”特征。

图5 主剖面多测道曲线Fig.5 The multi-channel graph of principal section

4 结论

本文将ANSYS软件应用在瞬变电磁场三维正演之中，介绍了利用ANSYS软件进行前处理、加载、单元分析、求解、后处理等流程中的参数设置方法。利用瞬变电磁均匀半空间解析解、层状地层模型的数字滤波解进行了ANSYS瞬变电磁场仿真的正确性验证，并进行了简单三维异常体的TEM响应试算。结果表明，ANSYS仿真结果与解析解和数字滤波解吻合、计算误差较小，三维仿真结果与瞬变电磁场理论相符。利用ANSYS进行瞬变电磁场三维正演具有较高的计算精度和速度，具备较强的实用性。

TheapplicationofANSYStoTEM3Dforwardmodeling