CN110008593B - 一种直流电阻率无单元单位分解法的单位分解方法 - Google Patents

Abstract

本发明提供了一种直流电阻率无单元单位分解法的单位分解方法，包括以下步骤，建立二维地电模型；利用单位分解积分将全局域积分转化为节点局部域积分；单位分解函数ψ_k(x)构造；直流电阻率无单元单位分解法计算，其特征在于，单位分解函数ψ_k(x)构造为：假设计算域Ω内N个节点的节点局部支持域Ω_k(k＝1,2,…,N)包含计算点x的节点数为n_k(n_k≥1)，这n_k个节点构成x的单位分函数域，即x的单位分函数域内包含n_k个节点，定义距离函数。本发明能够基于任意节点分布离散模型，利用一种距离函数作为单位分解函数，计算简单，对单位分解函数支持域节点数没有要求，使得直流电阻率无单元单位分解法求解稳定性和计算效率提高。

Description

一种直流电阻率无单元单位分解法的单位分解方法

技术领域

本发明涉及一种勘探地球物理领域的直流电阻率正演方法，特别涉及复杂地电模型的高精度、高灵活性和高适应性的基于单位分解积分的直流电阻率无单元正演方法。

背景技术

直流电阻率法勘探是地球物理勘探中的一种重要方法，被广泛应用于固体矿产资源勘探、水文地质勘察、环境治理与监测、工程地球物理勘查等领域。测量的视电阻率与地下介质的电阻率有着直接的关系，通过人工向地下供电，在地表或者井中观测视电阻率可以对地下电阻率异常体的分布进行判断。随着直流电阻率法勘探技术的发展，人们对复杂地电模型的高精度、高适应性和灵活性的正演方法的需求日益增长，无单元法是一种新兴的数值模拟方法(Belytschko,et al.,1994；Hadinia and Jafari,2015)，其仅需局部支持域内的节点信息，不依赖网格链接信息，摆脱了网格的约束而具有高灵活性和适应性的特点，同时由于采用高精度的插值方法其具有高精度的特点，因此，被广泛研究，本发明人也一直致力于这方面的研究，2017年在地球物理学报上与他人一起发表的《基于全局弱式无单元法直流电阻率正演模拟》，已在直流电阻率正演模拟中得到应用。

在进一步研究中发现，由于采用基于局部弱式的无单元法直流电阻率正演模拟的过程中，需要在计算域上进行背景单元剖分，仍然对单元有依赖，一定程度上损失了无单元的性质，同时，背景单元分布情况直接影响模拟精度，很难采用统一标准的背景单元剖分方法，获得针对不同模型和节点分布的最优化背景单元分布，往往采用密集的背景单元剖分适应不同模型和节点分布，以保证模拟精度。因此背景单元的使用，降低了无单元法的无单元属性，降低了无单元法对任意模型和任意节点分布的适应性。

为解决这一难题，本发明人又参与提出了《一种基于单位分解积分的直流电阻率无单元正演方法》CN201810444952.0，根据二维地电模型中介质电阻率的分布、电阻率异常体的几何形态和地形起伏形态建立计算域，利用不规则分布的节点离散地电模型；利用单位分解积分，将直流电阻率法的全局计算域积分方程转化为节点局部域积分方程；构造节点局部积分域，在每一个节点局部积分域内采用高斯积分计算积分式，获得计算域离散方程；求解离散方程获得节点电场场值，计算获得观测点的视电阻率参数。能够基于任意节点分布离散模型，利用单位分解积分将全局域积分转化为节点局部域积分，不再需要节点间连接信息和单元，摆脱了网格束缚，对任意复杂的地电模型适应性和灵活性强。

由于，在无单元单位分解法进行直流电阻率正演模拟的过程中，需要构造单位分解函数，单位分解函数和无单元法形函数为不同的函数集，两者本质上又相互独立，为了方便通常直接采用无单元法形函数本身作为单位分解函数。虽然该方案使得单位分解函数构造方便直接，但无单元法形函数计算对节点分布及支持域节点数要求相对较高，直接导致单位分解函数构造稳定性较差，为了提高稳定性需要扩大节点局部积分域，又会导致无单元法计算成本提高等问题。

