点源二维电场变分问题的无网格解法

引用本文

李俊杰, 严家斌. 点源二维电场变分问题的无网格解法[J]. 物探与化探, 2015,39(3): 606-609.
LI Jun-Jie, YAN Jia-Bin. A mesh-free method for the variational problem of the point source two-dimensional electric field[J]. Geophysical and Geochemical Exploration, 2015,39(3): 606-609.
Doi:10.11720/wtyht.2015.3.29 复制到剪切板

Permissions

《物探与化探》编辑部

点源二维电场变分问题的无网格解法

李俊杰¹, 严家斌²

1. 浙江省水利水电勘测设计院,浙江杭州 310002

2. 中南大学地球科学与信息物理学院,湖南长沙 410083

通讯作者: 严家斌(1969-),男,湖南常德人,副教授,从事电磁法数据处理及图像处理研究工作。Email:cspyry@csu.edu.cn

作者简介: 李俊杰(1989-),男,江西上饶人,硕士,助理工程师,从事电法勘探无网格化正演研究工作。 Email:838885421@qq.com

收稿日期: 2014-05-30

基金: 国家自然科学基金资助项目(40874055); 湖南省自然科学基金资助项目(07JJ5065)

摘要

无网格法作为一种新型数值方法,精度高、自适应分析容易,避免了复杂的网格生成过程,在计算力学领域应用广泛。尝试将无网格算法用于点源二维电场的计算,从点源二维变分问题出发代入移动最小二乘近似构造的形函数,推导了与之对应的无网格总体矩阵表达式并用含背景网格的高斯积分将其离散化;通过一个简单的二层模型算例验证了算法的正确性。

关键词: 无网格法; 点源; 移动最小二乘; 数值计算

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2015)03-0606-04

A mesh-free method for the variational problem of the point source two-dimensional electric field

LI Jun-Jie¹, YAN Jia-Bin²

1. Zhejiang Design Institute of Water Conservancy and Hydroelectric Power, Hangzhou 310002, China

2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China

Abstract

The mesh-free method, as a new numerical method which avoids the mesh generation and possesses the advantages of high precision and simple adaptive analysis, has been widely used in the field of computational mechanics. In this paper, the authors attempt to calculate the point source two-dimensional electric field by using mesh-free method. The moving least square approximation theory is introduced and stiffness matrix of mesh-free method is deduced into shape function structured by moving least squares approximation corresponding to the point source two-dimensional variational problem, and the expression is discretized by Gauss integral containing background grid. At last, the validity of the algorithm is verified by a simple calculation of two-dimensional model.

Keyword: mesh-free method; point source; moving least square; numerical calculation

Show Figures

点源二维电场正演 $\begin{matrix} ^{[13]} \end{matrix}$ 是研究直流电场分布规律的基础, 有限元法作为经典的网格数值算法适用于复杂物性分布和复杂边界形状的计算在地球物理电磁法领域应用广泛, 其最大缺陷在于求解高维问题时网格生成复杂。无网格法 $\begin{matrix} ^{[45]} \end{matrix}$ 通过局部支持域内的节点构造形函数及建立离散系统方程, 是一种基于节点的数值方法。该方法避免了网格生成的困难, 构造高阶形函数方便, 处理复杂模型较有限元便利, 作为有限元法的重要补充无网格法成功的解决了许多网格方法难以处理的问题 $\begin{matrix} ^{[68]} \end{matrix}$ , 现已成为计算力学领域的研究热点, 然而其在地球物理学电磁学领域的研究报导较少, 报导的文献主要研究了均匀节点分布下的无网格化二维电磁场响应 $\begin{matrix} ^{[913]} \end{matrix}$ , 点源二维电场的数值模拟至今未见报导。

笔者以点源二维电场问题为例介绍了无网格法的近似原理, 从点源二维变分问题出发结合无网格思想推导了总体矩阵的无网格化离散表达式, 并通过一个简单二层地质模型的计算验证了算法的有效性。

1 点源二维电场波数域变分问题

点源为三维源, 需先采用傅立叶变换将三维点源从空间域转化为波数域二维源, 当地下电性结构为二维时, 取走向为z轴, x轴与z轴垂直, y轴垂直向上, 求解域为矩形区域, 4个顶点依次以A、B、C、D逆时针编号, 波数域电位满足变分问题^[14]:

