1.1.3 有限單元法
有限單元法(Finite Element Methods, FEM)于20世紀(jì)60年代被提出, 最初用于結(jié)構(gòu)力學(xué)和應(yīng)力分析(齊基威次, 邱, 1967); 1965年, Winslow首先將有限單元法應(yīng)用于電氣工程問題; 其后, 在1969年Silvester將有限單元法推廣應(yīng)用到時諧電磁場問題; 1971年, Coggon發(fā)表了著名的關(guān)于電和電磁的有限單元法的論文, 拉開了有限單元法在電法勘探領(lǐng)域正演模擬的序幕, 他從電磁場總能量最小原理出發(fā), 實現(xiàn)了二維地電斷面有限單元法正演計算, 不過由于有限單元網(wǎng)格缺乏通用性, 計算精度和速度未能達到實用水平; 1976年Rodi發(fā)展了有限單元法的剖分方法, 采用矩形網(wǎng)格剖分, 以解決二維大地電磁測深正演問題, 使有限元向前發(fā)展了一大步; 1977年Rijo引入了一個通用性網(wǎng)格剖分方法, 使有限元正演的精度和速度得到了大幅度提高, 成為計算二維地電條件下電法/電磁法正演模擬的有效工具, 使有限單元法正式進入實用階段。20世紀(jì)70年代末, 有限單元法在甚低頻法正演(Kaikkonen, 1979)、 時間域電磁法的正演模擬(Jone 等, 1980)中得到應(yīng)用; Dey和Morrison(1979)將混合邊界條件引入到有限單元法的正演模擬, 使其得到進一步發(fā)展; 1978年, Pridmore等對有限單元法的三維地電斷面的電法/電磁法正演問題進行論述; Wannamaker等(1984, 1985, 1986)用有限單元法模擬大地電磁中的二維地形響應(yīng), 并且編制了求解二維大地電磁響應(yīng)的正演程序。Morgan(1990)詳細地研究了電磁散射中的有限單元和有限差分方法; Unsworth等(1993)發(fā)表了頻率域電流偶極源電磁場的2.5維有限單元法模擬; Albanese和Rubinacci(1998)利用有限單元法對三維渦流場進行了計算研究。Badea等(2001)研究了庫侖規(guī)范下的可控源音頻大地電磁法的有限元模擬; Pridmore等人(1981)研究了利用二次場進行三維模擬時的散度校正; Yuji Mitsuhata和Toshihiro Uchida(2004)研究了基于T-Ω的三維大地電磁有限元正演模擬; 為了解決三維大地電磁有限元法中的偽解問題, Kerry Key和Chester Weiss(2006)研究了二維大地電磁的非結(jié)構(gòu)化網(wǎng)格的自適應(yīng)有限單元法并對兩種后驗誤差的估計進行了比較。
在國內(nèi)的研究中, 我國著名的已故數(shù)學(xué)家馮康先生為有限元理論的建立做出了巨大的貢獻, 在20世紀(jì)50年代獨立于西方學(xué)者創(chuàng)立了數(shù)值求解偏微分方程的有限元方法, 形成了標(biāo)準(zhǔn)的算法形態(tài), 并先于西方建立了嚴(yán)密的理論體系。20世紀(jì)80年代初, 數(shù)學(xué)家李大潛發(fā)表專著《有限元法在電法測井中的應(yīng)用》, 有限單元法正式應(yīng)用到電法領(lǐng)域, 此后, 地球物理學(xué)家周熙襄、 鐘本善(1980, 1986)、 羅延鐘(1986, 1987, 2000, 2006)和徐世浙(1983, 1985, 1994a, 1994b, 1995, 1996, 1997)等對有限單元法進行了深入研究, 發(fā)表和出版了一系列有限單元法在電法中應(yīng)用的論文或?qū)V? 研究領(lǐng)域涉及直流電法、 電磁法領(lǐng)域, 網(wǎng)格剖分也由原來的簡單矩形剖分發(fā)展到三角單元和三角-矩形綜合剖分。周熙襄等在二維電阻率有限單元法正演中采用混合邊界條件等優(yōu)化措施, 提高了計算精度和速度, 對有限單元法的發(fā)展起到了推動作用; 從20世紀(jì)80年代開始, 國內(nèi)學(xué)者開始研究2.5維電磁場的數(shù)值模擬, 到目前, 已用二次場算法實現(xiàn)了二維地電構(gòu)造上諧變電偶極子電磁場的有限單元算法、 電偶源CSAMT法二維正演的有限單元算法、 主剖面上時間域和頻率域電偶源瞬變電磁場的2.5維有限單元正演模擬; 阮百堯等(2001, 2002)利用有限單元法實現(xiàn)了三維地電斷面的正演模擬并編制了相應(yīng)程序。黃臨平和戴世坤(2002)研究了復(fù)雜條件下3D 電磁場有限元計算方法, 推導(dǎo)出了地下變電導(dǎo)率σ條件下計算三維電磁場的有限元單元方程的解析表達式, 采用伽遼金方法推導(dǎo)出了散度校正有限元方程組, 為三維電磁場的數(shù)值計算提供了一條有效的新途徑。強建科(2006, 2007)深入研究了起伏地形下三維電阻率有限元正演模擬, 黃俊革(2002, 2003)實現(xiàn)了齊次邊界條件下三維地電斷面電阻率有限元的數(shù)值模擬。馮兵、 張繼鋒、 王永剛等人(2005, 2006)做了有限元方法探測油氣藏邊界的有限元數(shù)值模擬方面的研究。蔡軍濤, 趙國澤等人(2007)進行了復(fù)電阻率法的二維有限元數(shù)值模擬。
不同于FDM, FEM是基于電磁場的積分形式, 它是由電磁場的微分形式通過Green等定理變換而得, 通常也稱有限元法的解為微分形式的弱解, 并將定解區(qū)域剖分成規(guī)則或不規(guī)則的單元, 單元上的場由節(jié)點上的場值插值得到, 通過求泛函極值, 得到關(guān)于節(jié)點上電磁場值的線性方程組, AFEX=B, AFE為大型的對稱正定、 復(fù)數(shù)、 稀疏矩陣。FEM 并不一定要求模型能夠被剖分成規(guī)則單元, 如三角形與四面體單元(其被理論與實踐證明可以無限度精確地模擬地球物理模型), 因此, FEM能夠求解FDM不能夠求解的復(fù)雜地球物理模型。解方程組求出節(jié)點上的場值后, 由插值可得到整個定解區(qū)域上的電磁場。
FEM是依據(jù)變分原理或加權(quán)余值法導(dǎo)出與定解問題等價的積分弱解, FEM是邊值問題中數(shù)學(xué)理論最為完善的一種方法, 而且對單元形狀沒有要求。由于三角形和四面體單元可以足夠精確地剖分復(fù)雜的電磁模型, FEM能夠求解邊界復(fù)雜、 地形、 介質(zhì)變化大等FDM、 IEM難以求解的復(fù)雜模型, 且所形成的系數(shù)矩陣對稱正定稀疏、 容易求解、 收斂性好、 占用內(nèi)存少、 程序通用性強, 因此得到了廣泛的研究和應(yīng)用。國內(nèi)外許多學(xué)者對FEM做了大量的研究, 內(nèi)容包括FEM在直流電阻率、 大地電磁測深、 可控源電磁測深等方法中的應(yīng)用, 內(nèi)容涉及到2D、 2.5D及3D電磁場計算的基本理論、 單元及形狀函數(shù)、 吸收邊界條件、 線性方程組求解、 實際應(yīng)用等各個方面。但FEM單元的剖分常常沒有規(guī)律, 對無邊界區(qū)域的求解不方便。