Galerkin boundary element method for the forward problem of ERT

Flow Measurement and Instrumentation 21 (2010) 172–177

Contents lists available at ScienceDirect

Flow Measurement and Instrumentation journal homepage: www.elsevier.com/locate/flowmeasinst

Galerkin boundary element method for the forward problem of ERT Yaoyuan Xu, Feng Dong ∗ Tianjin Key Laboratory of Process Measurement and Control, School of Electrical Engineering and Automation, Tianjin University, Tianjin 300072, China

article

info

Article history: Received 19 June 2009 Received in revised form 15 December 2009 Accepted 22 December 2009 Keywords: Galerkin boundary element method Electrical resistance tomography Forward problem Complete electrode model

abstract Electrical resistance tomography (ERT) is a novel non-intrusive method that determines the crosssectional electric conductivity distribution. Currently, the ERT forward problem solvers are mainly based on a finite element method, which require a large number of mesh elements. A Galerkin boundary element method (GBEM) is presented for the ERT forward problem in both 2D and 3D. The complete electrode model is adopted. The contact resistivity of electrodes is measured by rectangular cells of different lengths, as well as the conductivity of the background medium. GBEM formulations solving the multi-domain problems are derived. Numerical results show that the method has a better performance on precision compared by both the finite element method and the collocation boundary element method. Crown Copyright © 2010 Published by Elsevier Ltd. All rights reserved.

1. Introduction Electrical resistance tomography (ERT) is a novel non-intrusive method that determines the conductivities of different materials and also their spatial distribution by applying currents to a pair of adjacent boundary electrodes and measuring the resulting voltages at the electrodes. Compared to computerized tomography (CT), which is a standard tool in medical diagnostics and nondestructive testing of materials, ERT has many advantages in cost, speed and safety. Examples of ERT applications include two-phase flow measuring [1,2], medical diagnosis [3,4], and geophysical prospecting [5,6]. A similar technique is electrical impedance tomography (EIT), which also considers permittivity. The two techniques have common governing equations and operation processes, so they share most of their algorithms. As a non-linear problem, ERT always adopts an iterative strategy to get a relatively acceptable result. Thus an accurate and fast forward solver is needed. Currently, ERT forward problem solvers are mainly based on the finite element method (FEM), for typical modeling of the problem see literature [7,8]. However, the FEM requires a large number of mesh elements to get an acceptable solution and the system matrix is always huge. An improvement for the FEM solver is reported in [9]. Brebbia [10] states that the FEM is inefficient and cumbersome to use in many engineering applications. The boundary element method (BEM) transforms the domain problem into a boundary problem. It does not need to discretize the entire domain which makes the process much easier, thus requires fewer unknowns and saves memory space and CPU time. Recently, there has been a growing number of articles addressed

on determining the location and shape of inclusions [11,12]. Hence the BEM is appreciated in solving both the forward and inverse problems. A 3-D boundary element based forward problem solver has been reported in [13], however, they did not present solutions based on a completed electrode model [4,14]. In this paper, the Galerkin Boundary Element Method (GBEM) is applied to solve the ERT forward problem based on a completed electrode model in 2D and 3D. The GBEM has the advantages of higher accuracy and theoretical error analysis compared to the collocation boundary element method (CBEM). Although the Galerkin approach requires additional calculations, the majority of the CPU resource cost is by solving the linear equations and the additional integrals can be ignored. GBEM formulations have been proposed in [15] for a 2D electrostatic field. The symmetric GBEM was investigated in elasticity in [16]. The methods are extended to multi-domain problems [17,18]. For methods to evaluate boundary integrals, see [19,20]. The major difference between the method in this paper and in the literature above lies in boundary conditions. The boundary conditions in the literature [16–20] are generally of the Dirichlet or Neumann type, however, for the ERT forward problem they are much more complex. Thus there is a need to derive GBEM formulations based on a completed electrode model. The paper is organized as follows. Section 2 introduces the mathematical model for multi-domain problems. Section 3 presents GBEM formulations based on the complete electrode model and organization to form the system matrix. In Section 4, experiments are carried out to illustrate the effectiveness of the GBEM which is compared by the FEM and CBEM. 2. The mathematical model