因此，有必要设计一种可保证数值稳定性和降低计算效率的单位分解函数用于直流电阻率无单元单位分解法正演。

发明内容

本发明旨在克服现有技术的不足，在CN201810444952.0的基础上，进一步提供一种直流电阻率无单元单位分解法的单位分解方法，利用该种单位分解函数实现直流电阻率的无单元法单位分解积分，降低对单位分解函数域节点数和分布要求，简化单位分解函数计算，使得直流电阻率无单元单位分解法求解稳定性和计算效率提高。

为实现上述目的，本发明提供了一种直流电阻率无单元单位分解法的单位分解方法，包括以下步骤：

步骤(1)、建立二维地电模型：

首先，根据二维地电模型中介质电阻率的分布、电阻率异常体的几何形态和地形起伏形态建立计算域，并设置好电极位置、观测装置和观测点位置；在计算域中将二维地电模型采用一组任意分布的节点进行离散，根据电阻率异常体的位置、几何形态和地形起伏形态以及电极位置布置节点；并根据正演模拟需求，在局部域任意加密节点，在场值变化不大或者远离场源电性不变的区域使用稀疏的节点分布；

步骤(2)、利用单位分解积分将全局域积分转化为节点局部域积分：

步骤(2.1)、采用第三类边界条件的2.5维直流电阻率边值问题对应的变分问题为式1)：

其中，Γ为边界符号，Γ_T为截断边界；σ为介质电导率；Ω为计算域，U为波数域电位，λ为波数；I₀为电流；δ为Kronecker delta函数；x为计算域内任意一点；A为场源点位置；

为梯度运算符；r_A为点源与截断边界上任意一点的直线距离，n为截断边界外法线单位向量，cos(r_A,n)为r_A与n的夹角余弦；K₀、K₁分别为第二类零阶、一阶修正贝塞尔函数；δ为变分符号；

步骤(2.2)、计算域Ω内设有N个节点，为每个节点k构造四边形的节点局部域积分Ω_k，k表示节点编号，k＝1,2,…,N，x_k为节点k的坐标：

步骤(2.3)、在每个节点k构造的节点局部域Ω_k内，建立满足如下条件的单位分解函数ψ_k(x)：

0≤ψ_k(x)≤1,x∈Ω_k 3)；

其中，Ω_k为节点k局部积分域；式2)～4)为单位分解性质，任意满足这些性质的函数均可作为单位分解函数；

步骤(2.4)、f(x)是定义在计算域Ω上的可积函数，则有如下等式成立：

式5)为单位分解积分；

步骤(2.5)、利用单位分解函数和单位分解积分将式1)中的第一式转化为节点局部域积分：

其中，

为截断边界Γ_T与节点k局部积分域Ω_k相交部分的子边界；

步骤(2.6)、式6)可首先通过变分原理获得对应的变分方程7)：

对于单个节点k，在其节点局部积分域Ω_k则有：

式7)可通过标准的无单元单位分解法进行离散，获得单个节点k对应的直流电阻率无单元单位分解法离散方程：

其中，设Ω_k内任意一点x的支持域内有n个节点，

和

为节点k局部积分域Ω_k的子系数矩阵，

和

均为n×n矩阵，F_k为节点k局部积分域Ω_k的子离散方程式9)的右端项，F_k ^T＝[f₁ f₂ … f_n]，U^T＝[U₁ U₂ … U_n]为支持域内n个节点对应的波数域电位向量，

式9)中：

式中，Φ^T＝[φ₁(x) φ₂(x) … φ_n(x)]为支持域内n个节点对应的形函数向量，Φ^T右上角的T表示转置；

步骤(3)、单位分解函数ψ_k(x)构造：

