A Discrete Resistance Network Based on a Multiresolution Grid for 3D Ground-Return Current Forward Modeling

Abstract

While the high-voltage direct current (HVDC) transmission system is in monopolar operation, it produces thousands of amperes of ground-return currents (GRCs). Accurate computation of the GRCs is essential for assessing safety implications for nearby industrial infrastructure. Current three-dimensional forward models of GRCs are typically constructed based on discrete differential equations, and their solving efficiency is constrained by the increased degrees of freedom resulting from the fine discretization grids in high-conductivity conductors and ground points. To address this issue, we present a new resistor network (RN) forward solver based on a multi-resolution grid approach. This solver utilizes an RN to avoid the massive degrees of freedom resulting from fine discretization of high-voltage conductors and enhances grid discretization efficiency near the surface grounding system through multi-resolution grids. We demonstrate, through multiple three-dimensional geoelectrical model cases, that the proposed method reduces the forward modeling misfit to 1% and possesses only 3‰ of the required discrete elements compared to traditional approaches. Furthermore, practical HVDC grid model analyses indicate the successful application of the proposed method for GRC analysis in complex geoelectric conditions.

1. Introduction

High-voltage direct current (HVDC) is suitable for projects involving the transmission of substantial power over extended distances. Due to the special erection and operation mode of HVDC transmission lines, when they are in a monopolar operation state, they inject high-power direct current (DC) into the ground, generating ground-return currents (GRCs) and earth surface potentials (ESPs) [1,2]. Unconsidered industrial stray current flowing into an alternating current (AC) power grid is harmful to power grid equipment, particularly if it leads to DC bias in transformer equipment. Therefore, accurate calculation of a substation’s neutral point current is crucial for the safe operation of AC/DC power grids [3,4].

Computing current and potential is a typical DC propagation problem [5] that is commonly solved via partial differential equation methods, represented by the moment method [6] and the finite element method (FEM) [7,8,9]. However, these differential methods struggle to discretize elongated, high-conductivity conductors, which can result in a massive system of differential equations with many degrees of freedom. Considering that these methods ignore the shunt effect of the overhead lines on the current and generate the solving misfit of the potential and neutral current, the resistance network (RN) method was developed [2,10]. Based on this, Kirchhoff’s resistance network was introduced, which constructs node current and voltage conservation equations for each differential unit [11]. This resistance-network-based differential discretization method avoids the problem of discretizing high-conductivity conductors encountered with traditional differential methods and is suitable for calculating problems involving long, straight conductors, steel casings, and ground electrical structures [11,12].

Building on Kirchhoff’s resistance network method, we utilized the node voltage (NV) principle [13] to construct a symmetric matrix system that includes both current conservation and voltage conservation equations. This system not only retains all unit information but also offers superior solving efficiency. However, the commonly used mesh for this system is relatively simple, primarily consisting of hexahedra. Using traditional structured grids to accurately represent rapid changes near the surface requires complex discretization, leading to increased computational costs due to the richness of the structured grid. The multiresolution (MR) grids presented in [14,15] strike a balance by employing primarily structured hexahedrons complemented by localized refinement based on nonconforming cells. This approach to grid construction enables adaptable resolution within a fundamentally structured grid, ensuring consistent mesh aspect ratios and enhancing the efficiency of grid discretization.

Following this introduction, Section 2 introduces the principles of the GRC. The proposed 3D RN based on MR grids is elaborated, including the processing of hanging nodes, boundary conditions, and grounding nodes, in Section 3. Section 4 contains numerical examples and analysis related to power grid applications, respectively. In Section 5, we draw our final conclusions regarding this work.

2. Related Works

2.1. Background

In this section, we introduce some basic concepts and principles of GRC and describe how the current during HVDC monopole operation propagates underground and leads to surface potential distribution. Additionally, we present the governing equations for this current propagation problem.

