Comparison between Dirichlet boundary condition and mixed boundary condition in resistivity tomography through finite-element simulation

2.1. Derivation of the finite-element equation of 2D geoelectric field with a point power source

Let A(x₀, y₀, z₀) be the coordinates of the point power source in an infinite medium whose resistivity is ρ(x, y, z), and I be the current intensity of the point power source. In the coordinate system, the x and y axes were arranged along the direction of the geological body while the z axis was vertically downward. Then, the potential function U of stable current field satisfies the following differential equation:

$\nabla \cdot \left( \sigma \cdot \nabla U \right)=-I\delta \left( r-{{r}_{A}} \right)$   (1)

where σ is the dielectric conductivity; δ is the Dirac equation; U is the potential at any point on the surface or in the medium; I is the supply current intensity; r is the field point radius vector; r_A is the source point radius vector.

Since

$\nabla \cdot \left[ \sigma \left( x,z \right)\nabla U\left( x,y,z \right) \right]=-I\delta \left( x-{{x}_{0}} \right)\delta \left( y-{{y}_{0}} \right)\delta \left( z-{{z}_{0}} \right)$    (2)

Applying Fourier transform to the basic differential equation of the steady current field along the y axis, we have:

$\phi \left( x,k,z \right)=\int_{0}^{\infty }{U\left( x,y,z \right)\cos \left( k\cdot y \right)dy}$   (3)

Transforming the potential U in space (x, y, z) to the generalized potential

$\nabla \cdot \left[ \sigma \left( x,z \right)\nabla \overline{U}\left( x,{{K}_{y}},z \right) \right]-K_{y}^{2}\sigma \left( x,z \right)\overline{U}\left( x,{{K}_{y}},z \right)=-\frac{1}{2}\delta \left( x-{{x}_{0}} \right)\delta \left( z-{{z}_{0}} \right)$   (4)

where K_yis the spatial wavenumber; δ is the Dirac function.

2.2. Delimited boundary condition

(1) Dirichlet boundary condition

The unique solution of equation (4) can be determined only under the definite condition of the boundary. In this paper, the Dirichlet boundary condition is selected, that is, the potential value φ is given on the boundary, because the boundaries must be far away from the field source and the non-uniform body. The Dirichlet boundary condition can be expressed as [21~30]:

$U\left( x,y,z \right)\left| _{\Gamma } \right.=\varphi \left( x,y,z \right)$   (5)

As the value of φ(x, y, z) is difficult to be obtained by formulas, one of the following two conditions is often adopted:

①

During the y direction Fourier transform of the U in equation (4), the cosine transform was applied with an integration interval between 0 and ∞, for potential U(x, y, z) is symmetric with respect to the plane xz, i.e. U(x, y, z)=U(x, -y, z).

$\widetilde{U}\left( x,k,z \right)=\int_{0}^{\infty }{U\left( x,y,z \right)\cos kydy}$     (6)

Considering the integral property of δ function, the δ function in equation (4) satisfies the following condition:

$\frac{1}{2}\delta \left( A \right)=\frac{1}{2}\delta \left( {{x}_{A}} \right)\delta \left( {{z}_{A}} \right)=\int_{0}^{\infty }{\delta \left( {{x}_{A}} \right)}\delta \left( {{y}_{A}} \right)\delta \left( {{z}_{A}} \right)\cos kydy$    (7)

Substituting equations (6) and (7) into equation (4), we have:

$\frac{\partial }{\partial x}\left( \sigma \frac{\partial \widetilde{U}}{\partial x} \right)+\frac{\partial }{\partial z}\left( \sigma \frac{\partial \widetilde{U}}{\partial z} \right)-{{k}^{2}}\sigma \widetilde{U}=-I\delta \left( A \right)$       (8)

For $\frac{\partial \varphi}{\partial n}=0$ on the boundary, the normal lies in the plane xz and has nothing to do with y. Thus, we have:

$F\left[ {{\left. \frac{\partial U}{\partial n} \right|}_{\Gamma }} \right]={{\left. \frac{\partial \varphi }{\partial n} \right|}_{\Gamma }}=0$   (9)