要计算10)～12)式，必须确定单位分解函数ψ_k(x)。假设计算域Ω内N个节点的节点局部支持域Ω_k(k＝1,2,…,N)包含计算点x的节点数为n_k(n_k≥1)，这n_k个节点构成x的单位分函数域，即x的单位分函数域内包含n_k个节点，定义如下距离函数：

作为单位分解函数ψ_k(x)；

式中，r_k为x与x_k之间的距离，r_i(i＝1,2,…,n_k)为x到其单位分函数域内n_k个节点的距离。显然，该距离函数满足单位分解性质式2)～4)，且计算简单，其对单位分解函数支持域节点数没有要求；

步骤(4)、直流电阻率无单元单位分解法计算：

通过标准的无单元单位分解法和利用单位分解函数式13)，获得单个节点k对应的直流电阻率无单元单位分解法离散方程式9)，将所有的节点对应的直流电阻率无单元单位分解法离散方程式9)的子系数矩阵

和

按照节点整体编号组装起来，获得总系数矩阵K：

同时，将所有的节点对应的直流电阻率无单元单位分解法离散方程式9)的子离散方程右端项F_k组装起来，获得离散方程右端项F：

结合Kronecker delta函数δ(A)的积分性质，当无单元法形函数具有Kroneckerdelta性质时，离散方程的右端项F的元素表达式为式16)，当无单元法形函数不具有Kronecker delta性质时，离散方程的右端项F的元素表达式为式17)：

其中，f_i表示右端项F的第i个元素。

最终获得计算域对应的离散方程：

KU＝F 18)。

本发明在CN201810444952.0的基础上创造性地提出了一套单位分解函数ψ_k(x)，采用式13)，即R_k(x)作为单位分解函数ψ_k(x)，不同于CN201810444952.0中直接将形函数Φ(x)作为单位分解函数ψ_k(x)。

在单位分解积分中用到两个函数集：形函数Φ(x)和单位分解函数Ψ(x)。单位分解函数和形函数本质上是两个不同的函数系，两者相互独立，任意具有单位分解性质的函数，即满足2)～4)式，均可作为单位分解函数，由单位分解积分可知，单位分解函数不会对积分计算精度产生影响。由于形函数具有单位分解性质，在求解域Ω内任意一点的形函数之和为1，即：

目前直流电阻率无单元单位分解法中通常将无单元法构造的形函数作为单位分解函数，即

ψ_k(x)＝φ_k(x),k＝1,2,…,N 20)

其中φ_k(x)(k＝1,2,…,N)为计算域任意一点x对节点k的形函数。然而，节点局部积分域Ω_k(k＝1,2,…,N)和x支持域本质上是相互独立的。由于无单元法形函数需要满足诸多性质，如紧支性、连续性、再生性等，当采用无单元法形函数作为单位分解函数时，需对节点局部积分域Ω_k包含x的n个节点构造无单元法形函数，为了保证形函数构造稳定，通常需要较多的支持域节点数，否则形函数可能构造失败，在节点任意分布情况下该种情况较容易发生。克服该问题的方法为使用较大的节点局部积分域Ω_k，但较大节点局部积分域Ω_k意味着需要更多的积分点计算积分以保证积分精度，增加了计算成本。另外，可将节点局部积分域Ω_k包含x的n个节点作为x的支持域，当采用无单元法形函数作为单位分解函数时，可避免单位分解函数的计算，但也存在上述问题，且x支持域依赖于节点局部积分域，使得x支持域可控性降低，通常会使支持域内节点数较多，导致支持域构造和形函数计算成本增加。

本发明提出的距离函数满足单位分解性质，计算简单，其对单位分解函数支持域节点数没有要求，任何一种情况下单位分解函数均能顺利构造，进而可将节点局部积分域减小，可解决将无单元法形函数作为单位分解函数时存在的问题，使得直流电阻率无单元单位分解法求解稳定性和计算效率提高，具有一定的实用价值和应用前景。