For initial cost considerations during the construction phase of the HVDC project, the system operated in a monopolar grounding configuration for several months until the activation of the second pole. Moreover, during HVDC monopole faults or maintenance periods, the system adopts asymmetric and monopolar operation modes in which the DC enters the earth through a circuit. Current enters the earth through one grounding electrode and exits from the grounding electrode at the opposite end. The generation of GRC occurs during HVDC monopolar operation, with DC ground currents (Figure 1) reaching up to thousands of amps [1].

When calculating a ground electrical model with an overhead power line above it (Figure 2), the shunting characteristics of the power line can alter the distribution of underground currents. This represents a typical field-circuit coupling problem. We calculate the resistance of the power line based on its length and cross-sectional area and integrate this resistance with the circuit representing the ground electrical structure to form a complete circuit for solving the field-circuit coupling problem. Since the computational model area is relatively small, this field-circuit coupling model can be discretized using the RN method. For larger research areas with significant longitudinal and transverse dimensions, this study employs a multi-resolution approach to establish a volumetric circuit model based on anisotropic meshes. This model allows for the flexible extension of the computational area, making it suitable for large-scale calculations over tens to hundreds of kilometers.

According to the law of conservation of charge, when the bulk charge density is 0, the differential form of the current continuity equation can be expressed as

$$\nabla \cdot j=0,$$

(1)

where the current density j can be divided into the known current density j₁ and the unknown current density j₂.

Considering the interplay among the current density, electric field E, potential φ, and conductivity σ, i.e., j = σE and E = ∇φ [16], we can obtain the following:

$$\nabla \cdot (\nabla \sigma \phi )=0,$$

(2)

2.2. GRC Calculation

Computing current and potential is a typical DC propagation problem [5] that is commonly solved via partial differential equation methods, including the moment method, the finite element method (FEM), the finite difference method, and the finite volume method. The simulation of current and potential distribution in DC propagation problems is commonly addressed using methods based on partial differential equations. These include the moment method [6], the finite element method (FEM) [7,8,9,17,18], the finite difference method [19], and the finite volume method [20]. Considering that these methods ignore the shunt effect of the overhead lines on the current and generate the solving misfit of the potential and neutral current, the resistance network (RN) method was developed [2,10]. It views the ground as a horizontal layered and a vertical layered medium to construct the equivalent resistance [21,22,23] and realizes the forward simulation of GRCs.

2.3. Discrete Grid

To address the inadequacy of the RN method in accurately simulating the propagation of GRC in a 3D geoelectric structure, a structured hexahedral unit-based RN was developed [2]. This type of structured grid (SG) is commonly used in numerical simulation methods such as the finite difference method [24,25] and FEM [26,27,28]. The grids are relatively simple and mostly hexahedral. However, achieving accurate representation of swift variations near the surface using conventional regular grids necessitates intricate discretization, resulting in escalated computational expenses attributed to the abundance of structured grids [29,30].

Due to flexible dispersion performance, the unstructured grid can dissect the target area in a targeted manner and avoid the transitional dissection and dispersion of non-target areas such as the structured grid. Unstructured grids, encompassing hexahedrons and having garnered significant attention in the past decade, find application in the finite volume method [31,32] and FEM [33,34]. A tetrahedron can be used to fit complex underground electrical structures, but the calculation accuracy is closely related to the mesh quality. The quality of the tetrahedrons is uncontrollable and can generate narrow elements, resulting in poor mesh quality [35,36].

Unstructured grids relying on hexahedrons pose several challenges, including the complexity of mesh setup, the risk of matrix ill-conditioning arising from unfavorable aspect ratios, and necessity for additional postprocessing efforts [37].

2.4. Contributions

We introduce a multiresolution grid to the node voltage network for a more efficient discretization of underground structures. This method offers the following advantages:

(1)

The proposed resistance network system does not require the continuity of regular elements to establish an equipotential surface for circulating currents on a face, and it maintains good symmetry, resulting in higher solving efficiency.

