E-Journal of Advanced Maintenance Vol.5-2 (2013) 113-120
Japan Society of Maintenology

On the Efficiency of Direct Sparse Solver for FEM-ECT Analysis
Shigeru HEMMI 1,*, Masayuki SARUHASHI 1 and Yasushi IKEGAMI 1
1

Itochu Techno-Solutions Corporation, Kasumigaseki Bldg, 3-2-5 Kasumigaseki, Chiyoda-ku, Tokyo
100-6080, Japan
ABSTRACTIn the analysis of time harmonic eddy current testing (ECT) by FEM, in which non-cored exciter/detector coils
are moved above the surface of metal object, the use of direct sparse solver has been studied. And it is foundpractically efficient because the factorized numerical matrix can be reused repeatedly in the calculation. The
validity of this method is discussed and a satisfactory result for an ECT model is shown.
KEYWORDS

ARTICLE INFORMATION

Eddy Current Testing, ECT, Finite Element Methods, FEM, Direct Sparse
Solver, MUMPS, Edge Element, Coulomb Gauge

Article history:
Received 9 November 2012
Accepted 10 June 2013

1.

Introduction

In Finite Element Method (FEM) model for the analysis of Eddy Current Testing (ECT),
especially for the case of time harmonic type, we have to treat sparse simultaneous equations
composed by complex numbers, that bring about a larger scale calculation. To solve such equations
efficiently, many authors tend to prefer iterative numerical solver, e.g., ICCG method[1,2,3,4].
However, in this paper, it is emphasized that direct sparse solver is also practically efficient, if the
problem size is not very large, and the moving coils do not involve ferromagnetic core, because thefactorized matrix for the first step in stepwise coil movement can be reused efficiently in thecalculation routines of following steps.

2.

Formulation and Direct Sparse Solver

2.1. Assumptions and Implementation Detail
In this note, only forward analysis is examined, not for inverse analysis. The phase convention oftime harmonic factor exp jt  is assumed where  is angular frequency. Four vector potentials
A,   are used as fundamental electromagnetic variables in SI units. The potential A isapproximated by edge-based functions[5,6,7], whereas  is approximated by nodal-based functions.
Only the lowest order shape functions are applied for both edge variables and nodal variables. Since
the electric displacement current is neglected, nodal  within space meshes are not included; electricfields in space region are assumed to be dependent quantities of other variables. For the numericalsamples herein, no boundary conditions applied except explicitly specified. Three types of FEM
element libraries are implemented in our code, i.e., hexahedra, penta(prism) and tetra. Tetra elements
are not used within this note because of its lower accuracy.
In our formulation, some aspects differ from other authors. Some differences excepting direct
solver usage, are summarized as follows: Total vector potentials are adopted; some authors are using reduced vector potentials[2].
 Since the space or vacuum is a type of electromagnetic substances, the space involving bothsurrounding region and crack region are also meshed. Contrary to ours, crack region is not meshedby some authors[4].
 Since meshing exciter/detector coils together with space and conductor as an integrated FEM meshis not an easy task, coil portions are not meshed as an integrated FEM mesh. Instead, equivalent
electric currents are applied to FEM edges as a right hand side of linear simultaneous equations.
 The authors prefer Coulomb gauge fixing[8], i.e. divA  0 (or rigorously   0 & divA  0 ), rather*Corresponding author, Tel: +81-3-6203-7450, E-mail: magna@ctc-g.co.jp
ISSN-1883-9894/10 © 2010 – JSM and the authors. All rights reserved.

113

S. Hemmi, et al./On the Efficiency of Direct Sparse Solver for FEM-ECT Analysis

than   0 gauge (or rigorously   0 & divA  0 gauge).
Although the   0 gauge is not popular in conventional books of electromagnetism, samples of  0 gauge are easily seen if the gauge transformation is applied to the problems in which Coulomb
gauge solutions are available. For example, consider an electric capacitor problem of static voltage
exerted. For this sample,   0 and A  0 is obtained by solving, e.g., Laplace equation, with properboundary conditions, and then if the gauge transformation is applied to this static problem so that it
satisfies   0 , then divA  0 is obvious in it. The   0 gauge is sometimes seen in other author’s
work of eddy current analysis, for example, see [2,9].
2.2. Ampere-Maxwell Equation in Weak Form and Singular Matrix of Edge Based Elements
The starting equation is Ampere-Maxwell equation expressed in term of four vector potentials
A,   . Then standard procedure of FEM is applied, i.e. multiplication with some weight functions
and integration over the entire domain. This is a weak form of Ampere-Maxwell equation.

 w  rot rotA    jA  grad dV   w  JExciteCoil

