一种自然电场三维多向映射耦合数值模拟方法与流程

本发明属于自然电场数值模拟方法技术领域，尤其涉及一种自然电场三维多向映射耦合数值模拟方法。

背景技术：

在地下污染或地下水监测中，地球物理方法以其成本低、效率高等明显优点，越来越多地被用来替代钻孔测试、定期采样分析等传统地下污染或地下水监测方法。

其中，自然电场法的野外观测尤其快捷便利，并且对地下水渗流、污染物扩散、离子迁徙、氧化还原反应等信号都非常敏感，特别适合地下污染、地下水的探测与监测。

有针对性地开展自然电场数值模拟工作，有助于反演解释，提高其在工程与环境污染监测、检测等方面的应用效果。对于自然电场的数值模拟，目前常规的做法是通过传统的有限元、有限差分或有限体积法等数值方法开展二维或三维的模拟工作。

这些方法各有优势，但却有共同的致命缺陷，即需要设定人工边界条件以满足计算精度。这样的限制条件使得研究区域必须剖分得足够大，且地电模型不能过于复杂，以保证电位在近似无穷远处尽可能衰减为零，以及人工边界的有效性。

同时对于常规的混合边界条件，边界积分矩阵需要累加到总刚度矩阵中，因此自然电场源的空间位置以及其数量的动态变化过程都会直接影响总刚度矩阵的计算，即总刚度矩阵需要随自然电场源的动态变化而进行多次计算。对于自然电场数值模拟来说，其自然电场源具有区域分布特征，且多为动态，常规数值模拟方法很难实现精确有效的模拟。

技术实现要素：

(一)要解决的技术问题

针对现有存在的技术问题，本发明提供一种自然电场三维多向映射耦合数值模拟方法，该方法在减少单元剖分的同时提高数值计算精度，减少总刚度矩阵的计算次数以提高数值模拟的计算效率，实现多源动态自然电场精确有效的数值模拟。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

一种自然电场三维多向映射耦合数值模拟方法，包括如下步骤：

s1、利用计算机中的工具软件，建立自然电场地电模型，所述模型的区域剖分包括有限单元和无限单元；

s2、通过理论基础，构建所述模型的边值问题；

s3、通过理论基础，对所述边值问题构造泛函并求其变分；

s4、利用工具软件，从变分问题出发，分别对所述有限单元和无限单元进行单元分析，求解得到有限单元刚度矩阵和无限单元刚度矩阵；

s5、利用工具软件，将有限单元刚度矩阵和无限单元刚度矩阵组装，得到总刚度矩阵；

s6、通过总刚度矩阵获得自然电场地电模型中各节点的电位值。

进一步地，在自然电场地电模型中，将所述有限单元和无限单元耦合在一起。

进一步地，所述无限单元包括单向映射无限单元、双向映射无限单元和三向映射无限单元。

进一步地，在所述步骤s4中，有限单元刚度矩阵的求解过程中仅采用形函数，无限单元刚度矩阵的求解过程中采用映射函数和形函数。

进一步地，在所述步骤s5中，将有限单元刚度矩阵和无限单元刚度矩阵按自然电场地电模型中的节点编号进行累加，得到总刚度矩阵。

(三)有益效果

本发明的有益效果是：

1、本发明提供的方法，将传统有限单元与三维多向映射无限单元在空间及物性上进行有效耦合之后，新耦合算法可实现电位逐渐衰减至无穷远的过程，能充分保证复杂模型的计算精度。

2、本发明提供的方法，大幅度减小有限单元的剖分区域，减小总自由度，提高计算效率，节约机时及运行内存。

3、本发明提供的方法，边界延伸至无穷远，不需要再考虑传统有限元法中的人工边界条件，总刚度矩阵与自然电场源再无联系，计算更加方便高效，且新提出的三种映射无限单元形函数进一步提高了计算精度。

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

附图说明

图1为本发明中复杂地形渗流模型；

图2为本发明中复杂地形渗流模型地形数据图；

图3为本发明中三维自然电场传统有限元模型；

图4为本发明中三维有限元与三维多向映射无限元的耦合示意图；

图5为本发明中一级六面体、二级四面体剖分方式示意图；

图6为本发明中四面体单元示意图；

图7a、7b为本发明中三维单向无限元映射过程；