(2)

The grid is refined at the source point and near the surface to enhance accuracy. This helps avoid redundant elements caused by lateral extension in regular grid refinement, reduces the degrees of freedom, and improves the solving efficiency of the resistance network system.

(3)

For hanging nodes in a multiresolution grid’s discretized resistance network, we compared three different interpolation methods and determined that the ghost fine interpolation resistance value method is suitable for MR resistance network discretization.

3. 3D RN Method

3.1. Stereo NV-RN on an SG

The integral form of Equation (2) can be expressed discretely as

$${\displaystyle \sum (I_1+I_2)}=0,$$

(3)

where I₁ is a vector composed of currents on all current inflow nodes, and I₂ is a vector composed of currents on all currents flowing out of nodes. This equation shows that the law of conservation of charge can fully derive Kirchhoff’s current law in electrical circuits under steady conditions. Considering that most HVDC grounding currents are DC or quasi-DC with very low frequency, when DC current propagates in the soil, it is usually affected only by the resistive effect of the soil. This resistance effect is equivalent to the blocking effect of the resistance element in the circuit on the current passing through it.

Taking node p as an example, it is specified that the current flows from node p into the other six subsidiary nodes, and the following can be obtained:

$$Y_{N\times N}{U}_{N\times 1}=I_{N\times 1},$$

(4)

$$\begin{array}{c}{Y}_{N\times N}=[Y_{pq}]\hspace{1em}{U}_{N\times 1}=[U_p]\hspace{1em}{I}_{N\times 1}={[i_p]}_{N\times 1}\\ p=1,2,\dots ,N\hspace{1em}q=1,2,\dots ,N\end{array},$$

(5)

here, p takes values from 1 to N, and q ranges from 1 to N, where N is the total number of cell nodes. In the node admittance matrix Y, the matrix elements are self-admittance when p = q, and the matrix elements are mutual admittance when p ≠ q. When the q-th node is not an adjacent subordinate node of the p-th node, the mutual admittance is 0. Therefore, the node admittance matrix exhibits sparsity and symmetry, with the width of the matrix remaining fixed based on a determined number of elements and the sorting method.

Obtaining the resistance between nodes is a necessary condition for using the node element method. Each regular hexahedral soil block is regarded as an independent unit, and the center of the unit is used as the potential sampling point. For a linear geoelectric case, the potential gradient of the grid boundary surface adopts the central difference format:

$$\nabla {\phi }_{r+1/2}={\displaystyle \frac{{\phi }_{r+1}-{\phi }_r}{h_{r+1/2}}}+\mathrm{O}(\mathsf{\Delta }{h}^2)\hspace{1em}r=i,j,k,$$

(6)

where h_r+1/2 = (h_r + h_r+1)/2 represents the average length of the grid edges of two adjacent cells in the axial direction.

3.2. MR Grid

MR grids offer an uncomplicated method for sparse data representation, utilizing fewer voxels for areas exhibiting minimal variability and a greater number of voxels for regions characterized by high variability. Unlike the structured SG grid, the MR grid has spatial recursive properties. The coarse grid is composed of multiple subgrids, as follows:

$$N_{grid}={\displaystyle \sum _{i=1}^{N_s}{N}_x^{C{L}^{(i)}}\times N_y^{C{L}^{(i)}}\times N_z^{C{L}^{(i)}}},$$

(7)

where N_grid denotes the total count of geoelectric cells. The coarseness level is represented by the coarseness parameter CL⁽ⁱ⁾, i = 1, 2, …, N_s. Here, CL is positively related to the grid spacing and can be assigned values of 1, 2, 3, and the like, and the minimum coarseness level is 1, represented by CL⁽¹⁾ = 1. The subgrids are composed of structured grids with dimensions N_x^CL(i) × N_y^CL(i) × N_z^CL(i), and the multilevel recursion is formed after subdivision. The difference in the coarseness level between adjacent coarseness level grids is unique, so the MR grid can be established by using a basic refined grid and the variables [CL⁽ⁱ⁾, Sp⁽ⁱ⁾], where Sp represents the subgrid spacing, and [CL, Sp] = [1, 5; 2, 10] provides a straightforward instance of an MR grid with a recursion depth of 2, as depicted in Figure 3.