dV

(1)



Here symbols  ,  stand for magnetic reluctivity and electric conductivity respectively. The symbol
J Excite Coil is current density supplied to the exciter coils and dV is volume of the differential elementand  is volumetric area of the entire domain involving both space and conductors.
For the weight function w , the two weight functions are applied simultaneously. One is w  A
that satisfies divA   0 and the other is w  grad   that satisfies rot grad    0 . The latter one
corresponds to the conservation law of electric current.
For a hexahedral edge element of rigid rectangular shape, the divergence in the approximatedfield is zero, because the edge shape functions themselves may have a zero divergence [7]. However,
if hexahedral shape is deformed and not rectangular shape, divergence free feature of fields may be an
approximation. Then the weight function can be expected to satisfy divA   0 .
Approximate divergence free feature of the edge shape functions is a fundamental basis ofdivergence free nature of numerical solution of entire domain. However, it is not sufficient condition
to get divergence free numerical solution. It is a well known fact that if gauge fixing is abandoned inedge element method, the matrix becomes singular one in the course of numerical process[4]. This
causes a difficulty upon iterative and also direct solvers. Therefore, for this code using edge-based
elements the diagonal components in assembled matrix are changed by adding very small values, inorder to avoid the occurrence of singular matrix.
2.3. Right Hand Side of Simultaneous Equations
Since the total vector potentials are used instead of reduced potentials in the formulation, thecalculation method for right hand side (r.h.s.) of equation (1) is necessary to be clarified. The
exciter/detector coils are divided into “meshes” but these “meshes” are not directly used in thesimultaneous equations of FEM. The r.h.s. of equation (1) is calculated through Biot-Savart
integration so that FEM edges have proper and equivalent exciting currents.
Paying attention to exciter coils system, the following equation is an identity based on
Biot-Savart law,

 1 J
r 
J ExciteCoil  rot H ExciteCoil  rot   ExciteCoil
dV 3r 4 


(2)where J ExciteCoil is the prescribed current density on exciter coils and H ExciteCoil is the corresponding
magnetic field.

114

E-Journal of Advanced Maintenance Vol.5-2 (2013) 113-120
Japan Society of Maintenology

When a discrete mesh structure of the model is made, the above can be rewritten as follows to
obtain an approximate value.

 1 J

r
dV
J edge  rot   ExciteCoilr3 4 


(3)

where J ExciteCoil is defined only for exciter coil region. However, in the case of non-FEM meshes,
any unknown variable are not assigned even where JExciteCoil resides. On the other hand, in case of J edge
it is defined on all edges of FEM-elements and unknown variables of simultaneous equations are
assigned. The calculated column of J edge from this equation (3) forms the r.h.s. of simultaneous linear
equations of electromagnetic FEM based on total vector potentials.
2.4. Electromotive Force or Impedance of Detector Coil
The voltage or electromotive force induced along the detector coil wire, is described as

Ve.m. f   j  A  dl

(4)

where A is the vector potential on the integration line and dl is the length of discrete element
along the integration line of detector coil.
The vector potentials A in equation (4) are calculated as follows. For each FEM elements, thevector potentials at the element centers are obtained by the formula

Ae 



Ai N i e(5)iElementEdges

where the subscript e stands for element center, i runs over all edges of the specific element, Ai
is the edge potential as the result of FEM computation, and N i e is the value of vector shape function
at the element center. The explicit form of vector shape function N i is given in [7]. Calculating thevector potential value at an arbitrary point, i.e., a point on the line integration line, is a problem ofinterpolating scattered data. In our case, interpolation is completed using above scattered element
center values by the formula:A   we A e /  we
e

(6)

e

where e runs over all FEM element center and we are positive scalar weights for interpolation. In
this code, we  1 / r 4 is chosen to enhance the contributions of closer elements much more, where r
is the distance of element center and integration point. There is no deeper inspection of our choice of
weight function; it was chosen because it is simple, robust and easy to code.
The change in induced voltage on detector coil is evaluated by the difference between values
obtained for plates with crack and without crack.
2.5. Direct Sparse Solver
In numerical ECT analysis, problems often tend to be a larger size. One reason is that the crack
size is so small that it requires finer mesh at the vicinity of cracks. The second reason is that
exciter/detector coils are also small and move around the metal plate to search flaws. Consequently
the mesh should be considerably fine for almost all regions involving metal plate and coil movement.
The use of sparse direct solver, in the context of three-dimensional FE modeling of Maxwell
equations, has been limited due to the requirement for larger memory size and longer CPU time.
Since CPU execution time of numerical factorization is proportional to DOF3, practical limit of theefficiency of direct sparse solver is a severe obstacle. That is one of the reasons why most authors
prefer iterative solvers [1,2,3,4,9]. Nevertheless, we choose direct solver here because of its