式中:▽是二维哈密顿算子, ▽= $\begin{matrix} \frac{\partial}{\partial x} \end{matrix}$ e_x+ $\begin{matrix} \frac{\partial}{\partial y} \end{matrix}$ e_y; σ 为电导率; k为波数; I为供电电流; δ 为单位脉冲函数, P为点源的位置; r为点源到边界Γ 上的距离; K₀与K₁分别为第二类零阶与一阶修正Bessel函数; cos(r, n)为矢径r与边界外法向n的夹角的余弦。

2 变分问题的无网格解法

无网格法利用位于积分点支持域内的场节点构造形函数, 用Galerkin法建立系统方程, 近似手段多采用移动最小二乘法, 其优点是能够通过权函数的适当选取来得到具有高阶连续性的形函数, 系统矩

阵表达式中包含的积分可采用含背景网格的高斯积分求解。图1为无网格法背景网格、支持域、积分点与场节点示意。

由于背景网格的积分常选用高斯积分法, 故积分点又称高斯积分点或高斯点, 对于节点均匀分布的情况支持域尺寸一般取横、纵向节点间距与支持域无量纲尺寸之积, 本文中支持域无量纲尺寸α =1.0。

	Figure Option View Download New Window
	图1 无网格法背景网格、支持域、高斯点与场节点示意

2.1 无网格离散系统方程的构造

将变分符号代入式(1), 代入MLS构造的的形函数 $\begin{matrix} ^{[46, 10]} \end{matrix}$ , 得

式中:Φ 为形函数矩阵; n为支持域内的节点数; U为支持域内n个节点的场向量。将式(3)、式(4)代入式(2), 有

式(5)采用支持域内n个场节点编号, 对于求解域内所有场节点还应有一个用于将局部节点矩阵组装成总体刚度矩阵的总体编号体系。当节点I和节点J位于不同支持域时, 其相应的被积函数为零, 因此可得到按求解域节点编号的总体编号体系方程

式(6)右端的积分代表点源对二维电场的贡献, δ 函数积分存在式(7)的数值解

$\begin{matrix} \int_{Ω} δ (P) dΩ = \int_{Ω} δ (x_{0}) δ (y_{0}) dΩ = \{\begin{matrix} \frac{1}{2} & P \in Ω \\ 0 & P \notin Ω \end{matrix} (7) \end{matrix}$

式中x₀与y₀为点源P所在位置的坐标, 将式(7)代入式(6)得

$\begin{matrix} \{\begin{matrix} δ U^{T} (KU - P) = 0 \\ K = \int_{Ω} [σ (\frac{\partial φ_{I}}{\partial x} \frac{\partial φ_{J}}{\partial x} + \frac{\partial φ_{I}}{\partial y} \frac{\partial φ_{J}}{\partial y}) + k^{2} σ φ_{I} φ_{J}] dΩ + \\ \int_{Γ_{\infty}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} \cos (r, n)] φ_{I} φ_{J} dΓ \\ P = {[0 \dots 0 \frac{I}{2} 0 \dots 0]}^{T} \end{matrix} (8) \end{matrix}$

由于δ U^T是任意的, 故式(8)成立的条件为

$\begin{matrix} KU = P 。 (9) \end{matrix}$

式(9)即为无网格法构造的系统方程。由于求解域为矩形, 具体计算时应先将K包含的无穷边界积分项展开成

$\begin{matrix} \begin{matrix} \int_{Γ_{\infty}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} \cos (r, n)] φ_{I} φ_{J} dΓ = \\ \int_{Γ_{AB}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} (- \frac{| x_{0} - x |}{r})] φ_{I} φ_{J} dΓ + \\ \int_{Γ_{BC}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} (- \frac{| y |}{r})] φ_{I} φ_{J} dΓ + \\ \int_{Γ_{CD}} [\frac{kσ K_{1} (kr)}{K_{0} (kr)} (- \frac{| x_{0} - x |}{r})] φ_{I} φ_{J} dΓ (10) \end{matrix} \end{matrix}$