Equation (5) can be transformed into:

$\widetilde{U}\text{=}F\left[ {{\left. \frac{\partial U}{\partial n} \right|}_{\Gamma }} \right]=\int_{0}^{\infty }{\frac{c}{\sqrt{{{r}^{2}}+{{y}^{2}}}}}\cos kydy=c{{K}_{0}}\left( kr \right)$    (10)

where r falls within the profile passing through point U(0,0,0) and perpendicular to y; Γ_∞ boundary is the distance from a point to U(0,0,0); K₀ is the zero-order modified Bessel function of the second kind. The differential quotient of K₀(x) is:

$\frac{d}{dx}{{K}_{0}}\left( x \right)=-{{K}_{1}}\left( x \right)$   (11)

K₁ is the is the first-order modified Bessel function of the second kind. Taking the derivative of equation (10) in the normal n, we have:

$\frac{\partial \widetilde{U}}{\partial n}=\frac{\partial \widetilde{U}}{\partial r}\frac{\partial r}{\partial n}=-ck{{K}_{1}}\left( kr \right)\cos \left( r,n \right)$    (12)

From equations (11) and (12), we have:

$\frac{\partial \widetilde{U}}{\partial n}+{{\left. k\frac{{{K}_{1}}\left( kr \right)}{{{K}_{0}}\left( kr \right)} \right|}_{\Gamma }}=0$  (13)

Then, the generic function can be established as:

$I\left( \widetilde{U} \right)\text{=}\int\limits_{\Omega }{\left[ \frac{\sigma }{2}{{\left( \nabla \widetilde{U} \right)}^{2}}+\frac{1}{2}{{k}^{2}}\sigma {{\widetilde{U}}^{2}}-Ix\left( A \right)\widetilde{U} \right]}d\Omega $   (14)

where the area Ω is surrounded by the surface boundary Γ_s and the infinity boundary Γ_∞. The corresponding variational problem can be expressed as:

$\begin{align}  & \delta I\left( \widetilde{U} \right)\text{=}\int\limits_{\Omega }{\left[ \sigma \nabla \widetilde{U}\cdot \nabla \delta \widetilde{U}+{{k}^{2}}\sigma \widetilde{U}\delta \widetilde{U}-I\delta \left( A \right)\delta \widetilde{U} \right]} \\ & =\int\limits_{\Omega }{\nabla \cdot \left( \sigma \widetilde{U}\delta \widetilde{U} \right)d\Omega +\int\limits_{\Omega }{\left[ -\nabla \cdot \left( \sigma \nabla \widetilde{U} \right)+{{k}^{2}}\sigma \widetilde{U}-I\delta \left( A \right) \right]}}\delta \widetilde{U}d\Omega  \\\end{align}$ (15)

Substituting equations (9) and (13) into equation (15), we have:

$\begin{align}  & F\left( \widetilde{U} \right)\text{=}\int_{\Omega }{\left[ \frac{\sigma }{2}{{\left( \nabla \widetilde{U} \right)}^{2}}+\frac{1}{2}{{k}^{2}}\sigma {{\widetilde{U}}^{2}}-I\delta \left( A \right)\widetilde{U} \right]}d\Omega  \\ & +\int_{\Gamma }{\sigma \frac{k{{K}_{1}}\left( kr \right)}{k{{K}_{0}}\left( kr \right)}\cos \left( r,n \right){{\widetilde{U}}^{2}}d\Gamma } \\\end{align}$   (16)

The $\tilde{U}$ of the split node can be obtained by solving the variational problem of equation (16) through finite-element method. The $\tilde{U}$ corresponds to a specific wave number k. The potential U can be obtained via inverse Fourier transform of a group of $\tilde{U}$ values corresponding to different k values.

(2) Mixed boundary condition

Mixed boundary condition can be expressed as:

${{\left. \left( \frac{\partial U}{\partial n}+AU \right) \right|}_{\Gamma }}=\varphi \left( x,y,z \right)$   (17)

where A is a known positive number.