图8a、8b为本发明中三维双向无限元映射过程；

图9a、9b为本发明中三维三向无限元映射过程；

图10a、10b分别为本发明中复杂地形渗流模型地表电位、x轴向主剖面地表电位的模拟结果；

图11为本发明中复杂地形渗流模型模拟结果(x轴向电位切片)；

图12a、12b分别为本发明中垃圾填埋场非渗漏、渗漏渗流模型；

图13a、13b分别为本发明中垃圾填埋场渗流模型的扩散路径、达西速度模型；

图14a～14h为本发明中垃圾填埋场渗流模型地表电位、x轴向主剖面地表电位的模拟结果；

图15a～15d为本发明中垃圾填埋场渗流模型模拟结果(x轴向电位切片)。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

实施例1

s1、如图1所示，利用计算机中的工具软件，构建复杂地形渗流模型。

复杂地形渗流模型中，有限元区域剖分为60×60×40个六面体网格单元。模型为两层，有限元区域尺度设置为90m×90m×110m，其中地表覆盖层厚度设为10m，电阻率取200ω·m，底部基岩电阻率取1000ω·m。地形数据如图2，依据地势变化设定3条优势流通道，渗流通道都位于x轴向主剖面，其中通道1、3向剖面中部倾斜，通道2保持直立，渗流通道横向尺度为0.5m×0.5m，电阻率为20ω·m。

s2、构建复杂地形渗流模型边值问题。如图3所示，为三维自然电场传统有限元模型，地下介质中的自然电场分布满足泊松方程为：

式中ω为研究区域，σ为电导率，u为电位，j为电流密度。在直流电法中，泊松方程右端的电源项为点电源；而一般自然电场源通过点源进行模拟，故而此处可通过三维多点源模型近似模拟自然电场源，其满足的泊松方程为：

式中ii为各点源电流，n为点源数量，δ(ai)为关于点源所在空间位置ai点的δ函数，其满足：

在地表γs上，电位法向导数为零：

在近似无穷远边界γ∞上，电位满足混合边界条件：

其中n为边界外法向单位向量，ri为各点源点到无穷远边界γ∞某点的距离。区域内部电导率边界为自然边界条件，不予考虑。

s3、对边值问题构造泛函并求其变分：

式中f(u)表示由边值问题构造的泛函，δf(u)为泛函的变分。

s4、从变分问题出发，分别对有限单元和无限单元进行单元分析，计算有限单元和无限单元的刚度矩阵。

在有限单元剖分区域中，先以六面体对有限单元区域进行规则剖分，然后将每个一级六面体二次剖分为5个二级四面体单元(如图4所示)，相应将地形数据加载到各个四面体单元中即可。

在四面体单元中(如图5所示)，空间坐标x、y、z及电位u的插值函数为：

上式中，分别表示空间坐标x、y、z，nk为形函数，且其中k＝i,j,l,m，ve为四面体体积。

其中，xi、yi、zi，xj、…、zm分别为四面体各节点的空间坐标，按右手法则循环计算得到aj、al、am、bj、bl、bm、cj、cl、cm、dj、dl、dm，其中k1依次取0、1、2、3，k2依次取1、2、3、4。

整理可得到四面体单元体积分系数矩阵为：

k1e＝(kij)(8)

其中kij为体积分系数矩阵元素，i、j＝1，…，4，对于边界积分项，近似cos(r,n)及r为常量，将其提到积分号外，进而整理得到边界积分为：

其中，符号s表示四面体单元边界面，i为各点源编号，k2e为边界积分系数矩阵，ue为单元节点电位向量。

有限元区域变分问题的积分由各子单元e上的积分系数矩阵ke扩展成由全体节点组成的有限单元刚度矩阵kf。在传统有限单元法中，ke＝k1e+k2e；在有限元-无限元耦合法中，ke＝k1e。

如图6所示，给出了三维有限元与三维多向映射无限元的耦合示意图。其中a区域表示有限单元剖分区域，位于研究区域中部，且其中一个面为地表，其余五个面与无限单元进行耦合；b区域表示单向映射无限元区域，位于四周以及底部；c区域表示双向映射无限元区域，位于四周以及底部共8条棱边；d区域表示三向映射无限元区域，位于底部边界的4个角点。

