一种基于自然单元的2.5维自然电场数值模拟方法与流程

本发明涉及自然电场数值模拟方法技术领域，尤其涉及一种基于自然单元的2.5维自然电场数值模拟方法。

背景技术：

随着我国经济的持续高速发展，在工业化、城镇化过程中的环境污染问题日益严重。即时检测即时控制污染源的扩散是保护土壤和地下水不受污染或少受污染的积极方法。与钻孔测试、定期采样分析等传统环境监测手段相比，自然电场法具有效率高、成本低、快速无损、无二次污染等优点，且对多种与污染扩散伴生的微生物活动信号、氧化还原反应以及渗流运动相当敏感，适用于地下污染物的检测以及长周期实时监测。

有针对性地开展自然电场数值模拟工作，有助于反演解释，提高自然电场法在工程与环境污染监测、检测等方面的应用效果。自然电场数值模拟所采用的常规数值方法有积分方程法、有限元法、有限差分法以及有限体积法等。这些数值方法各有优势，但也同时存在共同的缺陷，即需要构建网格单元以构造插值形函数，该过程不仅增加了前期处理工作量，也使得复杂模型数值模拟变得困难。而自然单元法节点布置灵活，形函数构造过程无需网格单元信息，前处理过程简单方便，模拟精度高，能有效实现自然电场等稳定电流场复杂模型的数值模拟。

技术实现要素：

本发明目的在于提供一种基于laplace插值构造自然单元形函数，实现一种基于自然单元的2.5维自然电场数值模拟方法。

为实现上述目的，本发明提供了一种基于自然单元的2.5维自然电场数值模拟方法，其具体步骤如下：

s1：建立二维自然电场地电模型，依据模型参数设置离散点集；

s2：构建2.5维自然电场地电模型边值问题，推导其基于自然单元的变分问题及其积分形式；

s3：构建基于laplace插值函数的自然单元法，推导2.5维自然单元形函数及其导数的基本方程；

s4：基于s2、s3推导的公式，计算各波数条件下的总刚度矩阵；

s5：处理边界条件及电源信息；

s6：求解大型稀疏方程组得到各波数条件下的波数域自然电位剖面，然后将各波数域自然电位剖面通过傅里叶反变换得到空间域下的自然电位剖面。

作为本发明的进一步方案：所述步骤s1中的所述二维自然电场地电模型的尺度设置为20m*30m(深度*横度)；所述离散点集设置有651个节点，地下水位位于2m深度处，上层电阻率为50欧姆米，下层电阻率为100欧姆米。

作为本发明的进一步方案：所述步骤s2中的所述构建2.5维自然电场地电模型边值问题的方法为采用第三类边界条件进行构建。

作为本发明的进一步方案：所述步骤s3中的自然单元法为一种基于voronoi图和delaunay三角化几何结构，以自然邻点插值为试函数的一种无网格数值方法；

作为本发明的进一步方案：所述步骤s4中对波数及其权值的选择具体如下：

波数k＝[0.004758,0.0407011,0.1408855,0.393225,1.088038]；

波数权值g＝[0.0099472,0.031619,0.0980327,0.2511531,0.7260814]。

所述总刚度矩阵的具体过程如下：

(1)对全域的各背景积分网格进行循环；

(2)对各背景积分网格中的高斯积分点进行循环；

(3)搜索各高斯积分点的自然邻点；

(4)计算各自然邻点相关于高斯积分点的形函数及导数；

(5)计算各高斯积分点的子刚度矩阵；

(6)将各子刚度矩阵加载到总刚度矩阵中，直至全域高斯积分点循环完毕。

作为本发明的进一步方案：所述步骤s6中各波数条件下的波数域自然电位剖面的具体计算过程如下：

(1)在各波数条件下，求解得到当前波数对应的总刚度矩阵；

(2)结合电源项，求解方程组，得到各波数条件下的波数域总自然电位剖面；

(3)通过傅里叶反变换将各波数域自然电位剖面转换为空间域中的总自然电位剖面。

应用本发明的技术方案，具有以下有益效果：

(1)本发明所提供的方法不依赖于网格单元，节点布置灵活，计算精度高，适用于自然电场复杂模型的高精度数值模拟。

(2)本发明提供的方法，通过基于laplace自然邻点插值的二维自然单元法，经傅里叶正反变换，实现波数域到空间域的2.5维自然电场数值模拟。

(3)本发明提供的方法，在模型边界处满足精确的线性插值，可以准确施加边界条件。

(4)本发明提供的方法，所有计算最终都可转化为简单的代数运算，计算效率高，编程简便。