∗

Corresponding author. Tel.: +86 22 27890276; fax: +86 22 27892055. E-mail address: [email protected] (F. Dong).

The forward problem of ERT is to evaluate the electric field distribution when the geometric distribution and physical properties

0955-5986/$ – see front matter Crown Copyright © 2010 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.flowmeasinst.2009.12.004

Y. Xu, F. Dong / Flow Measurement and Instrumentation 21 (2010) 172–177

173

where z and z 0 are the field point and source point; ε is a geometric coefficient relating to the free term of the singular kernel, for a smooth boundary ε = 0.5; u∗ and q∗ are the fundamental solution (also called the free space Green’s function) and its normal derivative, respectively. For a two dimensional problem, u∗ = − q∗ =

1

ln kz − z 0 k

2π


z − z0 · n ∂ u∗ =− ∂n 2π kz − z 0 k2

while in three dimensional coordinates, u∗ =

Fig. 1. An illustrative distribution of different media in cross-section.

of the domain are determined. Suppose an ERT sensor has L electrodes, and M objects located in sensing field as Fig. 1 shows, denote by Dm (m = 1, 2, . . . M ) the domains covered by these objects and assuming that Γm are their boundaries and σm are conductivities. D0 is the remaining region which indicates the homogeneous background with its conductivity σ0 . Let Γ0 = ∂ D0 , ΓEl (l = 1, 2, . . . , L) the partial boundary covered by the l-th metallic electrode, and ΓG the remaining insulating boundary. The electric potentials um in Dm and u0 in D0 , are governed by the Laplace equations,

 ∆ u0 = 0 ∆ um = 0

(1)

In order to take into account the contact resistance between the background medium and electrodes, the complete electrode model [4,14] is adopted.

Z ΓEl

on ΓEl , l = 1, 2, 3, . . . , L

(2)

σ0 q0 ds = Il

(3)

∂u

where q0 = ∂ n0 , n is the outward unit normal to the boundary which is supposed to walk in counterclockwise direction; Z is the contact resistivity, Il are the currents applied to the electrode, and Vl are the corresponding potentials. Due to the law of the conservation of charge, L X


z0 − z · n ∂ u∗ = . ∂n 4π kz − z 0 k3

Hence for a problem arising from the ERT, one may get the following boundary integral equation.

ε0 u0 +

=

Z

Il = 0.

(4)

u0 q∗ ds +

ΓG

L Z X

ΓEl

l =1

M Z X

Γm

m=1

in D0 in Dm , m = 1, 2, . . . M .

u0 + Z σ0 q0 = Vl

q∗ =

1 4π kz − z 0 k

q0 u∗ ds +

u0 q∗ ds +

m=1

L Z X

ΓEl

l =1

M Z X

Γm

u0 q∗ ds

q0 u∗ ds.

(9)

When dealing with problems where conductive objects are encountered, there is a need to generate an equation for each subdomain.

ε m um −

Z Γm

um q∗ ds = −

Z Γm

qm u∗ ds,

m = 1, 2, . . . , M .

(10)

Here the negative signs are ascribed to the unification of direction of the fluxes to away from the background domain. The conductive objects allow currents to flow through, thus the current densities are unknown. In the cause of computational efficiency, we attempt to avoid the calculation of current densities and reduce the multi-domain problem to a single-domain one. Since only the voltages on electrodes are measured and then used to image, the absence of these current densities does not hamper the process of solving the inverse problems. Furthermore, the fluxes on electrodes are expressed by Vi and u0 .

l=1

Because the solution of the Neumann boundary value problem above is unique up to an additive constant, a choice of ground is needed. L X

Vl = 0.

ε 0 u0 +

=

(5)

∂ u0 σ0 ∂n

Γm

∂ um = σm ∂n

(6)

Γm

u0 |Γm = um |Γm .

