Plasma and Fusion Research: Regular Articles

Volume 11, 2401043 (2016)

Simulation of Contactless Crack Detection in HTS Films:
Application of H-Matrix Method to Fast Matrix-Vector
Multiplication∗)
Teruou TAKAYAMA, Ayumu SAITOH1) , Atsushi KAMITANI1) and Hiroaki NAKAMURA2)
1)

Department of Informatics, Yamagata University, Yonezawa 992-8510, Japan
Graduate School of Science and Engineering, Yamagata University, Yonezawa 992-8510, Japan
2)
National Institute for Fusion Science, Toki 509-5292, Japan(Received 30 November 2015 / Accepted 3 March 2016)

Detectability of the multiple cracks in a high-temperature superconducting (HTS) film by a scanning permanent magnet method has been investigated numerically. To this end, a numerical code has been developed forperforming a shielding current analysis in an HTS film with cracks. Furthermore, the GMRES(k) method andthe H-matrix method are adopted in the code. By using the code, the scanning permanent magnet method can besimulated successfully. The results of the computations show that the x-coordinate of the crack endpoint cannot
be estimated with increasing the magnet radius. On the other hand, the y-coordinate of the crack endpoint can bedetected accurately.
c 2016 The Japan Society of Plasma Science and Nuclear Fusion Research

Keywords: crack detection, FEM, high-temperature superconductor, permanent magnet, simulation
DOI: 10.1585/pfr.11.2401043

1. Introduction
A superconductivity is a phenomenon characterized
by the zero electrical resistance when materials are cooled
to the critical temperature or less. As is well known, ahigh-temperature superconductors (HTSs) are used for developing various devices and systems such as nuclear fusion reactor, flywheel, MRI, and etc., and they are characterized by some parameters. In particular, a critical current
density jC is one of the most important parameters, and it
is important to measure the value of jC accurately.
As a contactless method and a high-speed method formeasuring a critical current density jC in an HTS film,
Hattori et al. proposed the scanning permanent magnet
method [1]. In the method, a permanent magnet is moved
along a HTS surface. They found that a spatial distributionof jC can be estimated from the distribution of an electromagnetic force acting on the film.
In order to simulate a scanning permanent magnet
method, a numerical code was developed for analyzing
the time evolution of a shielding current density in an
HTS film with multiple cracks [2, 3]. After an initialboundary-value problem of the shielding current density is
spatially discretized with the finite element method (FEM),
the resulting ordinary diﬀerential system is solved by using the backward Euler method and the Runge-Kutta (R-K)
method with an adaptive step-size algorithm. In particular,
it is found that the R-K method is a useful tool with inauthor’s e-mail: takayama@yz.yamagata-u.ac.jp∗) This article is based on the presentation at the 25th International Toki
Conference (ITC25).

creasing the number of nodes for the FEM. However, it isnecessary to consider the further high-speed of the shielding current analysis for a large-scale problem.
The purpose of the present study is to investigate thedetectability of the multiple cracks in the HTS film. Inorder to analyze a shielding current density, we apply theGMRES(k) method and H-matrix method to the linear systems obtained from the Backward Euler method.

2. Governing Equations
In a scanning permanent magnet method, a cylindrical magnet of radius R and height H is located just above
an HTS thin film, and a distance between a magnet bottom
and film surface is denoted by L. We assume a squareshaped HTS film of length a and thickness b, and thesquare cross-section of the film is denoted by Ω. If a crack
is contained in the film, Ω has not only the outer boundary
C0 but the inner boundaries C j ( j = 1, 2, · · · , m).
In the present study, we use the Cartesian coordinatesystem O : e x , ey , ez , where the z-axis is parallel to thethickness direction and the origin O is the centroid of thefilm. In terms of the coordinate system, the symmetry axisof the coil is determined as (x, y) = (xM , yM ). For characterizing the magnet flux density at (x, y, z) = (xM , yM , b/2)
for v = 0 mm/s. The movement of the magnet is assumed
as xM (t) = vt − a/2, where v is the magnitude of the scanning velocity.
As is well known, a shielding current density j in anHTS is closely related to the electric field E. The relation
can be written as2401043-1