Therefore, the boundary condition of equation (4) can be transformed into the following boundary value problem [31~36]:

$\left. \begin{matrix}   \frac{\partial }{\partial x}\left( \sigma \frac{\partial \overline{U}}{\partial x} \right)+\frac{\partial }{\partial z}\left( \sigma \frac{\partial \overline{U}}{\partial z} \right)-{{\lambda }^{2}}\sigma \overline{U}={{f}_{1}}  \\   {{\left. \frac{\partial \overline{U}}{\partial n} \right|}_{{{\Gamma }_{1}}}}=0  \\   {{\left. \left( \frac{\partial \overline{U}}{\partial n}+A\overline{U} \right) \right|}_{{{\Gamma }_{2}}}}=0  \\\end{matrix} \right\}$   (18)

where

${{f}_{1}}=-\frac{1}{2}\sum\limits_{k-1}^{n}{{{I}_{k}}\delta \left( x-{{x}_{k}},z-{{z}_{k}} \right)}$   (19)

$A\text{=}\lambda \frac{{{k}_{1}}\left( \lambda r \right)}{{{k}_{0}}\left( \lambda r \right)}\cos \theta \left( r\ne y \right)$   (20)

where Γ₁ is the ground boundary; Γ₂ is the remaining boundary; I_k is the current intensity of the k-th point current source; k₁(λr) is the first-order modified Bessel function; k₀(λr) is the zero-order modified Bessel function; θ is the angle between the radius vector r from the power point to the boundary point and the outer normal n of boundary.

Since the A in equation (20) is a constant, the following can be obtained through the Fourier transform of the U in equation (17):

$\varphi \left( x,y,z \right)=A{{k}_{0}}\left( \lambda r \right)$    (21)

where k₀(λr) is the zero-order modified Bessel function, and

$\frac{\partial \varphi }{\partial n}=-A\lambda {{k}_{1}}\left( \lambda r \right)\cos \theta $  (22)

where k₁(λr) is the first-order modified Bessel function; θ is the angle between the radius vector r and the outer normal n.

According to equation (21), we have:

$A\text{=}\frac{\varphi \left( x,y,z \right)}{{{k}_{0}}\left( \lambda r \right)}$  (23)

Substituting equation (23) into equation (22), we have:

$\frac{\partial \varphi }{\partial n}+\lambda \frac{k{}_{1}\left( \lambda \cdot r \right)}{{{k}_{0}}\left( \lambda \cdot r \right)}\varphi \cos \theta =0$ or $\frac{\partial \varphi }{\partial n}+\alpha \varphi \cos \theta =0$     (24)

Substituting the mixed boundary condition (24) after the variation of equation (17) into the condition, equation (18) can be converted into an equivalent variational problem [37~42]:

\[W\left( \overline{U} \right)=\iint\limits_{D}{\left\{ \begin{align}  & \sigma \left[ {{\left( \frac{\partial \overline{U}}{\partial x} \right)}^{2}}+{{\left( \frac{\partial \overline{U}}{\partial z} \right)}^{2}} \right]+ \\ & \sigma K_{y}^{2}{{\overline{U}}^{2}}-I\delta \left( x-{{x}_{0}} \right)\delta \left( z-{{z}_{0}} \right)\overline{U} \\\end{align} \right\}dxdz+\frac{1}{2}\int\limits_{T}{\sigma a\overline{U}ds}}\] (25)

where I is the current intensity; Γ is the boundary of D, the area of the 2D geoelectric profile; ds is the integral variable of the integrals along the boundary, that is:

$a={{k}_{y}}\frac{{{k}_{1}}\left( {{k}_{y}}r \right)}{{{k}_{0}}\left( {{k}_{y}}r \right)}\cos \theta $  (26)

where k₀ and k₁ are zero-order and first-order modified Bessel functions, respectively; θ is the angle between the radius vector r from the power point to the boundary point and the outer normal n of boundary. Then, the real potential U can be obtained through the following steps: discretizing the target area; solving the generalized potential $\bar{U}$ by the finite-element method; applying inverse Fourier transform to the generalized potential $\bar{U}$.