3.3. MR Stereo NV-RN

3.3.1. Hanging Nodes

The resistance value of the subgrids with the same CL is the same as that of the structured grid, and the mesh interface between adjacent CLs generates hanging nodes and hanging surfaces. Using Figure 4a as an example, the CL of cell 1 is larger than that of cells 2, 3, 4, and 5. Nodes n_l, n₂, n₃, and n₄ are so-called hanging nodes.

The hanging surface delineates the interface between coarse and fine meshes. In the process of electromagnetic flux numerical simulation based on MR grids, the fine grid surfels are obtained via the interpolation of coarse grid surfels, and fine mesh surfels are omitted from the system matrix. Using the same volumetric flux discretization as the structured grid, we obtain the following:

$$\int \partial V_{i,j,k}\left(\sigma (m)\nabla \phi \right)\cdot \overrightarrow{n}dS={\displaystyle \underset{S_{i,j,k}}{\int }(\sigma (m)\nabla \phi )}dS\approx {\displaystyle \sum _{r=i,j,k}{\displaystyle \sum _{\mathsf{\Delta }h=\pm \frac{1}{2}}{\left(\sigma (m)\nabla \phi \right)}_{r+\mathsf{\Delta }h}{w}_{r+\mathsf{\Delta }h}{S}_{r+\mathsf{\Delta }h}}},$$

(8)

where ω is the interpolation matrix.

When addressing the gradient operator at the junction of coarse and fine grids, using the x-direction refinement in Figure 4a as an example, the most natural method [38] is to approximate the gradient at the interface adjacent to the hierarchical grid as follows:

$${(\nabla \phi )}_{x2}\approx {\displaystyle \frac{{\phi }_2-{\phi }_1}{h_{2-1}}},$$

(9)

where h_2-1 is the distance between nodes in the x-axis of the coarse and refined grids. We use the grid spacing h_c of the coarse grid as the benchmark. Then, the finite difference of Equation (9) around the interface remains incomplete. Therefore, the formula may not even offer first-order accuracy.

Three types of grid spacing are defined in Figure 5: coarse spacing (CS), fine spacing (FS), and ghost fine (GF) spacing. The difference method is set to a coarse spacing (CS) grid (Figure 5a). We substitute Equation (9) into Equation (8), and the resistance coefficient between coarse and refined grids can be derived through surface interpolation.

3.3.2. Hanging Node Interpolation

We finally obtain the resistance relation between grids 1 and 2 as follows:

$$R_{i+1/2,j,k(1)}={\displaystyle \frac{3}{4}}\cdot {\displaystyle \frac{{\rho }_{i+1/2,j,k}\cdot h_c^x}{S_{i+1/2,j,k}}},$$

(10)

where ρ_i+1/2,j,k = W₁ρ_i,j,k + W₂ρ_i+1,j,k and W is the volume-weighted average operator. The 3D equivalent resistance connection is shown in Figure 6a.

The grid spacing size can be reduced to improve the difference accuracy. The resistance between the two nodes can be determined by connecting the resistance of the coarse and refined grids on either side of the interface. The resistance of the coarse grid is determined by the spacing size h_c, and for the refined grid, it is based on the spacing size h_f, designated as the FS grid (Figure 6b). The setting method for setting the 3D equivalent resistance, followed by the equation, is illustrated in Figure 6.