c 2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchPlasma and Fusion Research: Regular Articles
E = E(| j|)[ j/| j|].

Volume 11, 2401043 (2016)(1)

where u isu = μ0 ŴS (n) − (B(n+1) · ez  − B(n) · ez ).

As a function E( j), we assume the power law
E( j) = EC [ j/ jC ]N ,(2)where EC is the critical electric field and N is a constant.
By using the thin-layer approximation [4, 5], theshielding current density j can be expressed as
j = (2/b)∇S × ez ,(3)and the time evolution of the scalar function S (x, t) is governed by the following integro-diﬀerential equation:
μ0 ∂t (ŴS ) + (∇ × E) · ez = −∂t B · ez ,(4)where a bracket   is an average operator over the thickness of the HTS film. Also, Ŵ is defined by
Q(|x − x |)S (x , t)d2 x + (2/b)S (x, t),ŴS =Ω(5)and an explicit function Q(r) in (5) is expressed as follows:
Q(r) = −(πb2 )−1 [r−1 − (r2 + b2 )−1/2 ].

(6)

The initial and boundary conditions to (4) are assumed
as follows:
S = 0,∂S= 0,∂s h(E) ≡on

C0(7)

on

Ci

(8)

E · tdt = 0.

(9)

CiHere, s and n are an arclength along Ci and a normal unit
vector on Ci , respectively.

3. Numerical Methods
3.1Discretization

In the present study, we discretize the initialboundary-value problem of (4) in time and space using
the finite element method (FEM) and the Backward Eulermethod, respectively. A region Ω is divided into a number
Ne of square elements. On the other hand, a crack is given
by a line segment.
By using the Backward Euler method, the initialboundary-value problem of (4) is transformed into the nonlinear boundary-value problem as follows:⎧⎪= 0,in Ω(10a)G(S (n+1) )
⎪⎪
⎪
⎪
⎪
(n+1)
⎪
S= 0,on C0(10b)
⎪⎪
⎪
⎪
⎪
⎪
⎨∂S (n+1) /∂s= 0,on Ci(10c)(∗1) ⎪⎪⎪
⎪
(n+1)
⎪⎪
E· tds = φi ,(10d)⎪⎪⎪⎪Ci⎪⎪⎪⎩N(E(n+1) )
= 0.
(10e)
Here, a superscript (n) is a value at time t = tn (≡ nΔt). Inaddition, G(S(n+1) ) is defined byG(S (n+1) ) ≡ μ0 ŴS (n+1) + Δ(∇ × E(n+1) ) · ez − u, (11)(12)

Note that numerical solution contains the error due towhich a boundary condition (9) is not satisfied. In order to resolve this problem, we propose the virtual voltage
method [5]. In the method, we use a numerical boundary condition (10e), where N(E) is evaluated from h(E).
Therefore, the solution of the problem (*1) include not
only S (n+1) but also φi .
When the FEM and the Newton method are applied tothe problem (*1), we get simultaneous equations system:⎤⎡ ⎤ ⎡⎤⎡⎢⎢⎢ A(S) C F(φ)⎥⎥⎥ ⎢⎢⎢δS⎥⎥⎥ ⎢⎢⎢ g(S, φ)⎥⎥⎥⎥⎢ ⎥ ⎢⎥⎢⎢⎢⎢ C TOO ⎥⎥⎥⎥⎥ ⎢⎢⎢⎢ λ ⎥⎥⎥⎥ = ⎢⎢⎢⎢ 0 ⎥⎥⎥⎥⎥ ,(13)
⎢⎢⎣⎢
⎥
⎢
⎦
⎦
⎣
⎦
⎣
ODT (S) O
δφ−h∗ (S)at each time step of the Newton method. Here, S is then-vector, and the n is a number of nodes for the FEM. The(m−n)-by-n matrices A(S)(≡ W + J(S)) is determined fromŴ, where J(S) is the Jacobian matrix. The matrix C, F(φ),D(S), and the vector g(S, φ) are obtained by the cracks.
As a result, the initial-boundary problem of (4) is reduced
to the simultaneous equations system at each time step ofthe Newton method. The procedure for the system is asfollows:
1. By solving the simultaneous equations, we get thecorrection δS and δφ.
2. We update an approximate solutions by means ofS( j) ⇐ S( j−1) + δS and φ( j) ⇐ φ( j−1) + δφ.
Two procedures are repeated until the termination condition:max (||δS||/||S||, ||δφ||/||φ||) ≤ ε,(14)is satisfied. Here, ε is a constant, and || || denotes themaximum norm. Under the numerical method, we develop the numerical code for analyzing the time evolutionof the shielding current density in an HTS film with multiple cracks.