附图说明

图1为本发明中二维地电模型及其节点分布示意图；

图2为采用本发明的直流电阻率无单元单位分解法单位分解函数、独立构造单位分解函数域和支持域的模拟结果断面图；

图3为采用无单元法形函数作为单位分解函数、独立构造单位分解函数域和支持域的模拟结果断面图；

图4为采用无单元法形函数作为单位分解函数、单位分解函数域和支持域一致的模拟结果断面图；

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。

参见图1，本发明涉及的一种直流电阻率无单元单位分解法的单位分解方法，包括以下步骤：

步骤(1)、建立二维地电模型：

首先，根据二维地电模型中介质电阻率的分布、电阻率异常体的几何形态和地形起伏形态建立计算域1，并设置好电极位置、观测装置和观测点位置；在计算域中将二维地电模型采用一组任意分布的节点(图1中的附图标记6)进行离散，根据电阻率异常体的位置、几何形态和地形起伏形态以及电极位置布置节点；并根据正演模拟需求，在局部域任意加密节点，如地形变化区域、电极附近、电阻率异常体及其周围区域中的一个或多个区域加密节点，在场值变化不大或者远离场源电性不变的区域使用稀疏的节点分布。

步骤(2)、利用单位分解积分将全局域积分转化为节点局部域积分：

步骤(2.1)、采用第三类边界条件的2.5维直流电阻率边值问题对应的变分问题为式1)：

其中，Γ为边界符号，Γ_T为截断边界；σ为介质电导率；Ω为计算域，U为波数域电位，λ为波数；I₀为电流；δ为Kronecker delta函数；x为计算域内任意一点；A为场源点位置；

为梯度运算符；r_A为点源与截断边界上任意一点的直线距离，n为截断边界外法线单位向量，cos(r_A,n)为r_A与n的夹角余弦；K₀、K₁分别为第二类零阶、一阶修正贝塞尔函数；δ为变分符号；式1)中第一式为一个以U为函数变量的函数，第二式表示令函数F(U)的变分为0；

步骤(2.2)、计算域Ω内设有N个节点，为每个节点k构造四边形的节点局部域积分Ω_k，k表示节点编号，k＝1,2,…,N，x_k为节点k的坐标：根据计算域内的节点分布情况，构造节点局部积分域Ω_k：为了方便布置高斯点进行积分计算，采用矩形节点局部积分域；对于矩形域在二维情况下，使用x和z两个方向的尺寸参数qx和qz来确定范围，以节点k为中心，在全局节点内搜索节点k的节点局部积分域Ω_k包含的节点数N_q，通常N_q设置为6～24。

步骤(2.3)、在每个节点k构造的节点局部域Ω_k内，建立满足如下条件的单位分解函数ψ_k(x)：

0≤ψ_k(x)≤1,x∈Ω_k 3)；

其中，k表示节点编号，k＝1,2,…,N；Ω_k为节点k局部积分域；式2)～4)为单位分解性质；

步骤(2.4)、假设f(x)是定义在计算域Ω上的可积函数，则有如下等式成立：

式5)为单位分解积分；

步骤(2.5)、利用单位分解函数和单位分解积分将式1)中的第一式转化为节点局部域积分：

其中，

为截断边界Γ_T与节点k局部积分域Ω_k相交部分的子边界；

步骤(2.6)、对步骤(2)中式6)的电位函数U进行变分：

对式7)展开变分并进行整理可得：

对于单个节点k则有：

步骤(3)、计算单个节点直流电阻率单位分解法子离散方程：

为节点局部积分域内任意一点x构造支持域，假设支持域内有n个节点，使用该n个节点构造形函数(φ_i(x),i＝1,2,…,n)对x点处场值进行插值：

式中，Φ^T＝[φ₁(x) φ₂(x) … φ_n(x)]为形函数向量，U^T＝[U₁ U₂ … U_n]为支持域内节点的波数域电位向量；将式10)代入式9)可得：