$$R_{i,j,k}={\displaystyle \frac{1}{2}}\cdot {\displaystyle \frac{{\rho }_{i,j,k}\cdot h_c^x}{S_{i,j,k}}}\hspace{1em}{R}_{i+1,j,k(1)}={\displaystyle \frac{1}{2}}\cdot {\displaystyle \frac{{\rho }_{i+1,j,k}\cdot h_f^x}{S_{i+1,j,k}}},$$

(11)

After connecting two resistors in series, we obtain the following:

$$R_{i+1/2,j,k}=R_{i,j,k}+R_{i+1,j,k(1)},$$

(12)

By referring to the difference method in [39], a higher-order precision approximation is used, and the Taylor series expansion directly at φ₂ is not considered, but a ghost point is introduced, and the potential gradient is obtained through second-order precision interpolation. The interface gradient is approximated as follows:

$${(\nabla \phi )}_{x2}\approx {\displaystyle \frac{{\phi }_1+{\phi }_2-2{\phi }_3}{3h_f^x}},$$

(13)

This seems to be the approximate value of the first-order accuracy O(h), but it was shown in [39] that the method actually offers second-order accuracy O(h²). We substituted Equation (13) into Equation (8) and obtained the resistance coefficient between coarse and refined grids through surface interpolation (Figure 5c), which is referred to as the GF. Finally, we obtained the relationship coefficients at the hanging nodes of adjacent coarse and fine grids in the Cartesian coordinate system to determine the resistance coefficient connected between grids 1 and 2 as follows:

$$R_{i+1/2,j,k(1)}=6\cdot {\displaystyle \frac{{\rho }_{i+1/2,j,k}\cdot h_f^x}{S_{i+1/2,j,k}}}$$

(14)

3.3.3. Boundary Conditions

Our network automatically truncates at the ground–air boundary to ensure that current does not flow from underground into the air (Figure 7). For the underground boundary, we did not use a complete uniform Dirichlet boundary. Instead, we added a resistance with a value lower than that of the boundary inside, simulating the effect of an infinite range on the current in the computation area.

To avoid difficulties in determining resistances on co-planar edge boundaries, our computational cells use the cell-centered method rather than the vertex-centered method. This approach limits our ability to compute surface voltages, providing only the voltage at a depth of half the cell’s vertical dimension underground. However, since the neutral point of the substation is buried underground, this cell-centered approach does not affect the calculation of the neutral point’s potential. Moreover, we believe that when the computational area is sufficiently large, the final results calculated via the cell-centered and vertex-centered methods are virtually identical. For smaller computational models, under laboratory conditions, truncated boundaries can be used to prevent current from flowing out of the model’s boundary area.

4. Numerical Examples

This section presents two synthetic model examples to evaluate the adaptability of the proposed model in addressing the two characteristics of traditional ground return: wide propagation range and complex propagation medium results. The homogeneous model case is used to illustrate the process of constructing the MR grid, comparing the applicability of the RN system under broad-area conditions. It also evaluates the improvements in accuracy and computational cost provided by the proposed RN method. The anomaly model case is used to demonstrate the adaptability of the proposed RN to complex electrical structures.

The grid discretization is performed using the multiresolution grid shown in Figure 3. For hanging nodes, we compared three methods of handling hanging nodes (Figure 5), with the corresponding resistances in the resistance network calculated using Equations (10), (12) and (14). The resistance network system, as shown in Equation (5), is preprocessed using the ILU method and solved iteratively with the Biconjugate gradient stabilized (BiCGStab) solver. The modeling and solving process is executed using MATLAB 2023a.

4.1. Homogenous Model

The size of the homogeneous half-space model is designed to be 98 km × 98 km × 98 km, with a resistivity of 100 Ω·m. Both the current entering and exiting the ground are 6250 A. The coordinates of the current exit and entry points, which are the positive and negative electrodes of the current source, are (0, −20, 0) km and (0, 20, 0) km, respectively. The coordinates of the grounding points at both ends of the conductor are (−25, 0, −25) km and (25, 0, 25) km. The positions of the current source and measurement points are shown in Figure 8.