(7)

The (potential) boundary integral equation for the Laplace equation in some domain D with its boundary Γ can be derived by the Green–Gauss theorem [10]:

εu z +

Z

u (z ) q∗ z 0 , z ds =




Γ

Z

q (z ) u∗ z 0 , z ds




Γ

u0 q∗ ds +

L Z X

ΓEl

l=1

L Z X

u0 q∗ ds +

M X

Z cm

m=1

Γm

V l − u0 ∗ u ds Z σ0

ΓEl

u0 q∗ ds

(11)

where

ε0 =



ε0 + εm · κ ε0

z on Γm z on Γ0

cm = 1 − κ

κ = σm /σ0 . Rewrite Eq. (11) so that the unknown variables are on the left side, we get

3. Galerkin boundary element methods


0

ΓG

l=1

l=1

Finally, the current densities that flow into heterogeneous media and also the potentials are continuous.

Z

(8)

ε u0 + 0

+

Z

∗

ΓG

M X m=1

u0 q ds +

l=1

Z cm

L Z X

Γm

u0 q∗ ds −

 ΓEl

u0

∗

q +

u∗ Z σ0

 ds

Z L X Vl u∗ ds = 0. Z σ 0 ΓEl l=1

(12)

174

Y. Xu, F. Dong / Flow Measurement and Instrumentation 21 (2010) 172–177

In order to solve this problem, Eq. (12) should be discretized by expanding the boundary potential distribution in a set of base functions. u0 (z ) =

N X

αi φi (z )

(13)

i=1

where N is the number of nodes. And following [8], potentials on electrodes are expressed by V =

L−1 X

β l nl

(14)

l =1

where n1 = [1, − 1, 0, . . . , 0]T , n2 = [1, 0, −1, 0, . . . , 0]T et al. are RL×1 vectors. The choice of nl ensures that the condition (5) satisfied automatically. The potentials on electrodes are: V1 = PL−1 βl , Vl+1 = −βl , l = 1, 2, 3, . . . , L − 1. The classic collocation l boundary element method attempts to satisfy Eq. (12) at a discrete number of points, but ignores other points along the boundaries (they are not satisfied generally). In order to fulfil the requirement of error function to be zero on the entire boundary, the weighted residual technique may be used. When the weighting functions selected are the same as base functions in Eq. (13), the Galerkin formulation is fulfilled.

αi

Z

0

Γ0

ε 0 φi φj ds + αi

+ αi

L Z X

+ αi

M X

−

βl Z σ0

φj φj

Z

L

XZ Γ0

l =1

ΓG

φi q∗ dsds0

φi q +

Z

Γ0

m =1

φj

Z



ΓEl

Z cm

Γ0

Z

φj

Γ0

l=1

Z

Γm

ΓEl

u∗

∗

0

φi q∗ dsds0 nl u∗ dsds0 = 0.

(15)




C x=0

(16)




E x=I

(17)

where D(i, j) = −

Z

Z

ΓE1

φi ds −

i = 1, 2, 3, . . . , N , E(i, j) =



Z

+

!

Z ΓE (l+1)

i 6= j ej+1 Z

φi ds

(18)

j = 1, 2, 3, . . . , L − 1

 e1   , Z  e1



C . E

(21)

And linear equations to determine the coefficients can be written as 0 . I

(22)

According to ERT procedures, multiple injections of currents are required which signifies a family of I. Thus a QR decomposition using Householder reflections is adopted to solve Eq. (22). This method is good in numerical stability, the complexity of the algorithm for a square matrix with size n is O(n3 ). It is a pity that the system matrix A is non-symmetric, although many studies [16–18] successfully figured out the symmetric Galerkin boundary element methods suitable for solving problems governed by the Laplace equation, which suggests using the potential boundary equation on the Dirichlet boundary and the flux equation on the Neumann boundary. The key to achieve a symmetric GBEM is to get symmetric kernels. The complete electrode model in ERT contains integral expressions and also the Robin condition, which make it hard to get a symmetric system matrix. Anyway, the Galerkin method provides a powerful tool for the problem. 4. Experiments and results 4.1. Practical experiments