对式11)进行整理可得：

式12)中，如Φ^T、

U^T的右上角的T表示转置，黑体符号Φ表示形函数向量，黑体符号U表示支持域节点的波数域电位向量；

对式12)进行整理可得：

式13)中：

步骤(4)、单位分解函数ψ_k(x)构造：

要计算14)～16)式，必须确定单位分解函数ψ_k(x)。假设计算域Ω内N个节点的节点局部支持域Ω_k(k＝1,2,…,N)包含计算点x的节点数为n_k(n_k≥1)，这n_k个节点构成x的单位分函数域，即x的单位分函数域内包含n_k个节点，定义如下距离函数：

作为单位分解函数ψ_k(x)；

式中，r_k为x与x_k之间的距离，r_i(i＝1,2,…,n_k)为x到其单位分函数域内n_k个节点的距离。显然，该距离函数满足单位分解性质式2)～4)，且计算简单，其对单位分解函数支持域节点数没有要求；

步骤(5)、直流电阻率无单元单位分解法计算：

由于δU是任意的，故式13)成立的条件为：

其中，

和

为节点k局部积分域Ω_k的子系数矩阵，F_k为节点k局部积分域Ω_k的子离散方程式18)的右端项。

假设节点k局部积分域Ω_k中采用N_g个高斯点x_g＝(x_g,z_g)进行积分计算，高斯点对应的权重w_g和雅可比值J_g，g＝1,2,L,N_g；高斯点支持域内包含n个节点，则有：

其中，

B＝σλ²Φ(x_g)Φ^T(x_g) 22)；

将式21)～23)分别展开可得：

式24)～26)中φ_i为高斯点支持域内第i个节点对应的形函数，

表示形函数的一阶导数，i＝1,2,L,n；对于单个高斯点有：

和

为节点k局部积分域Ω_k内第g个高斯点对应的子系数矩阵。此时，高斯点对应的子离散方程表达式为：

其中，U_kg为节点k局部积分域Ω_k内第g个高斯点支持域节点的电位场值构成的电位向量，F_kg为节点k局部积分域Ω_k内第g个高斯点支持域节点构成的子离散方程的右端项；

式27)和式28)为高斯点对应的子系数矩阵，两个矩阵元素分别为：

式30)和式31)中，i和j表示高斯点支持域Ω_q中的节点局部编号，i＝1,2,L,n，j＝1,2,L,n；φ_i(x_g)表示高斯点x_g支持域内第i个节点的形函数，φ_j(x_g)表示高斯点x_g支持域内第j个节点的形函数，见φ_i(x),i＝1,2,…,n，只不过φ_i(x)中的x表示计算域中的任意一点或者表示节点局部积分域内的任意一点，而x_g表示高斯点，x包含了x_g。

将步骤(4)的距离函数式17)作为单位分解函数，即ψ_k(x_g)＝R_k(x_g)；此时式30)和式31)分别写为：

将所有的节点局部积分域内所有高斯点的子系数矩阵按照节点整体编号组装起来，获得总系数矩阵K：

将所有的节点局部积分域内所有高斯点的子离散方程右端项F_kg组装起来，获得离散方程右端项F：

结合Kronecker delta函数δ(A)的积分性质，当无单元法形函数具有Kroneckerdelta性质时，离散方程的右端项F的元素表达式为式37)，当无单元法形函数不具有Kronecker delta性质时，离散方程的右端项F的元素表达式为式38)：

式37)和38)中，f_i表示右端项F的第i个元素，详见式36)。

最终获得计算域对应的离散方程：

KU＝F 39)。

步骤(6)、计算视电阻率：

对步骤5中获得的计算域对应的离散方程进行求解，获得节点电场场值，再根据观测装置观测到的相关参数计算获得观测点的视电阻率参数。

以下为利用本发明一种基于单位分解积分的直流电阻率无单元正演方法计算一个复杂二维地电模型的高密度电法温纳装置视电阻率观测的实例。

