1 INTRODUCTION

The numerical simulation of direct current (DC) resistivity based on the partial differential equation is an active research topic in computational geophysics. At present, there are several numerical approaches available for DC resistivity forward simulation, including finite difference method (FDM), finite element method (FEM) and other methods. Among them, FEM has been widely accepted as the most suitable numerical approach for the simulation of DC resistivity method. At present, FEM is widely used in the 2D and 3D forward modeling of electrical prospecting (Xu, 1994; Ruan and Xiong, 2001; Ma and Liu, 2014). Unfortunately, FEM use the predefined mesh discretization, resulting in some disadvantages: (1) FEM has limitations for the problems of dramatic field changes; (2) FEM is not very suitable for complex shape or complex physical parameters distribution model; (3) Although the idea of domain division is really ingenious, meshing of the analysis domain can sometimes be very laborious and time consuming, such as when the number of mesh is large and especially discretized by unstructured mesh.

To avoid those undesired problems, an element-free method (EFM) is introduced to DC resistivity forward simulation. This element-free algorithm only requires the distribution of nodes instead of the mesh discretization. Therefore, EFM eliminates the limitations arising from the requirements of careful mesh discretization which are required in FEM. The moving least square (MLS) method was first used to construct the shape functions by Lancaster and Salkauskas (1981) to solve the continuity problem of shape function. EFM was proposed by Belytschko and Lu (1994) on the basis of the improvement of imposing boundary conditions and numerical integral processing problems. Since then, EFM has attracted the interest of scholars in various fields. EFM was first used to solve the problem of two dimensional electrostatic field of voltage transformer successfully in the field of electromagnetic field calculation by Cingoski and Miyamoto (1998). The basic boundary and dielectric interface problems of EFM were studied by Viana and Aparecida (1999), Herault and arechal (1999). Liu and Yang (1999) employed the EFM to calculate the electric field of parallel plate capacitor, the electric field of the long metal directly and the magnetic field of a current carrying conductor. Cingoski and Miyamoto (2000) tried to combine EFM and FEM in the analysis of electromagnetic problems and achieved good results. Liang and Zhi (2004) solved the steady-state and quasi-steady-state problems by EFM. Brighenti (2005) solved the problems of 3D fracture mechanics by EFM. Peng and Li (2011) used a complicated alterable EFM in the simulation of 2D elastoplastic problem, compared with the conventional EFM, the accuracy of simulation is higher. In the field of geophysics, EFM has not been widely used and less relevant literature was published. Jia and Hu (2006) used the EFM to simulate the seismic wave field. Wang (2007) studied the influence factors of EFM in the forward modeling of seismic wave field. Feng and Wang (2013) employed the MLS shape function in EFM for the 2D forward modeling of ground penetrating radar (GPR). Yan and Li (2014) also employed the MLS shape function in EFM for the 2D forward modeling of magnetotelluric (MT) method.

The radial point interpolation method (RPIM) is based on the radial base function (RBF). RBF has become an important tool for many current computing areas due to its excellent mathematical properties. For example, in the latest study in the field of geophysics, RBF was used for the local gravity field reconstruction (Wu et al., 2016). Unlike the MLS shape function, the RPIM shape function has some characteristics as follows: (1) The essential boundary conditions can be enforced by the direct method due to the Kronecker delta function property of RPIM, making the boundary processing simple in the EFM. (2) The slight change of the position or number of nodes in the support domain has little effect on the shape functions due to the higher stability of RPIM shape functions, which results in the strong adaptability of EFM to any irregularly distributed nodes. The study of EFM based on the RPIM shape functions in the field of DC resistivity forward modeling has not been reported. In this paper, the radial point interpolation is employed to construct the shape functions, the global weak form formulation is employed to derive the element-free equation and the integration is calculated by using a set of background integral cells. And then, the correctness and effectiveness of this algorithm are verified by comparing the results of EFM, FEM and analytic solutions. The results of this paper show that the EFM is a high-precision and flexible method for the DC resistivity forward modeling.

2 BOUNDARY VALUE PROBLEM AND VARIATION FORMULATION OF DC RESISTIV-ITY

Assuming that the electrical properties of underground media are two-dimensionally distributed, the elec-tric field generated by the point power supply is three-dimensionally distributed, in this case we call it 2.5-dimensional (2.5D) problem. For the conducted electric exploration, 2.5D problem is appropriate for the practi-cal problems of exploration, and the computational cost of 2.5D problem is much lower than that of 3D problem. So, the 2.5D problem is widely used in the conducted electric exploration. In the Cartesian coordinate system, the potential (U) of 2.5D stable electric field satisfied the partial differential Eq.(1):