式中r= $\begin{matrix} \sqrt[]{(x_{0} {- x)}^{2} + (y_{0} {- y)}^{2}} \end{matrix}$ (x、y为高斯点的横、纵坐标), 求解出波数域的电位U后, 就可利用表1给出的一组最优化波数^[14]将波数域电位转化为空间域电位, 转换表达式为

表1 最优化波数

2.2 求解域背景网格积分

K的表达式中包含对求解域Ω 与求解域边界Γ 的积分, 为计算这些积分, 须将求解域离散成一组背景网格(图1), 总体积分可表示成这些单元积分之和的形式, 每个单元的积分利用高斯积分法求解:

式中:n_c与n_cΓ分别为背景单元及边界单元的个数; n_g与n_gΓ分别为背景单元及边界单元内的高斯点数目, 本文中每个背景单元采用4(2× 2)个高斯点, 每个边界单元采用2个高斯点; ω _i为第i个高斯点X_Qi的权系数; $\begin{matrix} J_{ik}^{D} \end{matrix}$ 与 $\begin{matrix} J_{il}^{B} \end{matrix}$ 分别为背景单元k在积分点X_Qi处面积分及子边界l在积分点X_Qi处曲线积分的雅可比(Jacobian)矩阵。式(12)即为总体刚度矩阵的最终离散表达式, 它是带状、对称、稀疏的^[4]。

3 算法验证

点源二维问题数值计算精度一般受求解域大小、场节点分布方式及最优化波数的选取控制, 影响因素较为复杂。不同电极距的计算往往涉及节点的重新分布, 节点的变动有时会导致较差的计算效果。本文旨在介绍二维变分问题的无网格解法, 故先以单个电极距(r=1 m)的计算来验证算法的正确性, 影响无网格法求解精度的因素仍需进一步研究。

研究一个二层介质模型:第一层电阻率ρ ₁=1 Ω · m, 层厚h₁=2 m; 第二层电阻率ρ ₂=19 Ω · m。表2为该模型在偏移距r=1 m处的计算电位值。

表2 r=1 m时模型不同数值方法的电位计算结果

4 结论

(1) 将无网格算法应用于点源二维问题的正演计算, 简述了无网格法基本原理。

(2) 从点源二维变分问题出发, 代入无网格法构造的形函数, 推导了点源二维波数域电位的总体矩阵的表达式, 均匀半空间模型电位的数值计算验证了无网格算法的正确性, 其计算精度与有限元法双二次插值相当。

由于场源点是奇异点, 在点源附近随着r的增大电位衰减很快, 因此r较小时双二次插值的精度高于线性插值。表2显示无网格法计算结果与有限元双二次插值相当。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	汤井田, 肖晓, 杜华坤, 等. ANSYS 在直流电法正演中的应用[J]. 地球物理学进展, 2006, 21(3): 987-992. [本文引用:1]
[2]	汤井田, 王飞燕. 基于非结构化网格的2. 5D直流电阻率模拟[J]. 物探化探计算技术, 2008, 30(5): 413-418. [本文引用:1]
[3]	汤井田, 辛会翠, 王冉. 点电源下复杂角域地形影响及校正[J]. 吉林大学学报: 地球科学版, 2012, 42(1): 254-261. [本文引用:1]
[4]	Belytschko T, Lu Y Y, Gu L. Element-free galerkin methods[J]. International Journal for Numerical Methods in Engineering, 1994, 37(2), 229-256. [本文引用:1]
[5]	李俊杰, 严家斌. 无网格法进展及其在地球物理学中的应用[J]. 地球物理学进展, 2014, 29(1): 452-461. [本文引用:1]
[6]	Belytschko T, Lu Y Y, Gu L, Fracture and crack growth by element-free Galerkin methods[J]. Modelling and Simulation in Materials Science and Engineering, 1994, 2(3): 519-534. [本文引用:1]
[7]	Wang S L N. A large-deformation Galerkin SPH method for fracture[J]. Journal of Engineering Mathematics, 2011, 71(3): 305-318. [本文引用:1]
[8]	Liu H S, Xing Z W, Yang Y Y. Simulation of sheet metal forming process using reproducing kernel particle method[J]. International Journal for Numerical Methods in Biomedical Engineering, 2010, 26(11): 1462-1476. [本文引用:1]
[9]	冯德山, 王洪华, 戴前伟. 基于无单元Galerkin法探地雷达正演模拟[J]. 地球物理学报, 2013, 56(1): 298-308. [本文引用:1]
[10]	李俊杰, 严家斌. 无网格法进展及其在地球物理学中的应用[J]. 地球物理学进展, 2014, 29(1): 452-461. [本文引用:1]
[11]	李俊杰, 严家斌. 无网格点插值法大地电磁二维正演数值模拟[J]. 石油物探, 2014, 53(5): 617-626. [本文引用:1]
[12]	Wittke J, Tezkan B. Meshfree magnetotelluric modeling[J]. Geophysical Journal International, 2014, 198(2): 1255-1268. [本文引用:1]
[13]	李俊杰, 严家斌. 无网格局部Petrov-Galerkin法大地电磁二维正演[J]. 煤田地质与勘探, 2014, 42(6): 101-104. [本文引用:1]
[14]	徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 159-170. [本文引用:2]