如图1所示，建立一个宽200m(X:-100～100m)、高80m(Z:0～80m)的二维地电模型，采用任意不规则节点离散模型，该二维地电模型中设置11个电阻率异常体，第一层电阻率ρ₁＝40Ω·m，第二层电阻率ρ₂＝50Ω·m，第三层电阻率ρ₃＝60Ω·m，第四层电阻率ρ₄＝30Ω·m，第五层电阻率ρ₅＝50Ω·m，第六层电阻率ρ₆＝60Ω·m，第七层电阻率ρ₇＝70Ω·m，第八层电阻率ρ₈＝80Ω·m，第九层电阻率ρ₉＝10Ω·m，第十层电阻率ρ₁₀＝5Ω·m，第十一层电阻率ρ₁₁＝100Ω·m。在地表X:-58～58m范围内等间隔2m布置59根供电和观测电极，对二维地电模型进行高密度电法温纳装置视电阻率观测正演模拟。分别采用采用无单元法形函数作为单位分解函数和本发明的方法进行直流电阻率无单元单位分解法正演模拟。

首先，采用本发明方法进行正演模拟，即将距离函数作为单位分解函数，且独立构造单位分解函数域和支持域，即高斯点的支持域节点与计算单位分解函数是不一致的，支持域内采用距离其最近的8个节点，模拟结果如图2所示。其次，采用无单元形函数作为单位分解函数，且独立构造单位分解函数域和支持域，支持域内采用距离其最近的8个节点，模拟结果如图3所示。最后，采用无单元法形函数本身作为单位分解函数，将节点局部积分域包含高斯点的n个节点作为高斯点的支持域，即不独立构造单位分解函数域和支持域，高斯点的支持域即为节点局部积分域包含高斯点的n个节点，此时高斯点的支持域节点与计算单位分解函数是一致的，模拟结果如图4所示。将图2、图3和图4对比分析可知，两种方法的模拟结果基本一致，有效表明了本发明方法对直流电阻率正演模拟的有效性。

表1为上述是那种方法的计算时间，即不同单位分解函数和不同支持域类型计算时间。由表1分析可知，采用距离函数作为单位分解函数并独立构造支持域时无单元单位分解法花费的计算时间最少，而采用无单元法形函数本身作为单位分解函数时，无论支持域是否独立构造，均明显花费更多计算时间。其主要原因是为了保证形函数构造成功，使用了较大的节点局部积分域，以使得用于构造形函数的节点较多，避免节点过少(如少于4个节点)或节点分布极不合理时，容易导致形函数构造失败的情况，而为了保证积分精度，积分域越大需使用更多高斯点计算，使得计算成本增加。以上分析结果表明，本发明的方法对单位分解函数支持域内节点数要求低，可有效提高直流电阻率无单元单位分解法的计算效率，将其作为单位分解函数可有效克服将无单元法形函数作为单位分解函数时存在的问题。

表1不同单位分解函数构造方式计算时间

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种直流电阻率无单元单位分解法的单位分解方法，包括以下步骤，步骤(1)、建立二维地电模型；步骤(2)、利用单位分解积分将全局域积分转化为节点局部域积分；步骤(3)、单位分解函数ψ_k(x)构造；步骤(4)、直流电阻率无单元单位分解法计算，其特征在于，步骤(3)、单位分解函数ψ_k(x)构造为：假设计算域Ω内N个节点的节点局部支持域Ω_k(k＝1,2,…,N)包含计算点x的节点数为n_k(n_k≥1)，这n_k个节点构成x的单位分函数域，即x的单位分函数域内包含n_k个节点，定义如下距离函数：

作为单位分解函数ψ_k(x)；

式中，r_k为x与x_k之间的距离，r_i(i＝1,2,…,n_k)为x到其单位分函数域内n_k个节点的距离。

2.根据权利要求1所述的一种直流电阻率无单元单位分解法的单位分解方法，其特征在于，包括以下步骤：

步骤(1)、建立二维地电模型：

