智慶全，武軍杰，王興春，孫懷鳳， 楊毅，張杰，鄧曉紅，陳曉東,杜利明

(1.中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，河北 廊坊 065000; 2.自然資源部地球物理電磁法探測技術(shù)重點實驗室，河北 廊坊 065000; 3.山東大學(xué) 巖土與結(jié)構(gòu)工程研究中心，山東 濟南 250061; 4.中國冶金地質(zhì)總局 山東正元地質(zhì)勘查院，山東 濟南 250101)

0 引言

瞬變電磁法是應(yīng)用最廣泛的電磁探測方法之一，近年來在金屬礦勘查、煤田地質(zhì)勘探、隧道超前預(yù)報等領(lǐng)域[1-5]發(fā)揮了重要作用。瞬變電磁法一般采用不接地回線或接地導(dǎo)線作為發(fā)射源，供以穩(wěn)恒電流，觀測電流關(guān)斷后地下地質(zhì)體感生的二次電磁場，通過分析二次場的特征來推斷地下地質(zhì)體的電性、位置、深度等信息。利用瞬變電磁法進行勘探工作的關(guān)鍵技術(shù)在于根據(jù)瞬變電磁響應(yīng)規(guī)律推斷地質(zhì)體信息，而瞬變電磁正演是研究瞬變電磁響應(yīng)規(guī)律的有效途徑。隨著計算機技術(shù)的應(yīng)用，瞬變電磁正演方法有了快速發(fā)展[6-9]，其中，一維層狀模型的瞬變電磁正演通?；陬l譜法，采用線性數(shù)字濾波法進行快速計算，二維、2.5維的瞬變電磁正演主要采用有限單元法和有限差分法，三維瞬變電磁正演方法主要為有限單元法、積分方程法、時域有限差分法等。目前，一維瞬變電磁正演方法經(jīng)過多年發(fā)展已漸趨成熟，基于一維正演的一維反演技術(shù)是當(dāng)前進行瞬變電磁反演的主要方法，而二維、2.5維瞬變電磁方法由于其二維近似的條件常常難以滿足，在實際應(yīng)用中存在諸多限制；三維瞬變電磁正演方法相對一維、二維方法能更精確地模擬實際地質(zhì)情況，基于三維正演的解釋方法能取得更高的解釋精度，但一般需要花費很高的時間成本。

瞬變電磁的高精度、快速三維正演技術(shù)已成為當(dāng)前瞬變電磁正演研究的熱點，國內(nèi)外很多學(xué)者都針對這一問題展開了研究。目前，瞬變電磁的三維正演技術(shù)主要可分為頻時轉(zhuǎn)換法和直接時域法兩大類。頻時轉(zhuǎn)換法首先計算三維頻率域電磁響應(yīng)，然后采用余弦變換、G-S變換等頻時技術(shù)將頻率域響應(yīng)轉(zhuǎn)換到時間域，進而獲得瞬變電磁響應(yīng)。時域法一般通過時間離散來直接求取時域響應(yīng)，主要包括時域有限差分法、有限體積法、有限單元法等。時間域有限差分法基于交錯網(wǎng)格和顯式時間差分格式[10-12]，該方法具有數(shù)學(xué)概念直觀、表達簡單的特點，但是需要滿足穩(wěn)定性條件，其網(wǎng)格尺寸和時間步長均受到嚴格限制。時間域有限體積法和時間域有限元法通過采用后退歐拉時間離散格式對控制方程進行無條件穩(wěn)定的隱式時間離散[13-17]，極大地放寬了對時間步長的限制，而時間域有限元法可以采用非規(guī)則四面體網(wǎng)格進行空間離散，便于處理起伏地形和復(fù)雜地電結(jié)構(gòu)的情況。

隨著數(shù)值計算技術(shù)的發(fā)展，瞬變電磁三維正演的速度和精度均有了很大程度的提高，但仍然缺乏模型建立和剖分簡便、交互便利、便于推廣使用的成熟瞬變電磁三維正演軟件。在目前已有的商業(yè)軟件中，ANSYS是比較流行的大型有限元分析軟件之一[18]，其Multiphysics模塊可以分析解決多學(xué)科問題，如：結(jié)構(gòu)、力學(xué)、熱學(xué)、流體、電磁場以及任意兩種或兩種以上物理場之間的耦合問題，其低頻電磁學(xué)模塊所采用的控制方程與瞬變電磁法正演模擬的控制方程相同(均為擴散方程)[19-22]，其正演算法為時域有限元法，具備網(wǎng)格剖分操作簡單、便于處理大型復(fù)雜構(gòu)造問題、計算速度較快的特點，適合用于瞬變電磁法的三維正演模擬。鑒于此，本文嘗試利用ANSYS軟件進行瞬變電磁正演，并通過算例來分析ANSYS軟件的計算精度。

1 ANSYS電磁仿真理論基礎(chǔ)

電磁感應(yīng)現(xiàn)象可由麥克斯韋方程組進行描述。麥克斯韋方程組實際上由4個電磁學(xué)定律組成，即安培環(huán)路定律、法拉第電磁感應(yīng)定律、高斯電通定律和高斯磁通定律，其微分形式為：

