CN108873084B - 一种基于单位分解积分的直流电阻率无单元正演方法 - Google Patents

Abstract

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

Description

一种基于单位分解积分的直流电阻率无单元正演方法

技术领域

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

背景技术

直流电阻率法勘探是地球物理勘探中的一种重要方法，被广泛应用于固体矿产资源勘探、水文地质勘察、环境治理与监测、工程地球物理勘查等领域。测量的视电阻率与地下介质的电阻率有着直接的关系，通过人工向地下供电，在地表或者井中观测视电阻率可以对地下电阻率异常体的分布进行判断。随着直流电阻率法勘探技术的发展，人们对复杂地电模型的高精度、高适应性和灵活性的正演方法的需求日益增长，无单元法是一种新兴的数值模拟方法(Belytschko,et al.,1994；Hadinia and Jafari,2015)，其仅需局部支持域内的节点信息，不依赖网格链接信息，摆脱了网格的约束而具有高灵活性和适应性的特点，同时由于采用高精度的插值方法其具有高精度的特点，被广泛研究，目前无单元法在直流电阻率正演模拟中已获得了应用，例如麻昌英和柳建新等人2017年在地球物理学报上发表的《基于全局弱式无单元法直流电阻率正演模拟》，但在采用基于局部弱式的无单元法直流电阻率正演模拟的过程中，需要在计算域上进行背景单元剖分，仍然对单元有依赖，一定程度上损失了无单元的性质，同时背景单元分布情况直接影响模拟精度，很难采用统一标准的背景单元剖分方法，获得针对不同模型和节点分布的最优化背景单元分布，往往采用密集的背景单元剖分适应不同模型和节点分布，以保证模拟精度。因此背景单元的使用，降低了无单元法的无单元属性，降低了无单元法对任意模型和任意节点分布的适应性。

因此，有必要设计一种进一步降低对单元的依赖、无需背景单元剖分的、更高适应性和数值稳定性的直流电阻率无单元正演方法。

发明内容

本发明旨在克服现有技术的不足，提供一种基于单位分解积分的直流电阻率无单元正演方法，利用单位分解积分将全局域积分转化为节点局部域积分，在节点局部域内计算积分，不再需要背景单元剖分，进一步降低对单元的依赖。

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

步骤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(x)：

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

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

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

式5)为单位分解积分；

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

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

步骤3、构造节点局部积分域，计算离散方程：

步骤3.1、根据计算域内的节点分布情况，构造节点局部积分域；

步骤3.2、在每一个节点局部积分域内采用高斯积分计算积分式，首先为每一个高斯点构造支持域，然后使用支持域内包含的节点构造形函数对高斯点处场值进行插值，形成高斯点对应的子离散方程；

步骤3.3、将计算域内所有高斯点对应的子离散方程按照节点整体编号组装起来，获得计算域对应的离散方程；

步骤4、计算视电阻率：

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

进一步的，对所述步骤2中式6)的电位函数U进行变分：

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

对于单个节点k则有：

进一步的，为节点局部积分域内任意一点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)进行整理可得：

式13)中：

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

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

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

其中，

A₁＝σ▽Φ(x_g)·[▽Φ(x_g)]^T 20)；

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

将式20)～22)分别展开可得：

式23)～25)中φ_i为高斯点支持域内第i个节点对应的形函数，i＝1,2,…,n；对于单个高斯点有：

此时，高斯点对应的子离散方程表达式为：

其中，U_kg为高斯点支持域节点的电位场值构成的电位向量，F_kg为高斯点支持域节点构成的子离散方程的右端项；

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

其中，i和j表示高斯点支持域Ω_q中的节点局部编号，i＝1,2,…,n，j＝1,2,…,n；

采用节点局部积分域包含高斯点x_g的节点作为x_g的支持域，则计算单位分解函数ψ_k(x_g)使用的节点与构造形函数φ_i(x_g)使用的节点相同，由于无单元法形函数满足单位分解性质，即

因此采用无单元法形函数作为单位分解函数，即ψ_k(x_g)＝φ_k(x_g)，其中φ_k(x_g)为当前计算节点k对应的无单元法形函数；此时29)式和30)式分别写为：

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

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

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

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

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

KU＝F 38)。

相比于现有技术，本发明具有以下有益效果：

(1)、本发明的一种基于单位分解积分的直流电阻率无单元正演方法，利用单位分解积分将全局域积分转化为节点局部域积分，在节点局部域内计算积分，不再需要背景单元剖分(不再需要背景单元计算积分)，进一步降低对单元的依赖，仅需要输入模型节点信息，使得直流电阻率正演更加简便，对任意复杂的模型和节点分布具有很强的适应性，且计算精度高。同时，解决了采用背景单元积分时存在的问题，如在近场源区域相邻背景单元剖分尺寸大小不同时导致较大的积分精度差异，这种差异可明显影响模拟精度。

(2)、本发明的一种基于单位分解积分的直流电阻率无单元正演方法，无需节点连接信息以及背景单元剖分，可以对二维任意几何形态、地形起伏、电阻率参数分布复杂的地电模型进行直流电阻率法正演计算，相比基于节点连接信息剖分网格的有限单元法和基于背景单元积分的无单元法具有更强的灵活性，具有一定的实用价值和应用前景。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为本发明一种基于单位分解积分的直流电阻率无单元正演方法的计算示意图；

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

图3为采用本发明的基于单位分解积分的直流电阻率无单元正演方法的模拟结果断面图；

图4为采用基于局部弱式的无单元法直流电阻率正演模拟的模拟结果断面图；