首先，根据二维地电模型中介质电阻率的分布、电阻率异常体的几何形态和地形起伏形态建立计算域，并设置好电极位置、观测装置和观测点位置；在计算域中将二维地电模型采用一组任意分布的节点进行离散，根据电阻率异常体的位置、几何形态和地形起伏形态以及电极位置布置节点；并根据正演模拟需求，在局部域任意加密节点，在场值变化不大或者远离场源电性不变的区域使用稀疏的节点分布；

步骤(2)、利用单位分解积分将全局域积分转化为节点局部域积分：

步骤(2.1)、采用第三类边界条件的2.5维直流电阻率边值问题对应的变分问题为式1)：

其中，Γ为边界符号，Γ_T为截断边界；σ为介质电导率；Ω为计算域，U为波数域电位，λ为波数；I₀为电流；δ为Kronecker delta函数；x为计算域内任意一点；A为场源点位置；

为梯度运算符；r_A为点源与截断边界上任意一点的直线距离，n为截断边界外法线单位向量，cos(r_A,n)为r_A与n的夹角余弦；K₀、K₁分别为第二类零阶、一阶修正贝塞尔函数；δ为变分符号；

步骤(2.2)、计算域Ω内设有N个节点，为每个节点k构造四边形的节点局部域积分Ω_k，k表示节点编号，k＝1,2,…,N，x_k为节点k的坐标：

步骤(2.3)、在每个节点k构造的节点局部域Ω_k内，建立满足如下条件的单位分解函数ψ_k(x)：

0≤ψ_k(x)≤1,x∈Ω_k 3)；

其中，Ω_k为节点k局部积分域；式2)～4)为单位分解性质，任意满足这些性质的函数均可作为单位分解函数；

步骤(2.4)、f(x)是定义在计算域Ω上的可积函数，则有如下等式成立：

式5)为单位分解积分；

步骤(2.5)、利用单位分解函数和单位分解积分将式1)中的第一式转化为节点局部域积分：

其中，

为截断边界Γ_T与节点k局部积分域Ω_k相交部分的子边界；

步骤(2.6)、式6)可首先通过变分原理获得对应的变分方程7)：

对于单个节点k，在其节点局部积分域Ω_k则有：

式7)可通过标准的无单元单位分解法进行离散，获得单个节点k对应的直流电阻率无单元单位分解法离散方程：

其中，设Ω_k内任意一点x的支持域内有n个节点，

和

为节点k局部积分域Ω_k的子系数矩阵，

和

均为n×n矩阵，F_k为节点k局部积分域Ω_k的子离散方程式9)的右端项，F_k ^T＝[f₁ f₂ … f_n]，U^T＝[U₁ U₂ … U_n]为支持域内n个节点对应的波数域电位向量，

式9)中：

式中，Φ^T＝[φ₁(x) φ₂(x) … φ_n(x)]为支持域内n个节点对应的形函数向量，Φ^T右上角的T表示转置；

步骤(3)、单位分解函数ψ_k(x)构造；

步骤(4)、直流电阻率无单元单位分解法计算：

通过标准的无单元单位分解法和利用单位分解函数式13)，获得单个节点k对应的直流电阻率无单元单位分解法离散方程式9)，将所有的节点对应的直流电阻率无单元单位分解法离散方程式9)的子系数矩阵

和

按照节点整体编号组装起来，获得总系数矩阵K：

同时，将所有的节点对应的直流电阻率无单元单位分解法离散方程式9)的子离散方程右端项F_k组装起来，获得离散方程右端项F：

结合Kronecker delta函数δ(A)的积分性质，当无单元法形函数具有Kronecker delta性质时，离散方程的右端项F的元素表达式为式16)，当无单元法形函数不具有Kroneckerdelta性质时，离散方程的右端项F的元素表达式为式17)：

其中，f_i表示右端项F的第i个元素；

最终获得计算域对应的离散方程：

KU＝F 18)。

Images