(1)

式中：E、H分別為電磁場強度矢量，D、B分別為電位移矢量和磁感應(yīng)強度，ρ為體電荷密度。為了計算簡便，常常定義磁通量勢和電標量勢：

(2)

使得耦合方程(1)轉(zhuǎn)化為關(guān)于電場和磁場的兩個偏微分方程

(3)

式中：μ、ε分別為磁導(dǎo)率和介電常數(shù)。利用有限元方法求解上述偏微分方程，即可獲得磁矢勢和電勢的值，然后根據(jù)式(2)即可導(dǎo)出電磁場強度、渦流密度等物理量。

實際上，求解偏微分方程(3)時，必須要先給定初始條件和邊界條件，才能形成一個完整的初邊問題。一般地，利用ANSYS進行瞬變電磁問題模擬時，初始條件往往要通過一個靜態(tài)分析獲得；邊界條件一般使用Dirichlet邊界條件(即直接給定磁矢勢和電勢在邊界上的值)或Neumann邊界條件(即通量平行或垂直條件)。

2 ANSYS電磁仿真基本流程

ANSYS作為一種有限元仿真工具，求解電磁問題的一般思路就是依照前節(jié)給出的理論基礎(chǔ)，將所處理區(qū)域劃分為有限個單元，然后將磁矢勢和電勢所滿足的變分問題極小化以獲得節(jié)點上的磁矢勢和電勢值，繼而導(dǎo)出計算區(qū)域的其他待求量。其基本流程為：

1) 定義物理環(huán)境，主要包括單元類型選取、材料定義、坐標系選用等；

2) 進行幾何建模，并分配幾何體材料特性、單元類型等屬性，對求解區(qū)域按照一定規(guī)則進行剖分；

3) 根據(jù)問題的特征，施加邊界條件和載荷；

4) 求解和后處理，利用ANSYS求解器求解磁矢勢和電勢，然后利用其后處理器獲得與問題相關(guān)的其他量。

3 基于ANSYS的瞬變電磁三維正演算例

為分析和驗證利用ANSYS軟件進行瞬變電磁三維正演的正確性和精度，分別設(shè)計了均勻半空間模型、層狀地層模型和單個三維異常體模型進行了仿真模擬。

3.1 均勻半空間模型

3.1.1 模型描述

如圖1所示，仿真模型為回線源激發(fā)的均勻半空間，回線半徑為100 m，均勻半空間電阻率值為100 Ω·m。剖分區(qū)域分別由線圈區(qū)、空氣區(qū)和地層區(qū)三部分組成， 其外邊界為一個半徑為10 000 m的球面。在t<0時，回線中通以20 A的電流，t=0時刻，撤去發(fā)射電流，觀測此后的電磁場變化規(guī)律。

圖1 瞬變電磁模型示意Fig.1 The diagram of TEM model

3.1.2 幾何模型

幾何模型按照物理模型及各分析對象的特性建立，使用直角坐標系，其x-y平面與地層和空氣的分界面重合，z軸正方向定義為地層深部方向。線圈位于空氣層中，與地層—空氣層分界面相接。為方便電流密度載荷的施加和剖分，使用較粗的導(dǎo)線圍成發(fā)射線圈，其橫截面積為1 m2，作為真實線圈的近似。應(yīng)該指出的是，在距離線圈較近時(<5 m)，這種近似的仿真場和細導(dǎo)線產(chǎn)生的電磁場具有較大的偏差，而在較遠的區(qū)域，這種近似引起的誤差可以忽略。同時，對于這種具有明顯柱對稱性質(zhì)的問題，將其簡化為二維問題可大大減小仿真的計算量，此處考慮到進行三維仿真的普適性，仍將其作為三維問題進行分析。

3.1.3 有限元模型

1)單元類型。ANSYS軟件為三維電磁場分析提供了solid97、solid236、solid237等單元類型，其中solid236、solid237為矢量單元類型，基于矢量單元的有限元分析將由ANSYS軟件自動采用矢量有限元法進行分析?？紤]到基于標量有限單元法求解三維電磁場模擬問題時的“偽解”現(xiàn)象，本次分析使用矢量單元solid236，按照單元選項不同分為2類(表1)。

表1 單元類型

2)材料屬性。本次分析地層、線圈、空氣3種介質(zhì)。對應(yīng)這三種介質(zhì)，選取3種材料屬性，如表2所示。

表2 材料屬性

3)網(wǎng)格控制。網(wǎng)格剖分的質(zhì)量直接影響著計算結(jié)果的精度。Solid236單元為一階單元，其四面體結(jié)構(gòu)屬于常應(yīng)變單元，六面體結(jié)構(gòu)屬于梯度單元，當(dāng)剖分粒度相同時，六面體網(wǎng)格相比于四面體網(wǎng)格具有更高的計算精度[23]，因此應(yīng)盡量將計算區(qū)域剖分為六面體。本例中，線圈使用掃掠方式[24]生成規(guī)則網(wǎng)格，地層和空氣區(qū)域采用自由網(wǎng)格劃分，自動細化級別設(shè)為4。剖分完成后，共有549 036個單元，92 165個節(jié)點。