(5)本发明提供的方法，能推动自然电场等多源稳定电流场复杂模型数值模拟工作的开展，为自然电场法应用到污染监测、检测以及其他工程与环境地球物理问题中提供数值计算基础。

(6)本发明中步骤s1所建立的二维自然电场地电模型，该模型可近似模拟微生物降解地下有机污染物时，在地下水位附近发生的氧化还原反应所引起的自然电位异常。

(7)本发明中步骤s3中所采用的自然单元法具有无网格法和经典有限元方法的优点，又克服了两者的缺陷；其形函数满足插值性质，可以像有限元法一样直接施加本质边界条件，不存在基于移动最小二乘拟合的无网格方法不能直接施加本质边界条件的难题。

(8)本发明中步骤s3中的所述delaunay三角化具有最大最小角和空外接圆两个重要性质，且在二维情况下，delaunay三角形网格还具有所形成的三角形互不重叠和所形成的三角形可以覆盖整个平面的特点。

(9)本发明中，自然邻点插值函数选用laplace插值。laplace插值的优点主要体现在：a、直接利用节点的一阶voronoi结构的边长和点到voronoi边的距离计算插值基函数，避免了sibson插值所采用的waston面积叠加算法需要反复进行迭代的麻烦，计算量大大降低；b、除凸边界外，在凹边界上也满足精确的线性插值，可以准确施加边界条件。

(10)本发明中通过采用将各波数域自然电位剖面通过傅里叶反变换得到空间域下的自然电位剖面，其通过数值模拟结果表明，地质微生物降解地下有机污染物时所发生的氧化还原反应能在地表产生负的自然电位异常；同时表明本发明所提的数值模拟方法是有效的可行的。

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

附图说明

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

图1是本发明基于自然单元的2.5维自然电场实施例地电模型；

图2是本发明中节点a的voronoi单胞；

图3是本发明中7节点voronoi图；

图4是本发明中delaunay外接圆的示意图；

图5是本发明中delaunay三角化的示意图；

图6是本发明中7节点voronoi图中的x计算点的示意图；

图7是本发明中x计算点的自然邻点及其laplace插值元素的示意图；

图8是本发明中基于laplace插值的自然邻点插值函数示意图；

图9是本发明中自然邻点插值函数对x的导数的示意图；

图10是本发明中实施例地表自然电位异常的示意图。

具体实施方式

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

参见图1至图10所示，本发明提供的一种基于自然单元的2.5维自然电场数值模拟方法，其具体步骤如下：

s1：建立二维自然电场地电模型，依据模型参数设置离散点集(离散节点的布置参见图1所示)；

s2：构建2.5维自然电场地电模型边值问题，推导其基于自然单元的变分形式及其积分形式；

s3：构建基于laplace插值函数的自然单元法，推导2.5维自然单元形函数及其导数的基本方程；

s4：基于s2、s3推导的公式，计算各波数条件下的总刚度矩阵；

s5：处理边界条件及电源信息；

s6：求解大型稀疏方程组得到各波数条件下的波数域自然电位剖面，然后将各波数域自然电位剖面通过傅里叶反变换得到空间域下的自然电位剖面。

作为本发明的进一步方案：所述步骤s1中的所述二维自然电场地电模型的尺度设置为20m*30m(深度*横度)；所述离散点集设置有651个节点，地下水位位于2m深度处，上层电阻率为50欧姆米，下层电阻率为100欧姆米；

假定两处有微生物降解有机污染物现象，其中地下水位上部发生氧化反应，积累负电荷，地下水位下部发生还原反应，积累正电荷，通过离散正负点电源近似正负电荷的区域性分布，设定各点电源幅值为1ma。

作为本发明的进一步实施例：所述步骤s2中的所述构建2.5维自然电场地电模型边值问题的方法为采用第三类边界条件进行构建。

采用第三类边界条件进行构建2.5维自然电场地电模型，则其2.5维自然电场波数域总电位u的边值问题的表达示为：

其中：σ为介质电导率，i0为点源电流强度，δ(a)为与点源位置有关的δ函数，n为边界外法向单位矢量，为2维哈密顿算子、且k为波数，u为三维空间中的电位u沿y方向进行傅里叶变换的波数域总电位，k0(kr)、k1(kr)分别为零阶、一阶第二类修正贝塞尔函数，ω、γs、γ∞分别为经傅里叶变换后的研究区域、地表边界、无穷远边界，r为各点源点到外边界任意点的距离。

所述2.5维自然电场波数域总电位u的边值问题所对应的变分问题的表达示为：

δf(u)＝0(5)

对(4)式右端进行变分，则有：

假设某一计算点包含n个自然邻点，利用该n个节点构造自然单元形函数，采用该组形函数作为试函数近似计算点x＝(x,z)电位u，则有：