N ×(L−1)

where B is a R matrix, and C belongs to R . x = (α, β)T is a coefficient vector to be determined, where α = (α1 , α2 , α3 , . . . , αN ) and β = (β1 , β2 , β3 , . . . , βL−1 ). Since there are N + L−1 unknowns, more equations are required, which can be derived from Eqs. (3) and (4).

1

A=

B D

dsds

N ×N

D



 

Quadratic elements are adopted for 2D problems, the natural coordinate locations of nodes are (−0.5, 0, 0.5). Exact solutions of the element integral using quadratic elements have been presented in [21]. This work can help to evaluate the inner integrals in Eq. (15), however, it is very hard to derive further analytical expressions for the outer integrals. Thus gauss numerical integrations are performed for the outer integrals (using 6 gauss points). While in 3D coordinates, triangle elements are used to mesh the boundary, see Fig. 2. The integrations are fulfilled numerically (7 points) for quadratic elements case by the method in [22]. Then it is easy to express Eq. (15) in a matrix form, B

where ej is the length of the j-th electrode in 2D and the area in 3D. Similar formulations, which discretize Eqs. (3) and (4) for the FEM, can be found in [8]. Thus the system matrix takes the form

Ax =



Z σ0

Fig. 2. An illustration of triangle boundary mesh in 3D(coarse).

,

i, j = 1, 2, 3, . . . , L − 1

(19)

i=j

I = (I1 − I2 , I1 − I3 , . . . , I1 − IL )T

(20)

In order to illustrate the effectiveness of GBEM, a set of experiments is carried out. As Fig. 3 shows, 16 electrodes, each of which is 1 cm by width, 0.1 cm by thickness and 3 cm by height, are distributed on the boundary of a round area with its diameter 12.5 cm in equal intervals. The equipment is filled with Na2 SO4 solution. An adjacent mode is adopted. 5 mA 60 kHz AC current is applied to electrode 1 and electrode 16 is grounded. The conductivity of the solution and the contact resistance are measured using the rectangular cells shown in Fig. 4. The cells are (1.354, 1.590, 1.972, 2.210, 2.510, 2.852) cm in length. The resistances of the solution filled up different cells are as shown in Fig. 5. r0 and slope K can be calculated using the least squares error method. r0 is the sum of the contact resistances of a pair of electrodes, and the contact resistance of each electrode may be acquired by organizing and planning different pairs of measuring electrode pairs. Actually, the electrodes are made in the same conditions. Thus one may assume that the electrodes have

Y. Xu, F. Dong / Flow Measurement and Instrumentation 21 (2010) 172–177

175

Fig. 5. Resistance as a function of cell length.

Fig. 3. ERT sensor.

Fig. 6. Results by the GBEM compared with measurements.

Fig. 4. Rectangular cells in different lengths.

the same electrode contact resistance. The slope of the slope of the line is related to the conductivity.

σ =

1 AE K

(23)

where AE is the area of each electrode. Finally, the measured conductivity of the background material is 662 µS/cm. And the contact resistivity is 57.8 /cm2 . Thus one may use these parameters to simulate the boundary voltages. Fig. 6 compares the voltages on the electrodes obtained by practical experiments and numerical methods using GBEM (975 unknowns). As electrical properties of the electrodes may not completely agree with each other and noise exists in practical measurements, the measured voltages fluctuate within some range. Mean values (indicated in circles) are also shown in Fig. 6. The maximal absolute relative error is less than 5%. The numerical result could well simulate the practical measurements, benefitting from the complete electrode model. 4.2. Numerical experiments This part presents the performance of the GBEM as the mesh size grows, compared with the FEM and CBEM. The parameters