4)加載與求解。在線圈中加載20 A/m2的電流密度以實現(xiàn)發(fā)射電流的加載，所有外邊界節(jié)點上加載簡單狄利克雷邊界，即令A(yù)z=0，V=0。分析類型選擇Transient，求解方法選擇FULL?？紤]到發(fā)射線圈中在t<0時有穩(wěn)定電流，應(yīng)先將時間積分效應(yīng)開關(guān)關(guān)閉，進行一次靜態(tài)磁場分析以作為瞬態(tài)分析的初始條件，然后打開時間積分效應(yīng)開關(guān)，進行瞬態(tài)分析。自動分析采用自動時間步長，最大時間步為1 ms，最小為1 μs，初始時間步設(shè)為10 μs。

5) 后處理及分析

ANSYS直接求解的是磁矢勢和電勢，在求解完成后，可以在后處理模塊提取電磁場值、渦流等參數(shù)。為驗證ANSYS的模擬精度，首先提取中心點時間域磁場響應(yīng)，并與中心點磁場垂直分量的解析解進行對比。其解析解表達式為

(4)

對比結(jié)果如圖2所示。容易看出，ANSYS解和解析結(jié)果吻合較好，最大相對誤差約為3.5%，說明ANSYS模擬結(jié)果具有較高的計算精度。本次正演模擬在CPU核心數(shù)為4核，主頻為2.20 GHz的PC上運行，耗時約為30 min。

圖2 均勻半空間ANSYS計算磁場響應(yīng)和解析結(jié)果對比Fig.2 The comparison diagram of ANSYS and analytical solution of homogeneous half-space

3.2 層狀介質(zhì)模型

為分析ANSYS對于層狀大地瞬變電磁響應(yīng)的模擬精度，在均勻半空間中加入一個低阻層，使其成為H型大地模型，然后利用ANSYS進行模擬計算，將計算結(jié)果與數(shù)字濾波解進行對比。本算例的設(shè)計參數(shù)為：采用方形回線進行發(fā)射，回線邊長為100 m，發(fā)射電流為1A；H型地層模型的背景電阻率值取為100 Ω·m，低阻層電阻率值為10 Ω·m，深度范圍為200～300 m；觀測時間為10-5～10-1s，按照對數(shù)間隔設(shè)置21個時間道。圖3中給出了中心點的衰減電壓計算結(jié)果，并和數(shù)字濾波解進行了對比。計算結(jié)果表明，對于H型地層，ANSYS的模擬結(jié)果和數(shù)字濾波解吻合較好，最大相對誤差約為4.5%。

圖3 H型地層ANSYS計算衰減電壓和數(shù)字濾波結(jié)果對比Fig.3 The comparison diagram of ANSYS and analytical solution of H-type model

3.3 單個異常體模型

在均勻半空間中放置一個異常體，然后計算主剖面上的瞬變電磁響應(yīng)，以分析單個異常體的瞬變電磁響應(yīng)規(guī)律。如圖4所示，模型的設(shè)計參數(shù)為：采用方形回線進行發(fā)射，回線邊長為100 m，發(fā)射電流為1 A；大地模型的背景電阻率值取為100 Ω·m；異常體電阻率值為1 Ω·m，中心位置為(0 m, 0 m, 250 m)，規(guī)模為200 m×200 m×100 m；觀測時間為10-5～10-1s，按照對數(shù)間隔設(shè)置21個時間道；主剖面測線范圍為-200～200 m，測點間距50 m。

圖4 單個異常體模型示意Fig.4 The diagram of single anomalous body model

圖5為主剖面上的多測道曲線。從多測道圖中可以看出，當(dāng)測點通過低阻異常體正上方時，第10～15道的衰減電壓出現(xiàn)了“上凸”的特征。這是符合瞬變電磁場理論的：由于在瞬變電磁場傳播過程中，遇到低阻體將在其中感生出較強的渦流，使得瞬變電磁場的衰減速率減緩，在多測道曲線上即表現(xiàn)為“上凸”特征。

圖5 主剖面多測道曲線Fig.5 The multi-channel graph of principal section

4 結(jié)論

本文將ANSYS軟件應(yīng)用在瞬變電磁場三維正演之中，介紹了利用ANSYS軟件進行前處理、加載、單元分析、求解、后處理等流程中的參數(shù)設(shè)置方法。利用瞬變電磁均勻半空間解析解、層狀地層模型的數(shù)字濾波解進行了ANSYS瞬變電磁場仿真的正確性驗證，并進行了簡單三維異常體的TEM響應(yīng)試算。結(jié)果表明，ANSYS仿真結(jié)果與解析解和數(shù)字濾波解吻合、計算誤差較小，三維仿真結(jié)果與瞬變電磁場理論相符。利用ANSYS進行瞬變電磁場三維正演具有較高的計算精度和速度，具備較強的實用性。

TheapplicationofANSYStoTEM3Dforwardmodeling