3.2 High-speed method for shielding current
analysis
In this subsection, we propose a high-speed analysis
of a shielding current density in an HTS film. To this end,
we use the GMRES(k) method as a solver of simultaneous
equations (13), and we measure the CPU time for the simulation of a scanning permanent magnet method Here, aconvergence criterion of the GMRES(k) method is denotedby εG .
As is well known, a matrix-vector product is calculated at the each iteration in the GMRES(k) method. Inorder to calculate the matrix-vector product with a further
high-speed, we apply the H-matrix method [6] with theadaptive cross approximation [7–9] to the matrix W. In
the H-matrix method, the computation of the matrix-vector
product can be reduced by compressing the matrix.

2401043-2

Plasma and Fusion Research: Regular Articles

Volume 11, 2401043 (2016)

Fig. 2 Dependence of the CPU time on the restart coeﬃcient k
for the case with n = 2601.

Fig. 1 Residual history of the Newton method for the case with(a) xM = −9.9 mm and (b) xM = −4.3 mm. Here, k = 500
and n = 2601.

Throughout the present study, the physical and geometrical parameters are fixed as follows: a = 20 mm, b = 1µm, jC = 2 MA/cm2 , EC = 1 mV/m, N = 20, L = 0.5 mm,
H = 2 mm, BF = 0.3 T, v = 10 cm/s, and ε = 10−7 . In addition, we suppose the two cracks, C1 and C2 , are includedin the HTS film, and its shape is a line segment which isparallel to the x-axis. The crack distance d and the crack
size Lc are given by d = 4 mm and Lc = 1.6 mm.
Firstly, we investigate the influence of the Newton
method on the convergence determinant of GMRES(k)method. In Figs. 1, we show the residual history of theNewton method. We see from Fig. 1 (a) that, for εG = 10−1 ,the value of the residual norm almost becomes constant,whereas the Newton method converges in 11 and 4 iterations for εG = 10−2 and εG = 10−3 , respectively. On theother hand, we cannot obtain the solution at this step for
εG = 10−3 (see Fig. 1 (b)). In addition, for εG = 10−4 andεG = 10−5 , the Newton method converges in 4 iterations.
Consequently, it is necessary that the convergence determinant is εG ≤ 10−4 . In the following, the value of εG isfixed to εG = 10−5 .
Next, we investigate the influence of the restart coeﬃcient k on the CPU time. In Fig. 2, we show the dependenceof the CPU time on the restart coeﬃcient k. From this figure, the CPU time decreases with increasing the value ofk monotonously for 25 ≤ k ≤ 75, and it almost becomesa constant. Especially, the CPU time drastically increases
for k = 25. As a result, it is found that the CPU time decreases by using the large value of k. Hereafter, we use k =500 as the the value of the restart coeﬃcient for analyzingthe shielding current density.
We found that, by using the GMRES(k) method to asa linear-system solver, the time evolution of the shielding
current can be analyzed with high speed. In addition, thespeed of the GMRES(k) method can be accelerated by using the H-matrix method.