(1)

where σ = (x, z) is the conductivity of medium, U (λ, x, z) is the wave-field potential, λ is the wavenumber, I is the current, δ is Dirichlet function, A is the location of the point power supply, is the 2D Hamiltonian Operator.

The boundary value problem of the third kind of boundary conditions can be described as Eq.(2):

(2)

where Ω is the computational domain bounded by Γ_s and Γ_∞, cos(r, n) is the cosine of the angle between the vector r and the external normal n, K₀ and K₁ is the second kind of first-order and zero-order modified Bessel functions respectively. Boundary value problem (2) corresponds to the variational problem (3):

(3)

3 RPIM SHAPE FUNCTION

In the computational domain Ω, the RPIM approximation u(x) is expressed as Eq.(4).

(4)

where R_i(x) is the radial base function, n is the number of radial base functions, p_j (x) is the monomial of spatial coordinates x^T = [x, y], m is the number of polynomial basis functions, a_i and b_j are unknown coefficients. The multiquadric function is used as the radial base function. There are two shape parameters: α_c and q in the multiquadric function. By the experimental testing we noted that α_c about 1.0 and q about 1.03 in the DC resistivity simulation can achieve good results. d_c is the feature length between the calculating point x and the nodes of local support domain, it is usually taken as the average spacing of nodes in the support domain.

In the radial base function R_i(x), r_i is the distance between the calculating point x and the node x_i, in the case . In the formula (4), the general form of complete polynomial degree p is expressed as . To determine the coefficients a_i and b_j of formula (4), we construct the support domain of calculation point x. And then, we can get n nodes of the support domain. Let u(x) be equal to the field value of n nodes, we can get n equations as follows

(5)

where is field value vector. The coefficient matrix of radial basis functions is expressed as follows

The coefficient matrix of polynomial basis functions is expressed as follows

is the coefficient vector of radial basis functions. is the coefficient vector of polynomial basis functions. There are n + m variables in the Eq.(5), so m equations using the m constraints shown in Eq.(6) are added into the Eq.(5).

(6)

We get Eq.(7) by combining Eq.(5) and Eq.(6).

(7)

where and . Since R₀ is a symmetric matrix, so G is a symmetric matrix. Then, by solving Eq.(7) we can get

(8)

Substituting Eq.(8) into Eq.(4) we can get the Eq.(9).

(9)

where the RPIM shape functions are expressed as the Eq.(10).

(10)

where the RPIM shape functions of n nodes are expressed as the Eq.(11),

(11)

It should be noted that, there may be cases that the local nodal distribution is extremely unreasonable when the nodes are arbitrarily distributed, which may make the G^-1 in Eq.(10) irreversible. So, the extremely unreasonable nodal distribution should be avoided when generating nodes.

4 THE GLOBAL WEAK FORM EFM FORMULATION

As shown in Fig. 1, the global weak form element free method uses a set of rectangular background cells to cover the computational domain Ω. Then, the Gaussian integral points are arranged in the background cell, and the Gaussian integral method is employed to calculate the integration in the background cell. After the support domain Ω_q of each computation point (i.e., Gauss integration point) is defined, the RPIM shape functions can be constructed by using the nodes inside the support do-main. The integrations of all the background cells are summed to get the integration of the computational do-main Ω, and then the variation problem (3) is solved.

We can get Eq.(12) by expanding Eq.(3).

(12)

Assuming that there are n nodes inside the support do-main Ω_q. Then, the RPIM shape functions are constructed by using those nodes and used as the trial functions to approximate the wavenumber field potential of an arbitrary point x = (x, z).

Substituting into Eq.(12) we can get

(13)

The area integral equation of the first term of Eq.(13) is expressed in matrix form as Eq.(14):

(14)

The boundary integral equation of the last term of Eq.(13) is expressed in matrix form as Eq.(15):

(15)

Thus, Eq.(13) is expressed in matrix form as Eq.(16):

(16)

where

In the Eq.(16), because the δU is arbitrary, so we can get (K¹ + K²)U = F. In the background cell Ω_e, (K¹ + K²)U = F is written as , where

where K_e¹ and K_e² is the area and boundary coefficient sub-matrix of background cell Ω_e, respectively. They can be written as follows

(17)

(18)