图5为采用基于非结构化三角形的有限单元法直流电阻率正演模拟的模拟结果断面图；

图6为基于局部弱式的无单元法直流电阻率正演模拟的背景单元；

图7为有限单元法直流电阻率正演模拟的三角形单元；

其中，1、计算域，2、节点k，3、节点k局部积分域，4、高斯点，5、高斯点支持域，6、节点。

具体实施方式

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

参见图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(x)：

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

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

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

式5)为单位分解积分；

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

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

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

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

对于单个节点k则有：

步骤3、构造节点局部积分域，计算离散方程：

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

步骤3.2、为节点局部积分域内任意一点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、(▽Φ)^T、U^T的右上角的T表示转置，黑体符号Φ表示形函数向量，黑体符号U表示支持域节点的波数域电位向量；

对式12)进行整理可得：

式13)中：

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

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

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

其中，

A₁＝σ▽Φ(x_g)·[▽Φ(x_g)]^T 20)；

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

将式20)～22)分别展开可得：

式23)～25)中φ_i为高斯点支持域内第i个节点对应的形函数，i＝1,2,…,n；对于单个高斯点有：

此时，高斯点对应的子离散方程表达式为：

其中，U_kg为高斯点支持域节点的电位场值构成的电位向量，F_kg为高斯点支持域节点构成的子离散方程的右端项；

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

式29)和式30)中，i和j表示高斯点支持域Ω_q中的节点局部编号，i＝1,2,…,n，j＝1,2,…,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。

采用节点局部积分域包含高斯点x_g的节点作为x_g的支持域，则计算单位分解函数ψ_k(x_g)使用的节点与构造形函数φ_i(x_g)使用的节点相同，由于无单元法形函数满足单位分解性质，即

因此采用无单元法形函数作为单位分解函数，即ψ_k(x_g)＝φ_k(x_g)，其中φ_k(x_g)为当前计算节点k对应的无单元法形函数；此时式29)和式30)分别写为：

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

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

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

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

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

KU＝F 38)。

步骤4、计算视电阻率：

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

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

如图2所示，建立一个宽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根供电和观测电极，对二维地电模型进行高密度电法温纳装置视电阻率观测正演模拟。分别采用本发明的方法、基于局部弱式的无单元法和基于非结构化三角形的有限单元法进行正演模拟，采用基于非结构化三角形的有限单元法模拟时，由于其对边界与供电点要求较大距离，故在上述计算域基础上需扩边至宽2000m(X:-1000～1000m)，高1000m(Z:0～1000m)，以降低边界对模拟结果的影响。如图3至图5所示，为分别采用上述三种方法模拟的视电阻率模拟结果断面图，该三种方法的模拟结果基本一致，仅仅在少数局部区域有细微的差别，有效表明了本发明方法对直流电阻率正演模拟的有效性。

图6为现有技术中基于全局弱式的无单元法剖分的背景单元；图7为现有技术中基于非结构化三角形的有限单元法剖分的三角形单元。图6和图7表明该两种方法均需要较为复杂的单元剖分，基于局部弱式的无单元法中背景单元与节点关系不大，但为了提高正演稳定性，需将背景单元依据节点的分布情况进行指导剖分，节点密集的区域对应密集的背景单元分布，节点稀疏区域对应稀疏的背景单元剖分。有限单元法中三角形单元根据节点的连接信息剖分，单元定点即为节点，非结构化网格的剖分为图形学中重要的一门学科，非专业人员实现任意节点分布的非结构化网格剖分较困难。本发明方法无需单元剖分，可实现对任意复杂模型和节点分布的直流电阻率正演模拟。

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

Claims (3)

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(x)：

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

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

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

式5)为单位分解积分；

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

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

步骤3、构造节点局部积分域，计算离散方程：

步骤3.1、根据计算域内的节点分布情况，构造节点局部积分域；

步骤3.2、在每一个节点局部积分域内采用高斯积分计算积分式，首先为每一个高斯点构造支持域，然后使用支持域内包含的节点构造形函数对高斯点处场值进行插值，形成高斯点对应的子离散方程；

步骤3.3、将计算域内所有高斯点对应的子离散方程按照节点整体编号组装起来，获得计算域对应的离散方程；

步骤4、计算视电阻率：

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

2.根据权利要求1所述的基于单位分解积分的直流电阻率无单元正演方法，其特征在于，对所述步骤2中式6)的电位函数U进行变分：

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

对于单个节点k则有：

3.根据权利要求2所述的基于单位分解积分的直流电阻率无单元正演方法，其特征在于，为节点局部积分域内任意一点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)进行整理可得：

式13)中：

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

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

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

其中，

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

将式20)～22)分别展开可得：

式23)～25)中φ_i为高斯点支持域内第i个节点对应的形函数，i＝1,2,…,n；对于单个高斯点有：

此时，高斯点对应的子离散方程表达式为：

其中，U_kg为高斯点支持域节点的电位场值构成的电位向量，F_kg为高斯点支持域节点构成的子离散方程的右端项；

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

其中，i和j表示高斯点支持域Ω_q中的节点局部编号，i＝1,2,…,n，j＝1,2,…,n；

采用节点局部积分域包含高斯点x_g的节点作为x_g的支持域，则计算单位分解函数ψ_k(x_g)使用的节点与构造形函数φ_i(x_g)使用的节点相同，由于无单元法形函数满足单位分解性质，即

因此采用无单元法形函数作为单位分解函数，即ψ_k(x_g)＝φk(x_g)，其中φ_k(x_g)为当前计算节点k对应的无单元法形函数；此时29)式和30)式分别写为：

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

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

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

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

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

KU＝F 38)。