4. Detectability of Multiple Cracks
In this section, we numerically investigate the detectability of the multiple cracks by the scanning permanent magnet method. For this purpose, we adopt a defectparameter ds [5] as usual. By using the defect parameter,a relatively value can be obtained near the cracks, and its
value almost vanishes around the film edge.
Firstly, we investigate the distribution of the defect parameter ds . In Figs. 3, we show the contour maps of thedefect parameter ds for the two types of the magnet radius
R. From Fig. 3 (a), ds has a small value on both sides of
the two cracks, whereas the maximum value of ds is taken
near the center of the HTS film. Figure 3 (b) indicates that
we obtain a region of the small value by using the smallradius. These results mean the crack position can be estimated from the smallest value of the defect parameterds . In Figs. 3, we also show the smallest value of ds bya rhombus symbol in the first quadrant. Hereafter, the xycoordinates of this symbol are denoted by (x, y) = (xs , ys ).
The results of the computations show that we get (xs , ys ) =(2.3 mm, 2 mm) and (xs , ys ) = 0.67 mm, 2 mm for the casewith R = 3 mm and R = 0.5 mm, respectively.
Next, let us investigate the detect accuracy of the multiple cracks numerically. To this end, as a measure of theaccuracy, we define relative errors: e x ≡ |xc∗ − xs |/xc∗ and2401043-3

Plasma and Fusion Research: Regular Articles

Volume 11, 2401043 (2016)

Fig. 4 Dependence of the relative error e x on the magnet radiusR.

creases by the large value of the restart coeﬃcient k.
2. The x-coordinate of the crack endpoint cannot be
fairly evaluated with increasing the magnet radius.
On the other hand, the y-coordinate of the crack endpoint can be detected accurately.
In the future, we plan to compare the numerical results andthe experimental ones next year.

Acknowledgment

Fig. 3 Contour maps of the defect parameter ds for (a) R = 3mm and (b) R = 0.5 mm.

ey ≡ |yc − ys |/yc , where xc∗ is defined by xc∗ ≡ Lc /2. The
relative error ε x is calculated as a function of the magnet
radius R and is depicted in Fig. 4. We see from this figurethat the relative error e x decreases with R for R ≤ 1 mm,and its value monotonously increases with R for R ≥ 1mm. On the other hand, the value of ey becomes ey = 0%
for 0.5 mm ≤ R ≤ 3 mm, and the y-coordinate of the crack
position can be detected with a high accuracy.

5. Conclusion
Conclusions obtained in the present study are summarized as follows:
1. The Newton method does not converge at a certain
time step when the convergence determinant of theGMRES(k) method is greatly larger than that of theNewton method. Also, we found the CPU time de-This work was supported in part by Japan Society for
the Promotion of Science under a Grant-in-Aid for Scientific Research (C) (No. 26520204 and No. 15K05926)
A part of this work was also carried out with the supportand under the auspices of the NIFS Collaboration Researchprogram (NIFS15KNTS041, NIFS15KNXN297). Moreover, the numerical computations were carried out on Fujitsu PRIMEHPC FX100 of the LHD Numerical Analysis
Server in NIFS.
[1] K. Hattori, A. Saito, Y. Takano, T. Suzuki, H. Yamada, T.
Takayama, A. Kamitani and S. Ohshima, Physica C. 471,1033 (2011).
[2] T. Takayama and A. Kamitani, IEEE Trans. Appl. Supercond. 23, 9001107 (2013).
[3] T. Takayama and A. Kamitani, IEEE Trans. Appl. Supercond. 24, 9001505 (2014).
[4] A. Kamitani and S. Ohshima, IEICE Trans. Electron. E82-C,766 (1999).
[5] A. Kamitani, T. Takayama and S. Ikuno, Physica C. 494, 168(2013).
[6] M. Bebendorf, Numer. Math. 86, 565 (2000).
[7] S. Kurz, O. Rain and S. Rjasanow, IEEE Trans. Magn. 38,421 (2002).
[8] Y. Takahashi, C. Matsumoto and S. Wakao, IEEE Trans.
Magn. 43, 1277 (2007).
[9] V. Le-Van, B. Bannwarth, A. Carpentier. O, Chadebec. J. M.
Guichon and G. Meunier, IEEE Trans. Magn. 50, 7010904(2014).

2401043-4