According to the Kronecker delta function property of RPIM shape functions, the expression of the elements of right hand term F is the Eq.(19):

(19)

In the computational domain Ω, the RPIM shape functions are used as the trial functions to approximate the conductivity of an arbitrary point x = (x, z). Substituting into Eq.(17) and (18) we can get Eq.(20) and (21):

(20)

(21)

Of course, other methods can also be used to calculate the conductivity of an arbitrary point too. For example, we can utilize the neighbor nodes to calculate the conductivity of the calculation point. Eq.(20) corresponds to the body integration of the intersection area of the background cell Ω_e and computational domain Ω. Eq.(21) corresponds to the boundary integration of the intersection boundary of the background cell Ω_e and computational domain Ω. N Gaussian integration points x_r = (x_r, z_r), r = 1, 2, …, N are set in the background integral cell Ω_e, and their corresponding Gaussian weights are w_r, r = 1, 2, …, N. Table 1 lists the coordinate coefficients and weight coefficients of Gaussian integral points. J_r (r = 1, 2, …, N) are the area integral Jacobian ratio of local integral domain for Gaussian integral points. Using Gaussian integration in Eq.(20) we can get

(22)

Using Gaussian integration in Eq.(21) we can get

(23)

In this situation, J_r (r = 1, 2, …, N) is the Jacobian ratio of the curve integral at the Gaussian integration point x_r of the subboundary.

The global equations KU = F are obtained by assembling all sub-coefficient matrices K_e¹ and K_e. Wave-field potential of nodes in the computational do-main Ω is obtained by solving the system of equations. Then, the space domain potential of nodes is obtained by the inverse Fourier transform.

5 NUMERICAL EXAMPLES 5.1 The Point Source Potential in a Homogeneous Half-Space Model

The computational domain Ω of a homogeneous half-space model is established with horizontal range from –60 m to 60 m in length and vertical range from 0 to 60 m in depth. And both the horizontal and vertical distances between the two neighbouring nodes are set to be 2 m. As a result, there are 61 nodes in the horizontal direction and 31 nodes in the vertical direction. The total number of nodes is 1871. The resistivity of model medium is 100 Ω·m. The point source locates at the middle of ground surface where X = 0 m, and the point source electric current I=1 A. The background integral cells are the rectangular cells whose vertices coincide with the nodes. As a result, there are 60 cells in the horizontal direction and 30 cells in the vertical direction. The total number of cells is 1800. 64 Gauss integration points are used to compute the integration in each cell (i.e., 8 Gauss integration points are used in one-dimension). The coordinate coefficients and weight coefficients of Gauss integration points are shown in Table 1. The support domain G_q of each Gauss integral point in the background integral cell is set as the background integral cell itself. So, the RPIM shape functions are constructed using the four nodes corresponding to the vertices of background integral cell. The linear basis functions are used (i.e., m = 3) in the RPIM shape functions. In the homogeneous half-space model, the point source analytical potential of an arbitrary point (except for the location of point source) can be calculated by the formula , where I is the electric current, ρ is the resistivity of medium, r is the space distance between the calculated point and the point source. The simulation point source potentials u^h calculated by EFM and the point source analytical solutions u of the ground nodes are shown in Fig. 2. The results show that the EFM results are very consistent with the analytical solutions.

5.2 The Comparison of EFM and FEM

The above-mentioned point source potential in a homogeneous half-space model is simulated by EFM and FEM respectively. The two methods use the same distributed nodes. The background integral cells used in the EFM are the same as the rectangular cells used in the FEM. In each cell, the two methods utilize the four nodes of rectangular cell to construct the shape functions. The linear interpolation in rectangular cell is employed in the FEM. The results of two methods and the analytical solutions are showed in Fig. 3 and Table 2.

As shown in Fig. 3 and Table 2, the simulation point source potential values of two methods agree with the analytical solutions very well. It is shown that the global weak form element free method based on the RPIM shape function is reliable for the forward modeling of steady current field, and the correctness and validity of EFM program is verified. The results of two methods have some distortions in the vicinity of point source due to the rapid change of the electric field here. However, the distortion of EFM is significantly slighter than that of FEM. It is shown that the simulation accuracy of EFM is higher than that of FEM under the same condition near the point source. According to Table 2, when the distance r between the calculation point and the point source is relatively large (4 m ≤ r ≤30 m) the relative error of EFM forward results is always lower than that of FEM. It is shown that the simulation accuracy of EFM is higher than that of FEM under the same condition at relatively far distance from the point source. As the distance between the calculation point and the source point increases further, the simulation accuracy of EFM and FEM tends to be consistent.