2006

0.0

汤井田, 肖晓, 杜华坤, 等. ANSYS 在直流电法正演中的应用[J]. 地球物理学进展, 2006, 21(3): 987-992.

With point electric source as example, The paper introduces process of carying out DC method forward modeling with ANSYS. Bycomparing and analyzing ANSYS modeling result and theoretic value of typical model, It proves the validity of DC method forward modeing with ANSYS. ANSYS has the virtue for powerful analysis capability, quick and accurate forward modeing and convenient operation. It is important in appying ANSYS in geophysical modeling.

以点电源场为例，详细叙述了直流电法正演在ANSYS上的实现过程.通过正演模型ANSYS模拟解与理论值的对比分析，验证了应用ANSYS进行直流电法正演模拟的正确性.鉴于ANSYS强大的分析功能，快速、准确的模拟以及便于操作等优点，对ANSYS在地球物理正演中的应用的研究是有重大意义的.

2008

0.0

2012

0.0

汤井田, 辛会翠, 王冉. 点电源下复杂角域地形影响及校正[J]. 吉林大学学报: 地球科学版, 2012, 42(1): 254-261.

Angle-domain is the basic unit to describe rugged terrain. Apparent resistivity anomaly of complicated terrain can be drawn using the method of angle domain overlapping and then the topographic correction can be performed. We superpose the distorted apparent resistivity of the simple angle-domain terrain to get the distorted apparent resistivity of complicated angle-domain terrain. The method analyses the error of combination and superposition. It’s also compared with the result of ANSYS modeling. It’s proved that the apparent resistivity anomaly of complicated angle-domain terrain can be gained by combination and superposition when the terrain does not vary greatly. And finally we show some examples in which the topographic effect is obtained by the method of angle domain overlapping,and the topographic correction is conducted in terms of the terrain correction formulation. The comparison with the result of horizontal terrain shows that the comparison method can weaken the topographic influence，and extrude the anomaly of the local geological body.

Key Laboratory of Metallogemic Prediction of Nonferrous Metals,Ministry of Education/School of Geosciences and Info-Physics, Central South University,Changsha410083,China

角域是描绘起伏地形的基本单元,通过角域叠加可以获得复杂纯地形产生的视电阻率异常，进而可以进行地形校正。将简单角域视电阻率畸变进行组合叠加得到复杂角域地形的视电阻率畸变，分析了组合叠加误差，并与ANSYS软件数值模拟结果进行了对比。结果表明，在地形起伏不大的情况下，可以通过组合叠加获得复杂角域地形的视电阻率异常。最后给出了利用组合叠加法获得纯地形影响后应用比较法进行地形改正的模型算例；与水平地形时的结果对比表明，比较法可以削弱地形影响,突出局部地质体的异常。

1994

0.0

... 式(12)即为总体刚度矩阵的最终离散表达式,它是带状、对称、稀疏的^[4] ...

2014

0.0

李俊杰, 严家斌. 无网格法进展及其在地球物理学中的应用[J]. 地球物理学进展, 2014, 29(1): 452-461.