对于单向映射无限单元，其映射过程见附图7a、7b。单向映射无限元经过有限单元剖分区域的边界节点，沿边界面向外垂向映射至无穷远。单元映射由有限元边界面的四个节点1、2、3、4，分别经过节点5、6、7、8映射至无穷远，相应的将空间呈单向无限延伸的四棱柱映射为局部空间中边长为2的规则六面体。其中有限元边界面也是无限单元的映射起始面，面为无限单元中间节点面，可通过水平地表中心映射原点m到面的垂向距离来控制电位在无限单元中的衰减程度。无穷远的4个节点电位值为0，无需考虑。

设定ζ正向为映射无穷远方向，无限单元中任意点坐标映射公式可表示为：

其中，x、y、z为无限单元中任意点的空间坐标，x1、…、x8、y1、…、y8、z1、…、z8为无限单元节点空间坐标，ξ、η、ζ为无限单元任意点的局部坐标，mi为单向映射函数，其中i＝1、…、8，映射函数表达式为：

无限单元中任意点(x，y，z)的电位满足：

其中，u为无限单元任意点的电位值，ui为无限单元各节点电位值，i＝1、…、8，ni为本发明提出的单向映射无限元形函数，其表达式为：

对于双向映射无限单元，其映射过程见附图8a、8b。双向映射无限元在有限元区域棱边沿两个方向垂向映射至无穷远。单元映射由无限元边界面的四个节点1、2、3、4以及面的四个节点1、2、6、5，分别经过节点5、6、7、8及节点4、3、7、8映射至无穷远，相应的将空间呈双向无限延伸的六面体映射为局部空间中边长为2的规则六面体。其中有限元边界棱边也是无限单元的映射起始棱边，棱为无限单元中间节点棱边，面和面为无限单元中间节点面，可通过水平地表中心映射原点m到面及的垂向距离来控制电位在无限单元中的衰减程度。无穷远的8个节点电位值为0，无需考虑。

设定η、ζ正向为映射无穷远方向，无限单元中任意点坐标映射公式可表示为：

其中，x、y、z为无限单元中任意点的空间坐标，x1、…、x8、y1、…、y8、z1、…、z8为无限单元节点空间坐标，ξ、η、ζ为无限单元任意点的局部坐标，mi为双向映射函数，其中i＝1、…、8，映射函数表达式为：

无限单元中任意点(x，y，z)的电位满足：

其中，u为无限单元任意点的电位值，ui为无限单元各节点电位值，i＝1、…、8，ni为本发明提出的双向映射无限元形函数，其表达式为：

对于三向映射无限单元，其映射过程见附图9a、9b。三向映射无限元在有限元区域底部角点沿三个方向垂向映射至无穷远。单元映射由无限元边界面的四个节点1、2、3、4，面的四个节点1、2、6、5，面的四个节点1、4、8、5，分别经过节点5、6、7、8，节点4、3、7、8，节点2、3、7、6映射至无穷远，相应的将空间呈三向无限延伸的六面体映射为局部空间中边长为2的规则六面体。其中有限元底边界角点“1”也是无限单元的映射起始角点，角点“7”为无限单元中间角点，面及为无限单元中间节点面，可通过水平地表中心映射原点m到面及的垂向距离来控制电位在无限单元中的衰减程度。无穷远的12个节点电位值为0，无需考虑。

设定ξ、η、ζ正向为映射无穷远方向，无限单元中任意点坐标映射公式可表示为：

其中，x、y、z为无限单元中任意点的空间坐标，x1、…、x8、y1、…、y8、z1、…、z8为无限单元节点空间坐标，ξ、η、ζ为无限单元任意点的局部坐标，mi为三向映射函数，其中i＝1、…、8，映射函数表达式为：

无限单元中任意点(x，y，z)的电位满足：

其中，u为无限单元任意点的电位值，ui为无限单元各节点电位值，i＝1、…、8，ni为本发明提出的三向映射无限元形函数，其表达式为：

无限元单元刚度矩阵的求解过程中，jacobi变换矩阵表达式为：

偏导数中的可通过jacobi变换矩阵计算，有：

记fix、fiy、fiz分别为ξ、η、ζ的函数，整理得到无限单元刚度矩阵表达式为：