are the same as those in part 4.1 for a 2D simulation. In the 3D condition, the height of the cylinder is 10 cm. Because of complexities in the structure and the boundary conditions of the problem, it is hard to figure out an analytical solution. However, numerical results gradually approximate to it as the mesh grows. In order to illustrate how numerical methods approach the exact solution, the norm of potentials on electrodes is investigated. Fig. 7 shows the norms solved by different methods versus the size of unknowns (the x-axis for the FEM is on the top of the figure) in 2D. Fig. 8 shows the norms versus the size in 3D. It is obvious that the differences among the results provided by the three methods reduces as the mesh size grows in Figs. 7 and 8. It can be neglected when the mesh size is big enough, and thus the results could be regard as an ‘exact’ solution. Results show that the GBEM needs the smallest size of unknowns for a special error tolerance. To satisfy a requirement of the relative error of 1% in 2D, the size of unknowns should be at least 17 thousand for the FEM, and 399 for the CBEM, but 207 for the GBEM. While in 3D the required size is 70 thousand for the FEM, 13 thousand for the CBEM, but 4 thousand for the GBEM. The potentials on electrodes obtained in 3D are smaller than those in 2D. As the QR decomposition is adopted the complexity of which is O(n3 ), a decrease in the size of unknowns means an obvious saving in computer time. Fig. 9 shows the computer time of a QR decomposition in MATLAB. The method is quite effective for small size of unknowns but not very suitable for large scales. The required size of the GBEM is small, thus QR decomposition can be used.

176

Y. Xu, F. Dong / Flow Measurement and Instrumentation 21 (2010) 172–177

Fig. 10. Differential voltages on electrode pairs. Fig. 7. Norm of potentials on electrodes in 2D.

(160 elements on electrodes, 160 on ΓG , 30 on Γ1 and 40 on Γ2 ). The corresponding differential voltages are shown in Fig. 10. As the adjacent current injection mode is adopted, the boundary differential voltages on electrode pairs keep U shape distributions. When σ2 becomes higher, the mean conductivity of the crosssection becomes smaller and attracts current to travel through D2 . Thus the voltages on the electrodes are lower, and also the differential voltages considering that object 2 is close to the injection electrodes and vice versa. Potentials on D1 are approximately equipotential due to its high conductivity, which causes the least differential voltages on the nearby electrode pairs. 5. Conclusions

Fig. 8. Norm of potentials on electrodes in 3D.

ERT is a novel non-intrusive method that determines the electric conductivity distribution of a cross-section. Currently, ERT forward problem solvers are mainly based on the finite element method, which requires a large number of mesh elements. In this paper, a Galerkin boundary element method is presented for the ERT forward problem in multiply-connected domains with piecewise constant conductivity distributions in 2D and 3D. The complete electrode model is adopted. GBEM formulations solving the multi-domain problems are presented. A set of experiments are carried out to illustration the effectiveness of the GBEM. The values calculated agree with measurements with a maximal absolute relative error less than 5%. Results shows that the GBEM has considerable precision compared with the FEM with a huge number of mesh elements, while its mesh size is much smaller. Also the GBEM is superior to the CBEM in the sense of precision. Acknowledgements The author appreciates the support from National Natural Science Foundation Committee of China (No. 50776063), Natural Science Foundation of Tianjin, China (No. 08JCZDJC17700) and Program for New Century Excellent Talents in University of China (No. NCET-06-0230).

Fig. 9. Computer time of a QR decomposition.

4.3. Instance of the GBEM as an analyzing tool This part gives examples of solving multi-domain forward problems in 2D, where the GBEM acts as an analyzing tool. Suppose two objects are located on D1 and D2 in Fig. 1, the background conductivity σ0 and contact resistance are as measured in Section 4.1. σ1 is assumed to be 10 times as σ0 , while the multiple for σ2 is optional in (0, 0.5, 1, 10). The boundary is meshed by 490 elements