115

S. Hemmi, et al./On the Efficiency of Direct Sparse Solver for FEM-ECT Analysis

robustness for such intricate computation and growing capability of recent computer hardware making
the direct solver applicable within a reasonable size. Within our limited hardware, 2 million complex
DOF is close to the limit of our calculation at this moment.
As an efficient direct sparse solver, we chose an open source sparse direct solver, MUMPS [10].
MUMPS itself can work in MPI environment, which can be another measure expanding the model
scale, but the function is not utilized here.

3. Code Validation
Some sample problems are described here to validate the code. The first one (A.) is the problemof an infinitely long and axially symmetric solenoid coil as a sample having rigorous theoreticalsolution. The second one (B.) is a comparison between perfect wire ring and ring with a cut, both ofwhich are exerted by a low frequency magnetic field. This sample highlights a significance of the grad term in Coulomb gauge. And the last one (C.) is a benchmark on ECT, where an
impedance profile is obtained as the result and compared with that of a published experiment.
A. Infinite Axially Symmetric Solenoid
Here an infinite solenoid coil problem is solved. The problem is a mimic of the sample in
Adventure project [9]. As seen in Fig. 1, an infinite cylindrical conductor is surrounded with an
infinite solenoid coil. For this problem, theoretical result is given in the book [11]. In this case only
one segment region of the whole model can be the object of analysis due to axial symmetry, andmodeled in FEM as depicted in Fig. 2. Here, hexahedral and penta element types are used. Boundary
restraints are imposed upon two side faces such that A  n  0 and   0 where n is thenormal vector on the surface.

Fig. 1. Axisymmetric infinite solenoid problem.

Fig. 3. Imaginary part of eddy current densities

Fig. 2. FEM segment model of infinite solenoid problem

Fig. 4. Theory and FEM results of magnetic flux densities

116E-Journal of Advanced Maintenance Vol.5-2 (2013) 113-120
Japan Society of Maintenology

The prescribed current density in the current region is 50.0 [A/m2], the frequency is 60[Hz] andthe electric conductivity is 7.7  106 S / m . Vector plot of calculated eddy current densities ofimaginary part is given in Fig. 3. A comparison of the theory and FEM result is given in Fig. 4 formagnetic flux densities against r coordinate. As in Fig. 4, the agreement of the theory and FEM
result is very good, showing validity of the basic code.
B. A Wire Ring Connoted Versus A Wire Ring Disconnected
Although the infinite solenoid coil and conductor problem (A.) is a good example of code
validation, there is still a room to discuss and notice. For rigorous axially symmetric problems, grad disappears because  r ,1 , z    r , 2 , z  for 1   2 . As the role of  grad term issignificant in Coulomb gauge, we must look forward to another example which clearly shows
importance of  grad in eddy current analysis.
If A.C. magnetic fields of uniform magnitude are exerted upon a metallic wire ring so that thedirection of magnetic fields are perpendicular to the ring plane (see Fig. 5), circulating “eddy” current
is induced as an introductory example of induction problem. And when the wire is cut at a point
circulating current disappears. These two problems are modeled as the second example of validation.
In the calculation, the inner diameter of ring is 21[mm], the ring has rectangular cross section of
1[mm] x 1[mm]. For exerted magnetic field, frequency is 50[Hz] and amplitude is 10000[A/m]. The
electric conductivity of ring wire is 60.0  106 S / m .
The entire mesh is not shown in Fig. 5, the mesh of metallic ring shown in the figure areembedded in the entire space mesh. As the entire space mesh, elements are composed of hexahedral
and penta types, number of nodes is 67,002, number of edges 204,097 and complex DOF is 204,401.

Fig. 5. Metallic wire ring example.

Fig. 6 Circulating “eddy” currents if connected

Fig. 7b. Imaginary part of  if disconnected

Fig. 7a. Eddy currents disappears if disconnected

Fig. 6, Fig. 7a and Fig. 7b are the FEM results. Fig. 6 and Fig. 7a are vector plots of eddy current
densities for the connected ring (Fig. 6) and for the disconnected ring (Fig. 7a). The “eddy” current iscirculating in the connected case, while it disappears in the disconnected case. And evident  grad
occurs as shown in Fig. 7b for the disconnected case. On the other hand, for the connected case, the grad term does not occur due to axial symmetry. Then the figure of contour plot of  is trivial