4.1.1. Comparison of Solving Efficiency

The system matrices generated via the two RN methods are shown in Figure 9. The system matrix generated via K-RN is mainly composed of two diagonal matrices combined, which correspond to the behavior of simulating the propagation process of current using the Kirchhoff current and voltage rules, respectively. The proposed NV-RN generates a symmetric diagonal system matrix without losing any unit information.

Using BiCGStab as the solver and ILU as the preconditioner [40,41], the solution processes of the two methods are compared, as shown in Figure 10. When the ILU preconditioner is not used, both K-RN and NV-RN are difficult to solve; the residuals of K-RN do not decrease, and the residuals of NV-RN decrease slowly with iterations. With the ILU preconditioner, the performance of the GMRES solver is significantly improved. Processing K-RN with ILU required 1 h, 21 min, and 30 s, while processing NV-RN required only 4 min and 55 s. This indicates that the diagonal system matrix generated via NV-RN is more favorable for solving.

4.1.2. Mesh Discretization Efficiency

Two different mesh (

Table 1. Mesh generation of RN and FEM.

	NV-RN	KRN	MFEM-1	MFEM-2	MFEM-3	MFEM-4	MFEM-Bench
Mesh Type	Hexahedron				Tetrahedron
Ground NoE (1.0 × 106)	0.98	4.3	0.97	0.97	0.05	0.06	2.4
Line NoE			10,598	62,222	6996	61,191	117,481
Total NoE (1.0 × 106)	0.98	4.3	0.98	1.04	0.06	0.12	2.52
Run time (s)	23	127	666	1253	640	846	33,744

) based finite element methods were used to compare the accuracy of the RN method. M_FEM-1 and M_FEM-2 use regular hexahedral grids, while M_FEM-3, M_FEM-4, and M_FEM-5 use unstructured tetrahedral grids. Apart from the different grid types, these meshes also vary in density. Among them, M_FEM-5 uses an adaptive algorithm based on a posteriori error estimates to generate the mesh used for calculating the standard response values [42].

M_FEM-1 and M_FEM-2 employ regular hexahedral grids to discretize the underground structure, with a lateral element length of 1 km. In both M_FEM-1 and M_FEM-2, the number of underground discretization elements is 0.97 × 10⁶. In M_FEM-1, the number of mesh elements for discretizing the conductor is 10,598. In M_FEM-2, the conductor is more densely discretized, resulting in 62,222 mesh elements. The total number of mesh elements in M_FEM-1 and M_FEM-2 are 0.98 × 10⁶ and 1.04 × 10⁶, respectively. M_FEM-3 and M_FEM-4 use unstructured tetrahedral grids with varying densities to discretize the underground structure. In M_FEM-3, the total number of mesh elements is 0.06 × 10⁶, with 6996 elements for the conductor. In M_FEM-4, the total number of mesh elements is 0.12 × 10⁶, with 61,191 elements for the conductor. M_FEM-Bench represents the adaptive standard solution, where the underground structure is discretized using an unstructured tetrahedral grid and adaptive refinement at measurement points. The total number of mesh elements is 2.52 × 10⁶, with 117,481 elements for the conductor.

Comparing the solution results of K-RN and NV-RN, NV-RN achieved higher numerical accuracy with fewer elements (Figure 11). This is attributed to NV-RN’s strategy of constructing basic resistance network elements around the cell center. Unlike the K-RN method, NV-RN does not rely on the continuity of regular elements to establish an equipotential surface for a face’s circulating current. As a result, NV-RN can utilize MR grids or multi-level grids.

4.2. Abnormal Body Model

