Numerical Analysis of the Electrical Potential Calculation in a Transversely Isotropic Media

The electrical potential due to a point source of current supplied at the surface of a transversely isotropic medium is calculated using a finite element formulation. The finite and infinite elements are applied to calculate the potential for arbitrary electrical conductivity profiles. The accuracy of the scheme is checked against results obtainable using Chave's algorithm. Transversely isotropic | transverse conductivity | vertical conductivity | direct current | electrical potential | finite element | infinite element |


Introduction
The analytical solution for the electrical potential in transversely isotropic earth is represented as an integral equation form.Some authors calculated the electrical potential by solving the integral equation numerically.Sampaio [6] investigated cases in an isotropic medium.Stoyer and Wait [8] analyzed a two-layer model of isotropic medium, where the lower layer has an exponentially varying conductivity.Sato and Sampaio [7] considered a half-space whose resistivity varies as the power of an expression depending linearly on depth.However, their methods are limited to specific forms of conductivity.
In this paper, the electrical potential ) , ( z r φ in transversely isotropic media is chosen as a model.The potential resulting from a single current electrode located at the origin of a cylindrical coordinate system, which is axially symmetry satisfies a Boundary Value Problem (BVP) which consists of a partial differential equation: (ii) The electrical potential must approximate to zero at infinite distance, this satisfies where l σ is the conductivity in horizontal directions and v σ is the conductivity in vertical current flow.
A finite element scheme is developed to solve this BVP for any form of conductivity profile in a medium.The basic steps involved in solving this method are: (1) The formulation of a variational statement with an appropriate space of admissible function identified.
(2) The constructions of an approximation of the variational BVP on a finite dimensional subspace h H .
(3) The construction of a finite element mesh and piecewise-polynomial basis functions defined on the mesh.(4) Solving the system of equations.The first step of finite element method will be shown on the following section.

Variational Statement of the Boundary Value Problem
The equation (1.1) can be written as follows: In a cylindrical coordinate system which is axially symmetry Applying the divergence of F into equation (2.4), the BVP will consist of an equation: and the boundary conditions (1.2) and (1.3).The total weight residual error of equation (2.5) over a region Ω (see Figure 1) needs to be zero, that is: where n is an orthogonal vector to the boundary.
The integral on the boundary is calculated on three subregions of Ω ∂ , those are:

Z=depth
Figure 1: The domain Ω On the boundary 1 ∂Ω we have: (2.9)By adding the results (2.7), (2.8) and (2.9), the integral on the boundary becomes: Combining the equation (2.10) and (2.11), the total weight residual error over region Ω becomes: Hence the variational of the BVP can now be stated concisely in the following term: find where w(r, z) is a weight function and the space of admissible function

The Galerkin Approximation
A Galerkin approximation for solution φ is obtained by posing the variational problem on a finite dimensional subspace h H of admissible function and the weight function chosen as a function in h H . Specifically, we seek For the left hand side of (3.18) we have β are arbitrary, the equation (3.17) represents N equations to be satisfied by the i α defining h φ rather than the single equation.By choosing is substituted into equation (3.17) then the following equation is obtained: = is substituted into abov equation and we get: Now, we arrive at the system of N linear equation with N unknown j α as follows: The K is called the Stiffness matrix order (N x N), α is unknown and F is called load vector.Matrix K in equation (3.18) is invertible because is linearly independent, so the coefficient j α will be unique.
Therefore the Galerkin's approximation h φ of the solution φ is on the form: where j α is the solution of the linear system (3.18) and is the basis function of h H .

Finite Element Method
The FEM provides a general and systematic technique for construction the basis functions for the Galerkin's approximation.The procedure is started by divided Ω into two regions, those are a finite element mesh: and an infinite element mesh The finite element mesh is the set of master element shape functions (see [2]): For the region i Ω , the method as outlined in Zienkiewicz [9] is used for mapping the master element to infinite element with the following element shape functions: The integral (3.19) is evaluated by using The Gaussian Quadrature formula of order 3 in every element.The stiffness matrix K is constructed based on these results and the linear system (3.18) is solved to obtain the value of the electrical potential in every nodal point.

Interpolation Error
The actual solution φ of the boundary value problem will be interpolated by a finite element representation h φ which contains complete polynomial of degree 2. The interpolation error function can be expanded in a Taylor series about any interior point as follow (see [5] where 0 r is between r and r , 0 z is between z and z .Since h φ is interpolation of φ then E is zero at end point Assume that e h is the largest distance between two points in e Ω , then If this condition is applied and is assumed then we arrive on the following result: . So the interpolation error satisfies: , where 2 and h e is the largest distance between two points in e Ω .Similarly, , for h is the largest distance between two points in every element Ω or φ converge to h φ with rate of convergence h 2 .

The Accuracy of the Finite Element Approximation
The global matrix K in linear system (3.18) is a sparse matrix.If K is stored in a 2-dimension array, it is leading to an extravagant waste of computer storage and computer time.Therefore, the storage of K is modified by storing the non zero entries in K into an one dimensional array.If K is stored in a 2-dimension array, the code can only solve a problem in a mesh with 3600 nodes, but if K is stored in a one-dimension array, then the code can solve it for a mesh with 14400 nodes.By increased the number of nodes, the accuracy of the results will be improved.Tables 1, 2, 3 show the potential values and the enhancement of the relative error on the surface for distance 10 units, 20 units, 40 units and Table 4, 5, 6 show for distance 15 units, 30 units, and 60 units.The results are compared with the results obtained by using Chave's algorithms.The largest relative error in the range shown is less than 3% for the node located around the source.

Conclusion
The verification confirms the accuracy of the finite element code.Although in most practical problems only the potential on the surface of the earth is measurable, our scheme does give the potential at every node in the computational domain.Further research will investigate the application of the code to the inverse problem of determining various conductivity profiles based on surface measurements.
The vertical component of the current density must be zero at ground surface, this satisfies

7 )the boundary 2 ∂Ω
On , the boundary condition (1.1) is applied and the following result is obtained by choosing w integral in the domain Ω is obtained as follow: h w in(3.16) is substituted into equation (3.15), then the following result is obtained: mapping a nine-point isoparametric master element Ω with coordinate ) , ( η ξ in Cartesian and the origin (0, 0) is located at the center of Ω .The master element is transformed into the r z-plane by using the following coordinate transformation:

Table 1 .
The potential and the relative-error on the surface with AB=10 units and the conductivity

Table 2 .
The potential and the relative-error on the surface with AB=20 units and the conductivity

Table 3 .
The potential and the relative-error on the surface with AB=40 units and the conductivity

Table 4 .
The potential and the relative-error on the surface with AB=15 units and the conductivity

Table 5 .
The potential and the relative-error on the surface with AB=30 and the conductivity

Table 6 .
The potential and the relative-error on the surface with AB=60 units and the conductivity