117

S. Hemmi, et al./On the Efficiency of Direct Sparse Solver for FEM-ECT Analysis

and is not shown here.
For the connected case, a comparison of imaginary part of electric current density in a
representative element is given in Table 1. It shows good agreement with the theory.
Table 1. Comparison of imaginary part of current density in an element.
J imag [ A / m 2 ]Theory-1.214E6

FEM result-1.211E6

 almost cancels out with each other, resulting in
For the disconnected case,  grad and  Adisappearance of circulating eddy current. Such explicit cancellation is an obvious and significant
feature of Coulomb gauge. In ECT problems, crack or slit is always assumed to exist and suchcancellation mechanism should work in conductor regions especially in the vicinity of the crack or slit.
From this sample and the consideration, we can conclude that the magnitude of  grad term in
Coulomb gauge can’t be neglected in eddy current analysis of ECT problems.
C. An ECT Benchmark Problem of Single Pancake Type CoilAn experimental result of an ECT problem, in which single pancake type coil is placed above aflat metal plate involving a crack, is published in the literature [12] as a JSEM benchmark problem.
Fig. 8 illustrates the conditions of the problem. Among several cases given in [12], the case of arectangular slit residing in the metal surface is investigated here.
For the prescribed electric current, the frequency is 300[kHz], the amplitude is 8.0 [mA rms] andthe number of coil turns is 140, the electric conductivity of the plate (Inconel 600) is 1.0  106 S / m .

Fig. 8. An ECT problem of single pancake type coilFig. 9 FEM mesh and coil (essential part)Fig. 10. Experimental impedance change given in [12]Fig. 11. Calculated impedance change

Fig. 9 illustrates the meshed model, where only the essential part of model is shown for briefness.
For the entire mesh involving space mesh, the type of FEM elements is hexahedral only, number ofnodes is 603,032, number of edges is 1,782,142 and DOF is 1,866,994. Fig. 10 is the experimental
result referred to [12]. Fig. 11 is the FEM calculated result. Here only 1/4 area of xy plane iscalculated and the others are extrapolated based on the symmetry rule. As seen in Fig. 11, the shape of
impedance change profile is almost same to the experimental result of Fig. 10, although a small
difference is observed in peak values between two figures that is Z exp  0.121 while Z ours  0.160.

118

E-Journal of Advanced Maintenance Vol.5-2 (2013) 113-120
Japan Society of Maintenology

4. “TP01/A type/along crack” ECT problem – Numerical Experience
The Japan Society of Maintenology organizes round-robin tests and gathers non-destructive ECT
signals freely available for anybody on the website [13]. Since we are not on the stage to make an
inverse analysis and such a case of arbitrary shape crack made by tetrathionic acid is hard to model, it
is difficult to show some useful results here. However, as a starting point of future work and for anumerical experience, a following analysis is possible assuming rectangular crack with some specificsize. As an example, we chose “TP01/A type/along crack” ECT problem. Space, metal plate with slit
and exciter coils are modeled as an FEM mesh and the conditions are;Table 2. "TP01/A type/along crack" ECT problem
OS

Linux

RAM size

Nodes

657,018

Edges

24Gbytes Operation mode Out of core

Frequency[kHz]20,50,100

1,937,691

Element typeHexahedral

DOF(Complex)

2,080,071

The RAM size of 24G bytes was not enough to solve the problem with in-core mode. The
illustration of the model is given in Fig. 12 and a sample output of eddy current profile is shown in
Fig. 13. Electro-motive force with arbitrary scale is shown in Fig. 14.

Fig. 12. A wireframe view of entire FEM Model

Fig. 13 A snap shot of eddy current density

Fig. 14. Numerical result of electro-motive force in an arbitrary scale

The calculation is carried out for 20 coil positions as an example. Time consuming numericalfactorization of sparse matrix is done only at the first step position. For other steps of 2nd to 20th,
factorized matrix is reused repeatedly and calculation time is reduced as seen in Table 3. In Table 3,
65 minutes corresponds to numerical factorization while 38 minutes corresponds to Biot-Savart
integration and solve phase of simultaneous equations. As seen in Table 3, time consuming numericalfactorization is needed only at the first position of exciter coils.
The size of the problem in this section, is not extremely large, but it is near the end of the current
author’s hardware ability and is close to the realistic limit of our computation. If a bigger hardware isavailable for our work, our current practical limit, i.e., 2 million complex DOF, may be improved a