Meshless method is a new numerical method as an important complement and development of the grid method. This paper focuses on meshless approximation methods and the weighted residual methods to elaborate the principle of meshless method. This paper also review the advantages and disadvantages of meshless method and summarize its research progress and detailedly introduced the solving process of Galerkin weak form meshless method according to magnetotelluric two-dimensional problem. Finally, Application of element-free Galerkin method in geophysics is discussed.

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China

无网格法作为网格方法的重要补充和发展，是一新兴数值方法.本文围绕无网格法的近似方式及加权残量法这两个核心内容，详细阐述了无网格法的基本原理，评述了无网格法的优缺点，并对无网格法的研究进展作了综述；以大地电磁二维问题为例，详细介绍了Galerkin弱式无网格法的求解过程；最后介绍了无单元Galerkin法（EFGM）在地球物理学中的应用.

1994

0.0

2011

0.0

2010

0.0

2013

0.0

冯德山, 王洪华, 戴前伟. 基于无单元Galerkin法探地雷达正演模拟[J]. 地球物理学报, 2013, 56(1): 298-308.

Forward simulation of Ground Penetrating Radar (GPR) is conducted by the element-free Galerkin method (EFGM) which has advantages such as simple preprocessing, high precision and solution of the high-order continuous derivatives, so it only requires nodes without elements and adopts the moving least-squares (MLS) fitting field function. The wave equation that the GPR forward simulation needs to satisfy is derived starting with the Maxwell equation, and then the constructing method of the MLS shape function is presented in detail. To address the complicated processing of the imposed boundary conditions, a penalty factor is employed to process the imposed boundary condition because the EFGM cannot meet the interpolation condition. Meanwhile, to eliminate ultra-intense reflection from the truncated boundary when the EFGM is employed in the GPR forward simulation, the truncated boundary condition is used to transmit the reflected waves from the truncated boundary, so as to restrain these reflected waves. A Mat lab program for GPR forward simulation based on the EFGM is compiled and used in the forward simulation for the typical GPR geoelectric model. A comparison of the resultant profile from the simulation with the FEM forward-simulation profile based on linear interpolation shows that the EFGM simulation has better continuity and smoothness, and the consistency of the simulation results indicates the correctness and effectiveness of the EFGM used for GPR forward computation.

Feng De-Shan 1,2,3 ;Wang Hong-Hua 1,2,3 ;Dai Qian-Wei 1,2,3 ;

无单元Galerkin法采用滑动最小二乘法拟合场函数,只需节点无需单元,具有前处理简单、精度高、解高次连续等优点,被用于求解探地雷达(GPR)正问题.本文从Maxwell方程出发,推导了GPR正演需满足的波动方程;详细介绍了滑动最小二乘法形函数的构造方法.针对EFGM不满足插值条件导致强加边界条件的处理变复杂的特性,采用罚因子法对强加边界条件进行了处理;同时为了消除EFGM进行GPR正演模拟时来自截断边界处的超强反射,采用透射边界条件把GPR波在截断边界处的反射波透射出去,进而压制了来自截断边界处的反射波.然后,编制了EFGM的GPR正演模拟Matlab程序,应用该程序对典型GPR地电模型进行了正演模拟,并把该正演剖面图与基于线性插值FEM正演剖面图进行了对比,结果表明了EFGM用于GPR正演计算的正确性及有效性,并且在相同节点数条件下,EFGM比矩形剖分的FEM的精度要高,更有利于指导雷达剖面的数据解译.

2014

0.0

李俊杰, 严家斌. 无网格法进展及其在地球物理学中的应用[J]. 地球物理学进展, 2014, 29(1): 452-461.

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China

2014

0.0

2014

0.0

2014

0.0

1994

0.0

... 1 点源二维电场波数域变分问题点源为三维源,需先采用傅立叶变换将三维点源从空间域转化为波数域二维源,当地下电性结构为二维时,取走向为z轴,x轴与z轴垂直,y轴垂直向上,求解域为矩形区域,4个顶点依次以A、B、C、D逆时针编号,波数域电位满足变分问题^[14]: ...

... 式中r=(x0-x)2+(y0-y)2(x、y为高斯点的横、纵坐标),求解出波数域的电位U后,就可利用表1给出的一组最优化波数^[14]将波数域电位转化为空间域电位,转换表达式为 ...