5.3 The Comparison of Refining Nodes

We refine the surface nodes of the above-mentioned model. The first refining case: in the surface depth 0~4 m, the nodes are evenly distributed at 1 m intervals, below the depth 4 m the nodes are evenly distributed at 2 m intervals, in the horizontal and vertical directions. The second refining case: in the surface depth 0~4 m, the nodes are evenly distributed at 0.5 m intervals, in the depth 4~8 m the nodes are evenly distributed at 1m intervals, below the depth 8 m the nodes are evenly distributed at 2 m intervals, in the horizontal and vertical directions. The apparent resistivity of symmetrical quadrupole sounding configuration centered at X = 0 m is simulated by EFM. The results of unrefined surface node model is compared with the results of refined surface node model, as shown in Fig. 4.

5.4 Complex Geoelectric Models

As shown in Fig. 5, a low resistivity circular anomaly body with radius r = 5 m and center located at (x, z)=(0 m, 15 m) is set in the homogeneous half-space model. The resistivity of circular anomaly body is 1000 Ωm and the resistivity of background medium is 100 Ωm.

A computational domain Ω in horizontal range of –60~60 m and vertical range of 0~50 m of the low resistivity circular anomaly model is established. The horizontal and vertical distances of nodes are evenly spaced at 1 m. In order to improve the simulation ac-curacy of near-field source and terrain, we refine the nodes of near surface. The circularly distributed nodes are used to model the circular anomaly body, as shown in Fig. 5. The resistivity survey of symmetrical quadrupole sounding configuration for the low resistivity circular anomaly model is simulated by the EFM. The apparent resistivity of model is shown in Fig. 6. There is a significant low apparent resistivity anomaly at the location of low resistivity circular anomaly body.

A computational domain Ω with horizontal range from –60 m to 60 m and vertical range from –4 m to 50 m depth of a complex undulating terrain model is established. As presented in Fig. 7, the 4 m undulating ridge and valley relative to the horizontal ground surface are located at X = -20 m and x = 20 m, respectively. The resistivity of the background medium of model is 100 Ωm.

The distribution of nodes is shown in Fig. 7. We refine the nodes of near surface, and the parallel quadri-lateral background integral cell is used to simulate the undulating terrain. The resistivity survey of symmetrical quadrupole sounding configuration for the complex undulating terrain model is simulated by the EFM. The apparent resistivity of model is shown in Fig. 8. The apparent resistivity anomaly caused by the ridge and valley undulating terrain model is anti-symmetrical about X= 0 m. The ridge corresponds to the low apparent resistivity anomaly and the valley corresponds to the high apparent resistivity anomaly.

6 CONCLUSIONS

(1) In this paper, the forward modeling of typical geoelectric models was calculated by the global weak form EFM. The results of EFM are compared with the results of FEM and the analytical solutions. Then, the correctness and validity of global weak form EFM based on the RPIM shape function for DC resistivity forward simulation are verified by the results of comparison.

(2) Our results show that the high order approximation of RPIM shape function makes the simulation accuracy of EFM higher than that of FEM under the same conditions. In the FEM, we can only obtain the linear approximation in the case of the element having only 4 nodes. However, we can obtain high order approximation in the case of only four nodes inside the support domain in the EFM because the shape function of EFM has the characteristics of high order approximation. Especially where the field changes drastically, such as near the field source, the simulation accuracy of EFM is significantly higher than that of FEM.

(3) In the case of the distribution of nodes is reasonable, we can arbitrarily refine the nodes to improve the simulation accuracy and the refining of nodes is independent without changing the other parameters in the EFM. It is difficult to implement the refining of nodes in local domain for the FEM based on the structured cells. Although it is convenient to implement the refining of nodes in local domain for the FEM based on the unstructured cells, but the unstructured cells need to be regenerated. The EFM shows a higher flexibility than the FEM which makes it superior in the adaptive simulation. So, the adaptive element-free method can be considered as the further research content.

(4) The EFM requires the independent construction of shape functions for each computation point (i.e., Gaussian integration point). Usually, in order to improve the integral accuracy, a larger number of computation points are used. However, the calculation time of EFM is proportional to the number of computation points. Therefore, the computational efficiency of EFM is lower than that of FEM. The parallel computing can be used to improve the computational efficiency. So, the parallel computing element-free method can be considered as the further research content.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (41674080, 41674079, 41574123).