119

S. Hemmi, et al./On the Efficiency of Direct Sparse Solver for FEM-ECT Analysis

few factors. But we are guessing that improvement of several factors, i.e., 10 times or larger size, is
not within the reality at this moment since CPU execution time of numerical factorization is
proportional to DOF3.
Table 3. Elapsed machine time for moving 20 positions of exciter/detector coils
Calculation phase
Elapsed time
# of positions Total elapsed timeFirst position
103m (= 65m + 38m)
1
1h 43m
2nd~20th positions
38m (= 0m + 38m)
19
12h 02m
13h 45m

5. Conclusion
In spite of many authors in the field of low frequency electromagnetic FEM tend to prefer
iterative solvers for sparse simultaneous equations, the authors have paid attention to direct solver. It
is emphasized that direct sparse solver is enough efficient in FEM-ECT analysis as long as the size is
not very large (in short, problems having 2 million complex DOF are not extremely large within ourcurrent available hardwares). Especially for the case of moving exciter coils without ferromagnetic
core, reuse of a numerically factorized matrix has been proved much useful and can be a practicalmeans to reduce the computation time. This is an advantage of our approach compared withconventional works. While a disadvantage of direct solver is: that since CPU time of numericalfactorization is proportional to DOF3, long CPU time and large memory are required if problem size
getting larger.
It is also emphasized that the term  grad and the term  A in Coulomb gauge, have almost
close values and almost cancels with each other in the vicinity of crack and therefore neglect of term grad is not realistic in the framework of Coulomb gauge as discussed in the sample case of
disconnected wire ring.

Acknowledgement
We would like to express our sincere thanks to Mr. Osamu Kounosu for his valuable advice.

References[1] H. Fukutomi, T. Takagi, J. Tani, M. Hashimoto, J. Shimonoe and Y. Harada, “Numerical Evaluation of
ECT Impedance Signal due to Minute Cracks”, IEEE Transactions on Magnetics, vol. 33, no.2, pp
2123-2126, 1997.
[2] T. Takagi, H. Huang, H. Fukutomi and J. Tani, “Numerical Evaluation of Correlation between Crack Sizeand Eddy Current Testing Signal by a Very Fast Simulator”, IEEE Transactoins on Magnetics, vol. 34, no.5,
pp. 2581-2584, 1998.
[3] N. Yusa, Z. Chen, K. Miya, T. Uchimoto and T. Takagi, “Large-scale parallel computation for thereconstruction of natural stress corrosion cracks from eddy current testing signals”, NDT&E International
36 pp. 449-459, 2003.
[4] Y. Choua, L. Santandrea, Y. Le Bihan, and C. Marchand, Thin Crack Modeling in ECT With Combined
Potential Formulations, IEEE Transactoins on Magnetics, vol. 43, no. 4, pp 1789-1792, 2007.
[5] H. Whitney, “Geometric Integration Theory”, Princten, NJ: Princeton University Press 1957.
[6] J. C. Nedelec, “Mixed finite elements in R3”, Numer. Meth., vol. 35, 1980, 315-341.
[7] Jianming Jin, “The Finite Element Method in Electromagnetic”, JOHN WILEY & SONS, INC., 1993.
[8] As a review of “gauge fixing”, see the article of wikipedia: http://en.wikipedia.org/wiki/Gauge_fixing.
[9] Adventure project, “ADVETURE_Magnetic Non-Linear Magnetostatic / Time-Harmonic Eddy Current
Analysis with HDDM Version 1.3 User’s Manual”, 2007; See http://adventure.sys.t.u-tokyo.ac.jp/[10] P. R. Amestoy, I. S. Duff, J. Koster and J.-Y. L’Excellent, “A fully asynchronous multifrontal solverusing distributed dynamic scheduling”, SIAM Journal of Matrix Analysis and Applications, Vol 23, No 1,
pp 15-41 , 2001; See also http://graal.ens-lyon.fr/MUMPS/.
[11] For example, see Landau and Lifshitz, “Electrodynamics of Continuous Media”, Buttreworth Heinemann,
problem 3. Of page 207.
[12] J. Shimone, K. Maeda and Y. Harada, Impedance Measurement of Unique Shape Cracks by Pancake TypeCoil, JSAEM Studies in Applied Electromagnetics, 5, pp.12-17 ,1996.
[13] http://jsm.or.jp/jsm/at/scc/index_eng.htm

120