To verify the necessity of mesh refinement and the precision and effectiveness of MR mesh in 3D simulation, we consider the abnormal body model. We set the model symmetrically, as illustrated in the cross-sectional view of the main section at y = 0 m in Figure 12. The red cross denotes the point source with ±6250 A currents. The geological body has a surrounding resistivity of 100 Ω·m. The anomalous block area exhibits a distinct resistivity from the surrounding, measuring 1 Ω·m. The abnormal block’s dimensions are 11 km × 11 km × 5 km, with a burial depth of 1 km. SGs are regular hexahedral soil microelements of equal size, and the grid spacings are 90 m and 1000 m, denoted as SG₁ and SG₂, respectively. The MR grid is used to improve the grid resolution near the surface; the basic cell has an edge of 1000 m, and the fine grid spacing is 90 m, 200 m, and 330 m, labeled as MR₁, MR₂, and MR₃, respectively.

The gray part in Figure 12 represents the grid refinement on the surface near the abnormal body. The number of grid elements of each CL is shown in

Table 2. Grid subdivision strategies and solution efficiency for the 3D anomalous body model.

Grid Type	Range (X, km)	Depth (Z, m)	CL	Grid Spacing (km)	Grid Number	DoFs	Solving
							Iter	Time
MR1	−20.5–20.5	0–7	1	0.09	15661877	15816485	442	3 min 56 s
	−27.5–27.5	7–55	2	1	154608
MR2	−20.5–20.5	0–7	1	0.2	1470875	1625483	354	31 s
	−27.5–27.5	7–55	2	1	154608
MR3	−20.5–20.5	0–7	1	0.33	317709	472317	261	14 s
	−27.5–27.5	7–55	2	1	154608
SG1	−27.5–27.5	0–7		0.09	221445125	221445125	1654	20 h 3 min 29 s
SG2	−27.5–27.5	7–55	1	1	166375	166375	184	3.02 s

. We take the FEM calculation results as the benchmark solution, calculate the relative error (%) distribution of the MR grid and SG forward results, and conduct logarithmic transformation.

As shown in Figure 13a, the forward error of SG₂ above the anomalous body surface reaches 15% when the grid spacing is 1000 m, emphasizing the requirement for a finer grid in the target area. With the decrease in the grid spacing, the accuracy of the forward results above the anomalous body is improved by approximately 1.7%. However, the use of such a finely structured grid in the large area model increases the calculation time dramatically and places high demands on the CPU and memory.

The comparison of the forward results of MR₃, MR₂, and MR₁ is represented in Figure 13a–d. With the decrease in the grid step size, the relative error of ESP gradually decreases. When the grid spacing of MR₁ in the target area is reduced to 90 m, the forward accuracy above the anomalous body is less than 1.5%. Compared with SG₁, MR₁ ensures the accuracy of the target area and improves computational efficiency by 99.7%.

5. Conclusions

We propose a new MR grid-based RN forward calculation method. The purpose of this method is to improve the grid resolution of the RN method near the Earth’s surface to obtain higher computational efficiency and accuracy. In the implementation process, the MR grid is regarded as a space recursive grid. We addressed the equivalent resistance of three precisions at the hanging nodes and employed surfel interpolation to resolve the current exchange issue at the interface between the coarser and finer grids. The accuracy of the proposed algorithm was proven through comparisons with results from traditional structured grid-based computations in a typical layered-model scenario. On this basis, we used the MR grid to realize the efficient expansion of the electrical boundary. Compared with traditional proportional expansion, the computational efficiency of the MR grid was greatly improved under the premise of ensuring a calculation error of less than 1%. We showed through experiments that when there are abnormal bodies in the local superficial layer and the forward modeling accuracy is 1%, the calculation time of the proposed algorithm is reduced by 99.7% compared to that of the traditional RN. Given the practical necessity for higher accuracy in the surface potential at the substation junction, we employed the multi-resolution grid to improve the resolution in proximity to the substation and subsequently compared the results with FEM calculations. The experimental results show that compared with the traditional grid, the calculation results of the MR grid are closer to the real value.