2018, Vol. 61
2. 有色资源与地质灾害探查湖南省重点实验室, 长沙 410083
2. Non-ferrous Resources and Geologic Disasters Prospecting Emphases Laboratory of Hunan, Changsha 410083, China
无单元Galerkin法(Element-Free Galerkin Method, EFGM)是一种基于离散节点的数值模拟方法, 其利用局部支持域内的节点空间信息构造形函数, 通过Galerkin全局弱式将计算域上的弱式积分方程转化为离散系统方程, 具有模拟精度高, 预处理简单, 节点可根据需要任意布置, 适应性强等特点.起初, 随着滑动最小二乘(Moving Least Squares, MLS)被用于构造形函数, 使得无网格法(Meshless method)形函数的连续性问题而得以解决(Lancaster and Salkauskas, 1981), 在此基础上Belytschko等(1994a)改进了强加本质边界条件及数值积分处理等问题并提出EFGM, 随后EFGM在各个研究领域被广泛关注(Peng et al., 2011; Hajiazizi and Bastan, 2014; Hadinia and Jafari, 2015).EFGM已被证明是一种非常优秀的数值模拟方法, 在模拟精度、灵活性、便利性等方面均优于有限单元法(Finite Element Method, FEM), 不仅能应用到FEM的应用领域, 且可解决FEM由于网格的存在而难以处理的问题, 如场剧烈变化问题, 大变形问题, 复杂边界和形体问题等(Belytschko et al., 1994b; Li et al., 2012).在地球物理勘探领域方面, EFGM已经得到了一些研究和应用, 贾晓峰等(2006)对地震波场应用EFGM进行了正演模拟研究; 王月英等(2007)探讨了EFGM应用于地震波正演模拟中的各种影响因素; 冯德山等(2013)采用MLS构造形函数, 实现了探地雷达二维EFGM正演模拟; 严家斌和李俊杰(2014)也采用MLS构造形函数, 将EFGM应用于大地电磁法二维模拟, 计算了水平地形下大地电磁二维模型响应; 嵇艳鞠等(2016)利用EFGM计算了各向异性的2D大地电磁响应; 麻昌英等(2017)实现了直流电阻率EFGM正演模拟, 并与FEM进行了比较, 展现了EFGM的高精度和高灵活性特点.
然而, EFGM采用高斯积分计算积分, 由于EFGM的近似函数一般为有理函数而不是多项式, 为了保证模拟精度须使用较多的积分点和较高阶的高斯积分点(Dolbow and Belytschko, 1999)进行计算, 而EFGM计算量与积分点数成正比, 使得EFGM计算量大, 计算效率低.另外, 在直流电电阻率正演模拟中, 当采用第一类边界条件时, 计算域需要足够大, 此时在远离场源和异常体的计算域外围区域场值变化平缓, EFGM虽然对任意分布的节点具有很强的适应性, 但仍要求节点合理分布, 节点的不合理分布可能造成EFGM形函数构造失败, 因此在外围区域需要较多的节点覆盖以保证节点的合理分布, 这导致EFGM计算成本明显增加.虽然FEM由于单元网格的不可或缺而受到许多约束, 但FEM具良好的数值稳定性, 单元的形状及其分布变化情况不会造成模拟计算失败, 在某些特殊区域可充分利用该特点发挥FEM的性能, 例如在远离场源的场值变化平缓的区域, 可快速延伸单元扩大计算域以满足第一类边界条件, 快速延伸的单元虽然形状细长, 但FEM仍可顺利进行计算, 同时对模拟精度影响很小.针对上述直流电阻率EFGM正演中采用第一类边界条件时存在的问题, 本文尝试在外围区域使用FEM进行计算, 并尝试将EFGM与FEM直接耦合, 耦合法可将不同的数值模拟方法相结合, 发挥单一方法的优点的同时又能避免各自的缺点.EFGM与FEM的耦合已经得到了相关的研究和应用(Belytschko et al., 1995; Liu et al., 2006; Hadinia et al., 2016), 但在传统EFGM与FEM的耦合中, 由于传统EFGM基于MLS插值构造的形函数不具有Kronecker delta函数性质, 不能满足边界条件, 所以若将各自得到的代数方程直接耦合, 则会造成两个子域在耦合界面上场值不连续, 使得EFGM与FEM的耦合存在困难, 因此如何处理EFGM区域和FEM区域的耦合界面关系就成为EFGM与FEM耦合需要解决的首要问题, 在这方面已经展开了相关研究并取得了一定的成果(Kumar et al., 2014; Pathak et al., 2015; Samira et al., 2017), 但这些研究成果均为使用特殊处理手段实现EFGM与FEM耦合, 增加了算法复杂度.
本文从EFGM形函数方面进行研究, 尝试采用不同的形函数构造方法来解决EFGM与FEM耦合困难问题.地球物理研究领域中, 麻昌英等(2017)利用径向基点插值法(Radial Point Interpolation Method, RPIM)构造EFGM形函数, 获得了具有Kronecker delta函数性质的RPIM形函数; 李俊杰和严家斌(2015)初步展开了大地电磁法的EFGM与FEM耦合法的研究, 提出了有限元-径向基点插值法, 获得了良好的效果.RPIM以径向基函数(Radial Base Function, RBF)为基础, RBF具有优良数学特性使得RPIM形函数能满足EFGM形函数需要满足一些基本要求, 如适应性、稳定性、一致性、紧支性、相容性、Kronecker delta函数性质等(Liu and Gu, 2001).本文采用RPIM插值构造EFGM形函数解决EFGM与FEM耦合困难问题, 在计算域外围区域将EFGM与FEM直接耦合.同时, 针对采用第一类边界条件EFGM存在的问题, 基于矩形单元剖分的FEM, 采用快速扩展的单元网格覆盖远离场源的计算域外围区域, 以扩大计算域使得第一类边界条件得到满足.在此基础之上, 提出直流电阻率的无单元Galerkin-有限单元耦合法(EFG-FE).通过不同模拟方法对层状模型测深和均匀半空间点电源电位进行模拟对比, 验证了EFG-FE算法的有效性, 表明了EFG-FE可方便地通过外围FEM区域扩大计算域来满足第一类边界条件, 获得高精度的模拟结果.最后, 基于相同的EFGM计算域模型, 将EFG-FE与采用第三类边界条件的EFGM进行对比研究, 结果表明本文提出的EFG-FE相比于采用第三类边界条件的EFGM提高了计算效率, 具有更好的模拟性能.
1 直流电阻率边值问题采用第一类边界条件, 2.5维稳定电流场波数域总电位函数U(x, λ)满足的边值问题及其对应的变分问题分别表示为(阮百尧, 2001)：
|
(1) |
|
(2) |
式中, Ω为计算域, ΓS为地表边界, Γ∞为截断边界, σ=(x, z)为介质电导率, U(λ, x, z)为波数域电位, λ为波数, I0为电流, δ为Kronecker delta函数, x=[x, z]T为Ω上的任意一点, A为场源点位置, 
|
图 1 直流电阻率二维模型 Fig. 1 Sketch showing a 2D model of DC resistivity |
在2.5维直流电阻率EFGM正演模拟中, 为求解域Ω中任意点xT=[x, z]构造局部支持域Ωq, 假设Ωq中包含有n个节点{x1, x2, …, xn}, 它们对应的波数域电位值为{U1, U2, …, Un}, 则任意点的电位值波数域U(x)的RPIM近似可表示为
|
(3) |
式中, pj(x)为空间坐标x中的单项式, m为基函数个数, 完备的p阶单项式pj(x)的一般形式为pT(x)=[1x z x2 xz z2…xp zp].RT(x)=[R(x, x1) R(x, x2) … R(x, xn)]为径向基函数向量, R(x, xi)为径向基函数, 本文采用复合二次函数(4)式作为径向基函数，公式为
|
(4) |
式(4)包含αc和q两个形状参数, 通过实验αc取1.0左右, q取1.03左右在直流电阻率正演模拟中能取得良好的效果; dc为计算点x支持域中节点的平均间距, 可表示为
|
(5) |
式中
|
(6) |
其中径向基函数对应的系数矩阵为
|
(7) |
多项式基函数对应的系数矩阵为
|
(8) |
RPIM形函数为
|
(9) |
其中支持域中节点对应的RPIM形函数为
|
(10) |
为避免对G求逆, (9)式两边同时乘以G，即：
|
(11) |
变换可得到：
|
(12) |
利用标准线性方程组求解器求解(12)式可得到RPIM形函数.利用(12)式可计算RPIM形函数的一阶和二阶导数，公式为
|
(13) |
|
(14) |
更高阶导数同理可得.由(13)—(14)式可知, 为了获得RPIM形函数及其导数均需要独立求解方程, 这是导致EFGM计算效率较低的重要因素.
3 直流电阻率的无单元-有限单元耦合法如图 2所示, 假设计算域Ω被分割为EFGM区域Ω1和FEM区域Ω2, Ω=Ω1∪Ω2, Ω1中包含N1个EFGM节点, Ω2中包含N2个FEM节点.其中EFGM区域Ω1边界上的N3个节点为EFGM和FEM共用节点, 则计算域Ω中节点总数为N=N1+N2－N3.Ω1为模型的核心区域, 采用EFGM进行正演计算.Ω2为外围扩边区域, 采用FEM进行正演计算, 采用快速延伸网格扩大计算域, 以满足第一类边界条件.特别指出的是, 如图 2所示EFGM区域被FEM区域包围, 与截断边界Γ∞无交集, 因此EFGM区域内无需做边界计算.
|
图 2 无单元Galerkin-有限单元耦合 Fig. 2 Sketch showing the coupling of EFGM and FEM |
EFGM的求解过程如图 3所示.EFGM中首先使用一组总数为N1的节点离散EFGM区域Ω1, 然后使用一组背景单元将Ω1覆盖, 接着在每一个背景单元中布置高斯积分点计算积分, 将所有背景单元的积分求和得到Ω1上的积分.在EFGM计算过程中, 对每一个高斯积分点需要独立构造局部支持域Ωq, 并利用支持域内所包含的节点信息构造无单元法形函数.特别指出的是, 在EFGM中节点与背景单元是相互独立的, 节点用于模拟地电模型为EFGM的关键模型信息, 背景单元用于计算积分, 与FEM中的单元概念是完全不同的.在保证数值积分精度的情况下, 如果可以直接在Ω1上布置高斯积分点, 则可使得EFGM更加简便, 但通常这样很难控制积分计算精度, 因此引入背景单元便于布置高斯积分点进行积分计算.
|
图 3 无单元Galerkin法求解 Fig. 3 Sketch showing the solving process of EFGM |
采用Galerkin全局弱式法可导出(2)式的直流电阻率EFGM方程为(麻昌英等, 2017)：
|
(15) |
其中KEFG为N1×N1维的EFGM刚度系数矩阵, UEFG为EFGM区域节点波数域电位对应的1×N1维的列向量, F为1×N1维的EFGM方程右端项列向量.在背景积分单元Ωe内, 对于每一个积分点在其局部支持域Ωq内(15)式写为子方程组形式，即：
|
(16) |
其中：
|
(17) |
|
(18) |
Kq为局部支持域子刚度系数矩阵, 它的各个元素的表达式分别为
|
(19) |
其中i, j=1, 2, …, n, n为局部支持域内包含的节点总数.
(19)式表明EFGM刚度系数矩阵与FEM刚度系数矩阵一样, 均为对称稀疏矩阵.由于RPIM形函数具有Kronecker delta函数性质, 右端项Fq的各个元素表达式为
|
(20) |
将所有背景积分单元内的积分点的局部支持域子方程组(16)式组装起来可得到直流电阻率EFGM方程(15)式.
3.2 有限单元法如图 2所示, 在FEM区域Ω2使用矩形单元剖分, 假设矩形单元内电导率σ为常数(在计算中也可采用双线性插值获得任意点的电导率), 矩形单元内任意一点近似场值Uh(x)线性变换, 对Uh(x)采用双线性插值, 根据文献徐世浙(1994), 在矩形单元子域Ωe上直流电阻率变分问题(14)式可写为子方程组形式，即：
|
(21) |
式中Ke=Ke1+Ke2为FEM单元子刚度系数矩阵, 它们的表达式分别为
|
(22) |
|
(23) |
|
(24) |
其中
将FEM区域Ω2所有矩形形单元的子方程组(19)式组装得到直流电阻率FEM方程组为
|
(25) |
式中KFE为由全部矩形单元子域Ωe的Ke相加组成的N2×N2维的FEM刚度系数矩阵, UFE为FEM区域内所有矩形单元节点上的波数域电位子向量Ue组合成的1×N2维列向量, S为1×N2维的FEM方程右端项列向量.
3.3 无单元Galerkin法与有限单元法耦合在无单元Galerkin-有限单元耦合法中, 耦合界面Γint上两种方法需满足电位相容性条件(Gu and Liu, 2005; Liu et al., 2016)，公式为
|
(26) |
I表示耦合界面上两种方法的共同节点.若耦合界面为介质分界面, 两侧介质电阻率不同, 此时耦合界面上还需满足电位自然边界条件(徐世浙, 1994)，公式为
|
(27) |
传统的无单元法Galerkin法基于MLS形函数, 由于MLS形函数不具有Kronecker delta函数性质, 因此不具有相容性, 若将各自得到代数方程直接耦合, 则会造成两个子域在耦合界面上场值不连续, 产生较大误差.通常, 在传统EFGM与FEM耦合中在两种方法的子域之间设置连接区域，并引入斜坡函数等技术以克服相容性问题(Liu et al., 2006; Kumar et al., 2014; Pathak et al., 2015; Hosseini et al., 2017), 这些处理手段增加了耦合难度和算法的复杂性.本文采用具有Kronecker delta函数性质的RPIM形函数, 使得(26)式自然得到满足, 因此可以克服传统MLS形函数不具有Kronecker delta函数性质而产生的耦合困难.在耦合方法中, 相容性条件(26)式是相容性条件和自然边界条件两项要求中最重要的偶合条件, 其必须得到满足, 即使自然边界条件(27)式在一些耦合方法中不能完全满足, 也可得到良好的效果, 如果两种方法的形函数在耦合界面上满足Kronecker delta函数性质, 则这两种数值方法在耦合界面上可直接进行耦合(Gu and Liu, 2005).目前, Liu等(2016)将基于RPIM形函数的局部Petrov-Galerkin无单元-有限单元耦合法应用于电磁计算领域中, 展示了两种方法子域耦合界面上RPIM形函数与FEM形函数的相容性, 并采用两种方法直接耦合的方式简化了耦合过程, 并获得了良好的效果.
事实上, EFGM和FEM主要不同之处在于采用不同的插值方法, 但采用单元插值的FEM和采用局部支持域插值的EFGM均为分片插值, 它们的形函数均满足紧支性, 即FEM形函数在单元外为零, EFGM形函数在局部支持域外为零.因此, 当EFGM形函数满足Kronecker delta函数性质, 且局部支持域限定在EFGM区域内时, 则与FEM单元连接的EFGM局部支持域可视为FEM的一个连续插值单元(广义单元).贾亮等(2007)等基于广义单元的概念, 将FEM单元和EFGM局部支持域视为广义单元, 利用在每一个单元内假设的形函数来分片地表示求解域上待求的未知场函数, 将FEM和EFGM直接耦合, 获得了良好的结果.张琰等(2008)利用RPIM方法近似函数满足插值条件, 将无单元法与有限单元法直接耦合, 应用于应力分析获得了良好的效果.
基于上述分析结合直流电阻率地质模型, 本文在模型核心区域采用EFGM计算, 发挥EFGM的灵活性、适应性强和高效性等特点.选择耦合面两侧介质电阻率相同的区域分割FEM区域和EFGM区域, 使得计算域内部介质界面上满足自然边界条件(27)式, 同时将EFGM局部支持域限定在EFGM区域, 将两种方法子域直接耦合.因此, 将(15)和(21)式根据节点整体编号组装在一起得到直流电阻率EFG-FE的离散方程为
|
(28) |
式中K=KEFG+KFE, U=UEFG+UFE, P=F+S, K为N×N维的EFG-FE刚度系数矩阵.为方便两种方法的计算, 首先在两种方法子域中对节点分别采用局部编号.在FEM区域Ω2, 节点局部编号按常用的编号方法进行编号, 如图 4a所示, Nz表示垂向方向FEM节点数.在EFGM区域往往采用不规则分布的场节点, 节点局部编号方法如图 4b所示, 先按X方向从左到右排列, 再按Z方向由上到下排列, 事实上除了节点分布不规则外, 这样的编号方式与FEM区域采用的节点编号方式是一样的.最后, 在全局域上全部节点(包含FEM节点和EFGM节点)也按先X方向从左到右排列, 再按Z方向由上到下排列进行整体编号, 这有助于将刚度系数矩阵非零元素带宽减小, K与有限单元法形成的刚度系数矩阵相似, 它是稀疏对称正定矩阵, 使用定带宽的方式储存稀疏刚度矩阵(徐世哲, 1994), 减少零元素的储存.采用MKL库中开源求解器Pardiso求解方程组(28)式, 得到计算域Ω上N个节点的波数域电位, 通过反傅立叶变换可得到节点的空间域电位(潘克家等, 2012).
|
图 4 有限单元法和无单元法子域节点局部编号 (a) FEM区域; (b) EFGM区域. Fig. 4 Sketches showing the nodal numbering of the FEM subdomain and the EFGM subdomain (a) FEM region; (b) EFGM region. |
分别应用EFG-FE和基于矩形和线性插值的FEM对如图 5所示的三层地电模型进行对称四极测深模拟, 验证EFG-FE程序的正确性.如图 5所示, 第一层电阻率ρ1=100 Ωm, 厚度h1=5 m, 第二层电阻率ρ2=50 Ωm, 厚度h2=5 m, 第三层电阻率ρ3=200 Ωm.建立均匀半空间模型水平方向宽度200 m, 垂直方向高100 m的EFGM区域.在该区域内, 地表附近和地层界面附近采用稠密的网格分布, 离场源较远的区域使用稀疏的网格分布, 节点总数15725.当使用FEM计算时, FEM单元数15456, 由于使用矩形单元, 因此采用的节点分布没有采用无单元法计算时分布灵活.当采用无单元法计算时, 为了避免节点不合理分布导致EFGM形函数构造失败, 在EFM计算时的节点分布基础上适当改变了一些节点的分布.在外围区域使用快速延伸的FEM网格扩大计算域, 该区域使用1062个FEM单元, 1190个节点, 使得计算域水平宽度超过4 km, 垂直高度超过2 km.以地面X=0 m为测点, 模拟不同供电极距对称四极测深, 图 6给出了两种方法模拟的视电阻率曲线及其相对误差.两种方法的模拟结果与解析解均吻合得很好, 相对误差均在1%以内, EFG-FE结果平均绝对相对误差为0.20%, FEM结果平均绝对相对误差为0.37%, 两种方法均具有很好的模拟精度, 表明了EFG-FE算法的有效性.
|
图 5 三层地电模型 Fig. 5 Schematic of a three-layer geoelectric model |
|
图 6 三层模型对称四极测深视电阻率及其相对误差 Fig. 6 Apparent resistivity curves and corresponding relative errors from the test model |
如图 5a所示, 建立均匀半空间模型水平方向宽度120 m(X:-60~60 m), 垂直方向高60 m(Z:0~-60 m)的EFGM区域Ω1, 电阻率为100 Ωm, 点电源设置在地表X=0 m处, 电流强度1 A.模型节点分布如图 7a所示, 为了提高模拟精度同时兼顾计算成本, 利用EFGM的灵活性设置节点以点电源位置为中心向外逐渐由稠密逐渐稀疏分布, 节点主要集中分布在电场强烈变化的场源附近, 共3057个节点.如图 7b所示, 在Ω1外围建立FEM区域Ω2, FEM单元设置为向外成双倍扩展延伸8层, 水平范围X:-2014~2104 m, 垂直范围Z:0~-2104 m, 共790个节点, 702个单元.
|
图 7 点电源均匀半空间模型 (a) EFGM区域节点分布; (b) FEM区域单元分布. Fig. 7 Sketch showing a point-source homogeneous half-space model (a) Nodal distribution of EFGM region; (b) Element distribution of FEM region. |
首先, 对图 7所示的计算域Ω(Ω=Ω1+Ω2)采用第一类边界条件的EFGM(记为F-EFGM)进行点电源电位模拟, 此时在Ω2区域采用向外逐渐稀疏的节点离散, 以保证节点合理分布, 使得EFGM形函数构造顺利, 计算域总节点数为7712.然后, 对图 7a所示的节点模型和EFGM计算域Ω1, 采用第三类边界条件的EFGM(记为T-EFGM)(麻昌英等, 2017)进行点电源电位模拟.最后, 将图 7a和b所示的EFGM区域Ω1和FEM区域Ω2相结合的, 采用第一类边界条件的EFG-FE耦合法进行点电源电位模拟.在1~20 m范围内选取离点电源位置不同距离的35个地面节点的模拟电位进行对比, F-EFGM、T-EFGM和EFG-FE的模拟结果如图 8所示.如图 8所示, 三种方法的模拟结果与解析解均吻合得很好, 具有很高的模拟精度, 进一步表明了本文的EFG-FE算法的有效性.表 1列出了在同一台计算机上三种正演模拟方法分别进行10次模拟统计的平均计算时间和35个采样点的模拟电位平均绝对相对误差, 三种方法均得到了良好的模拟结果, F-EFGM由于采用了更多的节点而模拟精度更好.由表 1分析, F-EFGM中在Ω2中采用EFGM计算, 大大增加了EFGM计算区域, 使得EFGM计算时间大幅增加, 且为了保证EFGM计算稳定, 大量增加了节点, 增加了方程求解时间.表 1表明EFG-FE的计算时间相比于F-EFGM明显降低了, 表明本文提出的EFG-FE可有效解决采用第一类边界条件时由于扩大计算域导致EFGM计算量大量增加的问题.对于单个电极模拟, T-EFGM和EFG-FE计算时间很接近, 下文将展开对多电极供电观测的地电模型进行模拟, 进一步研究EFG-FE的模拟性能.
|
图 8 均匀半空间点电源模拟的电位 Fig. 8 Point-source potential in the homogeneous half-space model |
|
表 1 不同模拟方法点电源电位模拟的计算时间和平均相对误差 Table 1 Calculation time and average relative errors of different simulation methods for point source potentials |
本小节应用上节所用的T-EFGM和EFG-FE对复杂地形起伏模型进行视电阻率观测模拟.如图 9所示, 建立水平方向宽度80 m(X:-70~70 m)、垂直方向宽度40 m(Z:0~-40 m)的复杂地形起伏模型的EFGM区域Ω1, 在X=-14 m处为山脊最高处(最高为7.6 m), X=10处为山谷最低处(最高为3.4 m).在直流电阻率EFGM正演模拟中, 可在局部域任意加密节点提高模拟精度, 如在场源附近加密节点可显著提高近场源模拟精度(麻昌英等, 2017).因此, 如图 9所示, 针对该复杂地形起伏模型, 在电场场值变化剧烈的场源区域(即浅地表区域)和地形变化区域采用密集的节点分布以提高模拟精度, 同时利用EFGM的强适应性, 采用任意的节点分布以适应任意起伏的地形, 若有进一步加密节点需求, 可直接加密节点而无需其他处理, 相比于FEM中加入新的节点需要重新剖分单元网格, 这体现了EFGM的便利性和灵活性.随着与场源及地形距离的增加, 电场场值变化趋于平缓, 节点分布随之设置为逐渐稀疏, 以降低计算成本.EFGM区域Ω1共使用7749个节点.外围采用与上节相同形式的FEM单元网格扩展方式建立EFM区域Ω2(如图 7b所示), 水平范围超过4 km宽, 垂直范围超过2 km高, 共740个节点, 657个单元.复杂地形起伏模型介质电阻率为100 Ωm.在地表X:-58~58 m范围内等间隔2 m布置59根供电和观测电极, 对模型进行高密度电法温纳装置视电阻率观测正演模拟, 两种方法正演模拟结果及它们所花费的计算时间分别如图 10和表 2所示.
|
图 9 复杂起伏地形模型EFGM区域节点分布 Fig. 9 Sketch showing the nodal distribution of a complex undulating terrain model in EFGM area |
|
图 10 复杂起伏地形模型正演视电阻率观测拟断面 (a) T-EFGM；(b) EFG-FE. Fig. 10 Pseudo-cross sections of apparent resistivity of a complex undulating terrain simulated model (a) T- EFGM; (b) EFG-FE. |
|
|
表 2 复杂起伏地形模型视电阻率观测不同模拟方法正演的计算时间 Table 2 Calculation time of different simulation methods for apparent resistivity survey of the complex undulating terrain model |
将图 10a和b对比分析可知, 两种方法的模拟结果几乎一致, 进一步表明当采用相同的EFGM模型时, 本文的EFG-FE可获得与T-EFGM一致的模拟精度.表 2列出了两种方法对复杂地形起伏模型的正演计算时间, T-EFGM施加第三类界条件而增加了边界计算时间, 由于边界计算量与电极数、波数个数、边界段数、高斯点数等成正比, 因此在多电极供电观测条件下往往耗费较多计算时间.EFG-FE中由于FEM单元数很少, 且FEM计算效率高, 使得FEM计算时间很少, 但EFG-FE由于增加了FEM区域的节点需要花费更多方程求解时间, 然而EFG-FE中仅使用少量的FEM节点, 因此不会明显增加方程求解时间.根据表 2, 虽然EFG-FE需要花费更多的方程求解时间, 但T-EFGM需要较多的边界计算时间, 使得总体上EFG-FE比EFGM花费计算时间少.表明本文的EFG-FE相比于T-EFGM, 通过使用快速延伸的FEM网格扩展计算域, 仅少量增加节点的情况下使得第一类边界条件得到满足, 节省了边界计算时间, 从而提高计算效率.考虑到两者精度基本一致, 表明本文的EFG-FE具有更好的模拟性能.需要指出的是, 外围FEM单元剖分需适当选择, 若选择过于密集的FEM单元剖分, 可能造成FEM节点数量过多而导致EFG-FE方程求解时间明显增加, 反而使得EFG-FE比T-EFGM计算效率低.因此, EFG-FE中应使用合适的的网格剖分方式尽可能发挥FEM的高效性, 从而提高EFG-FE的模拟性能.
5 结论本文利用RPIM插值构造无单元法形函数, 在外围区域将EFGM与FEM直接耦合, 提出直流电阻率的无单元Galerkin-有限单元耦合法(EFG-FE).首先, EFG-FE实现了仅需节点离散的EFGM正演模拟, 使得任意复杂电模型可方便地使用任意分布的节点离散覆盖, 而无需网格剖分, 可根据模型和模拟需要任意加密或稀释节点, 利于使节点合理分布, 从而有助于提高模拟精度和控制计算成本.其次, EFG-FE利用FEM的数值稳定性和高效性, 基于矩形单元, 采用快速扩展的单元网格覆盖远离场源的计算域外围区域, 使用少量的FEM节点和单元网格的情况下实现计算域大范围扩大, 使得满足第一类边界条件, 且不造成EFG-FE计算成本大幅增加.最后, 通过多电极供电和观测的复杂模型视电阻率观测正演模拟对比, 结果表明选择适当的外围有限单元法区域及单元剖分, EFG-FE不仅可获得高精度的模拟结果, 且相比于采用第三类边界条件的EFGM计算效率高, 具有更好的模拟性能.虽然EFG-FE相比于传统的数值模拟方法(如FEM)计算效率仍然较低, 但本文的研究可为后续研究计算成本更低的高精度无单元模拟方法提供研究思路, 以及为直流电阻率正演方法的选择提供选项.
致谢特别感谢4位审稿专家中肯的修改意见，使得本文内容更加完善.
Belytschko T, Lu Y Y, Gu L. 1994a. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 37(1): 229-256. |
Belytschko T, Lu Y Y, Gu L. 1994b. Fracture and crack growth by element free Galerkin methods. Modelling and Simulation in Materials Science and Engineering, 2(3A): 519-534. DOI:10.1088/0965-0393/2/3A/007 |
Belytschko T, Organ D, Krongauz Y. 1995. A coupled finite element-element-free galerkin method. Computational Mechanics, 17(3): 186-195. DOI:10.1007/BF00364080 |
Dolbow J, Belytschko T. 1999. Numerical integration of the Galerkin weak form in meshfree methods. Computational Mechanics, 23(3): 219-230. DOI:10.1007/s004660050403 |
Feng D S, Wang H H, Dai Q W. 2013. Forward simulation of ground penetrating radar based on the element-free Galerkin method. Chinese Journal of Geophysics, 56(1): 298-308. DOI:10.6038/cjg20130131 |
Gu Y T, Liu G R. 2005. Meshless methods coupled with other numerical methods. Tsinghua Science & Technology, 10(1): 8-15. |
Hajiazizi M, Bastan P. 2014. The elastoplastic analysis of a tunnel using the EFG method:A comparison of the EFGM with FEM and FDM. Applied Mathematics and Computation, 234: 82-113. DOI:10.1016/j.amc.2014.02.024 |
Hadinia M, Jafari R. 2015. An element-free Galerkin forward solver for the complete-electrode model in electrical impedance tomography. Flow Measurement and Instrumentation, 45: 68-74. DOI:10.1016/j.flowmeasinst.2015.04.011 |
Hadinia M, Jafari R, Soleimani M. 2016. EIT image reconstruction based on a hybrid FE-EFG forward method and the complete-electrode model. Physiological Measurement, 37(6): 863-878. DOI:10.1088/0967-3334/37/6/863 |
Hosseini S, Malekan M, Pitangueira R, et al. 2017. imposition of dirichlet boundary conditions in element free galerkin method through an object-oriented implementation. Latin American Journal of Solids and Structures, 14(6): 1017-1039. |
Ji Y J, Huang T Z, Huang W Y, et al. 2016. 2D anisotropic magnetotelluric numerical simulation using meshfree method under undulating terrain. Chinese Journal of Geophysics, 59(12): 4483-4493. DOI:10.6038/cjg20161211 |
Jia L, Huang Q Q, Yin Z P. 2007. A new and better EFGM-FEM direct coupling method. Journal of Northwestern Polytechnical University, 25(3): 337-341. |
Jia X F, Hu T Y, Wang R Q. 2006. Wave equation modeling and imaging by using element-free method. Progress in Geophysics, 21(1): 11-17. |
Kumar S, Singh I V, Mishra B K. 2014. A multigrid coupled (FE-EFG) approach to simulate fatigue crack growth in heterogeneous materials. Theoretical and Applied Fracture Mechanics, 72: 121-135. DOI:10.1016/j.tafmec.2014.03.005 |
Lancaster P, Salkauskas K. 1981. Surfaces generated by moving least squares methods. Mathematics of Computation, 37(155): 141-158. DOI:10.1090/S0025-5718-1981-0616367-1 |
Li D M, Bai F N, Cheng Y M, et al. 2012. A novel complex variable element-free Galerkin method for two-dimensional large deformation problems. Computer Methods in Applied Mechanics and Engineering, 233/236: 1-10. DOI:10.1016/j.cma.2012.03.015 |
Li J J, Yan J B. 2015. Magnetotelluric two-dimensional forward by finite element-radial point interpolation method. The Chinese Journal of Nonferrous Metals, 25(5): 1314-1324. |
Liu G R, Gu Y T. 2001. A point interpolation method for two-dimensional solids. International Journal for Numerical Methods in Engineering, 50(4): 937-951. DOI:10.1002/(ISSN)1097-0207 |
Liu T X, Liu G, Xie Q, et al. 2006. An EFG-FE coupling method for microscale adhesive contacts. Journal of Tribology, 128(1): 40-48. DOI:10.1115/1.2114931 |
Liu Z H, Lu M, Li Q, et al. 2016. A direct coupling method of meshless local petrov-galerkin (MLPG) and finite element method (FEM). International Journal of Applied Electromagnetics and Mechanics, 51(1): 51-59. DOI:10.3233/JAE-150161 |
Ma C Y, Liu J X, Liu H F, et al. 2017. A global weak form element free method for direct current resistivity forward simulation. Chinese Journal of Geophysics, 60(2): 801-809. DOI:10.6038/cjg20170230 |
Pan K J, Tang J T, Hu H L, et al. 2012. Extrapolation cascadic multigrid method for 2.5D direct current resistivity modeling. Chinese Journal of Geophysics, 55(8): 2769-2778. DOI:10.6038/j.issn.0001-5733.2012.08.028 |
Pathak H, Singh A, Singh I V, et al. 2015. Three-dimensional stochastic quasi-static fatigue crack growth simulations using coupled FE-EFG approach. Computers & Structures, 160: 1-19. |
Peng M J, Li D M, Cheng Y M. 2011. The complex variable element-free Galerkin (CVEFG) method for elasto-plasticity problems. Engineering Structures, 33(1): 127-135. DOI:10.1016/j.engstruct.2010.09.025 |
Ruan B Y. 2001. 2-D electrical modeling due to a current point by FEM with variation of conductivity within each triangular element. Guangxi Sciences, 8(1): 1-3. |
Wang Y Y. 2007. Study of element-free galerkin method in the seismic forward modeling. Progress in Geophysics, 22(5): 1539-1544. |
Xu S Z. 1994. Finite Element Method in Geophysics. Beijing: Beijing Science Press.
|
Yan J B, Li J J. 2014. Magnetotelluric forward calculation by meshless method. Journal of Central South University (Science and Technology), 45(10): 3513-3520. |
Zhang Y, Wang J G, Zhang B Y. 2008. Coupled FEM and meshless radial point interpolation method. Journal of Tsinghua University (Science and Technology), 48(6): 951-954. |
Zhdanov M S, Cox L H, Endo M. 2015. Large-scale 3D inversion of airborne electromagnetic data based on the hybrid IE-FE method and the moving sensitivity domain approach. //International Workshop and Gravity, Electrical & Magnetic Methods and their Applications. Chenghu, China, 2015: 303-306.
|
冯德山, 王洪华, 戴前伟. 2013. 基于无单元Galerkin法探地雷达正演模拟. 地球物理学学报, 56(1): 298-308. DOI:10.6038/cjg20130131 |
嵇艳鞠, 黄廷哲, 黄婉玉, 等. 2016. 起伏地形下各向异性的2D大地电磁无网格法数值模拟. 地球物理学报, 59(12): 4483-4493. DOI:10.6038/cjg20161211 |
贾亮, 黄其青, 殷之平. 2007. 无网格-有限元直接耦合法. 西北工业大学学报, 25(3): 337-341. |
贾晓峰, 胡天跃, 王润秋. 2006. 无单元法用于地震波波动方程模拟与成像. 地球物理学进展, 21(1): 11-17. |
李俊杰, 严家斌. 2015. 大地电磁二维正演中的有限元-径向基点插值法. 中国有色金属学报, 25(5): 1314-1324. |
麻昌英, 柳建新, 刘海飞, 等. 2017. 基于全局弱式无单元法直流电阻率正演模拟. 地球物理学学报, 60(2): 801-809. DOI:10.6038/cjg20170230 |
潘克家, 汤井田, 胡宏伶, 等. 2012. 直流电阻率法2. 5维正演的外推瀑布式多重网格法. 地球物理学报, 55(8): 2769-2778. DOI:10.6038/j.issn.0001-5733.2012.08.028 |
阮百尧. 2001. 三角单元部分电导率分块连续变化点源二维电场有限元数值模拟. 广西科学, 8(1): 1-3. |
王月英. 2007. 地震波正演模拟中无单元Galerkin法初探. 地球物理学进展, 22(5): 1539-1544. |
徐世浙. 1994. 地球物理中的有限单元法. 北京: 科学出版社.
|
严家斌, 李俊杰. 2014. 无网格法在大地电磁正演计算中的应用. 中南大学学报(自然科学版), 45(10): 3513-3520. |
张琰, 王建国, 张丙印. 2008. 径向基点插值无网格法与有限元耦合法. 清华大学学报(自然科学版), 48(6): 951-954. |