其中，i、j＝1、…、8，kij为无限单元刚度矩阵元素。无限元区域变分问题的积分由各子单元的积分系数矩阵元素kij扩展成由全体节点组成的无限单元总刚度矩阵kif。

s5、将有限单元和无限单元的刚度矩阵按节点编号累加到耦合法的总刚度矩阵k相应位置，进而实现两种单元在空间及物性上的有效耦合，如下式：

其中，为耦合法的总刚度矩阵k中的刚度值，为有限单元刚度矩阵kf中相应位置的刚度值，为无限单元刚度矩阵kif中相应位置的刚度值，i、j为该节点在耦合法的总刚度矩阵中的位置标号。

s6、通过复杂地形渗流模型的自然电场源的空间位置计算其相应的节点编号，再将各点源的幅值按节点编号赋值到源向量p中。求解由耦合法的总刚度矩阵k、全域节点电位列向量u及源向量p组成的大型稀疏方程组：

ku＝p(26)

从而得到各节点的电位值。

再结合各节点的空间坐标信息，绘制自然电位分布图件。

如图10a、10b所示，分别为复杂地形渗流模型地表电位、x轴向主剖面地表电位的模拟结果。图11为复杂地形渗流模型x轴向切片电位的模拟结果。

实施例2

实施例2与实施例1的不同之处在于，实施例2中构建垃圾填埋场动态渗流模型。

垃圾填埋场动态渗流模型中，有限元区域剖分为60×60×40个六面体网格单元。有限元区域尺度为90m×90m×110m，填埋场位于模型中部区域，填埋空间为20m×20m×10m。地表覆盖层厚度为1m，电阻率200ω·m，围岩电阻率1000ω·m。防渗墙由混凝土及防渗膜等多层防渗材料组成，厚度设为1m，综合电阻率5000ω·m。填埋区电阻率50ω·m，渗滤液电阻率5ω·m。泄漏口位于x轴主剖面、防渗墙底板及x轴正向防渗墙面的三面交点处，且沿x轴正向渗漏，横向尺度0.5m×0.5m。在正常非渗漏情况下(如图12a所示)，假定渗滤液水头面平行地表，即填埋区内部横切面为等位面。而渗漏情况(如图12b所示)，渗滤液经由泄漏口向地下渗漏，相应流动电位随水头压差增大而逐渐升高，最终在无穷远处衰减至零。设定渗滤液扩散过程满足液体在饱和孔隙介质中适用的达西公式，其表达式为：

其中，v为达西速度，k为渗透率，φ为介质孔隙度，u为动黏度，p为压差，ρ为流体密度，g为重力加速度，z为高差。为简化模型，设定污染流可扩散的渗流路径为图13a所示，该路径位于x轴主剖面的相应空间位置。假定达西速度与压差及高差呈线性关系，将扩散速度模型简化为图13b所示，即可计算得到污染流水头与泄漏点的时空关系。

由于垃圾填埋场渗流属于动态源问题，因此相比实施例1，实施例2在s6前需增加一步：

更新所述步骤s6中的自然电场源的空间位置信息，重复步骤s6。

图14a～14h为垃圾填埋场动态渗流模型模拟结果，图14a为垃圾填埋场泄漏前地表电位等值线图，图14b为垃圾填埋场泄漏第一天地表电位等值线图，图14c为垃圾填埋场泄漏第二天地表电位等值线图，图14d为垃圾填埋场泄漏第三天地表电位等值线图，图14e为垃圾填埋场泄漏前x轴向主剖面地表电位曲线，图14f为垃圾填埋场泄漏第一天x轴向主剖面地表电位曲线，图14g为垃圾填埋场泄漏第二天x轴向主剖面地表电位曲线，图14h为垃圾填埋场泄漏第三天x轴向主剖面地表电位曲线。

图15a～15d分别为垃圾填埋场未泄漏、垃圾填埋场泄露第一天、垃圾填埋场泄露第二天和垃圾填埋场泄露第三天的x轴向电位切片数值模拟结果。

以上结合具体实施例描述了本发明的技术原理，这些描述只是为了解释本发明的原理，不能以任何方式解释为对本发明保护范围的限制。基于此处解释，本领域的技术人员不需要付出创造性的劳动即可联想到本发明的其它具体实施方式，这些方式都将落入本发明的保护范围之内。