References [1] Dickin F, Wang M. Electrical resistance tomography for process applications. Measurement Science & Technology 1996;7(3):247–60. [2] Lucas GP, Cory J, Waterfall RC, Loh WW, Dickin FJ. Measurement of the solids volume fraction and velocity distributions in solids-liquid flows using dual-plane electrical resistance tomography. Flow Measurement and Instrumentation 1999;10(4):249–58. [3] Barber CC, Brown BH. Imaging spatial distributions of resistivity using applied potential tomography. Electronics Letters 1983;19(22):933–5. [4] Cheng KS, Isaacson D, Newell JC, Gisser DG. Electrode models for electric current computed tomography. IEEE Transactions on Biomedical Engineering 1989;36(9):918–24.

Y. Xu, F. Dong / Flow Measurement and Instrumentation 21 (2010) 172–177 [5] Daily W, Owen E. Cross-borehole resistivity tomography. Geophysics 1991; 56(8):1228–35. [6] Koh CS, Kim MK, Jung HK, Hahn SY, Suh BS. Electric resistivity tomography for geophysical inverse problems. IEEE Transactions on Magnetics 1997;33(2): 1852–5. [7] Hua P, Eung JW, Webster JG, Tompkins WJ. Finite element modeling of electrode–skin contact impedance in electrical impedance tomography. IEEE Transactions on Biomedical Engineering 1993;40(4):335–43. [8] Vauhkonen PJ, Vauhkonen M, Kaipio JP. Three-dimensional electrical impedance tomography based on the Complete Electrode Model. IEEE Transactions on Biomedical Engineering 1999;46(9):1150–60. [9] Soleimani M, Powell CE, Polydorides N. Improving the forward solver for the complete electrode model in EIT using algebraic multigrid. IEEE Transactions on Medical Imaging 2005;24(4):577–83. [10] Brebbia CA, Dominguez J. Boundary Elements: An Introductory Course. 2nd ed. Southampton: WIT Press/Computational Mechanics Publications; 1992. [11] Aykroyd RG, Cattle BA. A boundary-element approach for the completeelectrode model of EIT illustrated using simulated and real data. Inverse Problems in Science and Engineering 2007;15(5):441–61. [12] Babaeizadeh S, Brooks DH. Electrical impedance tomography for piecewise constant domains using boundary element shape-based inverse solutions. IEEE Transactions on Medical Imaging 2007;26(5):637–47. [13] Babaeizadeh S, Brooks DH, Isaacson D. A 3-D boundary element solution to the forward problem of electrical impedance tomography. in: Proceedings of the 26th annual international conference of the IEEE EMBS. 2004, p. 960–3.

177

[14] Cheney M, Isaacson D, Newell JC. Electrical impedance tomography. SIAM Review 1989;41(1):85–101. [15] Beatovic D, Levin PL, Sadovic S, Hutnak R. A Galerkin formulation of the boundary element method for two-dimensional and axi-symmetric problems in electrostatics. IEEE Transactions on Electrical Insulation 1992;27(1): 135–43. [16] Sirtori S, Maier G, Novati G, Miccoli S. A Galerkin symmetric boundary element method in elasticity: Formulation and implementation. International Journal for Numerical Methods in Engineering 1992;35(2):255–82. [17] Gray LJ, Paulino GH. Symmetric Galerkin boundary integral formulation for interface and multi-zone problems. International Journal for Numerical Methods in Engineering 1997;40(16):3085–101. [18] Kallivokas LF, Juneja T, Bielak J. A symmetric Galerkin BEM variational framework for multi-domain interface problems. Computer Methods in Applied Mechanics and Engineering 2005;194(34–35):3607–36. [19] Aimi A, Diligent M. Numerical integration in 3D Galerkin BEM solution of HBIEs. Computational Mechanics 2002;28(3–4):233–49. [20] Gray LJ, Griffith BE. A faster Galerkin boundary integral algorithm. Communications in Numerical Methods in Engineering 1998;14(12):1109–17. [21] Fratantonio M, Rencis JJ. Exact boundary element integrations for two dimensional Laplace equation. Engineering Analysis with Boundary Elements 2000;24(4):325–42. [22] Hammer PC, Marlowe OP, Stroud AH. Numerical integration over simplexes and cones. Mathematical Tables and other Aids to Computation 1956;10: 130–7.