将(7)式代入(6)式中，则有：

对(8)式中的面积分项做变换，则有：

对(8)式中的边界积分项做变换，则有：

合并(9)、(10)式，则有：

δu^t[(k¹+k²)u-f]＝0(11)

其中：

由于δu是任意的，故(11)式成立的条件为：

(k¹+k²)u＝f(15)

即：

ku＝f(16)

其中：k＝k¹+k²为总系数矩阵。

若全局积分域用一组背景积分单元进行剖分，对单个积分单元有：

分别为背景积分单元ωe的面积分子系数矩阵、边界积分子系数矩阵、子右端项，分别为：

所述用于对全局积分域进行剖分的背景积分单元可为矩形网格、delaunay三角网等的其中一种或一种以上的组合。

若背景单元ωe中采用ng个高斯积分点xg＝(xg,zg)，g＝1,2,…,ng进行积分计算，那么，它们对应的权重wg，g＝1,2,…,ng和雅克比值jg，g＝1,2,…,ng，其高斯点支持域内包含n个自然邻点，则有：

其中：

b＝σk²φ(xg)φ^t(xg)(24)

展开(23)、(24)、(25)，则有：

其中：式(26)、(27)、(28)中的形函数及其导数，采用laplace自然邻点插值函数求取。

作为本发明的进一步实施例：所述步骤s3中的自然单元法为一种基于voronoi图和delaunay三角化几何结构，以自然邻点插值为试函数的一种无网格数值方法；

所述voronoi图及其对偶delaunay三角化(见图2、3、4、5)的概念来自于计算几何，是由一组不规则点定义的最基本的几何结构。

所述然邻点插值函数选用laplace插值。

对于两个voronoi单胞ti、tj，定义为集合ti的闭包；若d(xi,xj)≠0，则有：

所述laplace插值形函数定义为：

对于二维空间，所述laplace插值形函数的形式为：

其中，si(x,y)是与节点i关联的voronoi边的长度，hi(x,y)是插值点到节点i的voronoi边的垂直距离。

si(x,y)、hi(x,y)均为节点i的坐标(x,y)的函数；

由voronoi结构的定义可知，待求点的voronoi的结构顶点是由待求点x与相邻两自然邻节点组成的三角形外接圆圆心，si(x,y)实际是圆心间的距离；对于任意三角形δabc，设其顶点坐标为a(x1,y1)、b(x2,y2)、c(x3,y3)，则该三角形外接圆圆心v＝(vx,vy)的表达示为(设c为计算点，即积分点)：

其中：

d＝2[(x1-x3)(y2-y3)-(x2-x3)(y1-y3)](38)

α＝(x2+x1)(x2-x1)+(y2+y1)(y2-y1)(39)

d,x＝2(y1-y2)(40)

d,y＝2(x2-x1)(41)

设任意一条voronoi边的两个端点为v1＝(vx1,vy1)，v2＝(vx2,vy2)，则有：

由式(31)可得φi(x,y)关于点i坐标的一阶导数为：

其中αi关于坐标(x,y)的一阶导数为：

其中有：

基于laplace插值的自然邻点插值及其对x的导数如图8、9所示。

作为本发明的进一步实施例：所述步骤s4中对波数及其权值的选择具体如下：

k＝[0.004758,0.0407011,0.1408855,0.393225,1.088038]

g＝[0.0099472,0.031619,0.0980327,0.2511531,0.7260814]

所述总刚度矩阵的具体过程如下：

(1)对全域的各背景积分网格进行循环；

(2)对各背景积分网格中的高斯积分点进行循环；

(3)搜索各高斯积分点的自然邻点；

(4)计算各自然邻点相关于高斯积分点的形函数及导数；

(5)计算各高斯积分点的子刚度矩阵；

(6)将各子刚度矩阵加载到总刚度矩阵中，直至全域高斯积分点循环完毕。

作为本发明的进一步实施例：所述步骤s5中的边界条件经变分处理后，如式(19)所示，其具体计算过程如式(28)所示、电源项如式(20)所示，当其所形函数具有kroneckerdelta函数性质时，所述源项积分的表达示为：

作为本发明的进一步实施例：因地电模型为二维，而点源为三维源，故需要将空间域先转换为波数域进行计算，再反变换为空间域，一般称之为2.5维，其具体计算过程为：

(1)在各波数条件下，求解得到当前波数对应的总刚度矩阵；

(2)结合电源项，求解方程组，得到各波数条件下的波数域总自然电位剖面；

(3)通过傅里叶反变换将各波数域自然电位剖面转换为空间域中的总自然电位剖面，其计算结果如图10所示。

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